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

齐 念 叶继红

(东南大学混凝土及预应力混凝土结构教育部重点实验室，南京 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.根据牛顿第二定律，其运动控制方程为

式中，m，J分别为单元α的质量和转动惯量;r，ω分别为单元α的位置矢量和角速度矢量;分别为相邻单元j对单元α产生的接触力和接触力矩;t为时间.

1.1 单元间的接触本构模型

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

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

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

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

式中，Kn，Ks分别为弹簧法向刚度和切向刚度系数;ni(i=x，y，z)为 n 在3个坐标轴上的分量;分别为法向和切向相对位移增量，可通过接触点C处的速度计算得到.

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

式中，kn，ks分别为弹簧转动法向刚度和切向刚度系数;分别为法向和切向相对转角增量.

1.2 时间步长的确定

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

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

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

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

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

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

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

对于需要计算静力解的问题，可采用阻尼耗能或缓慢施加外荷载的方式来逼近结构的静止状态.本文是在运动方程式中添加一项阻尼力，利用阻尼产生的逆向运动来达到削减结构动力响应的目的.这里采用局部非线性黏滞阻尼形式［8］，其本质是在单元的合外力基础上，总体施加一个与其方向相反的作用力.单元α的阻尼力可表示为

式中，为阻尼系数，取值范围为0～1;F*为单元α的广义力;V*为单元α的广义速度.需要说明的是，此处的阻尼并非结构的真实阻尼，它实际上是一个虚拟的行为，其目的是为了获得结构的静力解而选择的一种耗能机制.

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

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

图2 DEM法的计算流程图

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

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

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

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

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

利用1个空间悬臂梁对该方法进行验证.梁长L=200 mm，其截面呈圆环形，外径φ=8 mm，壁厚 δ=0.5 mm，弹性模量 E=200 kN/mm2，剪切模量G=80 kN/mm2.梁端部承受集中力F=10 N，弯矩M=100 N·mm，扭矩T=2 kN·mm，材料为线弹性.悬臂梁DEM计算模型如图4所示.在模拟时将最左端圆球固定，因为是求解静力问题，取阻尼系数 =0.7，时间步长 Δt=1 ms.

图4 悬臂梁DEM计算模型

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

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

4 数值算例

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

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

如图5(a)所示，平面正方形刚架各构件的材料和尺寸均相同，相关几何及物理参数为:l=10 cm，S=0.5 cm2，I=0.0104 cm4，E=16 MN/cm2，材料为线弹性，P为外力.DEM计算模型如图5(b)所示，球半径取为1 cm，时间步长Δt=0.5 ms，阻尼系数 =0.7，材料密度取为单位1.

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

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

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

4.2 悬臂梁的大弯曲问题

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

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

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

图8 端部弯矩作用方式

工况1 对梁进行缓慢加载，如图8(a)所示.取Δt=1 ms，=0.7.在该弯矩作用下梁最终将弯成一个圆.由于加载速率很慢，同时又有阻尼的耗能作用，因此实际上是模拟结构的近似静力变形行为.表2为利用DEM法得到的梁自由端计算结果，与解析解［13］相比，两者最大误差为 0.6%，说明本文方法具有较高的精度.

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

工况2 梁的几何及材料参数与工况1相同，但梁端弯矩加载方式如图8(b)所示，阻尼系数 =0.1.利用DEM 法对其进行模拟，图9是悬臂梁分别于时间 t=9，10，11，13 s时的形状，发现梁最后卷曲成螺旋状，而非工况1中的圆形.这是因为此时的计算结果为动态解，与所采用的阻尼系数和加载速率有关，若保持最终加载力不变，随着时间的延续，由于阻尼耗能，动态解最终也会趋于静力解，螺旋状也将趋于圆状.

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

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

对于如图10所示的六角形空间刚架，结构的6个支座均为铰约束，各构件截面均为正方形，其几何物理参数的取值与文献［14］保持一致.利用DEM法进行模拟时，将构件离散成24个球单元，取 Δt=0.1 ms， =0.7，材料密度取为单位 1.

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

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

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

5 结论

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

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

3)DEM法本质上是求解结构的动力行为，对于需要计算静力解的问题则是通过增加阻尼项来进行模拟的.阻尼系数 的选取应根据分析问题的性质合理确定，对于静力解问题，综合考虑数值精度和计算效率，建议可取 =0.7.由于DEM 法是采用中心差分法求解运动方程，计算时间步长Δt的值一般取0.1～1 ms.所提的时间步长临界值估算方法，可供结构分析时借鉴.

References)

［1］Teh L H，Clarke M J.Co-rotational and lagrangian formulations for elastic three-dimensional beam finite element［J］.Journal of Constructional Steel Research，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］.Engineering Structures，2007，29(6):1189-1200.

［3］Thanh N L，Jean M B，Mohammed H.Efficient formulation for dynamics of co-rotational 2D beams［J］.Computational Mechanics，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］.Advances in Mechanics，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］.Pure and Applied Geophysics，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］.Structural Engineering and Mechanics，2010，36(3):1-19.

［8］Cundall P A.PFC user's manual(version 3.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］.Chinese Journal of Geotechnical Engineering，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］.International Journal for Numerical Methods in Engineering，1981，17(1):145-153.

［12］Wu Tongyue，Wang Renzuo，Wang Chungyue.Large deflection analysis of flexible planar frames［J］.Journal of the Chinese Institute of Engineers，2006，29(4):593-606.

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

［14］Shugyo M.Elastoplastic large deflection analysis of three-dimensional steel frames［J］.Journal of Structural Engineering，ASCE，2003，129(9):1259-1267.