基于离散元法的杆系结构几何非线性大变形分析

齐 念 叶继红

(东南大学混凝土及预应力混凝土结构教育部重点实验室,南京 210096)

在工程实际分析时,经常会遇到由几何非线性所引起的结构大变形问题,这类问题的特点是伴随大位移和大转动,因此必须考虑变形对平衡的影响.对于杆系结构几何非线性大变形问题,有限元法是目前最为常用的分析方法,学者们也提出了多种计算格式,如完全拉格朗日(TL)法、更新拉格朗日(UL)法和协同转动法等[1-3].但是,这些传统方法在处理几何非线性行为时,需要不断集成和修正切线刚度矩阵,经常会遇到刚度矩阵奇异或非线性方程迭代求解不易收敛等问题,需要对方法本身进行特殊处理,导致计算效率偏低.

离散元法(discrete element method,DEM)最早由美国学者Cundall等[4]提出,现已发展成为计算散体力学领域中一种新的数值方法.该方法直接应用牛顿第二定律,单元之间采用弹簧系统连接,接触力与接触位移之间的关系构成了DEM的接触本构模型.由于不要求满足位移连续和变形协调条件,因此它特别适用于各种非连续、非均匀以及结构大变形和失效破坏等复杂过程及其机理的研究[5].最初DEM法的研究对象主要是散土体材料,如今已在许多工程领域如岩体边坡滑动[6]、钢筋混凝土梁的断裂模拟[7]等方面取得了较好的应用.

DEM法在杆系结构中的应用还较少.本文提出将DEM法用于求解杆系结构(二维、三维)的几何非线性大变形问题.对3个经典算例的静力大变形问题及非线性动力行为进行了模拟,结果表明DEM法适宜于处理结构大变形问题.

1 DEM法

DEM法的基本思想是把研究对象离散成刚性单元的集合.任取一个单元α,设有n个单元与其相邻,作用在单元α上的外力为Fext,外力矩为Mext.根据牛顿第二定律,其运动控制方程为

(1)

1.1 单元间的接触本构模型

图1 相邻单元之间的接触模型

在接触平面内将接触力和接触力矩按矢量法则进行分解,得

(2)

将F和M分别向3个坐标轴进行投影,得

(3)

利用离散元法计算接触力和接触力矩时采用增量形式.在时间步长Δt内,由位移增量所引起的接触力增量为

(4)

接触力矩的计算过程与接触力的计算过程相似.在时间步长Δt内,球A与球B之间的相对转角增量所引起的接触力矩增量为

(5)

(6)

1.2 时间步长的确定

在求解运动方程时,DEM通常采用显示中心差分算法.中心差分法是条件稳定算法,其计算时间步长Δt必须小于由该问题所决定的某个临界值Δtcr,否则算法将不稳定.

临界值Δtcr的确定应同时满足结构物理步长和数学计算步长的要求.物理步长是指结构中弹性波来回一次所需要的时间,框架结构的物理步长约为1 ms[9].数学计算步长是指积分计算所需的时间步长,对中心差分法而言,解的稳定性条件为

(7)

式中,ωn为结构的最高阶固有振动频率;m为结构质量;k为结构刚度.

本文在确定数学计算步长时,按如下方法粗略估计临界值Δtcr:将式(7)中的k值取为1.1节中提到的4个弹簧刚度系数中的最大值,即

k=max{Kn,Ks,kn,ks}

(8)

最后,通过比较物理步长与数学计算步长的临界值,取二者较小值,以确定时间步长Δt的大小.

2 动力与静力解及大变形问题

与传统的结构静力分析不同,离散元法本质上是求解结构的动力问题,因为式(1)的解本身就是结构的动力反应.

(9)

分析传统结构时,小变形和大变形问题是需要进行严格区分的.前者的几何方程为线性的,而后者的几何方程则为非线性的.利用有限元法处理大变形问题时,一般是对单元切线刚度矩阵进行修正,采用增量迭代的方法求解结构的响应.而利用DEM法求解此类问题时,并不需要刻意区分是小变形还是大变形,因为在建立运动控制方程时并不涉及到几何方程,也不要求位移连续,无需组集刚度矩阵和迭代求解,这与有限元法有着本质的区别.因此,将DEM法用于结构分析时,可以采用统一的步骤分析结构的小变形和大变形问题,其计算流程图见图2.

3 弹簧接触刚度系数的确定

根据DEM法的基本理论,利用式(4)和(5)计算接触力和接触力矩增量时,需要事先给定弹簧接触刚度系数的值.本文的研究对象主要是连续体——框架结构,而传统的DEM法弹簧接触刚度系数计算公式是基于散粒体的,不适用于框架结构[10].基于简单梁理论,将通过弹簧连接的2个圆球比拟成1根梁(见图3).图3中,RA,RB分别为球A和球B的半径;L=RA+RB为梁的长度.设结构构件的截面积为S,截面惯性矩为I,扭转惯性矩为J,弹性模量和剪切模量分别为E和G.

图2 DEM法的计算流程图

图3 弹簧接触刚度的等价梁模型

根据结构力学中杆端力与杆端位移之间的关系,可以得到弹簧接触刚度计算公式.梁的轴向刚度系数Kn和切向刚度系数Ks分别为

(10)

2个单元之间发生相对转动产生转角Δθ,Δθ可分解成平行于接触面内的切向分量Δθs和垂直于接触面内的法向分量Δθn.切向分量Δθs所引起的接触力矩,与一端固定一端滑动的梁所产生的杆端弯矩等效.因此,转动法向刚度系数为

(11)

垂直于接触面内的转角法向分量Δθn引起的则是扭矩.若是平面问题,扭矩为0;对于一般空间问题,扭矩可根据截面扭转刚度与扭转角相乘得到.因此,转动切向刚度系数为

(12)

图4 悬臂梁DEM计算模型

悬臂梁的变形和内力计算结果见表1.由表可知,根据DEM法计算的悬臂梁变形及结构内力与理论解相比,最大误差仅为0.3%,说明本文确定弹簧接触刚度系数及Δt的方法是正确合理的.

表1 集中力和力矩作用下的悬臂梁计算结果

4 数值算例

根据DEM法的基本原理及计算流程,编制了空间框架结构几何非线性大变形分析程序.通过对文献中常被引用的3个典型算例进行模拟和分析,验证本文方法的正确性和适用性.

4.1 平面正方形刚架的大变形分析

图5 平面正方形刚架及离散元模型

图6为利用DEM法计算得到的荷载-变形曲线.与文献[11]中采用椭圆积分计算的结果进行对比,发现位移及转角的计算值均与文献解相吻合,最大误差不超过1%.

图6 平面正方形刚架对边受拉下的荷载-变形曲线

4.2 悬臂梁的大弯曲问题

如图7所示,在悬臂梁端部作用一集中弯矩M(t).弯矩随时间逐渐增加,梁将会由原来的静止状态弯成圆形或螺旋形.该典型算例常被用于结构大变形问题的分析中[12-13].梁的相关几何参数及物理参数采用文献[12]中的数据,弯矩作用方式采用如图8所示的线性渐增荷载,DEM模拟时将梁离散成9个半径为1.25 mm的圆球.

图7 受端弯矩作用的悬臂梁

下面分2种工况对该梁的大弯曲行为进行数值模拟.

图8 端部弯矩作用方式

表2 弯矩作用下悬臂梁自由端计算结果

图9 不同时刻悬臂梁的形状

4.3 六角形空间刚架的大变形分析

图10 六角形空间刚架(单位:mm)

对该结构进行几何非线性静力大变形分析.表3为利用DEM法得到的刚架顶点荷载-位移计算结果.由表可知,与文献[14]中利用有限元法分析的结果进行对比,两者最大误差不超过0.9%.

表3 六角形空间刚架顶点的荷载-位移结果

5 结论

1) 本文将DEM法推广应用于杆系结构(二维、三维)的几何非线性大变形问题.基于简单梁理论,推导了适用于杆系结构分析的弹簧接触刚度系数计算公式,并利用算例对其正确性进行了检验.

2) 基于本文推导的刚度系数计算公式,对杆系结构的几何非线性大变形问题进行分析.列出了3个典型数值算例,即2个平面框架和1个空间网格结构,分别对其静力和动力大变形行为进行了模拟,并将结果与其他计算方法的结果进行比较,两者吻合良好.此外,离散元法在处理几何非线性时无需组集刚度矩阵,也不用迭代求解非线性方程,故该方法适宜于处理杆系结构的大变形问题.

)

[1] Teh L H, Clarke M J. Co-rotational and lagrangian formulations for elastic three-dimensional beam finite element [J].JournalofConstructionalSteelResearch, 1998,48(3): 123-144.

[2] Yang Y B, Lin S P, Leu L J. Solution strategy and rigid element for nonlinear analysis of elastically structures based on updated Lagrangian formulation [J].EngineeringStructures, 2007,29(6): 1189-1200.

[3] Thanh N L, Jean M B, Mohammed H. Efficient formulation for dynamics of co-rotational 2D beams [J].ComputationalMechanics, 2011,48(2): 153-161.

[4] Cundall P A, Strack O D L. A discrete numerical model for granular assemblies [J].Geotechnique, 1979,29(1): 47-65.

[5] 刘凯欣,高凌天. 离散元法研究的评述[J]. 力学进展,2003,33(4):483-490.

Liu Kaixin, Gao Lingtian. A review of the discrete element method [J].AdvancesinMechanics, 2003,33(4): 483-490. (in Chinese)

[6] James F H, David S C, William S P. Simulation of unstable fault slip in Granite using a bonded-particle model [J].PureandAppliedGeophysics, 2002,159(1/2/3): 221-245.

[7] Azevedo N M, Lemos J V, Almeida J R. A discrete particle model for reinforced concrete fracture analysis [J].StructuralEngineeringandMechanics, 2010,36(3): 1-19.

[8] Cundall P A.PFCuser’smanual(version3.1) [M]. Minneapolis, MN, USA: Itasca Consulting Group, Inc., 2004.

[9] 喻莹. 基于有限质点法的空间钢结构连续倒塌破坏研究[D]. 杭州: 浙江大学建筑工程学院, 2010.

[10] 徐小敏,凌道盛,陈云敏，等. 基于线性接触模型的颗粒材料细宏观弹性常数相关关系研究[J]. 岩土工程学报，2010,32(7): 991-998.

Xu Xiaomin, Ling Daosheng, Chen Yunmin, et al. Correlation of microscopic and macroscopic elastic constants of granular materials based on linear contact model [J].ChineseJournalofGeotechnicalEngineering, 2010,32(7): 991-998. (in Chinese)

[11] Kjell M. Numerical results from large deflection beam and frame problems analysed by means of elliptic integrals [J].InternationalJournalforNumericalMethodsinEngineering, 1981,17(1): 145-153.

[12] Wu Tongyue, Wang Renzuo, Wang Chungyue. Large deflection analysis of flexible planar frames [J].JournaloftheChineseInstituteofEngineers, 2006,29(4): 593-606.

[13] 丁承先,段元锋,吴东岳. 向量式结构力学[M]. 北京:科学出版社,2012: 273-278.

[14] Shugyo M. Elastoplastic large deflection analysis of three-dimensional steel frames [J].JournalofStructuralEngineering,ASCE, 2003,129(9): 1259-1267.