第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]，通過這種方法獲得的參數不僅可以更好地滿足計算的要求，而且要比試驗方法省時省力?，F實世界是三維的，二維仿真無法真實地模擬物體之間相互作用的動態行為，為了更加形象地模擬一個過程或者系統的行為，應用離散單元法進行三維仿真的研究是可行的，而且十分必要。另外，通過仿真模型機具與土壤的動態行為，從而改進土壤機械的設計或制造，這種方法在實際應用中存在一定的誤差，所以應用三維離散單元法研究實際機具與土壤之間相互關系的三維動態仿真研究有重要的現實意義。
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期
張 銳等：離散單元法在土壤機械特性動態仿真中的應用進展
