海南熱區(qū)磚紅壤顆粒離散元仿真模型參數(shù)標(biāo)定

邢潔潔，張 銳，吳 鵬，張喜瑞，董學(xué)虎，陳 瑛，汝紹鋒

海南熱區(qū)磚紅壤顆粒離散元仿真模型參數(shù)標(biāo)定

邢潔潔1，張 銳1，吳 鵬1，張喜瑞1※，董學(xué)虎2，陳 瑛3，汝紹鋒1

（1. 海南大學(xué)機(jī)電工程學(xué)院，?？?570228；2. 中國熱帶農(nóng)業(yè)科學(xué)院農(nóng)業(yè)機(jī)械研究所，湛江 524091；3.加拿大曼尼托巴大學(xué)，溫尼伯 R3T 2N2）

為獲取適用于海南熱區(qū)磚紅壤與觸土部件相互作用的離散元仿真模型參數(shù)，該文利用EDEM仿真軟件中Hertz-Mindlin with JKR接觸模型對海南地區(qū)磚紅壤進(jìn)行相關(guān)模型參數(shù)標(biāo)定，首先以含水率為7.8%、實(shí)際休止角為32.35°的磚紅壤為研究對象，以物理試驗(yàn)獲取的磚紅壤-磚紅壤、磚紅壤-觸土部件間恢復(fù)系數(shù)、靜摩擦系數(shù)及滾動(dòng)摩擦系數(shù)為標(biāo)定對象，應(yīng)用Plackett-Burman 設(shè)計(jì)試驗(yàn)篩選出對休止角有顯著影響的參數(shù)：即磚紅壤-磚紅壤滾動(dòng)摩擦系數(shù)、JKR表面能、磚紅壤-28MnB5板滾動(dòng)摩擦系數(shù)、磚紅壤-磚紅壤恢復(fù)系數(shù)。進(jìn)一步以磚紅壤顆粒休止角為響應(yīng)值，基于Box-Behnken設(shè)計(jì)試驗(yàn)得到休止角與顯著性參數(shù)的二階回歸模型，并以實(shí)際休止角為目標(biāo)，針對顯著性參數(shù)進(jìn)行尋優(yōu)，得到最佳組合：磚紅壤-磚紅壤滾動(dòng)摩擦系數(shù)為0.07、磚紅壤-28MnB5板滾動(dòng)摩擦系數(shù)為0.24、磚紅壤-磚紅壤恢復(fù)系數(shù)為0.4、JKR表面能為4.31 J/m2。最后在標(biāo)定的參數(shù)下進(jìn)行休止角與破土阻力離散元仿真驗(yàn)證試驗(yàn)，對比計(jì)算得出仿真休止角與實(shí)測休止角相對誤差為0.62%，仿真破土阻力與實(shí)測破土阻力相對誤差為3.43%，并通過對比分析兩次試驗(yàn)中破土阻力變化曲線的擬合情況，得出兩曲線間的可決系數(shù)2=0.993 9，均方根誤差RMSE=0.017 7，該結(jié)果表明標(biāo)定所得相關(guān)參數(shù)可用作海南熱區(qū)磚紅壤離散元仿真。

磚紅壤；離散元模型；休止角；參數(shù)標(biāo)定；破土阻力

0 引 言

磚紅壤是中國最南端熱帶雨林或季雨林地區(qū)的地帶性土壤，相較于普通土壤有更強(qiáng)的膨脹性、流塑性及持水性，具有質(zhì)地黏重、土層深厚、比阻大等特點(diǎn)。對磚紅壤地進(jìn)行耕作作業(yè)時(shí)，會出現(xiàn)阻力大、效率低、觸土部件土壤粘附嚴(yán)重等問題[1-2]，嚴(yán)重影響生產(chǎn)效率。基于磚紅壤自身結(jié)構(gòu)的復(fù)雜性，研發(fā)適用于海南熱區(qū)磚紅壤地的高性能耕作裝備十分必要。在設(shè)計(jì)、優(yōu)化關(guān)鍵觸土部件時(shí)，往往需要解析磚紅壤-觸土部件的相互作用耦合機(jī)制，而實(shí)際土槽試驗(yàn)或田間試驗(yàn)僅能了解土壤宏觀運(yùn)動(dòng)情況，無法對土壤微觀運(yùn)動(dòng)規(guī)律進(jìn)行探究。目前國內(nèi)外已有專家學(xué)者利用有限元法探究耕作過程中土壤動(dòng)態(tài)擾動(dòng)行為[3-4]，然而有限元法是基于連續(xù)介質(zhì)的表達(dá)形式模擬土壤整體破壞行為，無法模擬土壤顆粒運(yùn)動(dòng)過程[5-6]。近幾年來，隨著農(nóng)業(yè)機(jī)械化水平的不斷提高，由Cundall等[7]提出的用以解決不連續(xù)介質(zhì)問題的離散元法逐漸被應(yīng)用在農(nóng)業(yè)工程領(lǐng)域，該方法可用于模擬顆粒間微觀和宏觀變形，在土壤-土壤、土壤-觸土部件間的接觸力學(xué)模型分析方面具有廣泛應(yīng)用前景[8-12]。

農(nóng)機(jī)具觸土部件在田間進(jìn)行耕作作業(yè)時(shí)，由于土壤類型的不同，土壤-觸土部件間的相互作用也會存在差異[13]。因此，對不同的土壤進(jìn)行離散元建模時(shí)，需要根據(jù)土壤的具體情況，選擇合適的土壤接觸模型并對相關(guān)模型參數(shù)進(jìn)行準(zhǔn)確性標(biāo)定[14]，以確保離散元仿真試驗(yàn)的真實(shí)性。Santos[15]通過旋轉(zhuǎn)滾筒法獲取了櫻桃果干顆粒的休止角，并標(biāo)定出顆粒間的接觸參數(shù)；Mustafa等[16]研究了土壤分別在有粘結(jié)力和無粘結(jié)力情況下顆粒塑性變形問題；王憲良等[17]驗(yàn)證了彈塑性接觸模型（Edinburgh elasto-plastic cohesion model，ECM）符合土壤本構(gòu)關(guān)系，并對土壤顆粒間的靜摩擦系數(shù)和滾動(dòng)摩擦系數(shù)進(jìn)行了標(biāo)定；石林榕等[18]建立了西北旱區(qū)農(nóng)田土壤顆粒接觸模型并對關(guān)鍵參數(shù)進(jìn)行離散元標(biāo)定；張銳等[19]提出一種標(biāo)定沙土顆粒仿真參數(shù)的方法，試驗(yàn)得出顆粒形狀對顆粒間靜摩擦系數(shù)影響較大；丁啟朔等[20]建立了濕粘水稻土的深松耕作離散元模型并對相關(guān)仿真參數(shù)進(jìn)行標(biāo)定。

目前，關(guān)于土壤建模參數(shù)標(biāo)定的研究主要是針對黏性較弱如沙土、旱土或巖石等土壤類型，對于黏性較大的磚紅壤土并不適用。本文針對磚紅壤離散元模型參數(shù)不清的問題，以海南熱區(qū)磚紅壤樣品為研究對象，選取EDEM軟件中“Hertz-Mindlin with JKR”（JKR）接觸模型，通過土壤休止角與破土阻力仿真對比試驗(yàn)，完成對磚紅壤土離散元模型參數(shù)的標(biāo)定，以期為后續(xù)探究磚紅壤-觸土部件耦合機(jī)理提供基礎(chǔ)數(shù)據(jù)。

1 材料與方法

1.1 接觸模型的選取

磚紅壤質(zhì)地厚重且存在粘附現(xiàn)象，經(jīng)典的Hertz-Mindlin接觸模型僅考慮彈性變形，不考慮顆粒模型間的粘結(jié)力，難以準(zhǔn)確模擬觸土部件在磚紅壤中的作業(yè)規(guī)律；Hertz-Mindlin with Bonding接觸模型雖然可用來模擬粘結(jié)顆粒，但該模型是利用有限大小的“膠粘劑”粘結(jié)，適用于模擬混凝土和巖石等較堅(jiān)硬介質(zhì)[21]；而Hertz-Mindlin with JKR接觸模型是一種建立在Hertz理論上的粘結(jié)性顆粒接觸模型，同時(shí)考慮到濕顆粒間粘結(jié)力對顆粒運(yùn)動(dòng)的影響，適用于模擬顆粒間因靜電、水分等原因發(fā)生明顯粘結(jié)和團(tuán)聚的物料，如農(nóng)作物和泥土等[22-23]。“Hertz-Mindlin with JKR”接觸模型中的法向彈性力可較好的表征顆粒間的粘彈性特征。

Hertz-Mindlin with JKR接觸模型的法向彈性力基于法向重疊量和表面能，即

式中JKR為JKR法向彈性力，N；為兩接觸顆粒之間法向重疊量，m；為兩接觸顆粒之間切向重疊量，m；為表面能，N/m；*為等效彈性模量，Pa；*為等效接觸半徑，m。等效彈性模量*與等效接觸半徑*定義為

式中1、2分別為兩接觸顆粒的彈性模量，Pa；1、2分別為兩接觸顆粒的泊松比；1、2分別為兩接觸顆粒的接觸半徑，m。

當(dāng)表面能=0時(shí)，JKR法向彈性力變?yōu)镠ertz-Mindlin法向力：

即使顆粒并不是直接接觸，Hertz-Mindlin with JKR接觸模型也提供吸引凝聚力[23]，顆粒間具有非零凝聚力的最大間隙為

式中c為顆粒間具有非零凝聚力時(shí)的法向最大間隙，m；c為顆粒間具有非零凝聚力時(shí)的切向最大間隙，m。

當(dāng)＞c時(shí)，顆粒間凝聚力變?yōu)?；當(dāng)顆粒并非實(shí)際接觸且間隔小于c時(shí)，凝聚力達(dá)到最大值：

式中cohesion為2個(gè)顆粒間凝聚力，N；*為等效接觸半徑，m。

海南熱區(qū)雨水較多，Hertz-Mindlin with JKR接觸模型也可模擬濕潤磚紅壤顆粒，此時(shí)將2個(gè)顆粒分開所需要的分離力取決于液體表面張力和濕潤角：

式中pullout為將2個(gè)顆粒分開所需的分離力，N；s為液體表面張力，N；為濕潤角，(°)；i、j為顆粒半徑，m。

1.2 模型參數(shù)的類型及標(biāo)定

模型參數(shù)包括材料本征參數(shù)、接觸參數(shù)和接觸模型參數(shù)，其中材料本征參數(shù)與接觸模型參數(shù)較為固定，可通過查閱資料得到。故本文主要針對磚紅壤模型接觸參數(shù)進(jìn)行試驗(yàn)標(biāo)定。

1.2.1 試驗(yàn)材料及本征參數(shù)

本次試驗(yàn)土壤樣品采自海南大學(xué)儋州市香蕉地試驗(yàn)田，地處19°11′～19°52′N，108°56′～109°45′E，測得采集土樣的含水率為7.8%，由于土樣不均勻，對土樣進(jìn)行風(fēng)干、碾碎、過篩處理后，根據(jù)土樣的pH值，加入適量氫氧化鈉溶液制成懸浮液并定容，根據(jù)文獻(xiàn)[24]所述方法計(jì)算各級顆粒的質(zhì)量，并計(jì)算其百分含量。測定出磚紅壤顆粒組成情況如表1所示。

表1 磚紅壤顆粒組成

磚紅壤顆粒材料本征參數(shù)為[25]：泊松比0.38；剪切模量1×106Pa；密度1 850 kg/m3；接觸模型參數(shù)主要為JKR表面能，將仿真顆粒堆積密度和仿真規(guī)模導(dǎo)入EDEM通用顆粒材料數(shù)據(jù)庫（Generic EDEM material model database, GEMM）中，得JKR表面能為3.8～10.2 J/m2。

1.2.2 磚紅壤休止角物理試驗(yàn)

顆粒狀物質(zhì)被傾倒于水平面時(shí)會形成錐體，錐體表面與水平面所成內(nèi)角即為休止角，是表征顆粒物料流動(dòng)、摩擦等特性的宏觀參數(shù)，與顆粒密度、顆粒表面積及顆粒摩擦系數(shù)有關(guān)，因此，休止角物理試驗(yàn)常被用作顆粒物料離散元參數(shù)標(biāo)定[26-27]。試驗(yàn)樣品土壤堆積密度為2.51 kg/m3，土壤含水率7.8%，采用文獻(xiàn)[28]所述方法進(jìn)行土壤休止角試驗(yàn)，試驗(yàn)結(jié)果如圖1所示，重復(fù)10次測定取磚紅壤實(shí)際休止角平均值為32.35°。

1.接料板 2.磚紅壤顆粒

1.2.3 恢復(fù)系數(shù)測定試驗(yàn)

在耕作作業(yè)過程中，土壤-土壤、土壤-觸土部件間會存在各種碰撞和擠壓，從而造成觸土部件不同程度的損壞?；謴?fù)系數(shù)是體現(xiàn)碰撞特性的力學(xué)參數(shù)，定義為碰撞后法向反彈速度與碰撞前法向靠近速度的比值[29]?；謴?fù)系數(shù)通常采用自由落體試驗(yàn)或平板碰撞試驗(yàn)測定，由于土壤顆粒大小不均勻，采用自由落體試驗(yàn)時(shí)，被測顆粒常發(fā)生旋轉(zhuǎn)或多點(diǎn)碰撞，碰撞后很難產(chǎn)生垂直方向上的彈回，測定結(jié)果會產(chǎn)生較大誤差。本試驗(yàn)中，基于運(yùn)動(dòng)學(xué)原理，采用圖2所示斜板碰撞試驗(yàn)裝置測定恢復(fù)系數(shù)。

注：h為土壤顆粒下落點(diǎn)距28MnB5板的高度，m；S1、S2和H1、H2分別為土壤顆粒兩次碰撞28MnB5板后在水平面內(nèi)、豎直平面內(nèi)的位移，m；28MnB5板傾角為45°；X、Y為坐標(biāo)軸；Vx、Vy分別為顆粒碰撞28MnB5板后沿X方向和Y方向的速度，m/s。

磚紅壤顆粒下落至28MnB5板處與板發(fā)生碰撞做斜拋運(yùn)動(dòng)，可分解為沿方向的勻速運(yùn)動(dòng)（速度V）和沿方向初速度為V，加速度為的勻變速運(yùn)動(dòng)，最終落在接料底板上，即：

式中S、H（=1, 2）分別為土壤顆粒碰撞28MnB5板后在水平面內(nèi)、豎直平面內(nèi)的位移，m。

改變接料底板的高度1，2，分別測得兩次磚紅壤顆粒碰撞后在接料底板上的水平位移量1，2，可得：

式中1、2和1、2分別為兩次土壤顆粒碰撞28MnB5板后在水平面內(nèi)、豎直平面內(nèi)的位移。

由恢復(fù)系數(shù)定義可求得恢復(fù)系數(shù)r得：

式中V為碰撞后法向反彈速度，m/s；′為碰撞前法向靠近速度，m/s；0為碰撞前顆粒速度，m/s。

磚紅壤-觸土部件間恢復(fù)系數(shù)通過磚紅壤顆粒撞擊28MnB5板測定；磚紅壤-磚紅壤間恢復(fù)系數(shù)可采用粘接劑將磚紅壤土顆粒均勻粘附在斜板表面，使其形成顆粒板，通過磚紅壤顆粒撞擊顆粒板測定。通過多次反復(fù)試驗(yàn)測定，最終得到：磚紅壤顆粒間恢復(fù)系數(shù)范圍為0.1～0.7，磚紅壤顆粒與28MnB5材料板間恢復(fù)系數(shù)范圍為0.05～0.7。

1.2.4 靜摩擦系數(shù)測定試驗(yàn)

在土壤實(shí)際耕作過程中，土壤-土壤、土壤-觸土部件之間大部分為多點(diǎn)接觸，產(chǎn)生的摩擦大多為靜摩擦；在離散元仿真試驗(yàn)中，兩顆粒間受到的切向合力通常是通過靜摩擦力來限制的，因此在采用離散元模擬土壤與觸土部件相互作用過程時(shí)，靜摩擦系數(shù)對模擬結(jié)果影響較大。

設(shè)計(jì)如圖3所示的斜板試驗(yàn)，通過改變斜板高度H，當(dāng)磚紅壤顆粒開始滑落時(shí)測定斜板傾角，進(jìn)而求得靜摩擦系數(shù)，反復(fù)試驗(yàn)測定20次，取平均值。磚紅壤顆粒剛開始滑落時(shí)在斜板上的平衡方程為

注：1為拉力，N；2為壓力，N；為重力，N；為摩擦力，N；為支持力，N；為靜摩擦系數(shù)；為斜板傾角，(°)；為底座長度，mm；H為斜面高度，mm。

Note:1is tension, N;2is pressure, N;is gravity, N;is friction, N;is support force, N;is angle between base and slope, (°);is the base length, mm;His the height of the slope, mm.

圖3 靜摩擦系數(shù)測定試驗(yàn)原理圖

Fig.3 Principle diagram of static friction coefficient test

磚紅壤顆粒在斜板上的靜摩擦系數(shù)為

式中1為拉力，N；2為壓力，N；為重力，N；為摩擦力，N；為支持力，N；為靜摩擦系數(shù)；為斜板傾角，(°)。

其中，磚紅壤-觸土部件間靜摩擦系數(shù)通過磚紅壤顆粒在28MnB5板滑動(dòng)試驗(yàn)測定；磚紅壤-磚紅壤間靜摩擦系數(shù)可通過磚紅壤顆粒與上述顆粒板滑動(dòng)試驗(yàn)測定。反復(fù)多次試驗(yàn)取平均值，最終得到磚紅壤-磚紅壤間靜摩擦系數(shù)為0.4～1.1，磚紅壤-觸土部件間靜摩擦系數(shù)為0.3～0.9。

1.2.5 滾動(dòng)摩擦系數(shù)測定試驗(yàn)

設(shè)計(jì)如圖4所示試驗(yàn)裝置測量滾動(dòng)摩擦系數(shù)，質(zhì)量為的磚紅壤顆粒在半徑為的導(dǎo)軌角處自由滾落，測定滾落高度、水平拋出高度及水平拋出距離，由微元法可知，當(dāng)顆粒轉(zhuǎn)動(dòng)一個(gè)較小的角度時(shí)，對導(dǎo)軌的正壓力近似不變，即：

式中為磚紅壤顆粒在整個(gè)運(yùn)動(dòng)過程中的摩擦力，N；為滾動(dòng)摩擦系數(shù)；為導(dǎo)軌傾角，(°)；為磚紅壤顆粒質(zhì)量，kg；為重力加速度，N/kg；為磚紅壤顆粒轉(zhuǎn)動(dòng)的較小的一個(gè)角度，(°)。

注：為導(dǎo)軌傾角，(°)；為導(dǎo)軌半徑，m；h為磚紅壤顆粒拋出點(diǎn)距地面高度，m；為水平拋出距離，m；為磚紅壤顆粒滾落高度，m。

Note:is angle between base and guide rail, (°);is guide radius, m;his height of latosol particles throwing point from the ground, m;is horizontal throw distance, m;is the rolling height of latosol particles, m.

圖4 滾動(dòng)摩擦系數(shù)測定試驗(yàn)原理圖

Fig.4 Principle diagram of rolling friction coefficient test

在此過程中，滑動(dòng)摩擦力所做的功1為

式中為導(dǎo)軌半徑，m。

磚紅壤顆粒在距地面h處水平拋出的動(dòng)能2為

式中為磚紅壤顆粒從導(dǎo)軌拋出時(shí)速度，m/s。

重力所做的功為

=（19）

式中為磚紅壤顆粒滾落高度，m；h為磚紅壤顆粒拋出點(diǎn)距地面高度，m。

由磚紅壤顆粒運(yùn)動(dòng)過程中的能量守恒，即12可求出滾動(dòng)摩擦系數(shù)，若要測量磚紅壤顆粒間滾動(dòng)摩擦系數(shù)，仍可采用上述粘接劑形成顆粒板，再將磚紅壤顆粒放置于顆粒板上自由滾落，多次重復(fù)試驗(yàn)測定求平均值得：磚紅壤-磚紅壤間滾動(dòng)摩擦系數(shù)為0.05～0.25，磚紅壤-觸土部件間滾動(dòng)摩擦系數(shù)為0.05～0.25。

1.3 參數(shù)標(biāo)定試驗(yàn)設(shè)計(jì)

綜上試驗(yàn)結(jié)果：磚紅壤顆粒之間的恢復(fù)系數(shù)范圍為0.1～0.7，靜摩擦系數(shù)范圍為0.4～1.1，滾動(dòng)摩擦系數(shù)范圍為0.05～0.25；磚紅壤土顆粒與觸土部件之間的恢復(fù)系數(shù)范圍為0.05～0.7，靜摩擦系數(shù)范圍為0.3～0.9，滾動(dòng)摩擦系數(shù)范圍為0.05～0.25；JKR表面能取值范圍為3.8～10.2 J/m2。應(yīng)用Plackett-Burman和Box-Behnken設(shè)計(jì)試驗(yàn)進(jìn)行參數(shù)標(biāo)定，其中試驗(yàn)因素選定為恢復(fù)系數(shù)、靜摩擦系數(shù)、滾動(dòng)摩擦系數(shù)和JKR表面能，試驗(yàn)指標(biāo)為休止角，確定各因素水平表如表2所示，共7個(gè)因素，各取一個(gè)高水平、一個(gè)低水平。

表2 仿真試驗(yàn)因素及水平

注：L為磚紅壤；M為28MnB5板。

Note: L is latosol; M is steel plate.

2 結(jié)果與分析

2.1 Plackett-Burman試驗(yàn)設(shè)計(jì)確定關(guān)鍵影響因素

Plackett-Burman是一種兩水平的試驗(yàn)設(shè)計(jì)方法，通過比較各個(gè)因子兩水平的差異與整體的差異來確定因子的顯著性。延用表2設(shè)置各因素水平，共進(jìn)行12次試驗(yàn)，試驗(yàn)方案及結(jié)果如表3所示，其中設(shè)置4個(gè)空白列用于誤差分析。

表3 Plackett-Burman Design 方案及結(jié)果

注：、、、為空白列。

Note:,,andindicate blank column.

對試驗(yàn)結(jié)果進(jìn)行分析，得到7個(gè)因素對休止角貢獻(xiàn)度如表4所示。由表4可知，7個(gè)因素對休止角均為正效應(yīng)，即休止角隨著該7個(gè)因素增大而增大，其中對休止角影響較為顯著的4個(gè)因素依次為：（磚紅壤-磚紅壤滾動(dòng)摩擦系數(shù)）、（表面能）、（磚紅壤-28MnB5板滾動(dòng)摩擦系數(shù)）和（磚紅壤-磚紅壤恢復(fù)系數(shù)），合計(jì)貢獻(xiàn)度為95.68%。

表4 Plackett-Burman Design試驗(yàn)結(jié)果分析

2.2 Box-Behnken試驗(yàn)設(shè)計(jì)

2.2.1 試驗(yàn)方案與結(jié)果

得到休止角顯著影響因素后，應(yīng)用Box-Behnken設(shè)計(jì)試驗(yàn)進(jìn)行響應(yīng)面分析并尋找最優(yōu)解，以（磚紅壤-磚紅壤滾動(dòng)摩擦系數(shù)）、（表面能）、（磚紅壤-28MnB5板滾動(dòng)摩擦系數(shù)）和（磚紅壤-磚紅壤恢復(fù)系數(shù)）為試驗(yàn)因素，休止角為試驗(yàn)指標(biāo)，延用表2各因素試驗(yàn)水平，共進(jìn)行29次試驗(yàn)，試驗(yàn)方案及結(jié)果如表5所示。

表5 Box-Behnken Design 方案及結(jié)果

2.2.2 休止角回歸模型分析

應(yīng)用Design-expert軟件對試驗(yàn)結(jié)果進(jìn)行多元回歸擬合分析，可得休止角回歸模型：

該回歸方程的決定系數(shù)2=0.991 1，表明有99.11%試驗(yàn)差異可用該模型解釋，與實(shí)際數(shù)據(jù)的擬合程度高；校正決定系數(shù)adj-2=0.982 2與2非常接近1，表明相關(guān)性較好；變異系數(shù)CV=0.60%，綜上表明該模型具有較好的可靠性，可以用于進(jìn)一步分析。對該模型試驗(yàn)結(jié)果進(jìn)行回歸方差分析，結(jié)果如表6所示。由表6可以看出，模型<0.000 1，說明該休止角回歸模型極顯著；在給定的試驗(yàn)因素水平范圍內(nèi)，（磚紅壤-磚紅壤恢復(fù)系數(shù)）、（磚紅壤-磚紅壤滾動(dòng)摩擦系數(shù)）、（磚紅壤-28MnB5板滾動(dòng)摩擦系數(shù)）、（JKR）對休止角影響極顯著；、、對休止角影響極顯著，、對休止角影響較為顯著，對休止角影響不顯著；2、2、2對休止角影響極顯著，2對休止角影響較為顯著。

表6 Box-Behnken Design 二次回歸模型方差分析

2.2.3 參數(shù)優(yōu)化與試驗(yàn)驗(yàn)證

在設(shè)置的試驗(yàn)因素水平范圍內(nèi)，基于Design-expert軟件的優(yōu)化模塊，以實(shí)際測得休止角32.35°為目標(biāo)對休止角回歸模型進(jìn)行尋優(yōu)，對得到的若干組解進(jìn)行休止角仿真驗(yàn)證，得到與物理試驗(yàn)形狀相近的一組最優(yōu)解，即（磚紅壤-磚紅壤恢復(fù)系數(shù)）為0.4，（磚紅壤-磚紅壤滾動(dòng)摩擦系數(shù)）為0.07，（磚紅壤-28MnB5板滾動(dòng)摩擦系數(shù)）為0.24，（JKR表面能）為4.31J/m2；在EDEM軟件中將、、、設(shè)置為上述最優(yōu)解，其余參數(shù)選取中間水平，即（磚紅壤-磚紅壤靜摩擦系數(shù)）0.75、（磚紅壤-28MnB5板恢復(fù)系數(shù)）0.38、（磚紅壤-28MnB5板靜摩擦系數(shù)）0.6，該優(yōu)化解下得到的休止角仿真結(jié)果與實(shí)際物理試驗(yàn)結(jié)果對比如圖5所示。

圖5 仿真試驗(yàn)與物理試驗(yàn)對比

在休止角測量試驗(yàn)中，散料堆會出現(xiàn)“凹形表面”和“凸形表面”等現(xiàn)象[30]，堆積邊界參差不齊，使得直接測量或采用EDEM軟件后處理量角器方法測量結(jié)果不夠準(zhǔn)確。為更精確測得休止角，本次試驗(yàn)借助Python的matplotlib模塊、pillow模塊及scikit learn的Linear regression模塊對試驗(yàn)圖像進(jìn)行處理并求解，處理過程如下：

1）讀取磚紅壤土堆單側(cè)圖像并進(jìn)行圖像的灰度及二值化處理（圖6a）；在處理完的圖片上選取坡度邊界，對獲取到的數(shù)據(jù)進(jìn)行最小二乘線性擬合，進(jìn)而確定擬合方程斜率（圖6b）；

2）按式（21）計(jì)算休止角：

式中ζ為磚紅壤土堆測量休止角，(°)；k為斜率。

由圖7可知，磚紅壤休止角單側(cè)輪廓采用最小二乘法擬合出的線性方程為=0.62+0.95，斜率=0.62，帶入式（21）中得休止角為31.79°，重復(fù)5次試驗(yàn)，得到磚紅壤休止角計(jì)算值分別為：31.79°、32.93°、31.96°、32.48°和32.53°，平均值為32.33°，與磚紅壤休止角實(shí)際測量值32.35°相對誤差為0.62%。

2.3 破土阻力對比試驗(yàn)

為進(jìn)一步核實(shí)標(biāo)定的模型參數(shù)的準(zhǔn)確性與合理性，設(shè)計(jì)破土阻力對比驗(yàn)證試驗(yàn)：在EDEM軟件中設(shè)置上述標(biāo)定最優(yōu)參數(shù)，并生成磚紅壤顆粒；通過EDEM軟件仿真模擬破土刃在磚紅壤中破土過程，并進(jìn)行實(shí)際試驗(yàn)，對比破土阻力變化的擬合情況。

2.3.1 EDEM仿真試驗(yàn)

在SolidWorks軟件中建立大小為350 mm×200 mm× 150 mm的土槽與破土刃三維模型并導(dǎo)入EDEM軟件中，待完全生成土壤顆粒后，讓破土刃以1 m/s2的加速度從靜止加速至1 m/s，之后繼續(xù)勻速前進(jìn)直至完全穿過磚紅壤土層，耕作深度為80 mm，仿真過程如圖7所示。待仿真結(jié)束后通過后處理模塊導(dǎo)出運(yùn)動(dòng)過程中阻力變化數(shù)據(jù)。

圖7 破土過程仿真試驗(yàn)圖

2.3.2 土槽試驗(yàn)

搭建如圖8所示試驗(yàn)裝置，取相同大小土槽置于試驗(yàn)臺上，土槽內(nèi)磚紅壤取自儋州市香蕉試驗(yàn)田，土槽上方安裝直線導(dǎo)軌。調(diào)整試驗(yàn)臺高度，使得測力儀與直線導(dǎo)軌位于同一水平線上，由伸縮電機(jī)帶動(dòng)測力儀和破土刃運(yùn)動(dòng)，控制破土刃沿直線導(dǎo)軌以1 m/s2的加速度從靜止加速至1 m/s，并勻速前進(jìn)穿過磚紅壤土層，與仿真試驗(yàn)中的運(yùn)動(dòng)一致，記錄運(yùn)動(dòng)過程中測力計(jì)讀數(shù)。

1.土槽 2.直線導(dǎo)軌 3.破土刃 4.測力儀 5.磚紅壤

2.3.3 阻力對比分析

繪制仿真運(yùn)動(dòng)過程與實(shí)際試驗(yàn)過程阻力變化曲線，如圖9所示，對比分析兩條曲線可知，在破土刃剛開始破土?xí)r所受阻力較大，隨著破土刃逐漸松開土壤，所受阻力逐漸趨于穩(wěn)定，對比可得兩阻力曲線形狀趨于一致；仿真試驗(yàn)所得穩(wěn)定破土阻力3.31 N，與實(shí)際測量所得穩(wěn)定破土阻力3.2 N相對誤差為3.43%，相對誤差較小；使用scikit learn中的metrics模塊分別計(jì)算模型的決定系數(shù)2與均方根誤差RMSE，得出兩曲線間的2=0.993 9，趨于1，說明兩曲線擬合狀況越好；RMSE=0.017 7，趨于0，說明兩曲線的離散程度越小。綜上，表明本文所標(biāo)定磚紅壤土顆粒仿真模型參數(shù)準(zhǔn)確可靠。

圖9 仿真與試驗(yàn)過程破土阻力對比分析

3 結(jié) 論

1）針對海南熱區(qū)特有的磚紅壤土，基于EDEM離散元仿真軟件，選用Hertz-Mindlin with JKR接觸模型對含水率為7.8%的磚紅壤樣品進(jìn)行離散元仿真，通過試驗(yàn)測定磚紅壤-磚紅壤、磚紅壤-觸土部件間恢復(fù)系數(shù)依次為0.1～0.7、0.05～0.7；磚紅壤-磚紅壤、磚紅壤-觸土部件間靜摩擦系數(shù)依次為0.4～1.1、0.3～0.9；磚紅壤-磚紅壤、磚紅壤-觸土部件間滾動(dòng)摩擦系數(shù)均為0.05～0.25。

2）應(yīng)用Plackett-Burman Design試驗(yàn)設(shè)計(jì)確定磚紅壤休止角關(guān)鍵影響因素包括：磚紅壤-磚紅壤恢復(fù)系數(shù)、磚紅壤-磚紅壤滾動(dòng)摩擦系數(shù)、磚紅壤-28MnB5板滾動(dòng)摩擦系數(shù)及JKR表面能。

3）應(yīng)用Box-Behnken試驗(yàn)設(shè)計(jì)確定4個(gè)顯著影響因素最優(yōu)值：磚紅壤-磚紅壤恢復(fù)系數(shù)0.40、磚紅壤-磚紅壤滾動(dòng)摩擦系數(shù)0.07、磚紅壤-28MnB5板滾動(dòng)摩擦系數(shù)0.24、JKR表面能4.31J/m2。以最優(yōu)參數(shù)組合得到的仿真試驗(yàn)休止角與實(shí)際測量休止角相對誤差為0.62%。

4）在EDEM軟件中仿真模擬破土刃在磚紅壤中破土過程，并進(jìn)行實(shí)際對比驗(yàn)證試驗(yàn)。兩次穩(wěn)定破土阻力相對誤差為3.43%，表明標(biāo)定磚紅壤離散元仿真模型參數(shù)準(zhǔn)確可靠。

本文研究結(jié)果對海南熱區(qū)磚紅壤地作業(yè)具有一定指導(dǎo)意義；但需指出的是，本文雖做了實(shí)際試驗(yàn)，結(jié)果較為理想，但基于磚紅壤土的特殊性和香蕉根系的復(fù)雜性，后期將開展更為深入的研究，本文為后期的實(shí)際土壤建模提供了基本參數(shù)支撐，可為進(jìn)一步研制適用于磚紅壤香蕉地作業(yè)條件下的觸土部件提供基礎(chǔ)數(shù)據(jù)。

[1] 張駿達(dá),孫向陽，劉艷，等. 海南省儋州市磚紅壤在中國土壤系統(tǒng)分類的歸屬研究[J]. 廣東農(nóng)業(yè)科學(xué)，2015，42(10)：46－52.

Zhang Junda, Sun Xiangyang, Liu Yan, et al. Attribution of laterite of Danzhou，Hainan province in Chinese soil taxonomy[J]. Guangdong Agricultural Sciences, 2015, 42(10): 46－52. (in Chinese with English abstract)

[2] 嚴(yán)曉麗，汪春，李明，等. 紅壤鑿形深松鏟的設(shè)計(jì)與試驗(yàn)[J]. 農(nóng)機(jī)化研究，2019，41(4)：135－142

Yan Xiaoli, Wang Chun, Li Ming, et al. Design and experiment on the chisel-shape deep-shovel in red loam[J]. Journal of Agricultural Mechanization Research, 2019, 41(4): 135－142. (in Chinese with English abstract)

[3] 徐中華，王建華. 有限元法分析土壤切削問題的研究進(jìn)展[J]. 農(nóng)業(yè)機(jī)械學(xué)報(bào)，2005，36(1)：134－137.

Xu Zhonghua, Wang Jianhua. Advances in finite element analysis of soil cutting[J]. Transactions of the Chinese Society for Agricultural Machinery, 2005，36(1): 134－137. (in Chinese with English abstract)

[4] Abo-Elnor M, Hamilton R, Boyle J T. Simulation of soil-blade interaction for sandy soil using advanced 3D finite element analysis[J]. Soil & Tillage Research, 2004, 75(1): 61－73.

[5] 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(3): 105－117.

[6] 張銳，李建橋，周長海，等. 推土板表面形態(tài)對土壤動(dòng)態(tài)行為影響的離散元模擬[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2007，23(9)：13－19.

Zhang Rui, Li Jianqiao, Zhou Changhai, et al. Simulation of dynamic behavior of soil ahead of the bulldozing plates with different surface configurations by discrete element method[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2007, 23(9): 13－19. ( in Chinese with English abstract)

[7] Cundall P A, Strack O D. A discrete numerical model for granular assemblies[J]. Geotechnique, 1979, 29(1): 47－65.

[8] 馬云海，馬圣勝，賈洪雷，等. 仿生波紋形開溝器減黏降阻性能測試與分析[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2014，30(5)：36－41.

Ma Yunhai, Ma Shengsheng, Jia Honglei, et al. Measurement and analysis on reducing adhesion and resistance of bionic ripple opener[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2014, 30(5): 36－41. (in Chinese with English abstract)

[9] 鄭侃，何進(jìn)，李洪文，等. 基于離散元深松土壤模型的折線破土刃深松鏟研究[J]. 農(nóng)業(yè)機(jī)械學(xué)報(bào)，2016，47(9)：62－72.

Zheng Kan, He Jin, Li Hongwen, et al . Research on polyline soil-breaking blade subsoiler based on subsoiling soil model using discrete element method[J]. Transactions of the Chinese Society for Agricultural Machinery, 2016, 47(9): 62－72. (in Chinese with English abstract)

[10] 鄧佳玉，胡軍，李慶達(dá)，等. 基于EDEM離散元法的深松鏟仿真與試驗(yàn)研究[J]. 中國農(nóng)機(jī)化學(xué)報(bào)，2016，37(4)：14－18.

Deng Jiayu, Hu Jun, Li Qingda, et al. Simulation and experimental study on the subsoiler based on EDEM discrete element method[J]. Journal of Chinese Agricultural Mechanization, 2016, 37(4): 14－18. (in Chinese with English abstract)

[11] 方會敏，姬長英，Chandio F A，等. 基于離散元法的旋耕過程土壤運(yùn)動(dòng)行為分析[J]. 農(nóng)業(yè)機(jī)械學(xué)報(bào)，2016，47(3)：22－28.

Fang Huimin, Ji Changying, Chandio F A, et al. Analysis of soil dynamic behavior during rotary tillage based on distinct element method[J]. Transactions of the Chinese Society for Agricultural Machinery, 2016, 47(3): 22－28. (in Chinese with English abstract)

[12] 于建群，付宏，李紅，等. 離散元法及其在農(nóng)業(yè)機(jī)械工作部件研究與設(shè)計(jì)中的應(yīng)用[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2005，21(5)：1－6.

Yu Jianqun, Fu Hong, Li Hong, et al. Application of discrete element method to research and design of working parts of agricultural machines[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2005，21(5): 1－6. (in Chinese with English abstract)

[13] Tong Jin, Mohammad M A, Zhang Jinbo, et al. DEM numerical simulation of abrasive wear characteristics of a bioinspired ridged surface[J]. Journal of Bionic Engineering, 2010，7(2): 175－181.

[14] 龔明. 材料物料特性參數(shù)標(biāo)定[R]. 北京：EDEM用戶大會報(bào)告，2013.

[15] Santos K G, Campos A V P, Oliveira O S, et al. Dem simulations of dynamic angle of repose of acerola residue : a parametric study using a response surface technique[J]. Blucher Chemical Engineering Proceedings, 2015, 1(2): 11326－11333.

[16] Mustafa U, Fielke J M, Saunders C. 3D DEM tillage simulation: validation of a hysteretic spring (plastic) contact model for a sweep tool operating in a cohesionless soil [J]. Soil & Tillage Research, 2014,144: 220－227.

[17] 王憲良，胡紅，王慶杰，等. 基于離散元的土壤模型參數(shù)標(biāo)定方法[J]. 農(nóng)業(yè)機(jī)械學(xué)報(bào)，2017，48(12)：78－85.

Wang Xianliang, Hu Hong, Wang Qingjie, et al. Calibration method of soil contact characteristic parameters based on DEM theory[J]. Transactions of the Chinese Society for Agricultural Machinery, 2017, 48(12): 78－85. (in Chinese with English abstract)

[18] 石林榕，趙武云，孫偉. 基于離散元的西北旱區(qū)農(nóng)田土壤顆粒接觸模型和參數(shù)標(biāo)定[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2017，33(21)：181－187.

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)

[19] 張銳，韓佃雷，吉巧麗，等. 離散元模擬中沙土參數(shù)標(biāo)定方法研究[J]. 農(nóng)業(yè)機(jī)械學(xué)報(bào)，2017，48(3)：49－56.

Zhang Rui, Han Dianlei, Ji Qiaoli, et al. Calibration methods of sandy soil parameters in simulation of discrete element method[J]. Transactions of the Chinese Society for Agricultural Machinery, 2017, 48(3): 49－56. (in Chinese with English abstract)

[20] 丁啟朔，任駿，Adam B E，等. 濕粘水稻土深松過程離散元分析[J]. 農(nóng)業(yè)機(jī)械學(xué)報(bào)，2017，48(3)：38－48.

Ding Qishuo, Ren Jun, Adam B E, et al. DEM analysis of subsoiling process in wet clayey paddy soil [J]. Transactions of the Chinese Society for Agricultural Machinery, 2017, 48(3): 38－48. (in Chinese with English abstract)

[21] 李俊偉，佟金，胡斌，等. 不同含水率黏重黑土與觸土部件互作的離散元仿真參數(shù)標(biāo)定[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2019，35(6)：130－140.

Li Junwei, Tong Jin, Hu Bin, et al. Calibration of parameters of interaction between clayey black soil with different moisure content and soil-engaging component in northeast China[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2019, 35(6): 130－140. (in Chinese with English abstract)

[22] Johnson K L, Kendall K, Roberts A D. Surface energy and the contact of elastic solids[J]. Proc R Soc London: A, 1971, 324(158): 301－313.

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

[24] 陳麗瓊. 比重計(jì)法測定土壤顆粒組成的研究[J]. 環(huán)境科學(xué)導(dǎo)刊，2010，29(4)：97－99.

Chen Liqiong. Research on structure of soil particle by hydrometer method[J]. Environmental Science Survey, 2010, 29(4): 97－99. (in Chinese with English abstract)

[25] Chen Y, Munkholm L J, Nyord T. A discrete element model for soil-sweep interaction in three different soils[J]. Soil & Tillage Research, 2013, 126: 34－41.

[26] 賈旭光，陳曦，李鑫. 不同粒度和堆載形態(tài)下散體瞬時(shí)自然安息角的實(shí)驗(yàn)研究[J]. 現(xiàn)代礦業(yè)，2015，549(1)：25－27.

Jia Xuguang, Chen Xi, Li Xin. Experiment research on instantaneous natural repose angle of granular slope under different granularity and stack form[J]. Modern Mining, 2015, 549(1): 25－27. (in Chinese with English abstract)

[27] 馮俊小，林佳，李十中，等. 秸稈固態(tài)發(fā)酵回轉(zhuǎn)筒內(nèi)顆粒混合狀態(tài)離散元參數(shù)標(biāo)定[J]. 農(nóng)業(yè)機(jī)械學(xué)報(bào)，2015，46(3)：208－213.

Feng Junxiao, Lin Jia, Li Shizhong, et al. Calibration of discrete element parameters of particle in rotary solid state fermenters[J]. Transactions of the Chinese Society for Agricultural Machinery, 2015, 46(3): 208－213. (in Chinese with English abstract)

[28] 朱震和. 表面活性劑－粉體和顆粒：休止角的測量[J]. 日用化學(xué)品科學(xué)，1981(4)：66－68.

[29] 葛藤，賈智宏，周克棟. 計(jì)算點(diǎn)接觸碰撞恢復(fù)系數(shù)的一種理論模型[J]. 機(jī)械設(shè)計(jì)與研究，2007(3)：14－15，22.

Ge Teng, Jia Zhihong, Zhou Kedong. A Theoretical model for the coefficient of restitution calculation of point impact[J].Machine Design & Research, 2007(3): 14－15, 22. (in Chinese with English abstract)

[30] Fr?czek J, Z?obecki A, Zemanek J. Assessment of angle of repose of granular plant material using computer image analysis[J]. Journal of Food Engineering, 2007, 83(1): 17－22.

Parameter calibration of discrete element simulation model for latosol particles in hot areas of Hainan Province

Xing Jiejie1, Zhang Rui1, Wu Peng1, Zhang Xirui1※, Dong Xuehu2, Chen Ying3, Ru Shaofeng1

(1.570228; 2.,,, 524091,; 3.,,322,)

Latosol is the zonal soil of tropical rainforest or seasonal rainforest in the southernmost part of China. Compared with ordinary soil, latosol has stronger expansibility, flow plasticity and water holding capacity. In the cultivation of latosol, there are many problems, such as high resistance, low efficiency, and serious soil adhesion of parts touching the soil. In recent years, with the continuous improvement of agricultural mechanization level, discrete element method is gradually applied in the field of agricultural engineering. But the existing research is mainly aimed at the soil with weak viscosity, such as sandy soil, dry soil or rock soil. These research results are not suitable for the latosol, and there are few researches on the parameter calibration of the discrete element simulation model of latosol. In this paper, the simulation contrast test of soil repose angle and soil breaking resistance, and calibrate the contact parameters of the discrete element model of latosol had been carried out with "Hertz-Mindlin with JKR" contact model in the EDEM software. The actual repose angle of latosol with moisture content of 7.8% was measured to be 32.35 °. The recovery coefficient, static friction coefficient and rolling friction coefficient between latosol, latosol and contact parts were measured by design test. Seven coefficients were selected as test factors by Design-expert software and factor level was set. The repose angle was taken as test index and design parameter calibration test was carried out. Using Plackett-Burman Design test to determine the key factors affecting the repose angle of latosol includerecovery coefficient of latosol-latosol, rolling friction coefficient of latosol-latosol, rolling friction coefficient of latosol-28MnB5 plate and surface energy. Using regressive analysis by Box-Behnken to fit the test results, and the repose angle regression model was obtained. The regression variance analysis of the model showed that the regression model was very significant, and the repose angle could be predicted according to the model.Taking the actual repose angle of latosol as the objective value, the regression model was optimized, and the optimal values of four significant factors were obtained: latosol-latosol recovery coefficient 0.40, latosol-latosol rolling friction coefficient 0.07, latosol-28MnB5 plate rolling friction coefficient 0.24, JKR surface energy 4.31 J/m2. The simulation experiment of the repose angle of EDEM was designed, and the image of the simulated repose angle was processed by python. The relative error between the simulated repose angle and the measured repose angle was 0.62%.The soil tank model was established to simulate the movement process of the breaking blade in the latosol, and the test was carried out to compare the fitting situation of the change of the breaking resistance. Comparion of soil breaking fitting resistance curve in numeral simulation and test, showed that the calibration of the parameters of the discrete element simulation model of latosol was accurate and reliable.

latosol; discrete element model; repose angle; parameter calibration; soil breaking resistance

2019-11-16

2020-03-08

國家自然科學(xué)基金（51865007）

邢潔潔，講師，主要從事熱帶智能農(nóng)業(yè)機(jī)械研究。Email：584731137@qq.com

張喜瑞，教授，主要從事熱帶智能農(nóng)業(yè)機(jī)械研究。Email：zhangxirui_999@sina.com

10.11975/j.issn.1002-6819.2020.05.018

S220.1

A

1002-6819(2020)-05-0158-09

邢潔潔，張 銳，吳 鵬，張喜瑞，董學(xué)虎，陳 瑛，汝紹鋒. 海南熱區(qū)磚紅壤顆粒離散元仿真模型參數(shù)標(biāo)定[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2020，36(5)：158－166. doi：10.11975/j.issn.1002-6819.2020.05.018 http://www.tcsae.org

Xing Jiejie, Zhang Rui, Wu Peng, Zhang Xirui, Dong Xuehu, Chen Ying, Ru Shaofeng. Parameter calibration of discrete element simulation model for latosol particles in hot areas of Hainan Province[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2020, 36(5): 158－166. (in Chinese with English abstract) doi：10.11975/j.issn.1002-6819.2020.05.018 http://www.tcsae.org