离散单元法在土壤机械特性动态仿真中的应用进展

第19卷第1期2003年1月农业工程学报T ran sacti on s of the CSA E V o l .19　N o.1Jan .　2003离散单元法在土壤机械特性动态仿真中的应用进展张　锐,李建桥,李因武(吉林大学)摘　要:简述了离散单元法的基本原理及其国内外发展现状,介绍了目前地面力学研究领域中离散单元法在土壤机械特性动态仿真中的应用情况,分析了离散单元法在土壤动态行为仿真中应用的可行性,指出离散单元法应用于土壤这样的多相不连续复合体系中,以离散单元的总体行为来描述土壤动态行为具有独特的优越性,提出了离散单元法在土壤动态特征研究中的发展趋势和近期需要解决的关键问题。

关键词:离散单元法;土壤动态仿真;力学模型中图分类号:O 242.26;S 152.9 文献标识码:A 文章编号:100226819(2003)0120016204收稿日期:2002204217　修订日期:2002206215基金项目:国家自然科学基金(50175045);教育部博士点基金(1999018504);教育部高等院校骨干教师基金([2000]6521)作者简介:张锐(1975-),男,博士研究生,长春市人民大街142号　吉林大学南岭校区生物与农业工程学院,1300251　引　言离散单元法(D istinct E lem en t M ethod ,简称D E M )是美国学者Cundall P .A .教授在1971年基于分子动力学原理首次提出并应用于分析岩石力学问题的一种不连续数值模拟方法[1]。

这种方法把介质看作由一系列离散的独立运动的粒子(单元)所组成,粒子的尺寸是细观的,其运动受经典运动方程控制,整个介质的变形和演化由各单元的运动和相互位置来描述。

1978年,Cundall P .A .和Strack O .D .L 开发出了二维圆形块体的B all 程序,用于研究颗粒介质的力学行为[2]。

自此离散单元法作为一种新兴的非连续体分析方法在岩石力学[3]、结构分析[4]、物料流动[5]及地面力学[6～8]等领域的数值模拟中得到了广泛应用。

1986年,Cundall P .A .与ITA SCA 咨询集团公司合作开发了三维离散单元程序3D EC [9],并且有了一些将三维离散单元法应用于工程实际中的尝试[10]。

由于三维问题数据结构和算法复杂,总体来说目前仍处于发展阶段。

离散单元法在国内起步相对较晚,王泳嘉等[11]于1986年将此法引入我国后,在采矿工程、岩土工程及水利水电工程等领域发展迅速。

但是,国内将离散单元法应用于土壤力学及土壤与机械相互作用关系等方面的研究目前鲜见报道。

作为多相混合体的土壤本身即是离散结构,被机械部分或植物根系切断或分离的土壤更是不连续的,因而将其作为离散颗粒的集合体用离散单元法进行分析,在土壤力学及土壤机械特性动态仿真的基础理论研究和应用技术开发方面具有十分重要的意义。

2　离散单元法的基本原理离散单元法是建立在最基本的牛顿第二定律基础之上的,具有牢固的理论依据。

这里以二维单元为基础,介绍离散单元法的基本原理和方法。

按离散体中单元的形状将其分为2类,一类为圆盘单元,另一类为多边形单元。

根据土壤颗粒的大致形状以及为了研究的方便,目前普遍将用于分析土壤的离散单元法中的单元视为圆盘形[12,13]。

本文主要针对日本学者的理论(D E M )对圆盘单元的离散单元法进行论述。

2.1　单元的接触力学模型每个单元都会从相邻接触单元获得触点压力。

力的大小是由单元的相对位移和相对速度来确定的。

在离散元仿真中,单元之间接触的弹性和非弹性性质用弹簧和阻尼器来表示。

弹簧代表单元的弹性,阻尼器代表单元的非弹性。

假设两个单元i 和j 之间存在法向弹性常数为k n 、切向弹性常数为k s 的一个弹簧,法向阻尼系数为Γn 、切向阻尼系数为Γs 的一个阻尼器,摩擦系数为Λ的一个滑块以及代表一个单元和其他单元之间没有拉力的非张力联接,单元之间的力学关系如图1所示[14,15]。

图1　离散元触点压力模型F ig .1　M odel of con tact fo rce in the D E M在离散元触点压力模型中,单元是靠与其相邻单元之间接触的“重叠”相互作用并产生触点压力的。

根据单元与单元之间的几何关系对单元的接触进行判断,即如果两个单元中心的距离比两个单元半径之和小,就认为这两个单元接触;反之,两单元不接触。

目前查找相邻单元的搜索方法有循环遍历搜索法、窗口法和区域法等[16],通常根据具体情况选择合适的搜索方法。

当两单元满足接触条件时,在t 时刻单元j 作用在单元i 上的法向弹性力增量∃e n 和切向弹性力增量∃e s ,法向非弹性力增量∃d n 和切向非弹性力增量∃d s 分别由图1所示的力学模型可以得到[17,18]:61∃e n=k n ∃u n,∃e s=k s ∃u s(1)∃d n=Γn ∃u n ∃t,∃d s=Γs ∃u s ∃t(2)式中　∃t——时间步长;∃u n——在∃t期间单元i在法向方向上的位移增量;∃u s——在∃t期间单元i在切向方向上的位移增量。

在t时刻,法向触点压力f n和切向触点压力f s分别为:f n=e n+d n,f s=e s+d s(3)式中　e n——t时刻的法向弹性力;e s——t时刻的切向弹性力;d n——t时刻的法向非弹性力;d s——t时刻的切向非弹性力。

由模型可知,法向触点压力不能为拉力,即当∃e n< 0时∃e n=∃d n=0(4)切向触点压力不能大于摩擦系数为Λ的摩擦力,即当f s> Λ ∃e s 时f s=sign(f s) Λ ∃e n(5)式中　Λ——单元之间的摩擦系数;sign(X)——符号函数。

2.2　单元的运动学模型在离散单元法中,每个单元的位置是在时间间隔∃t 中逐步确定的。

单元的位置和触点压力决定了这个单元的运动。

根据牛顿第二定律可得离散单元法的运动方程为: F=f+f body+f bound=m d2u d t2(6)M=m c+m bound=I dΞ d t(7)式中　F——单元所受到的合力;M——单元所受到的合力矩;f——单元之间的触点压力;f body——体积力;f bound——边界作用力;m c——触点压力的作用力矩; m bound——边界作用力矩;m——单元的质量;u——单元的位移;I——单元的转动惯量;Ξ——单元的角速度。

综合来说,离散单元法中所进行的计算是在单元上牛顿第二定律和力—位移接触定律应用的交替进行。

牛顿第二定律给出了作用在单元上的力所产生的这个单元的运动,力—位移接触定律计算了单元的位移所产生的触点压力。

3　土壤的离散单元仿真模型考察土壤和机械或植物根系之间的力学交互作用对设计和制造地面机械或评价耕作方法是非常重要的。

在涉及机械与土壤相互作用的研究领域,已经应用了许多连续介质理论,应用这些理论方法虽然可以描述土壤的基本力学性质,但是还不能准确地确定机械与土壤相互作用的力学关系和动态过程,这是因为连续介质理论较适合研究物体间的力学特性,无法描述在外力作用下土壤内部颗粒之间的相互作用以及土壤的剪切行为。

土壤由不同粒径的固体粒子堆砌而成,具有许多开式或封闭的孔隙和管道,其间存在空气或液体,它们在外力作用下,比土壤固体粒子更容易产生扩散迁移并影响着土体的整体结构和力学性能[19],所以假定土壤是连续体很难去仿真土壤颗粒的整体力学行为。

有限元(FE M)或边界元(B E M)经常被用来分析土壤和机械之间的相互关系[20],这些分析需要假设土壤是连续体,在不设置复杂边界条件的情况下,很难模拟土壤复杂的动态行为。

而将土壤作为离散颗粒的集合体更符合土壤的本质,用离散单元法去仿真伴随分离的土壤变形,将使在简单条件下模拟土壤的复杂行为成为可能,并且这种方法要比将土壤作为一个连续体来处理更加合理。

特别是如果在离散单元模型中能够考虑土壤颗粒的粘性特性的话,离散单元仿真计算结果会更接近实际情况。

土壤颗粒之间存在影响单元集合的粘附现象,传统的离散元力学模型难以描绘粘性壤土和农业生产常见的粘湿泥土的力学行为。

在传统的离散元模型中,非张力联接的影响是用式(4)来表示。

这个公式意味着两个接触单元之间的张力完全没有考虑,即此时法向触点压力f n为零。

但是在土壤与机械相互作用时,土壤的破碎实际上包含拉伸断裂[21](如图2所示),因此在土壤离散元模型中用式(8)代替式(4):当0<D t<D f时f n=e t+d n(8)式中　D t——两个单元边界之间的距离;D f——拉伸断裂的临界距离,用单元的直径与比例常数Εt之积来计算;e t——两单元边界距离D t和张力弹性常数k t产生的拉伸弹性力;d n——法向非弹性力。

因为这个条件,弹性力在压缩和拉伸方向都是有效的。

图2　单元之间的张力F ig.2　T en si on betw een tw o elem en ts另外,在传统的离散元模型中,切向方向的滑移摩擦限制了切向弹性力e s的范围,式(5)表明了切向方向的断裂标准。

在土壤离散元模型中,考虑到库仑定律,用式(9)代替式(5):当 e s >f c+e n tgΥ时f s=(f c+e n tgΥ) sign[e s](9)式中　Υ——土壤的内摩擦角;f c——土壤粘附系数c 和土壤内摩擦角Υ确定的力。

还有些学者为了考虑实际土壤颗粒之间的粘附性,使用系数C ad修正了传统的离散元模型[22],并且指出当式(10)成立时,两个单元之间就会产生法向张力。

r i+r j<D≤(1+C ad) (r i+r j)(10)式中　r i——单元i的半径;r j——单元j的半径; D——两单元中心之间的距离;C ad——系数。

71　第1期张　锐等:离散单元法在土壤机械特性动态仿真中的应用进展当两单元开始分离时,利用修正的离散元模型,使用式(11)可以计算法向张力f tens。

f tens=k ad (r i+r j-D)(11)式中　f tens——两单元之间的法向张力;k ad——单元之间粘附的弹性常数。

离散单元法在土壤动态行为仿真中的应用在国内尚鲜见研究报导,国外主要是日本在这方面已经有一定的发展,他们大多利用这种方法仿真机具作用下土壤的动态行为,据此设计并制造出更加适用的土壤机械。

这些学者的研究主要以二维为主,并且离散单元参数主要靠试验和经验公式获得[23,24]。

总体来看,该项技术尚处于初期阶段,无论在理论基础研究还是在技术应用领域仍很不成熟,由于如上所述的优越性,基于离散单元法仿真土壤动态行为存在很大的发展空间和良好的应用前景。

土壤本身是复杂的多相混合体,在建立土壤离散单元仿真模型时应该考虑土壤多方面的性质,即除了考虑土壤的粘附性以外,还应该考虑土壤的塑性特性,以及土壤内部水分、空气、粒径分布以及密度等对土壤动态行为的影响。