基于EDEM-Fluent耦合的氣爆松土效果仿真分析

陳金楚，奚小波，殷慧子，張翼夫，張寶峰，金亦富，張瑞宏

（揚(yáng)州大學(xué)機(jī)械工程學(xué)院，江蘇 揚(yáng)州 225127）

果園土壤合理有效深松有利于疏通土壤之間的間隙，降低水土流失，利于復(fù)合施肥作業(yè)，為肥料擴(kuò)散提供空間，提高農(nóng)作物產(chǎn)量[1-3]。傳統(tǒng)果園松土施肥作業(yè)是機(jī)器或人工開溝松土、施肥、覆土[4]，易出現(xiàn)傷害作物根系、破壞土壤結(jié)構(gòu)等情況[5]。氣爆松土技術(shù)不破壞土層結(jié)構(gòu)，通過輸送氣體對(duì)給定深度的土壤進(jìn)行爆氣疏松，并在土層中形成空隙，利于后續(xù)液肥擴(kuò)散[6]。

近年來，基于土體氣壓劈裂原理，孔德剛等提出氣壓式土壤深松方法[7]；增子利一等通過向果樹根部注入壓縮空氣實(shí)現(xiàn)快速松土和深層施肥[8]；張瑞宏等設(shè)計(jì)適用于液肥注施的氣動(dòng)深松施肥機(jī)[9-10]；奚小波等優(yōu)化氣爆松土參數(shù)[11-12]。目前氣爆松土一般通過試驗(yàn)分析指標(biāo)參數(shù)對(duì)土壤的大體擾動(dòng)效果，對(duì)土壤顆粒的具體擾動(dòng)效果缺乏理論研究。土壤顆粒運(yùn)動(dòng)的理論研究大多采用離散元方法[13]（Discrete element method,DEM），僅考慮土壤顆粒間產(chǎn)生的接觸力、摩擦力，忽略實(shí)際土壤間粘聚力和塑性變形[14-15]。Walton等針對(duì)塑性材料創(chuàng)建遲滯彈性接觸模型（Hysteretic spring contact model,HSCM）[16]，而后Janda等驗(yàn)證其準(zhǔn)確性[17-18]。Ucgul等將遲滯彈性接觸模型（HSCM）和法向粘聚模型（LCM）結(jié)合表征現(xiàn)實(shí)土壤顆粒之間的塑性變形和粘聚力[19-20]，使法向力隨顆粒間迭代量呈理想線性關(guān)系：當(dāng)外部施加力小于設(shè)定應(yīng)力值時(shí)，執(zhí)行彈性變形；當(dāng)外力超過設(shè)定應(yīng)力值時(shí)，發(fā)生塑性變形。考慮到土壤擠壓變形屬于非線性形變，Thankur等創(chuàng)建愛丁堡彈塑性附著力模型（Edinburgh elastoplastic adhesion model，EEPA）[21-22]，學(xué)者們根據(jù)土壤性質(zhì)將土壤結(jié)構(gòu)分為耕作層、犁底層和心土層，不斷優(yōu)化試驗(yàn)，驗(yàn)證EEPA模型合理性。

本文采用EDEM-Fluent耦合進(jìn)行氣爆松土仿真，研究不同爆氣流孔結(jié)構(gòu)及氣爆參數(shù)對(duì)不同深度土壤擾動(dòng)效果影響，擬合建立土壤擾動(dòng)軌跡線方程，通過試驗(yàn)驗(yàn)證仿真結(jié)果，為氣爆松土工藝參數(shù)制定及機(jī)具結(jié)構(gòu)設(shè)計(jì)提供理論依據(jù)。

1 氣爆松土總體結(jié)構(gòu)參數(shù)設(shè)計(jì)

如圖1所示，氣爆松土工作部件主要由進(jìn)氣接頭、氣缸和鉆桿3部分構(gòu)成。氣鏟標(biāo)準(zhǔn)驗(yàn)收氣壓為機(jī)械行業(yè)標(biāo)準(zhǔn)JB/T 8412-2016[23]規(guī)定0.63 MPa，氣鏟尾部鉆桿結(jié)構(gòu)重量為2.43 kg。根據(jù)標(biāo)準(zhǔn)JB/T 8412-2016中氣鏟各類型技術(shù)規(guī)格，考慮到鉆桿作業(yè)行程，選擇沖擊行程有優(yōu)勢(shì)、沖擊頻率低、運(yùn)行平穩(wěn)的C7型氣鏟。具體參數(shù)如下：機(jī)重7.4 kg，系統(tǒng)總重9.83 kg，機(jī)長445 mm，缸徑28 mm，沖擊頻率13 Hz，沖擊能17 J，耗氣量16 L·s-1，沖擊行程100.04 mm，氣管內(nèi)徑13 mm，工作氣壓0.63 MPa。

圖1 氣爆松土機(jī)整體結(jié)構(gòu)Fig.1 Overall structure of air burst soil loosening machine

考慮鉆桿松土?xí)r順暢程度和結(jié)構(gòu)剛度，盡可能使鉆桿外徑趨于最小值，以減少作業(yè)時(shí)在軸向上產(chǎn)生的合阻力，優(yōu)化后鉆桿壁厚5 mm、內(nèi)徑15 mm、錐角2a取60°。在流孔部分，孔心距離鉆桿平面l取20 mm，鉆桿長度L為950 mm[12]。流體分布器孔口結(jié)構(gòu)是影響液體分布器分布均勻性的因素之一[31]，本文在保持流孔體積不變情況下，根據(jù)截面形狀變化設(shè)計(jì)3種爆氣流孔結(jié)構(gòu)桿，其流孔數(shù)量均為4，如圖2所示。由圖2可知，3種爆氣流孔結(jié)構(gòu)桿分別為：①外擴(kuò)型，外口直徑大于內(nèi)口直徑，且成60°凹角，橫截面Φ1值為7 mm；②外縮型，外口直徑大于內(nèi)口直徑，且成60°凹角，橫截面Φ2值為7 mm；③直口型，外口直徑等于內(nèi)口直徑，橫截面Φ3值為9.9 mm。

圖2 不同流孔的鉆桿結(jié)構(gòu)參數(shù)Fig.2 Drill pipe structural parameters of different flow holes

2 離散元土壤模型建立

2.1 土壤顆粒模型

土壤顆粒組成通常為單顆粒核狀、雙顆粒條狀、多顆粒（≥3）團(tuán)簇狀等形式，如圖3所示。

圖3 土壤顆粒形狀Fig.3 Soil particle shape

在EDEM軟件中，基礎(chǔ)球顆粒數(shù)目越多，仿真時(shí)間越長，為減少顆粒數(shù)量，保證土壤體積，選擇相對(duì)較大粒徑的顆粒，考慮到果園土壤常年不翻耕，土層結(jié)構(gòu)緊實(shí)，故采用半徑3.5 mm、形式單一團(tuán)簇顆粒形狀，模擬粗粒組中粒徑較小的圓礫和角礫。犁底層顆粒間滾動(dòng)摩擦系數(shù)為0.21，靜摩擦系數(shù)為0.35[24]。顆粒仿真參數(shù)：密度2 600 kg·m-3，彈性模量1×106Pa，泊松比0.3，恢復(fù)系數(shù)0.6[25]。

2.2 土壤接觸模型

為更好體現(xiàn)土壤力-形變關(guān)系，本文采取土壤模型為EEPA，區(qū)別于傳統(tǒng)彈性接觸模型，考慮土壤顆粒間產(chǎn)生非線性粘聚和塑性變形[26-27]，可有效描述黏附顆粒應(yīng)力變化過程。顆粒間作用力為法向力，形變?yōu)閴毫π遁d后顆粒間保留重疊量。土壤之間力關(guān)系如圖4所示。

圖4 力關(guān)系原理[28-29]Fig.4 Schematic of force relations[28-29]

當(dāng)顆粒間法向力增加時(shí)，顆粒間重疊量δ相應(yīng)增加。接觸按照加載（Loading）曲線斜率加載；當(dāng)施加力消失，重疊量相應(yīng)降低，接觸按照卸載/重新加載（Unloading/reloading）曲線的斜率卸載。顆粒在分離前繼續(xù)施加力，按照k2加載。當(dāng)法向力小于0時(shí)，且趨于峰值f max時(shí)，顆粒間粘結(jié)力隨之消失。此后，法向力按照-kabh斜率恢復(fù)到顆粒間初始接觸力f0。EEPA模型相對(duì)于HSCM+LCM模型結(jié)合優(yōu)勢(shì)體現(xiàn)在圖4負(fù)半軸上。EEPA在EDEM軟件中設(shè)置參數(shù)如下：黏附力強(qiáng)度0 N，接觸黏附能50 J·m-2，接觸塑性比0.7，非線性曲線冪指數(shù)1.5，黏附分支功率5，切向剛度系數(shù)0.28571[30]。

2.3 土壤和鉆桿仿真模型

因氣爆松土作業(yè)對(duì)象為果園林木土壤，果園田間管理機(jī)械和收獲機(jī)械作業(yè)過程土壤壓實(shí)，且常年不翻耕，造成土質(zhì)板結(jié)，因此本文仿真可將其近似于土壤層結(jié)構(gòu)中的犁底層，在EDEM中生成的土壤模型設(shè)置厚度為500 mm。通過土壤顆粒模型搭建形成虛擬土體，尺寸為1 200 mm（長）×800 mm（寬）×500 mm（高）。犁底層土壤團(tuán)聚現(xiàn)象明顯，所以尺寸分布為固定值1，顆粒工廠產(chǎn)生700 000個(gè)團(tuán)簇體顆粒。生成的土壤模型如圖5所示。

圖5 土壤模型Fig.5 Soil model

采用Solidworks三維軟件創(chuàng)建對(duì)應(yīng)鉆桿幾何模型，材料選用304不銹鋼，材料密度ρ為7850kg·m-3，剪切模量G為1×1010Pa，泊松比ν為0.25，鉆桿和土壤之間恢復(fù)系數(shù)K為0.6。同犁底層之間靜摩擦因素為0.639；滾動(dòng)摩擦因素為0.13[16]。為更好滿足軟件間耦合時(shí)間步長匹配需要的整數(shù)倍關(guān)系，將時(shí)間步長設(shè)為0.00005 s，總體時(shí)間設(shè)為1 s，使EDEM中土壤顆粒生成后充分堆積接觸，單元尺寸（Cell size）設(shè)置為3Rmin生成最終土壤顆粒模型。3類流孔鉆桿與不同深度土壤層進(jìn)行接觸。

3 耦合仿真過程與結(jié)果分析

3.1 前處理過程

為簡(jiǎn)化模型結(jié)構(gòu)，節(jié)約仿真時(shí)間，本文僅對(duì)鉆桿部分進(jìn)行相關(guān)處理分析。將鉆桿三維模型導(dǎo)入Workbench軟件中。由圖2空腔作為鉆桿內(nèi)部流體域。上端面設(shè)為inlet，下端面氣孔4個(gè)面設(shè)為outlet，鉆桿圓柱內(nèi)壁設(shè)為wall。自動(dòng)劃分網(wǎng)格，確保fluent流體域網(wǎng)格尺寸大于顆粒體積。

3.2 耦合計(jì)算過程

在Fluent軟件中，選擇壓力出入口，根據(jù)實(shí)際工況設(shè)置出入口壓力。初始排氣壓力值分別為0.4、0.6、0.8 MPa，容積流量為16.667 L·s-1。不計(jì)工作過程中磨損等各因素。根據(jù)空壓機(jī)壓降公式：

式中，Δp為壓力降（Pa）；q為容積流量（L·s-1）；L為管子長度（m）；d為管道內(nèi)徑（mm）；p為排氣壓力（Pa）。最終計(jì)算得出口壓力值分別為397 436.46、

598 290.97、798 718.23 Pa。

因氣體與顆粒間作用僅發(fā)生在鉆桿流體域管道外部，只需對(duì)Fluent產(chǎn)生氣體與EDEM生成顆粒之間進(jìn)行單向計(jì)算。耦合是基于不考慮體積分?jǐn)?shù)影響的多相流耦合接口。在耦合時(shí)，接口以串行方式打開（并行時(shí)求解器進(jìn)程設(shè)為0）。將動(dòng)量交換源項(xiàng)的數(shù)學(xué)模型編譯進(jìn)edem_udf耦合文件，通過此文件加載連接EDEM和Fluent求解器。

湍流模型選擇Standardk-ε模型，依據(jù)本文局部固相體積保持比例占總體積比例，本文選用Lagrangian法耦合?？缦鄥?shù)交換中的阻力模型、熱傳遞模型均采用默認(rèn)設(shè)置，升力模型選擇Saff?man Lift和Magnus Lift體現(xiàn)速度梯度以及顆粒自旋所產(chǎn)生的升力。

3.3 不同流孔結(jié)構(gòu)的松土效果分析

通過Fluent設(shè)置的氣壓值、與EDEM生成土壤層不同接觸深度及3種流孔結(jié)構(gòu)鉆桿在氣爆壓力0.4 MPa、松土深度100 mm條件下仿真。時(shí)間步長為0.0005 s，EDEM端的10倍，步數(shù)為2000，耦合總體時(shí)間為兩者乘積值1 s，此時(shí)Residuals值開始收斂，適于本文分析。

圖6a、b、c是3種鉆桿出口氣流速度云圖，可見外擴(kuò)型最大氣流速度約為40 m·s-1，外縮型最大氣流速度約為50 m·s-1，直口型最大氣流速度約為60 m·s-1。由于重力作用，流孔出口速度云圖分布呈向下分布趨勢(shì)，外擴(kuò)型出口速度分布效果最差，外縮型最好，直口型次之。結(jié)合速度值和分布均勻性情況，最終確定直孔型流孔對(duì)土壤顆粒擾動(dòng)效果最佳。通過對(duì)3組數(shù)據(jù)生成的土壤顆粒從Y方向遠(yuǎn)離土壤模型中心70 mm、厚度150 mm的Slice處理，得出圖6d、e、f 3種鉆桿對(duì)土壤顆粒的擾動(dòng)效果。本文通過選取土壤顆粒運(yùn)動(dòng)速度＞0.045 m·s-1部分作為土壤擾動(dòng)活躍區(qū)域，通過EDEM中的速度調(diào)節(jié)及顏色劃分并去除邊界條件觀察，可見，隨出口氣流速度增大，土壤顆粒擾動(dòng)分布區(qū)域越大，顆粒輪廓線大體為拋物線，且拋物線頂點(diǎn)與氣爆深度接近。從土壤擾動(dòng)顆粒最外緣范圍看，外擴(kuò)型鉆桿為（-370 mm，400 mm），外縮型為（-375 mm，400 mm），直口型為（-385 mm，400 mm），可進(jìn)一步確定直口型為最佳鉆桿流孔結(jié)構(gòu)，其徑向擾動(dòng)區(qū)域范圍達(dá)到785 mm。

圖6 3種鉆桿出口氣流速度和土壤顆粒擾動(dòng)效果Fig.6 Three kinds of drill pipe outlet airflow velocity and soil particle disturbance effect

3.4 不同氣爆參數(shù)松土效果分析

上述仿真結(jié)果分析表明，直口型是符合作業(yè)能效最佳結(jié)構(gòu)形式。為進(jìn)一步研究氣爆壓力P與松土深度H對(duì)土壤擾動(dòng)效果影響，對(duì)9組P-H條件下土壤擾動(dòng)仿真分析，結(jié)果如圖7所示。可看出，當(dāng)H一致時(shí)，P越大，土壤顆粒擾動(dòng)區(qū)域也越大；當(dāng)P一致時(shí)，H不超過200 mm時(shí)，土壤顆粒擾動(dòng)區(qū)域隨H增大而增大，當(dāng)H為300 mm時(shí)，隨P變化，土壤顆粒擾動(dòng)區(qū)域分布差異不明顯，這是因?yàn)樯疃仍酱笸寥辣旧碇亓吭酱螅琍在一定范圍內(nèi)對(duì)其擾動(dòng)效果不明顯。

圖7 不同氣爆壓力P與松土深度H條件下的土壤擾動(dòng)分布Fig.7 Distribution of soil disturbance under different air burst pressure Pand loose soil depth H

利用Solidworks軟件進(jìn)行截圖采點(diǎn)，確定生成軌跡中代表性顆粒點(diǎn)坐標(biāo)，再用Origin軟件點(diǎn)狀圖分析擬合多項(xiàng)式方程，如圖8所示。

圖8 土壤顆粒速度分布的軌跡線Fig.8 Trajectory line of soil particle velocity distribution

擬合出的軌跡線方程為h=ar2+br+c型，由于理想條件下軌跡應(yīng)為與h軸對(duì)稱的拋物線，方程中b為0，同時(shí)為方便看出松土深度H最大擾動(dòng)深度，可將方程簡(jiǎn)化為h=ar2-H+c′型，最終得出9組土壤顆粒擾動(dòng)區(qū)域軌跡線方程：

從各組擬合曲線對(duì)應(yīng)的a值和c′值來看，當(dāng)作業(yè)深度為100mm時(shí)，土壤擾動(dòng)深度與作業(yè)深度一致；隨松土深度增大，顆粒擾動(dòng)最大深度下移，大于松土作業(yè)深度。

圖9為土壤表面土體擾動(dòng)全域分布情況。其擾動(dòng)范圍類似橢圓，長軸Δx在785~900 mm、短軸Δy在560~620 mm范圍內(nèi)變化。當(dāng)P一致時(shí)，H越大，Δx值變大，Δy值變??；當(dāng)H一致時(shí)，P越大，產(chǎn)生同樣變化趨勢(shì)。在土壤表面擾動(dòng)面積方面，通過橢圓面積公式計(jì)算得出，0.8 MPa-100 mm組合擾動(dòng)面積最大，0.6 MPa-100 mm組合次之，0.4 MPa-100 mm組合最小?？紤]到松土深度和實(shí)際作業(yè)時(shí)最佳效果，進(jìn)一步比較300 mm深度時(shí)擾動(dòng)范圍面積大小，得出0.6 MPa-300 mm為最合適參數(shù)組合。土壤表面土體擾動(dòng)范圍值見表1。

表1 土壤表面的土體擾動(dòng)范圍值Table 1 Disturbance range value of the soil surface

圖9 土壤表面的土體擾動(dòng)全域分布Fig.9 Global distribution map of soil disturbanceon soil surface

9組仿真數(shù)據(jù)的抬升高度Δh與不同氣爆壓力-松土深度組合之間的關(guān)系見圖10，可見，土壤模型整體抬升高度與深度呈反比，與氣壓值呈正比，隨氣壓值增大，其整體抬高變化量趨于平緩。

圖10 不同氣爆壓力P與松土深度H條件下土壤模型抬升高度Fig.10 Lifting height of thesoil model under the condition of different air burst pressure Pand loose soil depth H

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

4.1 試驗(yàn)?zāi)康呐c條件

為進(jìn)一步驗(yàn)證氣爆松土實(shí)際效果并判斷仿真結(jié)果準(zhǔn)確性，制作標(biāo)準(zhǔn)試驗(yàn)土體，其土塊質(zhì)地為壤土，緊實(shí)度為3.28 MPa，含水量為16.6%，室外溫度為18℃。試驗(yàn)裝置如圖11所示，其中空壓機(jī)型號(hào)為羅翔W-1.0/8，滿載氣壓為0.8 MPa，排氣量1 m3·min-1，功率7.5 kW，氣筒容量230L，轉(zhuǎn)速2 900 r·min-1，氣壓可在0.4~0.8 MPa范圍內(nèi)調(diào)節(jié)。試驗(yàn)設(shè)備還有TJS-100型土壤緊實(shí)度測(cè)定儀（測(cè)量深度0~450 mm，精度±0.1%），TZS-I型土壤水分測(cè)定儀（含水率0~100%，相對(duì)誤差≤3%）、鋼板尺（量程50 mm，精度1 mm）、卷尺（量程5 m，精度1 mm）等各類用于測(cè)量標(biāo)記的工具。

圖11 氣爆松土試驗(yàn)裝置實(shí)物Fig.11 Picture of gas explosion subsoiling experimental equipment

4.2 試驗(yàn)方法

氣爆松土試驗(yàn)時(shí)，氣體在土壤層由內(nèi)向外擴(kuò)散，抬起土體，最終疏松的土體實(shí)現(xiàn)整體脫離，為更好觀察氣爆后土體擾動(dòng)情況，處理土體剖面：距氣爆中心點(diǎn)50 mm處標(biāo)記，用鏟子進(jìn)行人工剖面。本試驗(yàn)參數(shù)選用具有代表性的氣爆壓力0.4 MPa、松土深度100 mm和氣爆壓力0.8 MPa、松土深度300 mm，測(cè)量土壤在橫縱方向上的擾動(dòng)范圍，對(duì)比分析產(chǎn)生的裂隙采點(diǎn)標(biāo)記坐標(biāo)及仿真結(jié)果。

4.3 試驗(yàn)結(jié)果與分析

氣爆后的土體內(nèi)部產(chǎn)生裂隙，通過圖12試驗(yàn)結(jié)果與仿真結(jié)果比較發(fā)現(xiàn)，0.4 MPa-100 mm組合的Δx試驗(yàn)值為730 mm，Δy試驗(yàn)值為540 mm，與仿真值誤差均小于7.6%；0.8 MPa-300 mm組合的Δx試驗(yàn)值為820 mm，Δy試驗(yàn)值為515 mm，與仿真值誤差均小于9.8%。

圖12 氣爆松土試驗(yàn)土體截面形貌Fig.12 Crosssection morphology of soil in air burst test

圖13為氣爆松土試驗(yàn)裂隙軌跡線的擬合曲線，可以發(fā)現(xiàn)，該擬合曲線亦是拋物線，與仿真結(jié)論一致，同時(shí)0.4 MPa-100 mm組合的軌跡線c′為-5.57 mm，對(duì)應(yīng)仿真值1.678 mm；0.8 MPa-300 mm組合軌跡線的c′為-136.21 mm，對(duì)應(yīng)仿真值-125.664 mm，二者相對(duì)誤差較小，說明仿真結(jié)果合理性與準(zhǔn)確性。試驗(yàn)與仿真結(jié)果存在誤差，因?yàn)樵囼?yàn)土體是根據(jù)果園土壤的緊實(shí)度通過人工堆積制作而成，與仿真模型存在一定偏差。

圖13 氣爆松土試驗(yàn)裂隙軌跡線的擬合曲線Fig.13 Fitting curveof crack track linein air burst soil loosening test

5 結(jié)論

a.本文采用EDEM-Fluent耦合作氣爆松土仿真，研究不同氣爆流孔結(jié)構(gòu)及氣爆參數(shù)對(duì)土壤擾動(dòng)效果的影響，建立土壤擾動(dòng)軌跡線方程，并通過試驗(yàn)驗(yàn)證仿真結(jié)果。

b.直孔型流孔對(duì)土壤顆粒擾動(dòng)效果最佳，在氣爆壓力P為0.4 MPa、松土作業(yè)深度H為100 mm時(shí)最大氣流速度達(dá)60 m·s-1，徑向擾動(dòng)范圍達(dá)785 mm。

c.不同氣爆參數(shù)對(duì)土壤擾動(dòng)效果的仿真結(jié)果表明，當(dāng)H一致時(shí)，土壤擾動(dòng)區(qū)域隨P增大而增大；當(dāng)P一致時(shí)，H＜200 mm時(shí)，土壤擾動(dòng)區(qū)域隨著H增大而增大，而當(dāng)H＞200 mm時(shí)，土壤擾動(dòng)區(qū)域隨H變化不明顯。同時(shí)，H＜100 mm時(shí)，土壤擾動(dòng)最大深度與H幾乎一致；H＞100 mm時(shí)，土壤擾動(dòng)最大深度大于H。土體抬升高度與H呈反比，與P呈正比，隨著P持續(xù)增大，抬升高度趨于平緩。

d.氣爆松土試驗(yàn)結(jié)果顯示，氣爆能使土體產(chǎn)生裂隙并擴(kuò)散，其擬合曲線為拋物線，與仿真結(jié)論一致，試驗(yàn)土體擾動(dòng)范圍與仿真值的誤差較小，仿真結(jié)果合理可靠。