石林榕，趙武云，孫 偉

?

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

石林榕，趙武云※，孫 偉

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

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

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

0 引 言

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

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

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

1 離散元土壤模型建立

1.1 接觸模型整合

試驗(yàn)土壤為沙壤土，既有散粒體物料特性又有一定的壓縮性。采用的土壤接觸模型為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為分別為顆粒和的屈服強(qiáng)度，Pa；屈服強(qiáng)度的默認(rèn)值（Pa）可以用楊氏模量和顆粒半徑預(yù)估。

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

卸載剛度為

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

每個(gè)時(shí)間步殘余重疊量由式（5）更新。

切向接觸力定義為

法向和切向阻尼力分別為

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

為制備不同含水率土壤進(jìn)行堆積角試驗(yàn)和直接剪切試驗(yàn)，通過(guò)式（18）制備。

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

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

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

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

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

1.2.2 接觸力學(xué)參數(shù)

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

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

圖2 試驗(yàn)中6種含水率土壤堆積角

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

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

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

表2 仿真試驗(yàn)方案及結(jié)果

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

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

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

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

表3 預(yù)測(cè)的土壤靜摩擦系數(shù)和動(dòng)摩擦系數(shù)

1.2.3 參數(shù)校核

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

圖4 直剪仿真試驗(yàn)

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

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

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

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

圖5 土壤貫入試驗(yàn)工具及壓力-沉陷曲線

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

2 驗(yàn)證試驗(yàn)

2.1 鴨嘴插入土壤模型仿真

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

2.2 結(jié)果與分析

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

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

注：Re為相對(duì)誤差。

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

3 結(jié) 論

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

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

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

[2] 趙成剛. 土力學(xué)的現(xiàn)狀及其數(shù)值分析方法中某些問(wèn)題的討論[J]. 巖土力學(xué)，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] 張銳，李建橋，李因武，等. 部件復(fù)雜表面影響土壤擾動(dòng)行為的離散元宏細(xì)觀分析[J]. 吉林大學(xué)學(xué)報(bào)：工學(xué)版，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]. 中國(guó)農(nóng)業(yè)大學(xué)學(xué)報(bào)，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] 胡國(guó)明. 顆粒系統(tǒng)的離散元素法分析仿真[M]. 武漢：武漢理工大學(xué)出版社，2010.

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

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

[19] 鄧時(shí)琴.土壤礦質(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] 董曉麗，李文利. 土力學(xué)與基礎(chǔ)工程[M]. 北京：清華大學(xué)出版社，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ù)字化設(shè)計(jì)方法研究[D]. 長(zhǎng)春：吉林大學(xué)，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] 李艷潔，吳騰，林劍輝，等. 基于離散元法的貫入圓錐對(duì)沙土顆粒運(yùn)動(dòng)特性分析[J]. 農(nóng)業(yè)工程學(xué)報(bào)，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ǎng)沙典型粉質(zhì)粘土的三軸試驗(yàn)及本構(gòu)模型研究[D]. 長(zhǎng)沙：中南大學(xué)，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] 楊菊. 基于擴(kuò)展有限元法和離散元法的土壤-耕具接觸研究[D]. 昆明：昆明理工大學(xué)，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] 孫紅月. 土力學(xué)實(shí)驗(yàn)指導(dǎo)[M]. 北京：中國(guó)水利水電出版社，2010：50－79.

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

[29] 南京水利科學(xué)研究院土工研究所. 土工試驗(yàn)技術(shù)手冊(cè)[M]. 北京：人民交通出版社，2007：33－49.

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

石林榕，趙武云，孫 偉. 基于離散元的西北旱區(qū)農(nóng)田土壤顆粒接觸模型和參數(shù)標(biāo)定[J]. 農(nóng)業(yè)工程學(xué)報(bào)，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è)科研專項(xiàng)（2015-03105）；甘肅省高等學(xué)?？蒲许?xiàng)目（2016B-039）；甘肅農(nóng)業(yè)大學(xué)盛彤笙科技創(chuàng)新基金（GSAU-SLS-1623）

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

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

中國(guó)農(nóng)業(yè)工程學(xué)會(huì)高級(jí)會(huì)員：石林榕（E041200842S）