齐 念

(中冶华天南京工程技术有限公司，江苏 南京 210019)

颗粒离散元法(Discrete Element Method，DEM)是20世纪70年代由美国学者Cundall教授[1]首先提出的，该方法最初应用于土体、岩石等材料的力学性能分析。其基本思想是把研究对象离散为刚性单元的集合，单元与单元之间通过弹簧连接，接触力与接触位移之间的关系构成了DEM的接触本构模型。各个刚性单元满足牛顿运动方程，DEM法分析时允许单元发生相对运动，不要求必须满足位移连续和变形协调条件，尤其适合非线性问题的研究。目前DEM法已广泛用于计算散体力学领域[2,3]，如粉末的加工、研磨，药品的运输等。本文通过理论推导和数值分析，将DEM法推广用于框架结构几何非线性大变形问题。然后通过一个经典算例的大变形行为进行模拟和分析并与其他方法的计算结果进行比较。结果表明，DEM方法适用于框架结构大变形问题的分析。

1 颗粒DEM法基本原理

1.1 单元之间的接触力分析

颗粒离散元，以颗粒为基本研究对象，将物体离散作为具有代表性的数个单元，利用颗粒流模型构建物体的力学性质。以平面问题为例，假定单元形状为圆盘，单元之间的接触采用柔性接触，即允许接触处产生重叠，图1为任意两个相邻圆A和圆B间的接触模型，C点为接触中心。

当圆A与圆B之间产生接触时，接触重叠量Un可表示为：

Un=R[A]+R[B]-D

(1)

其中，R[A]，R[B]，D分别为圆A、圆B的半径及两圆之间的中心距。

在离散元方法中，相邻圆之间的相互接触作用力是通过接触弹簧来实现[4]。如图1所示，在接触平面内，接触力F可分解成垂直于接触面的法向力Fn和平行于接触面的切向力Fs，即：

F=Fn+Fs=Fn×n+Fs×s

(2)

其中，Fn，Fs分别为法向力和切向力的大小；n，s分别为接触平面的法向和切向单位矢量。

此外，接触点C处因两圆盘发生相对转动还会引起大小为Mz的接触力矩，这与散粒体离散元中的接触模型不同。

将接触力分别向坐标轴进行投影，得：

Fi=Fn×ni+Fs×si(i=1,2)

(3)

根据力系平移定理，将Fi移至各单元的中心处，有：

(4)

(5)

1.2 接触力和接触力矩的计算

接触力和接触力矩的计算一般是通过力与位移之间的关系来表达。对于法向接触力Fn：

ΔFs=-KnΔUn

(6)

切向接触力Fs，用增量的形式进行计算：

ΔFs=-KsΔUs

(7)

ΔUs=VsΔt

(8)

其中，Kn,Ks分别为弹簧法向刚度和切向刚度系数;ΔUs为相对切向位移增量;Δt为计算时步;Vs为接触点C处的切向速度。

对于接触力矩M3的计算，也采用增量形式来表达：

(9)

其中，Kθ为弹簧转动刚度;Δθ为相对转角增量;ω[B],ω[A]分别为圆A和圆B的转动角速度。

然后，按以下方式进行接触作用力的更迭计算：接触形成时，对Fn,Fs,M3初始化为零，然后进行下一时步的计算，将由相对位移引起的接触力增量、相对转角引起的接触力矩增量累加到现值上，即：

(10)

1.3 运动控制方程

(11)

对于运动方程的求解，DEM通常采用显示积分算法，本文采取中心差分格式。

2 算例验证

平面正方形刚架大变形分析。

如图2所示，一平面正方形刚架在对边受拉力作用下发生大位移大转动行为，采用DEM法模拟并与文献[5]中采用椭圆积分计算的结果进行对比。刚架各构件的材料和尺寸均相同，相关几何及物理参数为：L=10 cm，A=0.5 cm2，I=0.010 4 cm4，E=1.6×107N/cm2，材料为线弹性。DEM建模时，将刚架各边构件均离散成11个半径为1.0 cm的圆盘。计算时步Δt=1.0×10-3s，阻尼系数ζ取0.6。

表1为用DEM法计算得到的刚架在对边受拉情形下的荷载—变形结果，定义荷载无量纲因子：λ=PL2/EI。与文献中KJELL采用椭圆积分分析进行对比发现，两种方法计算的结果相当接近，最大误差不超过0.9%，表明本文方法具有较高的数值精度。

表1 平面正方形刚架对边受拉下的荷载—变形结果

图3为不同荷载大小作用下平面正方形刚架的最终变形图。从图3中看出，荷载值越大，结构变形越显著，用DEM分析时很容易获得结构的最终变形，因为它计算时会自动考虑几何非线性的影响，而不去区分是大变形还是小变形，因此求得的解就是它真实的状态。

3 结语

本文将散体力学中应用广泛的颗粒DEM法推广用于求解平面框架结构几何非线性大变形问题。算例分析结果表明：本文方法计算精度较高，对框架结构几何非线性大变形(大位移、大转动)问题能很好的处理。

DEM法是以牛顿第二定律为力学基础，在求解结构大变形问题时，不用组集单元的刚度矩阵，不需要迭代，流程相对简单。另外，该方法不要求满足位移连续和变形协调条件，突破了传统非线性有限元在结构非线性问题中遇到的困难，为今后求解结构大变形这类问题提供了一种新的研究思路和方法。

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

[2] 徐 泳，孙其诚，张 凌.颗粒离散元法研究进展[J].力学进展，2003，33(2)：251-260.

[3] 唐志平.激光辐照下充压柱壳失效的三维离散元模拟[J].爆炸与冲击，2001，21(1)：1-7.

[4] Itasca Consulting Group Inc.,Particle Flow Code in 2-Dimensions(PFC2D),Version 3.10,Minneapolis,Minnesota,2004.

[5] KJELL MATTIASSON.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.