基于EDEM/FEM耦合算法平壟過程仿真

喬永昌，弋景剛，袁永偉

（河北農業(yè)大學 機電工程學院，河北 保定 071001）

根據(jù)農業(yè)部580個蔬菜重點縣信息監(jiān)測數(shù)據(jù)，2017年我國蔬菜累計種植面積1 313.25萬a［2］，而蔬菜總需求量在2020年將達到5.9億t［1］。溫室起壟栽培工藝具有增加土壤透氣性、提高水肥利用率、促進蔬菜根系生長、提高地表溫度等優(yōu)點［3-4］，所以溫室多以起壟種植為主。目前溫室整機機械主要以微耕機和大田旋耕機為主［5-7］，但是起壟種植模式使得地表高度差相對較大，旋耕機或微耕機在耕作時會因地表高低不平造成耕深一致性差、旋耕刀受力不均以及耕作后地表平整度差等問題，因此提高耕整土地耕深一致性、增加耕作后地表平整度以及降低旋耕刀波動受力是推廣平壟作業(yè)的主要問題。常見的平壟方式主要是用推土板將壟推平，該種方式的缺點是消耗功率大，碎土質量差，平整效果差，且一次作業(yè)不能達到土地的待播狀態(tài)。在平壟作業(yè)過程中，利用理論計算的方法很難準確分析出土壤對平壟部件的作用力，從而對平壟部件的強度計算及部件選型帶來困難。離散元法是通過建立土壤顆粒模型，模擬仿真觀察土壤顆粒運動情況，得出土壤顆粒間及土壤顆粒與工作部件之間的相互作用，從而對工作部件進行受力分析［8-9］。

本研究旨在基于EDEM/FEM耦合算法模擬平壟過程土壤對平壟部件產生的作用力及土壤粒子運動情況，對平壟部件工作參數(shù)選擇并對關鍵部件進行強度校核，為平壟裝置的設計提供可靠的依據(jù)。

1 整體結構設計

平壟作業(yè)是為了解決起壟種植土壤結構一次性旋耕作業(yè)不能達到旋耕作業(yè)耕深與作業(yè)后地表平整度要求的實際問題，平壟作業(yè)的主要形式與旋耕作業(yè)相同，都是靠刀具切削土壤實現(xiàn)松土的目的，但平壟作業(yè)與旋耕作業(yè)的區(qū)別在于平壟作業(yè)既要松土又要輸土，具體來說就是要將一部分壟上土壤輸送到壟溝處，實現(xiàn)耕作地面初步平整。平壟裝置主要結構組成如圖1所示。

圖1 平壟部件結構示意圖Fig.1 Diagram of ridge component structure

彎刀在刀軸上的布置形式不僅影響平壟作業(yè)質量、工作穩(wěn)定性和功率消耗，而且也對刀軸的受力、機具所受阻力矩等影響較大，因此在同一作業(yè)區(qū)域內布置4把彎刀按照左彎刀、右彎刀周向間隔90°排列，同側5組彎刀在平壟刀軸上按照4頭螺旋線（螺旋線頭數(shù)與同組內彎刀數(shù)量一致）布置，傳動裝置左右兩側彎刀對稱式布置，該種布置方式較單頭螺旋線方式工作時扭轉力矩小，刀軸受力更為均勻。作業(yè)時傳動裝置左右兩側對稱位置處左彎刀或右彎刀同時入土，減少刀軸所受軸向壓力，保證整機左右兩側受力均勻一致。與彎刀交叉間隔式布置的帶狀螺旋葉片可將碎土輸送至壟兩側的壟溝內，實現(xiàn)平壟作業(yè)。由圖1可以知，在四頭螺旋式彎刀與螺旋葉片組合布置的平壟裝置結構中，帶狀螺旋葉片只輸碎土而不切土，且?guī)盥菪~片與平壟刀軸之間具有一定間隙，可有效防止平壟過程的壅土現(xiàn)象。

刀軸的左右兩半軸上布置彎刀數(shù)Z1=Z/2=20把，其中Z為平壟部件彎刀總數(shù)，相鄰兩入土彎刀之間的轉角為：

根據(jù)上述計算，將刀軸上的40把彎刀按照4頭螺旋線、相鄰入土轉角18°進行布置，如圖2所示。

圖2 彎刀軸上布置平面展開圖Fig.2 Planar expansion diagram of bending cutter arrangement on axis

2 離散元仿真模擬

應用離散元EDEM 2.7軟件對彎刀與帶狀螺旋葉片組合布置的平壟部件進行耕作過程仿真模擬，既可以直觀地觀察平壟作業(yè)后土壤顆粒的運動狀況和作業(yè)質量，又可以得到平壟過程工作部件的受力狀況，為后期工作部件的受力分析提供可靠的依據(jù)。

2.1 離散元土壤接觸模型與參數(shù)確定

為了準確模擬出平壟部件的作業(yè)性能，需要建立適合的土壤接觸模型并確定合理的仿真參數(shù)。接觸模型的本質是固體顆粒在準靜態(tài)條件下的關于接觸力學的彈塑性分析結果，不同的接觸模型對計算固體顆粒的力與力矩產生很大影響，因此，為了提高模擬結果的可靠性，針對不同作業(yè)對象，需要搭建正確的接觸模型［10］。在離散元法顆粒運動過程中，顆粒i受到重力mig和接觸力Fi兩者的合力。其中，接觸力Fi包括顆粒之間或平壟部件與顆粒之間的法向碰撞接觸力Fn,ij、切向碰撞力Fτ,ij、法向阻尼、切向阻尼和法向結合力Fcon,ij［11］。

本實驗組作業(yè)區(qū)在河北省保定市，該地區(qū)土壤屬于輕壤土（顆粒表面的粘附力小，具有一定的壓縮性），因此土壤間接觸模型為Hysteretic Spring和Linear Cohesion，其中 Hysteretic Spring 接觸模型能夠很好地模擬土壤顆粒塑性變形，而Linear Cohesion接觸模型可以很好地模擬法向粘聚力［12-13］。通過查閱相關文獻［14-17］和實際測量得到主要仿真參數(shù)見表1。

表1 仿真參數(shù)表Table 1 Simulation parameter table

2.2 EDEM平壟仿真試驗設計

平壟部件左右兩側對稱式布置，其工作狀況基本一致，為了減少仿真軟件計算量與和縮短運行時間，現(xiàn)對平壟部件中一側刀輥進行仿真模擬，根據(jù)上一節(jié)對壟結構尺寸的分析，利用EDEM軟件分別建立壟層和壟底層的土壤模型，由于采用單刀輥仿真作業(yè)，故建立半條壟，壟層尺寸為（壟高×壟頂長×壟底長×壟長）150 mm×400 mm×500 mm×5 300 mm（土壤顆粒數(shù)為56 800），為方便觀察土壤粒子的運動情況，將壟底層適當加寬，則壟底層的尺寸為（長×寬×高）5 300 mm×1 200 mm×60 mm（土壤顆粒數(shù)為84 500）。將Creo三維軟件建立的平壟部件模型導入EDEM軟件中，仿真土槽模型與平壟部件模型如圖3所示。

圖3 仿真試驗結構圖Fig.3 Structural diagram of simulation test

仿真試驗以作業(yè)后耕作地表與非耕作地表高度差和輸土距離為主要評價指標，以平壟作業(yè)的前進速度、平壟刀軸轉速、螺旋葉片螺距為影響因素，設計1個3因素3水平的正交試驗進行平壟作業(yè)仿真模擬。通過查閱農業(yè)機械設計手冊及參照相關旋耕機作業(yè)參數(shù)，初步確定前進速度為0.42、0.69、0.97 m/s，平壟刀軸轉速120、180、240 r/min，螺旋葉片螺距120、160、180 mm。試驗過程每種因素組合重復3次，仿真時間均為16 s，時間步長1.4×10-5s。

2.3 平壟模擬結果分析

根據(jù)仿真試驗設計過程，每種組合因素試驗中每間隔0.5 m取土槽模型橫斷面進行測量；取耕作地表與非耕作地表高度差和螺旋葉片輸土距離2個參數(shù)值，結果如表2所示。

表2 仿真結果Table 2 Simulation result

續(xù)表：

由表2中極差分析可知，3種因素對耕作地表與非耕作地表高度差影響程度由高到低為：螺旋葉片螺距、平壟刀軸轉速、前進速度；3種因素對輸土距離影響程度由高到低為：平壟刀軸轉速、螺旋葉片螺距、前進速度。由2種評價指標中的均值，優(yōu)選出平壟部件選用作業(yè)參數(shù)為：前進速度0.97 m/s，平壟刀軸轉速180 r/min，螺旋葉片螺距160 mm，并以該參數(shù)組合在此進行試驗，其仿真模擬過程如圖4所示。

圖4 平壟部件仿真過程圖Fig.4 Simulation process diagram of ridge-flat parts

在圖4中由左至右劃分4個作業(yè)后區(qū)域，其中1號區(qū)域為未耕作區(qū)，由于平壟部件減速箱是中置結構，導致平壟作業(yè)時產生漏耕現(xiàn)象，因此該區(qū)域土壤顆粒仍保持原有狀態(tài)；2號區(qū)域為平壟作業(yè)區(qū)域，該作業(yè)區(qū)域內對壟底層土壤顆粒擾動量較小，作業(yè)后壟層剩余土壤厚度很均勻；3號區(qū)域是由平壟刀具、螺旋葉片與擋土板共同作用將土壤顆粒輸送而形成的，該區(qū)域內土壤顆粒堆積較多，但整體高度較1號區(qū)域內低且較為平坦；4號區(qū)域內壟層土壤顆粒較少，說明螺旋葉片輸土過程對非作業(yè)區(qū)域內的影響較小，不會影響已經平整后地表的精度。對仿真模擬過的土壤模型進行實際測量，其結果為：耕作地表與非耕作地表高度差為16 mm，平均輸土距離103 mm，滿足平壟作業(yè)要求，因此將該組參數(shù)確定為平壟部件作業(yè)參數(shù)。

3 耦合仿真分析

3.1 Hypermesh前處理

將導入到EDEM軟件中的平壟部件中帶狀螺旋葉片模型（.stp）導入到HyperMesh軟件中，首先以拓撲狀態(tài)并進行渲染的方式來檢查模型表面是否存在缺失或重疊，根據(jù)檢查結果對幾何模型進行清理修復［18-19］，最終成為閉合的幾何體；然后進行網(wǎng)格劃分，網(wǎng)格的類型和大小都會對仿真分析結果產生很大影響，在該軟件中以六面體網(wǎng)格的收斂性最好且計算精度高［20］，故對螺旋葉片劃分大小為5 mm的六面體網(wǎng)格，并對劃分完成的網(wǎng)格進行質量檢查，保證劃分網(wǎng)格時的網(wǎng)格變形控制在不影響仿真分析質量的范圍內。最終對帶狀螺旋葉片模型劃分的網(wǎng)格模型得到104 121個節(jié)點和50 302個單元，網(wǎng)格劃分模型如圖5所示。

圖5 網(wǎng)格劃分模型Fig.5 Mesh generation model

3.2 Workbench靜力學分析

靜力學分析之前，需要將前期仿真分析得到的數(shù)據(jù)模型導入到ANSYS Workbench中，在EDEM軟件中將仿真得到的土壤對螺旋葉片產生的應力結果以.axdt格式導出，然后再導入至ANSYS Workbench中，在HyperMesh中將劃分好網(wǎng)格模型的螺旋葉片以.inp文件格式導出后也再導入至ANSYS Workbench中，至此完成軟件間的數(shù)據(jù)交換。

工作過程的載荷主要由彎刀切削土壤時土壤對刀具的切削阻力和帶狀螺旋葉片排土過程受到的土壤阻力兩部分組成，本研究僅針對帶狀螺旋葉片進行分析。在靜力學分析時，將土壤對帶狀螺旋葉片部件產生的阻力設置成壓力，且該壓力大小數(shù)值與EDEM仿真過程導出的應力數(shù)據(jù)一致，將固定約束施加在帶狀螺旋葉片與刀軸連接位置處。最后，在ANSYS Workbench后處理模塊中進行求解分析，得到帶狀螺旋葉片的總變形云圖、等效應力云圖和等效應變云圖如圖6、圖7和圖8所示。

圖6 總變形云圖Fig.6 Total Deformation Nephogram

圖7 等效應變云圖Fig.7 Equivalent strain nephogram

圖8 等效應力云圖Fig.8 Equivalent stress nephogram

由帶狀螺旋葉片的靜力學分析結果可以得出：

（1）平壟作業(yè)過程中，帶狀螺旋葉片最大總體變形產生在其螺旋外邊緣處，且靠近平壟部件中間位置的帶狀螺旋葉片外邊緣處變形量最大。帶狀螺旋葉片輸土過程越靠近其外側輸土量越大，且?guī)盥菪~片外邊緣處距離其固定端最遠，所以造成該處變形量最大，最大變形為1.46 mm，在帶狀螺旋葉片使用許用范圍內。

（2）帶狀螺旋葉片最大應變和最大應力均產生在其與刀軸連接位置處，且在兩者固定連接位置處變形與應力最嚴重，最大應變?yōu)?.16×10-4mm，最大應力為122.22 MPa，等效應力、應變值均在帶狀螺旋葉片許用值范圍內，即帶狀螺旋葉片的剛度和強度滿足平壟作業(yè)的使用要求。

4 結論

（1）本文以EDEM、Hypermesh以及ANSYS Workbench三軟件聯(lián)合仿真的方法對平壟部件進行平壟過程仿真模擬以及關鍵部件的后處理分析，根據(jù)河北保定地區(qū)土壤類型選取土壤顆粒接觸類型并建立了離散元仿真結構模型。

（2）建立離散元仿真正交試驗，根據(jù)仿真模擬作業(yè)后耕作地表與非耕作地表高度差和輸土距離2個評價指標確定平壟部件最優(yōu)工作參數(shù)為：前進速度0.97 m/s，刀軸轉速180 r/min，帶狀螺旋葉片螺距160 mm，并在該參數(shù)下再次進行平壟模擬，其作業(yè)結果為耕作地表與非耕作地表高度差為16 mm，平均輸土距離103 mm，滿足平壟作業(yè)要求。

（3）根據(jù)有限元后處理分析得出，帶狀螺旋葉片最大變形為1.46 mm，最大應變?yōu)?.16×10-4mm，最大應力為122.22 MPa，其等效應力、應變值均在帶狀螺旋葉片許用值范圍內，表明平壟部件設計合理。