第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号 吉林大学南岭校区生物与农业工程学院,130025
1 引 言
离散单元法(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]:
6
1
∃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]。总体来看，该项技术尚处于初期阶段，无论在理论基础研究还是在技术应用领域仍很不成熟，由于如上所述的优越性，基于离散单元法仿真土壤动态行为存在很大的发展空间和良好的应用前景。
土壤本身是复杂的多相混合体，在建立土壤离散单元仿真模型时应该考虑土壤多方面的性质，即除了考虑土壤的粘附性以外，还应该考虑土壤的塑性特性，以及土壤内部水分、空气、粒径分布以及密度等对土壤动态行为的影响。离散单元法中的参数除了通过经验公式和试验的方法获得外，还可以通过理论推导或采用半经验公式的方法进行计算[25]，通过这种方法获得的参数不仅可以更好地满足计算的要求，而且要比试验方法省时省力。现实世界是三维的，二维仿真无法真实地模拟物体之间相互作用的动态行为，为了更加形象地模拟一个过程或者系统的行为，应用离散单元法进行三维仿真的研究是可行的，而且十分必要。另外，通过仿真模型机具与土壤的动态行为，从而改进土壤机械的设计或制造，这种方法在实际应用中存在一定的误差，所以应用三维离散单元法研究实际机具与土壤之间相互关系的三维动态仿真研究有重要的现实意义。
4 结 语
离散单元法在国内外的发展很快，主要应用在岩石力学、采矿工程及松散物料流动等方面的研究中，
相应的算法和软硬件已基本成熟，将粘附力的理论引入传统离散单元法可以实现其在地面机器系统中土壤动态行为的模拟。但离散单元法在土壤动态行为仿真中的应用还处于发展阶段，只有少数几个国家(主要是日本)的一些学者在这方面进行了一些探索性研究。由于土壤本身的性质十分复杂，土壤力学模型很难建立，仿真参数也不易确定，所以目前这方面的研究主要是在二维条件下，并且用离散单元法进行试验机具作用下的土壤动态行为仿真，同时仿真参数也主要靠经验来确定。国内将离散单元法应用于分析机械和土壤之间动态相互关系的研究很少，目前尚鲜见有相关文章发表。因此，离散单元法在土壤动态行为仿真中的应用有待进一步研究。
离散单元法也有其局限性：离散单元法是为分析离散颗粒的集合而建立的，它的计算需要花费大量的时间[8]；离散单元法的参数对于计算结果的影响十分复杂，目前主要通过经验公式和试验的方法确定参数，所以常存在很大的误差[13]；对单元间的力学关系进行的假设和近似与实际有差距，累加后误差较大，目前只能靠实际试验进行修正。离散单元法将所分析的物体看作是离散颗粒的集合体，这正适合于分析象土壤这样的物质在触土部件作用下的动态行为。离散单元法克服了有限元和边界元计算时需要复杂边界条件的缺点，使仿真某些复杂行为简单化。另外，利用离散单元法可以代替某些费时费力的试验，是一种研究机械作用下的土壤动态行为的有效方法和工具。由于土壤本身的复杂性，通过研究合理地确定离散单元参数以及建立合适的离散单元仿真力学模型，对于正确地模拟地面机器系统中土壤的动态行为十分必要。
[参 考 文 献]
[1] Cundall P A.A compu ter model fo r si m u lating p rogressive
large scale movem en ts in b locky system[A].In:M u ller L ed.P roceedings of Sympo sium of the In ternati onal Society of Rock M echan ics[C].Ro tterdam:A.A.Balkem a,1971
(1):8～12.
[2] Cundall P A,Strack O D L.A discrete num erical m ethod
fo r granu lar asm b lies[J].Geo techn ique,1979,29(1):47～65.
[3] 鲁 军，张楚汉，王光伦等.岩体动静力稳定分析的三维离
散单元法[J].清华大学学报,1996,36(10):98～104.
[4] 张楚汉.结构——地基动力相互作用问题[A].结构与介质
相互作用理论及其应用[M].南京：河海大学出版社,1993. [5] 徐 泳,Kafu i K D,T ho rn ton C.用颗粒离散元法模拟料仓
卸料过程[J].农业工程学报,1999,15(3):65～69.
[6] H iroak i T anaka,Ko ji Inooku,Yu ji N agasak i,et al.
Si m u lati on of so il loo n ing at sub su rface tillage u sing a vib rating type sub so iler by m ean s of the distinct elem en t m ethod[A].P roceedings of8th Eu ropean Conference of the In ternati onal Society of T errain2V eh icle System s[C],2000:
32～38.
[7] T anaka H,M omo tsu M,O ida A,et al.Con structi on of
model of m echan ical so il behavi o r by m ean s of distinct elem en t m ethod[J].J Kan sai B ranch of JSAM79,1996. [8] Fu jii H,O ida A,N akash i m a H,et al.A nalysis of lunar
terrain2w heel system in teracti on by D E M[A].P roceedings of the6th A sia2Pacific Conference of the In ternati onal Society of T errain2V eh icle System s[C],2001:129～136. [9] ITA SCA Con su lting Group.3D EC—32D D iscrete E lem en t
Code.1987.
[10] Shen B,Stephan sson O.Cyclic loading characteristics of
j o in ts and rock b ridges in a j o in ted rock speci m en[A].
P roceedings of the In ternati onal Conference on Rock Jo in ts
[C],1990:725～729.
[11] 王泳嘉.离散单元法——一种适用于节离岩石力学分析的
数值方法[A].第一届全国岩石力学数值计算及模型试验
讨论论文集[C],1986:32～37.
[12] H iroak i T anaka,Ko ji Inooku,O sam u Sum ikaw a,et al.
Si m u lati on of so il behavi o r at sub so iling by the distinct
elem en t m ethod[A].P roceedings of the6th A sia2Pacific
Conference of the In ternati onal Society of T errain2V eh icle
81农业工程学报2003年
System s [C ],2001:194
～200.[13] M asato sh i M omozu ,A k ira O ida ,H ito sh i N akash i m a .
Si m u lati on of shear box tex t by the distinct elem en t m ethod [A ].
P roceedings of
the
6th A sia 2Pacific
Conference of the In ternati onal Society of T errain 2V eh icle
System s [C ],2001:181
～188.[14] H iro sh i N akash i m a ,A k ira O ida .
Si m u lati on of so il 2tire in teracti on by a coup led distinct elem en t 2fin ite elem en t m ethod [A ].
P roceedings of
the
6th A sia 2Pacific
Conference of the In ternati onal Society of T errain 2V eh icle
System s [C ],2001:59
～63.[15] T aku Ibuk i ,A k ira O ida .
Si m u lati on to analyze the
in teracti on betw een so il and a tire lug by distinct elem en t m ethod [A ].P roceedings of the 8th Eu ropean Conference of the In ternati onal Society of T errain 2V eh icle System s [C ],2000:87～94.
[16] 陈龙斌，胡晓军，唐志平.离散元数值模拟中查找邻居元关
系的改进算法[J ].计算力学学报,2000,17(4):497～499.
[17] M asato sh iM omo tsu ,H iroak i T anaka ,A k ira O ida ,et al
.Si m u lati on of so il defo rm ati on and stress variati on at th ree ax ial comp ressi on by the distinct elem en t m ethod [A ].P roceedings of the 12th In ternati onal Conference of the In ternati onal Society of T errain 2V eh icle System s [C ],1996:21～28.
[18] H iroak i T anaka ,M asato sh i M omozu ,A k ira O ida ,et al .
Si m u lati on of so il defo rm ati on and resistance at bar penetrati on by the distinct elem en t m ethod [J ].Jou rnal of
T erram echan ics ,2000,37:41
～56.[19] B rady N C ,W eil R R .
E lem en ts of the natu re and
p roperties of so ils [M ].N ew Jery ,P ren tice 2H all ,Inc .2000.
[20] H iromm a T ,O h ta Y ,Kataoka T .A nalysis of the so il
defo rm ati on beneath a w heel by fin ite elem en t m ethod [J ].J Japane Society of A gricu ltu ralM ach inery ,1994,56(6):3～10.
[21] H iroak i T anaka ,Ko ji Inooku ,et al .N um erical analysis of
so il loo n ing in sub su rface tillage by a vib rating type sub so iler by m ean s of the distinct elem en t m ethod [A ].P roceedings of the 13th In ternati onal Conference of the In ternati onal Society of T errain 2V eh icle System s [C ],
1999:791～798.
[22] M omozu M ,O ida A ,Yam azak i M ,et al .Si m u lati on of
so il loo n ing p rocess by pendu lum type b lade by m ean s of modified distinct elem en t m ethod [A ].P roceedings of 13th In ternati onal Conference of the In ternati onal Society of
T errain 2V eh icle System s [C ],1999:71
～78.[23] O ida A ,H i Schw anghart ,O hkubo S ,et al .Si m u lati on of
so il defo rm ati on under a track shoe by the D E M [A ].P roceedings of the 7th Eu ropean Conference of the In ternati onal Society of T errain 2V eh icle System s [C ],1997:155～162.
[24] O hkubo S ,O ida A ,Yam azak iM .A pp licati on of D E M to
find the in teracti on betw een tire and so il [A ].P roceedings of the 5th A sia 2Pacific Conference of the In ternati onal
Society of T errain 2V eh icle System s [C ],1998:175
～183.[25] 邢纪波，余良群，张瑞丰等.离散单元法的计算参数和求解
方法选择[J ].计算力学学报,1999,16(1):47～51.
D evelop m en t of si m ulation on m echan ical dynam ic behav ior
of so il by d isti nct elem en t m ethod
Zha ng Rui ,L i J ia nq ia o ,L i Yinw u
(T he K ey L abora tory f or T erra in 2M ach ine B ion ics E ng ineering ,
M in istry of E d uca tion ,J ilin U n iversity ,Chang chun 130025,Ch ina )
Abstract :T he developm en t and basic theo ry of distinct elem en t m ethod (D E M )w ere in troduced in b rief .T he details of app licati on on si m u lating the m echan ical dynam ic behavi o r of
so il by D E M in the rearch field of terram echan ics w ere p ren ted .T he feasib ility of the app licati on of D E M to describe so il dynam ic p rop erty w as analyzed .T he sp ecial p ri o rity of describ ing so il dynam ic behavi o r u sing the w ho le theo ry of distinct elem en ts in app lying D E M to the discon tinuou s m u lti p ha m u lti p le system such as so il w as po in ted ou t .T he trends of the developm en t of D E M in the so il dynam ic p rop erty and the key p rob lem s that m ay ari from its u in the near fu tu re w ere p u t fo rw ard .
Key words :distinct elem en t m ethod (D E M );so il dynam ic si m u lati on ;m echan ical m odel
9
1 第1期
张 锐等：离散单元法在土壤机械特性动态仿真中的应用进展
