廣西濱海地區(qū)甘蔗耕種土壤離散元仿真模型的參數(shù)標(biāo)定

陶利民，石艷琪，黃柳益，林圣強(qiáng)，劉友建，李霄楠

（北部灣大學(xué) 機(jī)械與船舶海洋工程學(xué)院，廣西 欽州 535011）

廣西濱海地區(qū)日照充足、雨水豐盈，宜種甘蔗，海岸帶土壤資源具有很好的農(nóng)業(yè)開發(fā)應(yīng)用前景。甘蔗屬主要農(nóng)業(yè)經(jīng)濟(jì)作物，研發(fā)適用性強(qiáng)、更具推廣潛力的甘蔗種植機(jī)械，可促進(jìn)糖料蔗生產(chǎn)、加快臨海大工業(yè)發(fā)展進(jìn)程。開溝、覆土作為甘蔗精準(zhǔn)播種的重要環(huán)節(jié)，面臨著耕種土壤回流、黏附嚴(yán)重、作業(yè)阻力大等難題。為更切合實(shí)際的種植環(huán)境，針對(duì)廣西濱海地區(qū)甘蔗耕種土壤的特性，應(yīng)用離散元法對(duì)土壤離散元仿真模型的參數(shù)進(jìn)行了標(biāo)定［1］。

離散元仿真模型參數(shù)的精確度取決于顆粒建模參數(shù)的準(zhǔn)確性，主要分為本征物理參數(shù)和材料接觸參數(shù)。離散元法可直觀地顯現(xiàn)土壤的力學(xué)特性和動(dòng)態(tài)行為、可用于仿真分析散粒體為主的農(nóng)業(yè)物料與耕種機(jī)械間的相互作用關(guān)系［2］。國內(nèi)外專家基于離散元法對(duì)不同地區(qū)土壤的特性進(jìn)行了深入研究，為數(shù)值模擬仿真和耦合探究奠定了前期基礎(chǔ)。

Ucgul等［3］整合了non-linear elastic與hysteretic spring接觸模型，考慮土壤塑性形變對(duì)有、無黏結(jié)力的沙灘沙增大粒徑后進(jìn)行參數(shù)標(biāo)定。Chen等［4］分析粒徑對(duì)土壤流動(dòng)性的影響過程，對(duì)沙壤土模型進(jìn)行有效驗(yàn)證。邢潔潔等［5］針對(duì)海南地區(qū)磚紅土壤離散元模型，通過仿真和土槽試驗(yàn)驗(yàn)證破土阻力的一致性，選用Hertz-Mindlin with JKR接觸模型，確保參數(shù)可靠。石林榕等［6］選取6種含水率的西北旱區(qū)土壤驗(yàn)證延遲彈性模型Hysteretic Spring Contact Mode（HSCM）和線性內(nèi)聚力模型Liner Cohesion Mode（LCM）的可行性，以入土阻力為指標(biāo)進(jìn)行鴨嘴插入土壤和仿真試驗(yàn)驗(yàn)證了參數(shù)有效性。向偉等［7］對(duì)南方黏壤土參數(shù)進(jìn)行了標(biāo)定，成穴和土槽對(duì)比試驗(yàn)的結(jié)果表明標(biāo)定準(zhǔn)確可行。

鑒于土壤分布地區(qū)的差異性、土壤本身結(jié)構(gòu)的復(fù)雜性，上述研究并不適宜本區(qū)域甘蔗種植的開溝覆土作業(yè)所需。本文以廣西濱海區(qū)域的甘蔗耕種土壤為研究對(duì)象，選取“Hertz-Mindlin with JKR”接觸模型，采用實(shí)測(cè)、仿真方法對(duì)相關(guān)接觸參數(shù)進(jìn)行標(biāo)定，以便合理地設(shè)計(jì)開溝覆土裝置結(jié)構(gòu)，使各項(xiàng)性能滿足實(shí)際工況的要求；為分析土壤動(dòng)態(tài)擾動(dòng)行為和覆土結(jié)果提供基礎(chǔ)數(shù)據(jù)。依據(jù)土壤性質(zhì)、接觸特性，對(duì)土壤顆粒間、土壤與開溝、覆土裝置間進(jìn)行準(zhǔn)確的參數(shù)標(biāo)定顯得十分必要。

1 材料與方法

采集土壤樣本，測(cè)量其密度、含水率，并標(biāo)定土壤物理參數(shù)。隨后依據(jù)土壤質(zhì)地選擇土壤顆粒的接觸模型，并測(cè)量土壤堆積角，為后續(xù)參數(shù)標(biāo)定提供數(shù)據(jù)參考［8-9］。

1.1 材料基本參數(shù)

試驗(yàn)土壤取自廣西欽州市欽南區(qū)久隆鎮(zhèn)甘蔗地（108°46′11″N，22°4′56″E），在所選取5個(gè)測(cè)試點(diǎn)取得0~30 cm土壤層間的土壤基本數(shù)據(jù)。用烘干法將被測(cè)土壤內(nèi)置于實(shí)驗(yàn)室恒溫烘箱105 ℃條件下烘至恒重，測(cè)得含水率16.3%；采用環(huán)刀（100 cm3）與電子天秤（精度0.01 g）測(cè)得密度為1790 kg/m3。土壤和Q235號(hào)鋼的其他參數(shù)參照文獻(xiàn)［10-11］選?。ū?）。

表1 離散元模型仿真參數(shù)表

1.2 土壤顆粒間接觸模型選取

廣西濱海地區(qū)土壤以紅壤、赤紅壤和磚紅壤為主［8］。該地區(qū)甘蔗耕種的土壤土質(zhì)疏松、透氣性好、根系易吸收水分。由于不同類型土壤之間的黏結(jié)性、黏附性和土壤塑性的不同，對(duì)觸土部件產(chǎn)生的影響也不同［9］。故接觸模型的選擇應(yīng)考慮具有較高含水率、黏結(jié)性較強(qiáng)的土壤對(duì)耕種質(zhì)量的動(dòng)態(tài)影響。Hertz-Mindlin with Bonding黏結(jié)接觸模型主要針對(duì)黏性土，通常用來模擬破碎工況，Hysteretic Spring延遲彈性模型描述了粒子的塑性變形，應(yīng)用于可壓縮材料。Hertz-Mindlin with JKR Cohension模型賦予土壤顆粒間的吸引凝聚力，考慮了濕顆粒之間的黏結(jié)力，故本研究選用其作為土壤離散元仿真的接觸力學(xué)模型，更貼合仿真效果。

1.3 土壤堆積角物理試驗(yàn)

土壤堆積角試驗(yàn)用于評(píng)估松散土的流動(dòng)性和顆粒間摩擦。取含水率為16.3%的土壤，搭建實(shí)測(cè)試驗(yàn)平臺(tái)如圖1所示。將漏斗（進(jìn)料口徑170 mm、漏斗出料口徑17 mm）固定于鐵架臺(tái)一側(cè)，為取得良好的堆積效果，多次預(yù)試驗(yàn)后，取漏斗出料口徑與接料臺(tái)間距為300 mm。開始試驗(yàn)之前，塞住漏斗口，放入被測(cè)土壤后，再打開漏斗口使其自然順暢流出，待土壤靜置后，從試驗(yàn)臺(tái)不同方位平行角度拍攝靜態(tài)堆積角，取其平均值。

圖1 實(shí)測(cè)試驗(yàn)平臺(tái)

在堆積角的測(cè)量試驗(yàn)中，土壤堆積邊緣參差不一，很難使用常規(guī)角度軟件實(shí)現(xiàn)精準(zhǔn)測(cè)量。本次測(cè)量應(yīng)用Matlab軟件進(jìn)行圖像邊緣檢測(cè)處理，后利用高斯分布擬合來測(cè)量堆積角［12］。具體步驟為：先將原始圖像轉(zhuǎn)換成灰度圖；并通過Otsu大津算法，將圖像的興趣區(qū)域進(jìn)行有效分割，以此作為后續(xù)特征提取的目標(biāo)對(duì)象；隨后進(jìn)行邊緣檢測(cè)，選取了比Canny算子效果更好、更符合試驗(yàn)要求的Zerocross算子；最后，在Origin軟件中對(duì)導(dǎo)入的加深輪廓圖像進(jìn)行相對(duì)應(yīng)輪廓點(diǎn)的提取，采用Gauss分布對(duì)圖像輪廓進(jìn)行擬合并取拐點(diǎn)處的切線與X軸的夾角為土壤堆積角。堆積角的測(cè)量均值為35.8°，過程如圖2所示。

圖2 土壤堆積角測(cè)量過程

2 離散元參數(shù)標(biāo)定

2.1 土壤接觸參數(shù)標(biāo)定

將上述所測(cè)的密度、堆積角和仿真規(guī)模輸入GEMM（顆粒材料模型數(shù)據(jù)庫）中，以堆積角為試驗(yàn)指標(biāo)，進(jìn)行EDEM仿真標(biāo)定試驗(yàn)。并以實(shí)測(cè)堆積角度為標(biāo)準(zhǔn)進(jìn)行優(yōu)化，從而確定土壤間的最優(yōu)參數(shù)組合。

2.1.1 參數(shù)標(biāo)定試驗(yàn) 在UG中建立堆積角仿真試驗(yàn)平臺(tái)后導(dǎo)入EDEM中，因顆粒工廠通常被顆粒粒徑限制，為合理、高效地生成所需土壤顆粒，本文在仿真中設(shè)置半徑10 mm的默認(rèn)球形顆粒為供試土壤進(jìn)行仿真［13］。仿真試驗(yàn)平臺(tái)如圖3所示。將仿真規(guī)模、顆粒密度和堆積角等輸入顆粒材料模型數(shù)據(jù)庫GEMM（顆粒材料數(shù)據(jù)庫）中。當(dāng)JKR表面能太高、恢復(fù)系數(shù)太小、靜摩擦系數(shù)太大時(shí)，土壤易因黏重而塌陷偏移［14］，因此，參數(shù)范圍取值：恢復(fù)系數(shù)為0.35~0.75；靜摩擦系數(shù)為0.32~0.80；滾動(dòng)摩擦系數(shù)為0.05~0.15；JKR表面能為4~15 J/m2。依據(jù)上述參考值，應(yīng)用Design-Expert軟件的Box-Behnken優(yōu)化方法設(shè)計(jì)4因素3水平共29組參數(shù)組合，選取堆積角為評(píng)價(jià)指標(biāo)進(jìn)行標(biāo)定試驗(yàn)，利用CAD測(cè)量堆積角的角度。試驗(yàn)因素及水平如表2所示，試驗(yàn)結(jié)果見表3。

表2 土壤堆積角仿真試驗(yàn)因素及水平

圖3 仿真試驗(yàn)平臺(tái)

對(duì)表3中數(shù)據(jù)進(jìn)行多元回歸擬合分析，建立堆積角和變量x1、x2、x3、x4的二階回歸模型：

表3 土壤堆積角仿真試驗(yàn)設(shè)計(jì)及結(jié)果

回歸方程（1）的決定系數(shù)R2=0.9979，校正決定系數(shù)Adj R2=0.9957，均接近于1，說明回歸方程的擬合程度極佳，根據(jù)模型的方差分析結(jié)果可知（表4）：模型P＜0.0001，表明該回歸模型極顯著；失擬項(xiàng)P＞0.05，失擬不顯著顯示該模型擬合程度高且可靠。在4個(gè)因素給定范圍內(nèi)：JKR表面能x1、土壤恢復(fù)系數(shù)x2、土壤滾動(dòng)摩擦系數(shù)x4對(duì)堆積角的影響極其顯著，土壤靜摩擦系數(shù)x3較為顯著；在一次交互作用項(xiàng)中，x1x2、x2x4的P值＜0.05，表明此交互項(xiàng)對(duì)堆積角的影響尤為顯著，x1x3、x1x4略微顯著；試驗(yàn)因素的二次方項(xiàng)中，x22對(duì)堆積角的影響很顯著，x32與x42均顯著，x12不顯著。

表4 堆積角回歸模型的方差分析

2.1.2 參數(shù)優(yōu)化及驗(yàn)證 應(yīng)用Design-Expert軟件中的優(yōu)化模塊Numerical，以含水率為16.3%的土壤堆積角值35.8°為目標(biāo)對(duì)堆積角回歸模型進(jìn)行尋優(yōu)，得到最優(yōu)解：土壤顆粒之間的JKR表面能8.89 J/m2，恢復(fù)系數(shù)0.51、靜摩擦因數(shù)0.72、滾動(dòng)摩擦因數(shù)0.08。該最優(yōu)解下的堆積角仿真值為36.06°，與實(shí)測(cè)值（35.8°）相對(duì)誤差為0.73%。仿真試驗(yàn)與實(shí)測(cè)試驗(yàn)的形態(tài)對(duì)比如圖4所示。

圖4 土壤堆積角物理與仿真試驗(yàn)堆型

2.2 土壤與Q235號(hào)鋼接觸參數(shù)測(cè)定

對(duì)于土壤與觸土部件材料Q235號(hào)鋼間的接觸參數(shù)測(cè)定：基于斜面滑動(dòng)法原理進(jìn)行不同含水率下靜摩擦系數(shù)的測(cè)定與分析，并采用斜板滾動(dòng)試驗(yàn)對(duì)滾動(dòng)摩擦系數(shù)進(jìn)行測(cè)定。

2.2.1 靜摩擦系數(shù)測(cè)定驗(yàn)證 土壤與Q235號(hào)鋼之間靜摩擦系數(shù)的測(cè)定選用斜面滑動(dòng)法［15］，因廣西濱海土壤質(zhì)地黏重、易吸濕，相比旱土、沙壤土等的顆粒特性更繁雜，其靜摩擦特性也會(huì)隨含水量發(fā)生變動(dòng)。為探究其關(guān)系，遂配置不同含水率（10.15%、10.62%、12.03%、13.01%、14.60%、16.30%）的1 cm3的土塊進(jìn)行試驗(yàn)。含水率為16.3%的土壤與Q235鋼靜摩擦系數(shù)μ1為0.453。

基于不同含水率下的土壤靜摩擦系數(shù)變化圖可知（圖5）：含水率在10.0%~16.5%范圍變化時(shí)，靜摩擦系數(shù)在范圍在0.40~0.46之間。隨著含水率的增大，土壤隨黏著力的出現(xiàn)，摩擦系數(shù)也隨傾角的增大而增大。當(dāng)超過一定程度后，隨著含水率的增加，靜摩擦便會(huì)減少產(chǎn)生黏附現(xiàn)象，試驗(yàn)中未達(dá)到臨界值［16］。結(jié)合毗鄰北部灣的廣西濱海地區(qū)甘蔗耕種土壤環(huán)境現(xiàn)狀，為降低濕潤多雨氣候?qū)μ镩g耕種環(huán)境的影響，應(yīng)避免耕種裝置受土壤摩擦阻力過大而導(dǎo)致作業(yè)質(zhì)量變差。

圖5 不同含水率下的靜摩擦系數(shù)與傾角的變化

2.2.2 滾動(dòng)摩擦因素測(cè)定 土壤接觸Q235鋼板表面并相對(duì)其做純滾動(dòng)運(yùn)動(dòng)或產(chǎn)生純滾動(dòng)趨勢(shì)時(shí)，即產(chǎn)生滾動(dòng)摩擦力。參照斜板滾動(dòng)試驗(yàn)進(jìn)行滾動(dòng)摩擦系數(shù)的測(cè)定［17］，搭建的斜板實(shí)測(cè)試驗(yàn)臺(tái)由斜板和數(shù)顯液晶傾角儀、測(cè)量直尺3部分組成，如圖6a。將含水率16.3%的5 mm小球放置第一塊斜板中央起始位置，以初速度為0 m/s沿斜面自由平穩(wěn)滑落至第二塊斜板直至靜止，測(cè)量其水平滾動(dòng)距離S，重復(fù)10次操作并記錄，采用式（2）計(jì)算，得出S滾動(dòng)距離的平均值21.3 mm，滾動(dòng)摩擦系數(shù)μ2為0.158。

式（2）中：H2為小球下落前與水平面距離（mm）；L為小球在斜面的滾動(dòng)距離（mm）；θ2為斜面鋼板與水平面鋼板的夾角（°）；S為小球在水平面鋼板滾動(dòng)的距離（mm）。

仿真試驗(yàn)臺(tái)如圖6b所示，為消除土壤之間未產(chǎn)生接觸對(duì)試驗(yàn)結(jié)果造成的誤差影響，將土壤間接觸系數(shù)均設(shè)置為0。在預(yù)仿真試驗(yàn)結(jié)果的基礎(chǔ)上，保證仿真與實(shí)測(cè)試驗(yàn)臺(tái)的關(guān)鍵參數(shù)相統(tǒng)一，取彈性恢復(fù)系數(shù)為0.5，測(cè)定此條件下的仿真滾動(dòng)距離，共進(jìn)行10組試驗(yàn)，每組重復(fù)3次后取均值，利用EDEM自帶的Ruler工具測(cè)定滾動(dòng)距離，繪成折線圖，由圖7可知：仿真距離與實(shí)測(cè)距離的變化趨勢(shì)基本相同，根據(jù)土壤滾動(dòng)距離y和滾動(dòng)摩擦因數(shù)x的關(guān)系，對(duì)仿真值進(jìn)行曲線擬合得方程：

圖6 滾動(dòng)摩擦實(shí)測(cè)和仿真試驗(yàn)臺(tái)

圖7 不同滾動(dòng)距離下滾動(dòng)摩擦因數(shù)的實(shí)測(cè)值與仿真值

式（3）的決定系數(shù)R2=0.9604，接近于1，證明其擬合的可靠性和準(zhǔn)確度高。實(shí)測(cè)滾動(dòng)距離均值y=21.3 mm代入二次多項(xiàng)式擬合曲線方程（3）得x為20.68 mm，仿真中驗(yàn)證與實(shí)值的相對(duì)誤差為2.91%，驗(yàn)證了仿真標(biāo)定的準(zhǔn)確性，遂取滾動(dòng)摩擦因數(shù)μ2為0.158。

3 結(jié)果與分析

由于不同土層含水率存在差異性，土壤的平均含水率最大不超過30%［18］，使用上述所測(cè)得的土壤參數(shù)，構(gòu)建不同含水率（1.0%、10.0%、16.3%、30.0%）的土層。在0.3 m/s的前進(jìn)速度下，通過EDEM軟件進(jìn)行開溝器仿真試驗(yàn)，探究不同含水率下工作阻力的變化。由表5可知，開溝器在4種不同含水率的土壤下工作阻力范圍為653.63~ 1387.12 N，當(dāng)土壤含水率由1.0%增大至16.3%時(shí)，開溝器破土力隨土壤的彈性模量與剪切模量的下降而減?。?9-20］；而由16.3%增大至30.0%時(shí)，土壤顆粒之間的黏附力增大，結(jié)塊的大顆粒對(duì)開溝器造成沖擊，開溝工作阻力也隨之增大。

表5 不同含水率的開溝工作阻力

從微觀層面上分析表明：應(yīng)用EDEM中的Clipping功能添加剪裁平面，后勾選show legend顯示法向力的梯度表。由土壤顆粒受力云圖可知（圖8）：圖8a中開溝器前端紅色顆粒分布密集，表明其受到的阻力最大；由圖8b和圖8c得知，隨著含水率由1.0%增大至16.3%的過程中，紅色顆粒明顯減少，表明隨著含水率的增大工作阻力減小至適宜值。但伴隨含水率從16.3%增大至30.0%，增多的紅色顆粒反映出工作阻力因土壤黏性增強(qiáng)又持續(xù)增大（圖8d）。因此，在種植甘蔗時(shí)，應(yīng)盡量避免干旱或雨天時(shí)節(jié)，選擇含水率適中的土壤進(jìn)行耕種，使開溝覆土等作業(yè)滿足工況需求。

圖8 不同含水率下土壤顆粒受力情況

4 結(jié)論

確定廣西濱海地區(qū)甘蔗耕種土壤的本證參數(shù)及接觸模型，借助Matlab圖像處理技術(shù)，選擇Zerocross算子對(duì)土壤輪廓邊緣進(jìn)行處理，得土壤堆積角測(cè)量值為35.8°

通過4因素3水平正交組合試驗(yàn)，對(duì)土壤間的接觸參數(shù)進(jìn)行標(biāo)定，以實(shí)測(cè)堆積角的角度35.8°為目標(biāo)尋優(yōu)得土壤間接觸參數(shù)最佳組合：土壤—土壤恢復(fù)系數(shù)、土壤—土壤靜摩擦系數(shù)、土壤—土壤滾動(dòng)摩擦系數(shù)、表面能參數(shù)分別為：0.51、0.72、0.08、8.89 J/m2。通過實(shí)際試驗(yàn)和模擬仿真相結(jié)合的方法得出土壤與Q235號(hào)鋼的接觸參數(shù)：靜摩擦系數(shù)0.453、滾動(dòng)摩擦系數(shù)0.158。

通過開溝器仿真試驗(yàn)，驗(yàn)證了適宜含水率土壤對(duì)開溝作業(yè)工作阻力的影響。為減少土壤摩擦阻力在設(shè)計(jì)耕種部件時(shí)應(yīng)注重土壤黏附問題，土壤含水量較大時(shí)進(jìn)行松土以改善土壤透氣性。本文所測(cè)得的土壤模型參數(shù)在后續(xù)基于離散元法的甘蔗種植機(jī)開溝覆土作業(yè)過程仿真和裝置的設(shè)計(jì)優(yōu)化發(fā)揮可觀后效。