基于離散元的西北旱區(qū)農(nóng)田土壤顆粒接觸模型和參數(shù)標定

石林榕，趙武云，孫 偉

?

基于離散元的西北旱區(qū)農(nóng)田土壤顆粒接觸模型和參數(shù)標定

石林榕，趙武云※，孫 偉

（甘肅農(nóng)業(yè)大學機電工程學院，蘭州 730070）

為了解決利用離散元法模擬土壤作業(yè)過程在預測農(nóng)具阻力和土壤動態(tài)運動時存在失真等問題，整合延遲彈性模型（hysteretic spring contact model，HSCM）和線性內(nèi)聚力模型（liner cohesion model，LCM）優(yōu)勢建立西北旱區(qū)農(nóng)田土壤模型，以不同參數(shù)（靜摩擦系數(shù)、動摩擦系數(shù)和內(nèi)聚強度）組合下仿真得到的土壤仿真堆積角為響應值，基于Box-Behnken試驗法建立回歸模型，并根據(jù)該回歸模型進行了參數(shù)預測并驗證，對17組土壤仿真堆積角方差分析表明：靜摩擦系數(shù)、動摩擦系數(shù)、動摩擦系數(shù)和抗剪強度的交互項、動摩擦系數(shù)的二次項對仿真堆積角的影響極顯著；靜摩擦系數(shù)和動摩擦系數(shù)的交互項、靜摩擦系數(shù)的二次項對仿真堆積角的影響顯著。使用預測的參數(shù)進行6種不同含水率土壤直接剪切仿真和試驗對比可知，當含水率為1%～20%時，仿真與試驗間的抗剪強度相對誤差為1.18%～9.31%，仿真與試驗間的內(nèi)摩擦角相對誤差為0.55%～4.07%。對仿真和試驗鴨嘴插入阻力數(shù)據(jù)進行分析可知，仿真與試驗曲線在入土距離處于0～50 mm期間時，但仿真入土阻力曲線波動較大，仿真和試驗阻力走勢基本一致，玉米直插穴播最深處50 mm處的仿真和試驗入土阻力相對誤差為0.928%，可利用此時的入土阻力分析直插鴨嘴結構對強度的影響。

土壤；離散元法；摩擦系數(shù)；參數(shù)；標定

0 引 言

在農(nóng)業(yè)精耕細作、精密播種、高效收獲等作業(yè)過程及減阻降耗研究中，土壤力學是設計和優(yōu)化觸土部件優(yōu)先研究的問題[1-2]。多種方法已經(jīng)用于研究農(nóng)具－土壤互相作用過程的研究。傳統(tǒng)經(jīng)驗方法在分析土壤運動行為時只涉及土壤應力失效問題[3]。Asaf等[4]提出的基于Terzaghi’s被動土壤壓力理論的數(shù)值模擬方法--離散元法（discrete element method，DEM）和流體動力學（computational fluid dynamics，CFD）不僅可以計算模型作業(yè)阻力，還能模擬土壤破碎過程，但FEM和CFD只能對連續(xù)介質(zhì)的土壤進行整體分析[5]。農(nóng)田作業(yè)過程中存在土層的分離、混合和裂縫及土壤顆粒的流動，呈現(xiàn)出散粒體的特性，DEM可為研究土壤復雜的動態(tài)行為提供了新的思路和方法。

內(nèi)摩擦力、內(nèi)聚力和黏附力是影響田間作業(yè)阻力的關鍵土壤參數(shù)。內(nèi)摩擦力是土顆?；瑒犹匦灾弧?nèi)聚力為相鄰顆粒間粘結的分子力，由土顆粒間存在的液橋引起。黏附力為不同材料互相吸附的分子力?，F(xiàn)有基于DEM研究土壤文獻大體分為3類[6]：1）土顆粒間的液橋形成的內(nèi)聚力和黏附力，介于其方程的復雜性及顆粒的不規(guī)則，目前還未應用于土壤-農(nóng)具互作仿真模型；2）基于宏觀物理力學定律建立的涉及內(nèi)聚力和黏附力土壤模型，通過匹配試驗結果進行參數(shù)校準，其中并行黏結模型（parallel bond contact model，PBCM）多被用于建立具有內(nèi)聚力的土壤-農(nóng)具仿真模型[7]，PBCM模型中顆粒間由粘結鍵連接來限制顆粒法向和切向力，當法向和切向力大于預先設定的力值，粘結鍵會被刪除，其后顆粒間通過Hertz-Mindlin模型（Hertz Mindlin contact model）進行計算。石林榕等[8]借助Flunet-UDF將土壤單軸壓縮有限元模型網(wǎng)格信息導出，利用PBCM生成含水率為10.5%的不飽和土壤單軸壓縮模型，仿真結果顯示模型軸向壓縮力-應變曲線與試驗曲線基本一致。然而HMCM模型在預測土壤-農(nóng)具互作模型時會產(chǎn)生不真實的垂直阻力[9-10]；3）內(nèi)聚力和黏結力被添加到法向和切向力中[11]。然而，以上3類模型中均未涉及到土壤的塑性行為。Walton等[12]針對塑性材料創(chuàng)建延遲彈性模型（hysteretic spring contact model，HSCM）。當外力未達到預先設定的應力強度時，顆粒間執(zhí)行線彈性方程；當顆粒接觸面的總應力超過設定的應力時，顆粒間執(zhí)行塑性方程，即顆粒間總應力超過屈服點后模型有塑性變形行為。Ucgul等[13-14]將線彈性內(nèi)聚力和黏結力整合進HSCM模型中，通過直剪試驗得到不同含水率的剪切強度，土壤貫入試驗得到HSCM屈服點，仿真參數(shù)靜摩擦系數(shù)和動摩擦系數(shù)通過反復調(diào)整匹配土壤試驗結果，仿真結果顯示改進的HSCM能準確預測不同含水率、作業(yè)深度、垂直壓力、牽引速度深松犁的水平和垂直阻力，為基于DEM優(yōu)化農(nóng)具提供新的思路。

基于以上研究，為了解決利用DEM模擬土壤作業(yè)過程在預測農(nóng)具作業(yè)阻力和土壤動態(tài)運動時存在失真等問題，考慮顆粒間水對土顆粒流動特性的影響，將線性內(nèi)聚力模型（liner cohesion model，LCM）整合到HSCM模型中建立6種不同含水率西北旱區(qū)農(nóng)田土壤模型，并對該模型所需參數(shù)進行標定，以期為西北旱區(qū)農(nóng)田土壤離散元仿真參數(shù)標定提供參考。

1 離散元土壤模型建立

1.1 接觸模型整合

試驗土壤為沙壤土，既有散粒體物料特性又有一定的壓縮性。采用的土壤接觸模型為HSCM，并在其法向增加LCM，如圖1a所示。土壤的彈塑性如圖1b所示。

注：ra, rb分別為顆粒a和b的半徑，m；,分別為法向接觸力和阻尼力，N；,分別為切向接觸力和阻尼力，N；μ為摩擦系數(shù)；Uabn為顆粒間的法向重疊量，m；U0為顆粒間的殘余重疊量，m；K1, K2分別為加載和卸載剛度，N×m-1。下同。

在HSCM中法向接觸力定義為

式中r等效半徑，m；Y和Y為分別為顆粒和的屈服強度，Pa；屈服強度的默認值（Pa）可以用楊氏模量和顆粒半徑預估。

式中若＜100 Pa，可用0.003代替[12]。

卸載剛度為

式中為顆粒間的恢復系數(shù)。

每個時間步殘余重疊量由式（5）更新。

切向接觸力定義為

法向和切向阻尼力分別為

將內(nèi)聚力施加到法向力中。

式中F和F分別為顆粒與顆粒的內(nèi)聚力、顆粒與接觸物間的黏結力，N。

式中r顆粒接觸半徑，m。

式中E為等效楊氏模量，Pa。

式中E和E分別為顆粒和的楊氏模量，v,v分別為顆粒和的泊松比。

1.2 仿真模型參數(shù)確定

根據(jù)HSCM及LCM確定土壤基本物理參數(shù)、接觸力學參數(shù)和接觸模型參數(shù)?；疚锢韰?shù)包括顆粒粒徑分布、土粒密度、剪切模量、泊松比、容重。由于水的作用，土壤顆粒群內(nèi)形成的液橋使土壤顆粒之間產(chǎn)生內(nèi)聚力，本文通過土壤剪切試驗獲取不同含水率土壤的抗剪強度和內(nèi)摩擦角，LCM模型參數(shù)能量密度分別等于6種不同含水率條件下的土壤內(nèi)聚強度，HSCM模型參數(shù)屈服強度通過土壤貫入試驗獲得，顆粒間的阻尼系數(shù)和剛度系數(shù)采用EDEM軟件默認值，分別為0.05、0.95[15-17]。接觸力學參數(shù)包括恢復系數(shù)、靜摩擦系數(shù)和動摩擦系數(shù)。

1.2.1 土壤基本物理參數(shù)

通過標準土壤篩篩分試驗可知80%的顆粒尺寸介于0.25和5 mm之間。由于在大多數(shù)土壤中，礦物質(zhì)顆粒約占土壤固相物質(zhì)質(zhì)量的95%以上，離散元軟件模擬土壤顆粒需要用到土粒密度，取土粒密度值為2 600 kg/m3[18-19]。土壤彈性模量試驗測量困難，通過查閱文獻[20-21]來獲得土壤的彈性模量。

土壤剪切模量表征著土壤抵抗切應變的能力[21]，可由彈性模量和泊松比近似確定，如式（15）所示。

式中為彈性模量，Pa；為泊松比。

土壤的泊松比是指土壤在無側限條件受壓時，側膨脹應變與豎向壓縮應變的比值，通常按公式估算泊松比的值，根據(jù)材料力學中廣義胡克定律可推導求得土壤側壓力系數(shù)與泊松比的相互關系[22]，同時查閱相關文獻和資料，泊松比取值范圍為0.3～0.42[23-25]。泊松比通過式（16）計算。

式中0為側壓力系數(shù)。

土壤側壓力系數(shù)是指土壤在側限條件下受壓時，側向壓力與豎向壓力的比值，根據(jù)經(jīng)驗公式由土壤的內(nèi)摩擦角由式（17）近似求得[23]。

式中為內(nèi)摩擦角，（°）。

為制備不同含水率土壤進行堆積角試驗和直接剪切試驗，通過式（18）制備。

式中1為制備試樣所需要的加水量，kg；0為基土質(zhì)量，kg；0為基土含水率，%；1為制樣要求的含水率，%。

將土壤存放于塑料袋中24 h，使其水分均勻，用環(huán)刀切取土壤并稱質(zhì)量，土壤容重可由式1=/計算得到，通過式（19）可得不同含水率土壤的孔隙率。

式中e為土壤孔隙率，%；1為土壤容重，kg/m3；為土顆粒密度，kg/m3。

內(nèi)摩擦角、抗剪強度等力學參數(shù)對仿真結果準確性起著非常重要的作用。本文采用快剪試驗測定6種含水率在3種垂直壓力（70.22、151、256 kPa）下的的土壤的抗剪強度，根據(jù)土壤的剪切強度與垂直壓力關系獲取土壤的抗剪強度、內(nèi)摩擦角等力學參數(shù)[26-29]。試驗儀器：南京土壤儀器廠/ZJ型應變控制直剪儀，1000N-S型拉壓力傳感器，精度0.001 g的電子天平。將剪切速率控制在0.8 mm/min左右，3～5 min內(nèi)剪破。通過計算獲得不同含水率土壤的體積密度、孔隙率、泊松比和剪切模量，如表1所示。

表1 6種不同含水率土壤試驗數(shù)據(jù)

1.2.2 接觸力學參數(shù)

仿真顆粒的大小受限于計算能力，仿真時間是隨著生成的顆粒數(shù)量的增加而變長，建立與實際土壤顆粒大小的顆粒模型較難。試驗發(fā)現(xiàn)80%的顆粒尺寸介于0.25和5 mm之間。通過查閱資料可知，肥沃土壤的結構形態(tài)由土壤黏結成粒狀和小團塊狀，大體呈球形?；谝陨显颍瑫r為大幅縮短仿真時間，本文采用半徑為5 mm球體作為土壤基體，選擇球徑倍率分布為0.95～1.05。為使仿真結果接近試驗結果，需對半徑5 mm仿真顆粒的接觸參數(shù)進行校核。

在進行參數(shù)校核之前需測定6種含水率土壤的堆積角。采用無底的有機玻璃圓柱筒（內(nèi)徑54 mm′高度270 mm）裝填土壤至一定高度后，緩慢抬起圓柱筒，待土壤坡面穩(wěn)定后，垂直拍照，導入CAD軟件中標注土壤堆積角。不同含水率土壤堆積角如圖2所示。由圖2可知，隨土壤含水率由1%±1%增大至20%±1%，土壤堆積角由34.14°降至21.2°。

圖2 試驗中6種含水率土壤堆積角

通過仿真采用整合LCM和HSCM模型的土壤顆粒的堆積試驗，對堆積后形成的傾角與測得的內(nèi)摩擦角進行對比。顆粒半徑大小為5 mm，顆粒生成質(zhì)量通過前文的土壤容重進行計算，時間步長為1×10-5s，在一個圓柱體內(nèi)生成顆粒，完畢后撤掉圓柱體，讓顆粒自然堆積。如圖3為不同靜摩擦系數(shù)、動摩擦系數(shù)和抗剪強度條件下的仿真土壤堆積角。由圖3可知，隨靜摩擦系數(shù)由0.23增加至0.53，動摩擦系數(shù)和抗剪強度在相等條件下仿真土壤堆積角增大；隨動摩擦系數(shù)由0.16增加至0.30，靜摩擦系數(shù)和抗剪強度在相等條件下仿真土壤堆積角增大；隨抗剪強度由3.47 kPa增加至18.58 kPa，靜摩擦系數(shù)和動摩擦系數(shù)在相等條件下仿真土壤堆積角增大。

注：A0.23 B0.23 C3.47代表靜摩擦系數(shù)0.23-動摩擦系數(shù)0.23-抗剪強度3.47 kPa。

為后續(xù)研究方便，依據(jù)中心試驗設計方法，建立土壤仿真模型接觸參數(shù)預測模型。采用Box-Behnken設計模型能用較少的試驗次數(shù)進行全面的分析研究。以靜摩擦系數(shù)、動摩擦系數(shù)和抗剪強度為仿真試驗因素，仿真土壤堆積角為響應值，對各試驗因素水平進行編碼，依響應面法試驗安排共實施17組仿真試驗見表2。

表2 仿真試驗方案及結果

應用Design-Expert.8.06軟件對仿真試驗結果進行二次回歸分析，結果表明，靜摩擦系數(shù)、動摩擦系數(shù)、動摩擦系數(shù)和抗剪強度的交互項、動摩擦系數(shù)的二次項對仿真堆積角的影響極顯著（＜0.05）；靜摩擦系數(shù)和動摩擦系數(shù)的交互項、靜摩擦系數(shù)的二次項對仿真堆積角的影響顯著（0. 05＜＜0.1）；抗剪強度、靜摩擦系數(shù)和抗剪強度的交互項和抗剪強度的二次項對仿真堆積角影響不顯著（＞0.1）。剔除不顯著項，最終獲得仿真堆積角回歸方程為

=-30.84-146.67+478.02+324.29-9.62+164.612-849.742（20）

式中為靜摩擦系數(shù)的水平值；為動摩擦系數(shù)的水平值；為抗剪強度的水平值；為堆積角，(°)。

以實際土壤堆積角（圖2）為優(yōu)化目標，將快剪試驗得到的抗剪強度作為已知值，根據(jù)回歸模型（1）預測相應的靜摩擦系數(shù)和動摩擦系數(shù)，結果如表3所示。

表3 預測的土壤靜摩擦系數(shù)和動摩擦系數(shù)

1.2.3 參數(shù)校核

通過直剪仿真試驗對預測的土壤的靜摩擦系數(shù)和動摩擦系數(shù)進行校核。直剪仿真試驗如圖4所示。根據(jù)直剪試驗原理，建立剪切仿真試驗裝置，分上盒、下盒和壓力板，上盒、下盒外徑為62 mm，高10 mm，下盒上口開、下口閉，上盒上下口均開。通過對壓力板添加較小垂直向下速度來添加土壤模型垂直載荷。然后給下剪切盒一定的水平速度其大小為1 mm/s，下盒水平位移量達到6 mm時停止仿真，在此期間取峰值作為剪切強度，無峰值取4 mm處值。

圖4 直剪仿真試驗

將6種含水率條件下土壤快速直剪仿真與試驗得到的內(nèi)摩擦角和抗剪強度進行對比。由表4可知，當含水率為1%～20%時，仿真與試驗間的抗剪強度相對誤差為1.18%～9.31%，仿真與試驗間的內(nèi)摩擦角相對誤差為0.55%～4.07%，可以認為土壤快速直剪仿真結果與試驗結果近似一致。

表4 6種含水率的仿真土壤抗剪強度與摩擦角

為了獲取含水率為1%的土壤屈服強度，通過自制的土壤貫入試驗裝置獲取壤貫入試驗壓力-沉陷曲線，土壤貫入工具如圖5a所示。試驗中土壤貫入試驗工具以 1 mm/s速度貫入土壤。得到的土壤貫入試驗壓力-沉陷曲線如圖5b所示。如圖所示，在軸方向0～12 mm區(qū)間，由于傳感器精度問題，未有應力值；隨著工具繼續(xù)深入，當位移量處于軸方向12～26 mm區(qū)間時應力與位移量基本趨于直線關系；當位移量處于軸方向26～40 mm區(qū)間時，應力與位移量呈下彎走勢；當位移量繼續(xù)增加，土壤發(fā)生失效，應力極速下降，此時的土壤受到的應力值為0.23 MPa，最后趨于平緩。因此，HSCM接觸模型的屈服強度參數(shù)為0.23 MPa。

最終獲得的土壤離散元法仿真模型所需參數(shù)中，土壤含水率、泊松比、剪切模量見表1，土壤靜摩擦和動摩擦系數(shù)見表3，其他參數(shù)值見表5。

圖5 土壤貫入試驗工具及壓力-沉陷曲線

表5 離散元法仿真的參數(shù)

2 驗證試驗

2.1 鴨嘴插入土壤模型仿真

直插穴播工作過程中鴨嘴與土壤的接觸、垂直受力是鴨嘴設計重點。本文采用的直插鴨嘴驗證模型為韓國產(chǎn)提攜直插式玉米播種器，如圖6a所示。鴨嘴模型通過Solidworks軟件建模，鴨嘴模型具體尺寸如圖6b所示。將鴨嘴模型導入EDEM 中并設置材料參數(shù)與運動屬性參數(shù)，使用土壤含水率為1%對應土壤參數(shù)建立土壤模型，土壤顆粒半徑大小為5 mm，球徑倍率分布為0.95～1.05，仿真土槽模型外徑100 mm、高100 mm，未仿真之前調(diào)整鴨嘴尖距土槽上面2 mm，仿真時鴨嘴垂直插入土槽速率為1 mm/s，時間步4.17′10-6s，仿真總時間為50 s，使用壓板壓實土壤模型使其孔隙率為49.23%，壓板速率為1 mm/s。試驗前，調(diào)整鴨嘴處于圓形土槽中間位置，并使鴨嘴尖離2 mm，設置鴨嘴垂直插入土壤速率為 1 mm/s，當鴨嘴尖插入土壤5 s時停止運動，將傳感器鴨嘴入土阻力導入電腦進行處理。試驗裝置為HDV-1K電動雙柱拉壓力測試臺。其他條件與試驗一致。仿真時測量鴨嘴受到土壤模型的垂直阻力，當鴨嘴進入土壤模型50 mm處停止鴨嘴運動。建立的鴨嘴-土壤模型如圖6c所示。

2.2 結果與分析

通過比較仿真與試驗中鴨嘴插入土壤過程入土阻力來驗證所標定土壤參數(shù)的可靠性。由于優(yōu)化鴨嘴結構尺寸一般參考鴨嘴入土阻力的最大值，通過比較不同結構尺寸鴨嘴的入土阻力來確定最優(yōu)結構，因此，本文重點比較在同等條件下仿真和試驗垂直入土阻力的誤差，兼參考曲線走勢進行仿真模型和參數(shù)驗證。鴨嘴插入土壤試驗裝置如圖7a所示。

1. 固定鴨嘴 2. 活動鴨嘴 3. 固定座 4. 排種輪 5. 種管 6. 把手 7. 支撐座 8. 彎頂桿 9. 彈簧 10. 導桿

注：Re為相對誤差。

對仿真和試驗鴨嘴插入阻力數(shù)據(jù)進行分析（圖7b）可知，仿真與試驗曲線在入土距離處于0～50 mm期間時，仿真和試驗阻力走勢相差不大，但仿真入土阻力曲線波動較大，分析原因可知，由于采用半徑為5 mm的土顆粒，與實際土顆粒相差較大，雖然能大幅度減少仿真時間，但不可避免會產(chǎn)生入土阻力波動較大的問題。玉米直插穴播最深處50 mm處的仿真和試驗入土阻力相對誤差為0.928%，此時的仿真入土阻力最大，可利用最大阻力值分析直插鴨嘴結構對其強度的影響。

3 結 論

1）提出將HSCM和LCM整合于一體模擬西北農(nóng)田土壤，利用Box-Behnken試驗法建立以靜摩擦系數(shù)、動摩擦系數(shù)和抗剪強度為因素、以土壤仿真堆積角為指標的數(shù)學回歸模型，并借助土壤直接剪切仿真試驗對靜摩擦和動摩擦系數(shù)參數(shù)進行校核，校核結果發(fā)現(xiàn)：6種不同含水率土壤直接剪切仿真和試驗中，當含水率為1%～20%時，仿真與試驗間的抗剪強度相對誤差為1.18%～9.31%，仿真與試驗間的內(nèi)摩擦角相對誤差為0.55%～4.07%，表明HSCM和 LCM組合模型是可行的。

2）通過鴨嘴插穴土壤仿真與實際對比檢驗2種土壤接觸模型的適應性以及預測的靜摩擦系數(shù)和動摩擦系數(shù)的可靠性，結果表明：仿真與試驗曲線在入土距離處于0～50 mm期間時，但仿真入土阻力曲線波動較大，仿真和試驗阻力走勢基本一致，玉米直插穴播最深處50 mm處的仿真和試驗入土阻力相對誤差為0.928%，此時的仿真入土阻力最大，可利用最大阻力值分析直插鴨嘴結構對其強度的影響。

[1] Zadeh S R A. Modelling of Energy Requirements by a Narrow Tillage Tool[D]. Saskatoon: University of Saskatchewan, 2006.

[2] 趙成剛. 土力學的現(xiàn)狀及其數(shù)值分析方法中某些問題的討論[J]. 巖土力學，2006，27(8)：1361－1364.Zhao Chenggang. Discussion on state-of-art of soil mechanics and some problems with applications of numerical method to geotechnical engineering[J]. Rock and Soil Mechanics, 2006, 27(8): 1361－1364. (in Chinese with English abstract)

[3] Godwin R J, O’Dogherty M J, Saunders C, et al. A force prediction model for mould board ploughs incorporating the effects of soil characteristic properties, plough geometric factors and ploughing speed[J]. Biosystems Engineering, 2007, 97(1): 117－129.

[4] Asaf Z, Rubinstein D, Shmulevich I. Evaluation of link-track performances using DEM[J]. Journal of Terramechanics, 2006, 43(2): 141－161.

[5] Karmakar S, Ashrafizadeh S R, Kushwaha R L. Experimental validation of computational fluid dynamics modeling for narrow tillage tool draft[J]. Journal of Terramechanics, 2009, 46(6): 277－283.

[6] Tsuji T, Nakagawa, Matsumoto N, et al. 3D DEM simulation of cohesive soil pushing behavior by bulldozer blade[J]. Journal of Terramechanics, 2011, 49(1): 37－47.

[7] 張銳，李建橋，李因武，等. 部件復雜表面影響土壤擾動行為的離散元宏細觀分析[J]. 吉林大學學報：工學版，2009，39(5)：1218－1223. Zhang Rui, Li Jianqiao, Li Yinwu, et al. DEM macroscopic and mesoscopic analysis in disturbed behavior of soil acted by part with complex surface[J]. Journal of Jilin University: Engineering and Technology Edition, 2009, 39(5): 1218－1223. (in Chinese with English abstract)

[8] 石林榕，吳建民，趙武云，等. 基于離散單元法農(nóng)田土壤單軸壓縮模型的建立及參數(shù)校核[J]. 中國農(nóng)業(yè)大學學報，2015，20(4)：174－182. Shi Linrong, Wu Jianmin, Zhao Wuyun, et al. Establishment and parameter verification of farmland soil model in uniaxial compression based on discrete element method[J]. Journal of China Agricultural University, 2015, 20(4): 174－182. (in Chinese with English abstract)

[9] Ucgul M, Fielke J M, Saunders C. Three-dimensional discrete element modelling of tillage: Determination of a suitable contact model and parameters for a cohesionless soil[J]. Biosyst Eng, 2014, 121(2): 105－117.

[10] Ucgul M, Fielke J M, Saunders C. 3D DEM tillage simulation: validation for a sweep tool for a cohesionless soil[J]. Soil Till Res, 2014, 144(4): 220－227.

[11] Momozu M, Oida A, Yamazaki M, et al. Simulation of a soil loosening process by means of the modified distinct element method[J]. Journal of Terramechanics, 2002, 39(4): 207－220.

[12] Walton O R, Braun R L. Stress calculations for assemblies of inelastic spheres in uniform shear[J]. Acta Mechanica, 1986, 63(1/2/3/4): 73－86.

[13] Ucgul M, Fielke J M, Saunders C. Three-dimensional discrete element modelling (DEM) of tillage: accounting for soil cohesion and adhesion. Biosyst Eng, 2015, 129(5): 298－306.

[14] Ucgul M, Fielke J M, Saunders C. Three-dimensional discrete element modelling of tillage: determination of a suitable contact model and parameters for a cohesionless soil[J]. Biosystems Engineering, 2014, 121(2): 105－117.

[15] DEM Solution Ltd. EDEM 2.5 theory reference guide[R/OL]. 2014-12-05[2015-06-07].http://www.docin.com/p- 980174717.html

[16] 胡國明. 顆粒系統(tǒng)的離散元素法分析仿真[M]. 武漢：武漢理工大學出版社，2010.

[17] 王國強. 離散單元法及其在EDEM上的實踐[M]. 西安：西北工業(yè)大學出版社，2010.

[18] 邵明安，王全九，黃明斌. 土壤物理學[M]. 北京：高等教育出版社，2006：20－51.

[19] 鄧時琴.土壤礦質(zhì)顆粒及土壤質(zhì)地[J].土壤，1983(2)：77－80. Deng Shiqin. Soil mineral particles and soil texture[J]. Soil, 1983(2): 77－80. (in Chinese with English abstract)

[20] 董曉麗，李文利. 土力學與基礎工程[M]. 北京：清華大學出版社，2009：17－42.

[21] Academia. Some useful numbers for rock sand soils [R/OL].

2014-12-28 [2015-11-17]. http://www.academia.edu/4056287/ Some_Useful_Numbers_for_rocks_and_soils. 2015.

[22] 錢立彬.基于離散元法的開溝器的數(shù)字化設計方法研究[D]. 長春：吉林大學，2009：15－36.

Qian Libin. Digital Design Method of Furrow Openers Based On DEM[D]. Changchun: Jilin University, 2009: 15－36.

[23] Cheng Y P, Nakata Y, Bolton M D. Discrete element simulation of crushable soil[J]. Geotechnique, 2003, 53(07): 633－641.

[24] 李艷潔，吳騰，林劍輝，等. 基于離散元法的貫入圓錐對沙土顆粒運動特性分析[J]. 農(nóng)業(yè)工程學報，2012，28(24)：55－61.

Li Yanjie, Wu Teng, Lin Jianhui, et al. Influence of penetrating cone on motion characteristics of sandy soil particle using discrete element method[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2012, 28(24): 55－61. (in Chinese with English abstract)

[25] 文建鵬.長沙典型粉質(zhì)粘土的三軸試驗及本構模型研究[D]. 長沙：中南大學，2004. Wen Jianpeng. Triaxial Test and Constitutive Model on the Silty Clay in Changsha Area[D]. Changsha: Central South University, 2004. (in Chinese with English abstract)

[26] 楊菊. 基于擴展有限元法和離散元法的土壤-耕具接觸研究[D]. 昆明：昆明理工大學，2014. Yang Ju. Analysis of Soil-Tool Interaction Using Extended Finite Element Method and Discrete Element Method[D]. Kunming: Kunming University of Science and Technology, 2014. (in Chinese with English abstract)

[27] 孫紅月. 土力學實驗指導[M]. 北京：中國水利水電出版社，2010：50－79.

[28] 董曉麗，李文利. 土力學與基礎工程[M]. 北京：清華大學出版社，2009：110－135.

[29] 南京水利科學研究院土工研究所. 土工試驗技術手冊[M]. 北京：人民交通出版社，2007：33－49.

[30] Das B M. Advanced Soil Mechanics[M]. New York: Taylor & Francis, 1997: 14－20.

石林榕，趙武云，孫 偉. 基于離散元的西北旱區(qū)農(nóng)田土壤顆粒接觸模型和參數(shù)標定[J]. 農(nóng)業(yè)工程學報，2017，33(21)：181－187. doi：10.11975/j.issn.1002-6819.2017.21.021 http://www.tcsae.org

Shi Linrong, Zhao Wuyun, Sun Wei. Parameter calibration of soil particles contact model of farmland soil in northwest arid region based on discrete element method[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2017, 33(21): 181－187. (in Chinese with English abstract) doi：10.11975/j.issn.1002-6819.2017.21.021 http://www.tcsae.org

Parameter calibration of soil particles contact model of farmland soil in northwest arid region based on discrete element method

Shi Linrong, Zhao Wuyun※, Sun Wei

(,730070,)

In order to solve simulation distortion problems of soil resistance and dynamic movement with discrete element method, the hysteretic spring contact model (HSCM) and the liner cohesion model (LCM) were used to simulate farmland soils in order to investigate their soil plasticity behavior and dynamics, which could be used for predicting the resistance of farm tools. To achieve better simulation results, soil basic physical parameters, contact mechanical parameters and contact model parameters of HSCM and LCM were needed. Soil basic physical parameters included particle size distribution, soil particle density, elasticity modulus, shear modulus, poisson's ratio and bulk density. Among them, the elasticity modulus was obtained by the related literature. Soil contact mechanical parameters included cohesion between soil particles formed due to the action of the water liquid bridge, cohesive strength and internal friction angle obtained from soil shear test with 6 kinds of different soil moisture. Contact model parameters were composed of coefficient of restitution, static friction coefficient and dynamic friction coefficient. The soil particles were 5 mm so that the simulation time could be shortened. Because soil simulation particles were larger than actual ones, additionally in order to achieve better simulation results, a simulation model of soil contact parameter prediction model was established according to central experimental design method (Box-Behnken) used for a comprehensive analysis of fewer trial times. Taking static friction coefficient, dynamic friction coefficient and shear strength as test factors, and soil accumulation angle as response value, 17 sets of simulation tests were carried out according to response surface method test arrangement. The results showed that the influence of static friction coefficient, dynamic friction coefficient and the interaction between dynamic friction coefficient and cohesion strength on repose angle were significant. Interaction between static friction coefficient and dynamic friction coefficient, quadratic term of static friction coefficient had a significant effect on repose angle. Taking repose angle of soil with moisture of 1%, 4%, 8%, 12%, 16% and 20% as the response value, fixing the cohesion strength, we predicted the static friction coefficient and dynamic friction coefficient using the established regression model. Then, the results were validated by soil shear test. The results showed that the relative error of shear strength was 1.18%-9.31%, and the relative error of internal friction angle was 0.55%-4.07%. In order to verify the feasibility of the proposed soil contact model, a test by using duckbill opening hole test device in soil was carried out and the results showed that when the displacement was 0-50 mm the deviation of simulation and test was big. But when the displacemtn was 50 mm, the relative error was 0.928%, which could be used for design of duckbill structure in the future.

soils; discrete element method; friction coefficient; parameters; calibration

10.11975/j.issn.1002-6819.2017.21.021

O347.7;S22

A

1002-6819(2017)-21-0181-07

2017-03-05

2017-10-10

農(nóng)業(yè)部公益性行業(yè)科研專項（2015-03105）；甘肅省高等學?？蒲许椖浚?016B-039）；甘肅農(nóng)業(yè)大學盛彤笙科技創(chuàng)新基金（GSAU-SLS-1623）

石林榕，助教，主要從事西北旱區(qū)農(nóng)業(yè)工程技術與裝備研究。Email：getmao@qq.com

※通信作者：趙武云，教授，博士生導師，主要從事農(nóng)業(yè)工程技術與裝備研究。Email：zhaowy@gsau.edu.cn

中國農(nóng)業(yè)工程學會高級會員：石林榕（E041200842S）