濕粘水稻土深松過程離散元分析

丁啟朔 任 駿 BELAL Eisa Adam,2 趙吉坤 葛雙洋 李 楊

(1.南京農業大學工學院， 南京 210031； 2.扎林蓋大學農業工程系， 扎林蓋 06)

濕粘水稻土深松過程離散元分析

丁啟朔1任 駿1BELAL Eisa Adam1,2趙吉坤1葛雙洋1李 楊1

(1.南京農業大學工學院， 南京 210031； 2.扎林蓋大學農業工程系， 扎林蓋 06)

濕粘水稻土的深耕松對作業機具設計與作業參數的要求不同于旱作制，為探究水稻土條件下土壤深松擾動過程及其內在相互作用機理，本文結合田間實測土壤物理參數，借助離散元EDEM軟件，建立適用于粘性水稻土的深松耕作離散元模型。利用粘性水稻土模型對水稻土的機械深松耕作過程進行離散元模擬仿真，并結合田間試驗結果對機具耕作阻力、土壤宏觀擾動進行對比驗證，進一步從仿真的角度揭示深松土壤擾動的微觀過程、土壤失效破碎的微觀機理。結果表明，該模型下深松耕作阻力平均誤差為6.63%，土壤擾動的起壟寬度平均誤差為4.39%，起壟高度平均誤差為19.22%；DEM仿真對土壤的微觀擾動過程分析進一步論證了宏觀試驗測試結果及假說的正確，并且能夠從土顆粒接觸力學層面表達土壤擾動的邊界生成、土體內部破碎等過程；結合DEM仿真角度提出的反映土體破碎程度指標——斷裂系數，與傳統試驗的指標參數碎土系數對比兩者誤差為3.46%，該指標更有利于對土壤破碎過程微觀機理的表達。

深松； 土壤擾動； 離散元； 水稻土

引言

稻麥輪作制是長江中下游地區的主導高產種植模式，然而該模式的機械化耕作問題日趨嚴重，因稻麥輪作周年干濕交替、稻季淹水以及長期淺旋代深耕翻造成土壤硬化、孔隙率低、耕層淺薄等問題日趨嚴重[1-5]，嚴重限制了稻麥作物的高產創建和品質提升，水稻土的合理機械化耕作已成為制約南方稻作制發展的瓶頸問題。為改善水稻土的結構性，近年來人們開始探索將北方已成熟應用的深耕深松技術應用于南方耕作的可能性[6-9]。不過相關研究表明現有的深松技術尚難于滿足不同土壤條件下的深松要求[10]，因此對于特定土壤類型的深松性能差異性研究就顯得尤為重要。黃玉祥等[11-12]在西北地區塿土條件進行了系統研究，分析了該土壤條件的深松擾動行為且提供了作業效果的方法；鄭侃等[13]針對華北平原壤土區的特定土壤類型設計了一種新型的折線刃深松鏟。鑒于水稻土土壤結構的復雜性和特殊性，適用于北方旱作區的深松機具及理論并不完全適合南方稻麥輪作區域[9]。因此須對南方水稻土條件下深松過程深入研究。

現有研究指出，土壤耕作(尤其是深耕松)的機理非常復雜，機械擾動過程受土壤空間變異、耕作機械動力學、土壤破碎等眾多復雜因素影響，僅僅采用傳統的試驗方法難以深刻理解深松土壤的擾動行為[11]，因此需要進一步從內在機制上探討與機械耕作密切相關的土壤物理響應及其表現[14]。但鑒于宏觀測試方法的局限性，其定量化測試與描述并不容易實現，因此，近年來更強調借助離散元等數值模擬手段解析土壤的機械耕作過程[11,13,15-16]。實踐表明，數值模擬方法的確為深刻闡釋深耕松等機械耕作過程提供了有效的解決途徑[10]，鄭侃等[13]、李博等[17]、鄧佳玉等[18]對不同工況下的深松鏟耕作力學行為進行了研究，顯示應用離散元法分析深松鏟耕作過程的力學行為及效應的可行性和有效性，計算仿真已經成為揭示田間試驗耕作力學過程的優勢技術手段。

本文針對水稻土耕作力學行為的特殊性及其定量表達的困難性，以高粘結性為特征的水稻土離散元參數模型對深松土壤擾動過程進行機理分析，進一步結合田間試驗驗證，探討濕粘性水稻土機械化深松的土壤失效機理及深松機具的力學行為與土壤結構效應，以期從微觀擾動過程及破碎機理的仿真角度實現與田間宏觀測試結果的配合，為南方耕作機具的研究提供借鑒及方法學的支撐。

1 田間試驗

1.1 試驗地點

試驗地點位于南京市六合區八百橋試驗田(118°59′E、31°98′N)，地勢平坦，土壤類型為壤質粘土，該地區常年稻麥輪作。土壤pH值為7.6，土壤砂粒、壤粒、黏粒、有機物質量分數分別為24.06%、40.39%、35.55%和2.27%，液限為41.94%，塑限為25.18%。0～40 cm土層土壤容重、含水率、總孔隙度及硬度參數見表1。試驗前進行人工地表清茬。

表1 0～40 cm土層基礎物理參數Tab.1 Basic physical properties in 0～40 cm soil layer

1.2 試驗方法

為保證田間試驗的可靠性，本文在土壤耕作原位綜合試驗臺[19](圖1)上進行原位深松耕作試驗。試驗臺內部臺車上安裝2個LKL-101/1T型拉壓力傳感器用于測量深松鏟的水平耕作阻力。利用Excel 軟件進行數據分析，獲取耕作阻力波動曲線(圖2a),對信號進行濾波處理后(圖2b)截取波動穩定段并求其平均值，作為相應耕作深度下的水平耕作阻力。

圖1 土壤耕作原位綜合測試平臺Fig.1 Schematic view of in-situ experimental platform1.牽引電動機 2.控制柜 3.鏈條 4.信號拖鏈 5.臺車機構6.升降電動機 7.臺車導軌 8.電源拖鏈 9.信號終端計算機10.移動導軌

南方水稻土的適宜深耕深松深度尚難明確，本文進行5個耕作深度(100、150、200、250、300 mm)的深松深耕試驗，一方面可以完整反映耕深變化的耕作效應，另一方面則可觀察水稻土條件深松鏟不同部位的觸土性能差異，耕作速度為0.1 m/s。一次深松作業完成之后，隨機選取2 m耕作區域作為測量區域。沿著裂紋方向將測試區域內的不同尺度的土垡緩慢取出放于海綿墊上，自然風干后進行篩分處理，土垡尺度過大(大于64 mm)且難于篩分時采用手工測量。小尺度區間(32～64 mm、16～32 mm、8～16 mm、0～8 mm)的土壤破碎體使用篩分分級，分級尺度區間的上、下限算術平均值作為該尺度區間的平均土塊徑[20]，計算結果依次為48、24、12 mm，小于8 mm土壤粒徑默認為4 mm。上述試驗操作及測試方法見文獻[9]，所有試驗均重復3次。

圖2 耕作阻力波動曲線Fig.2 Wave curves of tillage draught

2 EDEM仿真分析

2.1 深松鏟建模

圖3 試驗用深松鏟Fig.3 Subsoiler used in field test and its structural parameters

本試驗選用鑿型深松鏟進行深松耕作，深松鏟實物圖如圖3a所示，鏟柄為圓弧形(仿制JB/T 9788—1999《深松鏟和深松鏟柄》所規定的圓弧形深松鏟柄)，其結構簡圖如圖3b所示。其中鏟柄高度600 mm，鏟柄寬度60 mm，鏟尖長度160 mm，鏟柄曲率半徑300 mm，入土角為22°。為保證仿真分析的準確性，對試驗用深松鏟(圖3)按照1∶1的比例，采用Pro/E建立深松鏟的三維結構模型，并保存為.igs格式導入EDEM軟件(圖4)。

圖4 EDEM仿真模型Fig.4 Simulation model of EDEM

2.2 土壤模型建立

本研究選用球形顆粒模擬田間土壤顆粒，離散元仿真中顆粒尺寸減小會導致仿真耗時呈幾何級數增加，故仿真中的顆粒尺寸受計算時間和存儲空間限制，總是比真實土壤尺寸大[15-16]，本文選定土壤顆粒半徑為8 mm，處于可接受范圍內。結合表1田間實測土壤狀態，建立兩種不同物理參數(耕層土壤、底層土壤)的土壤參數模型模擬田間水稻土分層結構，其中0～150 mm為耕層土壤，150～400 mm為底層土壤。

2.2.1 土壤接觸模型

接觸模型是離散元法的重要基礎，其實質是準靜態下顆粒固體的接觸力學彈塑性分析結果[21]。接觸模型的分析計算直接決定了顆粒受力和力矩的大小，對不同的仿真對象須建立不同的接觸模型，用以提高仿真結果的準確性。試驗地區為壤質粘土，顆粒表面粘附力較大(粘粒35.5%)，土顆粒接觸模型選定為顆粒之間存在粘結作用的Hertz-Mindlin with Bonding模型，該模型可以很好地表達粘性土壤的破碎過程[15-16,22-23]。

該模型接觸情況下，粘結作用既能傳遞作用力又能傳遞力矩，因此粘結可看作一系列具有恒定法向和切向剛度的彈性彈簧，這些彈簧均勻分布在兩顆粒之間的球形區域內。粘結作用產生后，外界給予土顆粒的力Fn、Fτ和力矩Mn、Mτ隨著時間步通過下面公式[24-25]從零開始增加。

δFn=-vnSnAδt

(1)

δFτ=-vτSτAδt

(2)

δMn=-ωnSnJδt

(3)

(4)

其中

式中A——接觸區域面積Rb——粘結半徑J——截面極慣性矩Sn、Sτ——粘結顆粒法向、切向剛度

vn、vτ——顆粒運動速度的法向、切向分量

ωn、ωτ——顆粒角速度的法向、切向分量

δt——時間步長

當外界作用力超過某個定義值時，粘結就被破壞。定義法向和切向臨界應力為

(5)

當外界的作用力破壞了土壤顆粒之間的粘結后，土顆粒之間的相互作用將不再受到粘結作用的影響。

2.2.2 土壤模型參數

在離散元EDEM軟件中，模型參數主要分為材料參數和接觸參數，材料參數包括剪切模量、密度、泊松比，接觸參數主要包括材料間的摩擦因數、恢復系數。材料參數和接觸參數可以通過實地、實驗測量及查閱已有本地區相同土壤類型的相關文獻得到，其中接觸參數通過查閱文獻[13,16,26-27]可得：土壤-土壤及土壤-深松鏟恢復系數均為0.6、土壤-土壤靜摩擦因數0.5、土壤-深松鏟靜摩擦因數0.6、耕層土壤間動摩擦因數0.6、耕層土壤-深松鏟動摩擦因數0.35、底層土壤間動摩擦因數0.24、底層土壤-深松鏟動摩擦因數0.13。本文選用的針對粘性土的Hertz-Mindlin with Bonding土顆粒接觸模型本身有5個微觀參數，分別是：粘結法向剛度Sn、粘結切向剛度Sτ、粘結法向臨界應力σmax、粘結切向臨界應力τmax及顆粒粘結半徑Rb。其中，前面4個參數反映顆粒之間的粘性大小，顆粒粘結半徑則可以反映濕顆粒含水率的高低，這5個參數綜合可以體現濕粘水稻土的高含水率、高粘性的特點。粘結半徑Rb可通過公式計算出，在顆粒半徑一定的情況下，濕顆粒的粘結半徑可根據材料的密度、含水率計算得到[27]。其他微觀參數由于不便通過試驗測出，故本文選擇虛擬仿真標定的方法得出。參照文獻[28-29]的方法，為了減少標定參數的個數，可以取Sn=Sτ，σmax=τmax，已有的研究表明該模型下顆粒行為對粘結剛度參數不敏感，因而可取粘結剛度為5×107N/m3[30]。臨界應力的取值關系到粘結強度，進而決定耕作過程中土壤的破碎程度及耕作阻力，因此本文在其他參數確定的情況下，以田間實測耕作阻力為標準校核臨界應力的取值，建立適合本試驗過程的深松鏟-土壤離散元耕作模型。具體方法是在EDEM軟件中建立如圖4所示的模型，以顆粒接觸模型中臨界應力為變量(其他參數如前文所述)進行深松耕作的仿真試驗，得出不同臨界應力下土壤對深松鏟的水平耕作阻力，然后參照田間試驗相同工況下的水平耕作阻力進行誤差分析，確定適宜的參數取值。

考慮到土壤模型的分層結構，分別以150 mm、300 mm為耕作深度校核耕層土壤和底層土壤的粘結強度。應先校核耕層土壤的粘結強度，得到確切的耕層土壤顆粒間臨界應力取值后，再以此為基礎校核底層土壤顆粒粘結的臨界應力值。田間試驗下測得150 mm、300 mm耕深下的耕作阻力分別為1.47 kN和5.21 kN。代表顆粒粘結強度的臨界應力可以以宏觀土壤剪切、壓縮破壞強度為參考[27]，田間取樣后做單軸壓縮試驗得耕層、底層土壤最大法向應力為105Pa級別，以此為基礎設定一系列代表粘結強度的臨界應力，進行仿真試驗得到相應的耕作阻力并與田間試驗真實阻力進行比較。仿真標定試驗結果如表2所示，最終確定耕層、底層土壤顆粒粘結強度的臨界應力為3×105Pa和5×105Pa，由于模型土壤顆粒與真實土壤顆粒粒徑的差異，模型臨界應力取值比真實土壤破壞強度稍大。最終仿真設計的所有參數如表3所示。

表2 土壤粘結臨界應力校核參數Tab.2 Calibration process and effects of particle critical stiffness on tillage draught

表3 離散元模型基本參數Tab.3 Basic parameters of discrete element model

2.3 EDEM耕作模型

本文采用離散元EDEM軟件進行整個深松耕作過程的仿真。為滿足深松作業的要求，在模型中建立長×寬×高為1 000 mm×800 mm×400 mm的虛擬土槽[11]，深松鏟模型導入后位于土體的一端，結合田間試驗深松鏟耕作速度為0.1 m/s，耕作深度可在模型建立時自主設定。結合田間實測土壤狀態，土體模型分為上下兩層。上層耕作層深度區間為0～150 mm，生成土顆粒按隨機排列，土壤孔隙度約為50%；下層犁底層深度區間為150～400 mm，按簡單立方晶格排列生成土顆粒，土壤孔隙度約為40%。仿真模型中共生成土壤顆粒79 824顆，仿真中涉及到所有參數如表3所示。對應表3在EDEM前處理器模塊依次進行接觸力學模型、土壤顆粒模型、幾何模型和顆粒工廠等的設置。仿真試驗為單因素試驗，通過改變耕作深度來分析不同情況下對土壤擾動情況的影響。在EDEM求解器模塊對仿真時步、仿真時間、數據保存間隔時間及網格大小等進行設置。仿真開始時生成土壤顆粒，待顆粒沉降穩定后深松鏟開始運動，仿真時間總共為12 s，0～2 s為土壤顆粒生成及沉降穩定階段，2～12 s為深松鏟運動階段。仿真結束后在EDEM后處理工具模塊進行仿真結果數據的導出和分析。

3 結果與分析

3.1 耕作阻力與土壤擾動對比分析

圖5a反映了不同耕作深度下田間實測水平耕作阻力與仿真水平耕作阻力的變化趨勢。其中耕作深度為150 mm時誤差最小，為0.68%，耕作深度200 mm時誤差最大，為11.05%，5個耕作深度下耕作阻力的平均誤差為6.63%。由圖5a可知耕作阻力隨著耕作深度的增加而增大，當耕作深度在200 mm以內時，耕作阻力平緩增加；當耕作深度大于200 mm時，耕作阻力急劇增加。

圖5 不同耕作深度下實測與仿真結果對比Fig.5 Comparisons of simulated and measured results under different tillage depths

圖5b反映了不同耕作深度下田間試驗土壤宏觀擾動狀態與仿真下相對應值的對比。田間試驗條件下起壟寬度隨著耕作深度先增大后減小，在耕作深度200 mm時達到最大值，為422 mm，耕作深度100 mm時達到最小值，為264 mm；起壟高度隨著耕作深度變化趨勢相對平緩，總體也呈現出先增大后減小，在耕作深度200 mm時達到最大值(163 mm)，耕作深度300 mm時為最小值(103 mm)。仿真環境下起壟高度變化趨勢與田間試驗結果高度吻合，各耕作深度下起壟寬度平均誤差為4.39%；起壟寬度變化總體趨勢一致，但在耕作深度150 mm及300 mm時誤差相對較大，分別為32.03%和39.81%，總體顯示5個耕作深度的起壟高度平均誤差為19.22%。

圖6是仿真環境下不同耕作深度的深松鏟垂直耕作阻力變化情況。其變化規律和水平耕作阻力變化一致，當耕作深度在200 mm以內時，增幅平緩；當耕作深度大于200 mm時，垂直耕作阻力急劇陡增。

圖6 仿真環境下不同耕作深度的垂直耕作阻力 變化曲線Fig.6 Variation curve of vertical tillage draught in simulation under different tillage depths

本文以田間實測耕作阻力為校核標準，結合田間實測土壤物理參數建立水稻土離散元耕作模型。以深松阻力行為作校核顯示仿真耕作阻力與田間實測阻力誤差相對較小，在可接受范圍內，表明DEM模型預測的匹配與適用性。但是鑒于田間試驗條件地表雖經清茬處理仍然存在根茬和秸稈殘留，而且雖然考慮到水稻土在深度方向的分層結構但并未描述其他方向的土壤力學空間變異，因此在土壤宏觀擾動狀態方面仿真與實測結果差異較大，尤其是起壟高度在某些點出現較大誤差。總體上看，仿真結果與田間試驗結果在耕作力學及土壤擾動兩方面的效應基本一致，表明針對粘性水稻土的離散元耕作模型的土壤接觸模型選定、仿真參數設定以及雙層土壤模型建立的準確性。

田間試驗和仿真結果都表明無論是耕作阻力(包括垂直耕作阻力)還是土壤宏觀擾動，200 mm耕作深度是水稻土耕作深度的特征臨界點，該臨界深度對合理選擇適宜的深松深度具有一定理論參考價值[9]，耕作的深度一定要因地制宜，尤其對于深松這種高能耗作業，相近工況下華北平原壤土區深松耕作水平阻力為1.80 kN[13]、東北平原黑土區土壤耕作阻力為3.86 kN[31]，然而在本文水稻土條件下，300 mm耕作深度時深松耕作牽引阻力已達到5.21 kN，這充分反映出濕粘水稻土深松作業與旱作土的差異性。因此，不同土壤條件的作業參數各不相同，且耕作深度等作業參數也需要進一步配合農藝要求。

耕層土壤剪切模量小、顆粒粘結強度低，底層土壤剪切模量大、顆粒粘結強度高，耕作阻力變化趨勢如圖5a所示。土壤宏觀擾動參數隨著耕作深度先增大后減小，在200 mm耕作深度時宏觀擾動參數達到最大，這顯然是土壤自身屬性和深松鏟特點共同作用的結果。當耕作深度較淺(100 mm)時，深松鏟的觸土效應主要由鏟尖和少部分曲柄段提供，此時耕作阻力小，對土壤的擾動也小；當耕作深度逐步增大到200 mm時，深松鏟鏟尖部分和曲柄弧段完全作用于耕層土壤和部分底層土壤，耕作阻力逐漸增大，土壤擾動量增大，曲柄弧形結構的抬土功能得到充分發揮，此時地表起壟參數達到最大值；當耕作深度繼續加到300 mm時，深松鏟作用于耕層土壤和底層土壤，由于底層土壤剪切模量大、粘結強度高，使得耕作阻力急劇增大，深松鏟對土壤的擾動能力被大大削弱，地表起壟參數降低。由此可見，粘性水稻土的DEM模型仿真為土壤宏觀擾動過程提供了機理性解釋，DEM仿真也是實現不同耕作效應的耕作機具設計的重要方法學支撐。

3.2 深松土壤擾動微觀機理

3.2.1 土壤擾動微觀過程分析

在仿真環境下選取耕作方向上中間段厚度為100 mm的局部土體，以300 mm耕作深度為例分析深松鏟對土壤的微觀作用過程。如圖7所示，分別為耕作方向及其側向在不同時刻的截圖，圖8為土體內部粘結結構變化情況，圖9為不同耕作深度下土體內部粘結斷裂總數隨時間的變化曲線。

圖7 不同時刻深松鏟對土壤擾動仿真結果Fig.7 Simulation result of soil disturbance of subsoiler at different time

圖8 關鍵時刻土體內部粘結圖示Fig.8 Diagrams of soil internal particle bond at important moment

圖9 各耕作深度下粘結斷裂總數隨時間的變化曲線Fig.9 Variation curves of total broken bond under different tillage depths

在3.3 s時刻深松鏟鏟尖即將進入土體(圖7)，此時在牽引力的作用下鏟尖部分通過擠壓土壤“鑿”進土體，此時(圖9)粘結斷裂數逐步增加，土壤開始產生破碎；在5.4 s時刻深松鏟曲柄刃段作用于土體，此時由于曲柄段的弧形結構使得作用區域內的土體被向上前方“抬”升(圖7)，此時刻粘結圖示(圖8)深松鏟作用下土體內部以抗剪強度為極限逐步形成深松擾動邊界，擾動區域內土壤在深松鏟作用下被抬升，達到松土效果；在6.6 s時刻深松鏟直柄端作用于土體(圖7)，此時隨著深松鏟的前進直柄端將前一階段抬升起來的土體進行“剪”切，在深松鏟兩側地表形成土垡及其裂紋，粘結圖示(圖8)表層下方土壤在直柄的剪切作用下進一步破碎，達到碎土效果；在5.4～6.6 s時段內(圖9)粘結斷裂總數劇增，這是由于深松鏟曲柄段和直柄端對土壤的抬升、剪切作用土體內部發生著劇烈的剪切失效使土壤充分破碎；在7.5 s時刻深松鏟完全離開土體(圖7)，此時土體會形成深松后的溝槽，溝槽兩側土壤會有一部分回“落”，回填溝槽，土壤粘結斷裂總數趨于一個穩定值(圖9)，也說明深松鏟對這一局部區域土體擾動完成。

傳統田間試驗條件下，傳統的試驗方法不便于深入理解認識深松鏟-土壤的作用機理。已有的研究[10-11,32]表明鑿型深松鏟向前上方擠壓土壤，隨著深松鏟的不斷前進土壤產生剪切裂紋，破裂線從鏟尖延伸到地表(圖8)并發生破碎。本文從微觀的DEM仿真角度進一步解釋和論證了宏觀試驗測試結果及假說的正確。文獻[11]在西北地區塿土土壤類型下對深松土壤擾動過程進行了試驗與仿真研究，著重分析了深松鏟對不同位置土壤的擾動機理，對不同層土壤的運動狀態進行了對比分析。但是其并未結合深松鏟自身結構特點分析機具-土壤作用過程。本文DEM仿真能夠更深一步從微觀層面充分結合深松鏟的結構特點論證從鏟尖到鏟柄對土壤作用的各個細觀過程，進一步明確了深松鏟各部分結構的觸土功能特點。同時結合仿真模型從土顆粒接觸力學層面嘗試表達了機具對土壤擾動的邊界生成、土體內部破碎等過程，更加深入明確了深松鏟對濕粘型水稻土的機械擾動過程，為優化深松鏟結構設計及更好適應南方水稻土深耕松提供了理論參考。

圖10 代表性土壤顆粒位置Fig.10 Location of representative soil particles

3.2.2 土壤擾動位移微觀分析

圖11 代表性顆粒不同方向位移隨時間的變化曲線Fig.11 Changing curves of displacement of representative particles with time in different directions

為進一步定量分析深松對土壤的擾動過程，描述土壤在深松擾動過程的位置變化至關重要。參照文獻[15]的方法，選取代表性位置土壤顆粒(即深松鏟耕作路徑上距地表0、100、200、300 mm處，如圖10中橙色顆粒)分析其在深松鏟擾動過程中X、Y、Z方向的位移變化趨勢。由圖11可知，在耕作方向(X)上0、300 mm處的表層和底層土顆粒運動位移較大，300 mm處顆粒較其他位置先被擾動；在側向(Y)上0 mm處表層土顆粒側向擾動位移明顯大于其他位置，但總體位移變化不明顯；在縱向(Z)上，開始擾動3.3～5.4 s階段在深松鏟的鏟尖和曲柄刃的抬升作用下各位置處土顆粒都向上運動，5.4～7.5 s時段深松鏟直柄端作用于這部分土壤，表層0、100 mm處土顆粒繼續向上運動，而較底層200、300 mm處土顆粒則向下運動。

對3個方向的位移分析表明，深松鏟觸土過程的土壤主要在耕作方向和縱向上運動，側向擾動不明顯。這與陳坤等[10]、 郭新榮[32]的研究結果一致，深松過程中土壤主要是向上前方運動。但是結合前文對擾動過程的定性和此處縱向位移定量分析表明，深松擾動過程中土壤有一個明顯的抬升然后再回落的過程，黃玉祥等[11]的仿真研究也發現類似的現象。這說明深松過程土壤在縱向上先被深松鏟抬起，隨著深松鏟的不斷前進，土壤再向下回落，因此土壤不只單單的向上前方運動，而是耕作方向上向前運動而縱向上先向上抬升再向下回落的運動軌跡(圖11c)。除此之外，粘性水稻土在整個過程中伴隨剪切失效進而發生破碎，達到松碎土壤、調理耕層的效果。

3.2.3 土壤破碎程度對比分析

南方濕粘型水稻土的機械耕作特征顯著不同于旱作制，古老東方文明有別于西方的一個典型的特征是稻作制的觸土凸面設計強化了濕粘型水稻土的碎土效應[33]。然而鑒于深松鏟并不具備犁的翻耕效應，因此強化深松鏟在濕粘型水稻土條件的碎土效果就愈發成為現代保護性耕作理論的重要課題。耕后土壤破碎體的尺度分析是土壤宏觀力學結構研究的基本內容[34]，傳統的篩分法可以利用碎土系數從宏觀上反映耕后土壤的破碎程度[35]，但這些都是基于耕后結果的后評價。DEM仿真的優勢之一是有可能從觸土過程中定量土壤的破碎程度。基于這一考慮，本文從仿真模擬手段在微觀角度定義一個反映土壤破碎程度的指標：斷裂系數I，其計算公式為

(6)

式中I——斷裂系數Nb——擾動區域土壤顆粒粘結斷裂數Ni——擾動區域土壤顆粒完整粘結數

擾動區域的邊界劃定為:以耕作方向上耕作穩定段中間的100 mm為區域長度，以耕作后起壟寬度作為區域寬度，耕作深度為區域深度。仿真結束后，提取每個耕作深度下擾動區域內的數據如表4所示。

表4 擾動區域粘結斷裂數據Tab.4 Bond broken data in disturbed zone

圖12是田間試驗后測得不同耕作深度下土壤破碎體尺度分布情況，由圖12可知本文水稻土深松耕作后，土壤破碎體粒徑主要分布在100～400 mm區間(質量百分比大于50%)。耕作后土壤破碎體尺度分布能有效反映土壤的破碎程度，參照文獻[35]的方法以小于200 mm的土塊質量占總碎土質量的百分比定義為碎土系數。當耕作深度在200 mm以內時，粒徑小于200 mm的土壤破碎體所占質量百分比處于一個相對較低的狀態，分別為31.16%、28.62%、21.62%；當耕作深度大于200 mm時則迅速增加，耕作深度300 mm時達到最大值，為61.19%。圖13為田間實測不同耕作深度下碎土系數和仿真情況下土壤顆粒粘結斷裂系數的變化趨勢，兩者變化趨勢基本一致，不同耕作深度下的兩者平均誤差為3.46%。仿真情況下斷裂系數在耕作深度為200 mm以內時，基本保持在30%左右；當耕作深度大于200 mm時，土壤破碎率開始增大，耕作深度300 mm時達到最大值52.52%。

圖12 田間實測土壤破碎體徑級分布Fig.12 Relative size distribution of measured soil clods

圖13 土壤破碎程度定量對比Fig.13 Comparison of measured and simulated soil crushing degrees

以上結果表明，田間試驗下耕作深度較小時(200 mm以內)，深松耕作產生的土壤破碎體主要是粒徑大于200 mm的大土垡結構(破碎土壤質量的60%以上)，土體整體破碎程度不高，此時深松擾動的效果是使被擾動土體破而不碎；仿真模擬環境下的結果顯示，此時土體內部土壤顆粒粘結斷裂系數相對較低，說明土顆粒之間的粘結斷裂不充分、不完全，呈現出的結果就是土體破裂產生大土垡結構，土壤松動、破而不碎。隨著耕作深度加大(超過200 mm)，土壤顆粒粘結斷裂系數增大，此時對比田間試驗土壤破碎體尺度分布結果也能印證這一點，耕作深度250、300 mm時碎土系數分別達到50.07%和61.19%，其中粒徑小于100 mm小尺度碎塊所占比例也分別達到23.02%、23.93%，整體破碎程度明顯高于其他耕作深度情況。

以上對比分析說明，仿真環境下模型土壤顆粒粘結斷裂系數指標也可以定量表達田間試驗下土壤的破碎程度，其大小及變化趨勢與傳統試驗下的碎土系數基本一致。這一指標也更有利于對深松土壤破碎過程微觀機理的表達和理解。旱作制的深松鏟優化設計強調觸土曲面優化[36-37]、自激或受迫振

動[38-39]、仿生減阻[40-41]等措施，例如文獻[13]針對華北平原壤土區不同土層特點設計新型折線破土刃深松鏟，通過優化深松鏟觸土曲面強化破土碎土過程從而達到節能降阻的目的。其旨在通過強化觸土過程達到減小耕作阻力的目的，并未過多關注耕后土壤破碎效果。然而南方水稻土的機械化深耕松更需要解決有效碎土及高效耕作層擾動的問題，本文針對水稻土耕作力學特性的離散元建模及仿真顯示可以從深松鏟的工作參數(耕作深度、耕作速度等)優化入手提高濕粘型土壤的深耕松效應。

4 結論

(1)基于田間實測的南方稻麥輪作區域水稻土耕層土壤、底層土壤的物理性質差異，選用帶有粘結作用的離散元接觸模型，結合田間試驗及虛擬仿真標定的方法建立了離散元雙層深松土壤模型。其中耕作阻力平均誤差6.63%，土壤擾動壟寬平均誤差4.39%、壟高平均誤差19.22%。該模型的建立可為特定區域土壤耕作機具的研究提供理論參考模型。

(2)對深松土壤的微觀擾動過程進行了分析，進一步明確了深松鏟與土壤相互作用機理。從仿真模型粘結斷裂的角度表達了深松擾動邊界的產生過程，土壤產生破碎的過程。對土壤運動位移的分析表明，土壤在耕作方向、垂直方向位移較大，側向土壤擾動不明顯。為南方因地制宜開展合理深松技術，例如新型深松鏟的優化設計、合理的深松鏟布局方式、最佳耕作深度的選擇等提供了科學的決策依據。

(3)應用離散元EDEM軟件，提出反映土壤破碎程度的指標“斷裂系數”。本研究中當耕作深度在200 mm以內時，斷裂系數穩定在30%左右；耕作深度大于200 mm時，斷裂系數逐漸增大到52.52%。與傳統田間試驗下“碎土系數”指標無論是在變化趨勢，還是相對誤差都在可接受范圍內。斷裂系數可以作為衡量耕后土壤破碎程度的指標，斷裂系數越大土壤破碎程度越高，反之亦然。

1 王恩姮，柴亞凡，陳祥偉. 大機械作業對黑土區耕地土壤結構性特征的影響[J]. 應用生態學報，2008，19(2): 351-356. WANG Enheng, CHAI Yafan, CHEN Xiangwei. Effects of heavy machinery operation the structural characters of cultivated soils in black soil region of Northeast China[J]. Chinese Journal of Applied Ecology, 2008, 19(2): 351-356. (in Chinese)

2 孫忠英，李寶筏. 農業機器行走裝置對土壤壓實作用的研究[J]. 農業機械學報，1998，29(3): 172-174. SUN Zhongying, LI Baofa. Research on the effect of soil compaction caused by the device of running agricultural machinery[J]. Transactions of the Chinese Society for Agricultural Machinery, 1998, 29(3): 172-174. (in Chinese)

3 雷友，曹國鑫，牛新勝，等. 土壤深耕對冬小麥根系在土壤剖面分布的影響[J]. 現代農業科技，2011(8): 272-273. LEI You, CAO Guoxin, NIU Xinsheng, et al. Influence of deep tillage on distribution of wheat root system in soil profile[J]. Modern Agricultural Sciences and Technology, 2011(8): 272-273. (in Chinese)

4 楊建昌，杜永，劉輝. 長江下游稻麥周年超高產栽培途徑與技術[J]. 中國農業科學，2008，41(6): 1611-1621. YANG Jianchang, DU Yong, LIU Hui. Cultivation approaches and techniques for annual super-high-yielding of rice and wheat in the lower reaches of Yangtze River[J]. Scientia Agricultura Sinica, 2008, 41(6): 1611-1621. (in Chinese)

5 何進，李洪文，高煥文. 中國北方保護性耕作條件下深松效應與經濟效益研究[J]. 農業工程學報，2006，22(10): 62-67. HE Jin, LI Hongwen, GAO Huanwen. Subsoiling effect and economic benefit under conservation tillage mode in Northern China[J]. Transactions of the CSAE, 2006, 22(10): 62-67. (in Chinese)

6 孫少華，張文斌. 機械化深松技術及蘇州地區適應性研究[J]. 安徽農業科學，2012，40(20): 10701-10702. SUN Shaohua, ZHANG Wenbin. Research of subsoiling techniques and their adaptability in Suzhou[J]. Journal of Anhui Agricultural Sciences, 2012, 40(20): 10701-10702. (in Chinese)

7 張文斌，孫少華，李志剛. 蘇州地區機械化深松技術試驗研究與分析[J]. 農業裝備技術，2012，38(5): 20-23.

8 張軍，丁啟朔，丁為民，等. 分層旋耕對水稻土耕層物理性狀的影響[J/OL]. 南京農業大學學報，2015, 38(6): 1016-1022.http:∥nauxb.njau.edu.cn/oa/DArticle.aspx?type=view&id=201506022.DOI:10.7685/j.issn.1000-2030.2015.06.022. ZHANG Jun, DING Qishuo, DING Weimin, et al. Effect of stratified rotary tillage on paddy soil physical properties[J/OL]. Journal of Nanjing Agricultural University, 2015, 38(6): 1016-1022. (in Chinese)

9 丁啟朔，葛雙洋，任駿，等. 水稻土的深松阻力及土壤擾動效果研究[J/OL]. 農業機械學報，2017，48(1):47-56. http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?file_no=20170107&flag=1. DOI：10.6041/j.issn.1000-1298.2017.01.007. DING Qishuo, GE Shuangyang, REN Jun, et al. Characteristics of subsoiler traction and soil disturbance in paddy soil[J/OL]. Transactions of the Chinese Society for Agricultural Machinery, 2017,48(1):47-56. (in Chinese)

10 陳坤，胡曉麗，趙新子，等. 國內外深松鏟研究現狀與展望[J]. 農業與技術，2010，30(3): 30-34.

11 黃玉祥，杭程光，苑夢嬋，等. 深松土壤擾動行為的離散元仿真與試驗[J/OL]. 農業機械學報，2016,47(7):80-88. http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?file_no=20160712&flag=1. DOI:10.6041/j.issn.1000-1298.2016.07.012. HUANG Yuxiang，HANG Chengguang，YUAN Mengchan,et al. Discrete element simulation and experiment on disturbance behavior of subsoiling[J/OL].Transactions of the Chinese Society for Agricultural Machinery, 2016,47(7): 80-88. (in Chinese)

12 黃玉祥，杭程光，李偉，等. 深松作業效果試驗及評價方法研究[J]. 西北農林科技大學學報:自然科學版，2015，43(11):228-234. HUANG Yuxiang，HANG Chengguang，LI Wei, et al.Subsoiling test and evaluation methodology of tillage quality[J].Journal of Northwest A&F University: Natural Science Edition,2015,43(11):228-234.(in Chinese)

13 鄭侃，何進，李洪文，等. 基于離散元深松土壤模型的折線破土刃深松鏟研究[J/OL]. 農業機械學報，2016，47(9):62-72. http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?file_no=20160910&flag=1. DOI:10.6041/j.issn.1000-1298.2016.09.010. 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/OL].Transactions of the Chinese Society for Agricultural Machinery,2016,47(9): 62-72.(in Chinese)

14 丁啟朔，丁為民. 現代土壤機械耕作研究的綜述[J]. 土壤通報， 2006，37(1): 149-153. DING Qishuo, DING Weimin. Modern soil tillage research: a review[J]. Chinese Journal of Soil Science, 2006, 37(1): 149-153. (in Chinese)

15 方會敏，姬長英，FARMAN Ali Chandio，等. 基于離散元法的旋耕過程土壤運動行為分析[J/OL]. 農業機械學報，2016，47(3): 22-28.http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?file_no=20160304&flag=1. DOI:10.6041/j.issn.1000-1298.2016.03.004. FANG Huimin，JI Changying，FARMAN Ali Chandio, et al. Analysis of soil dynamic behavior during rotary tillage based on distinct element method[J/OL].Transactions of the Chinese Society for Agricultural Machinery, 2016,47(3): 22-28.(in Chinese)

16 方會敏，姬長英，Ahmed Ali Tagar，等. 秸稈-土壤-旋耕刀系統中秸稈位移仿真分析[J/OL]. 農業機械學報，2016，47(1): 60-67. http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?file_no=20160109&flag=1. DOI:10.6041/j.issn.1000-1298.2016.01.009. FANG Huimin,JI Changying,AHMED Ali Tagar,et al. Simulation analysis of straw movement in straw-soil-rotary blade system [J/OL].Transactions of the Chinese Society for Agricultural Machinery,2016,47(1): 60-67.(in Chinese)

17 李博，劉凡一，陳軍，等. 深松鏟耕作阻力影響因素的離散元法仿真分析[J]. 農機化研究， 2015(2): 71-74. LI Bo, LIU Fanyi, CHEN Jun, et al. DEM Simulation on analyzing affecting factors of the soil resistance applied on the subsoiler [J].Journal of Agricultural Mechanization Research,2015(2): 71-74.(in Chinese)

18 鄧佳玉，胡軍，李慶達，等. 基于EDEM離散元法的深松鏟仿真與試驗研究[J]. 中國農機化學報，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)

19 楊艷山，丁啟朔，丁為民，等. 田間原位綜合耕作試驗臺設計與應用[J/OL]. 農業機械學報，2016， 47(1): 68-74. http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?flag=1&file_no=20160110&journal_id=jcsam. DOI:10.6041/j.issn.1000-1298.2016.01.010. YANG Yanshan, DING Qishuo, DING Weimin, et al. Design and application of multi-purpose in-situ tillage tool testing platform[J/OL]. Transactions of the Chinese Society for Agricultural Machinery, 2016, 47(1): 68-74. (in Chinese)

20 WATTS C W, DEXTER A R. Traffic and seasonal influences on the energy required for cultivation and on the subsequent tilth[J]. Soil & Tillage Research,1994, 31(4): 303-322.

21 胡建平，周春健，侯沖，等. 磁吸板式排種器充種性能離散元仿真[J/OL]. 農業機械學報，2014， 45(2): 94-98.http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?flag=1&file_no=20140216&journal_id=jcsam. DOI：10.6041/j.issn.1000-1298.2014.02.016. HU Jianping, ZHOU Chunjian,HOU Chong, et al. Simulation analysis of seed-filling performance of magnetic plate seed-metering device by discrete element method[J/OL]. Transactions of the Chinese Society for Agricultural Machinery,2014, 45(2): 94-98.(in Chinese)

22 張紹軍. 3種深松鏟對深松作業效果影響研究[J]. 農業科技與裝備，2012(6): 23-25. ZHANG Shaojun.Experimental analysis on the subsoiling effect of three subsoilers[J]. Agricultural Science & Technology and Equipment,2012(6): 23-25.(in Chinese)

24 EDEM. EDEM user guide[M]. Edinbuegh: EDEM-Solutions, 2010.

25 POTYONDY D O, CUNDALL P A. A bonded-particle model for rock[J]. International Journal of Rock Mechanics & Mining Sciences,2004, 41(8): 1329-1364.

26 AHMAD F. Performance evaluation of disc type furrow openers in no-till paddy soil environment[D]. Nanjing：Nanjing Agricultural University, 2015.

27 張永良. 基于離散元的逆轉旋耕施肥播種機拋土性能仿真及試驗研究[D].鎮江：江蘇大學, 2012.

28 ASAF Z, RUBINSTEIN D, SHMULEVICH I. Determination of discrete element model parameters required for soil tillage[J]. Soil & Tillage Research, 2007, 92(1-2): 227-242.

29 LINDE J V D. Discrete element modeling of a vibratory subsoiler[D]. Stellenbosch: Stellenbosch University, 2007.

30 MAK J, CHEN Y, SADEK M A. Determining parameters of a discrete element model for soil-tool interaction[J]. Soil & Tillage Research,2011, 118(5): 117-122.

31 齊關宇, 劉林, 趙艷忠, 等. 深松鏟入土深度及鏟形對耕作阻力影響[J]. 農機化研究，2015(11):161-165. QI Guanyu, LIU Lin, ZHAO Yanzhong, et al. Effects of subsoiler’s penetrating depth and spade shape on traction resistance[J]. Journal of Agricultural Mechanization Research, 2015(11):161-165. (in Chinese)

32 郭新榮. 土壤深松部件松土機理的分析[J]. 農機化研究，2005(3): 78-80. GUO Xinrong. Comparative study on mechanism of mutil-subsoiler deep loosening the soil[J]. Journal of Agricultural Mechanization Research,2005(3): 78-80.(in Chinese)

33 沈昌蒲，王秋華，溫錦濤. 從犁的工作部件看中國和西歐設計思想的差異[J]. 農業機械學報， 2007，38(6): 178-183. SHEN Changpu,WANG Qiuhua,WEN Jintao. Difference of design idea of plow share between Chinese and West European[J]. Transactions of the Chinese Society for Agricultural Machinery,2007, 38(6): 178-183.(in Chinese)

34 丁啟朔，丁為民，潘根興，等. 土壤宏觀力學結構與精準耕作[J]. 中國農業科學，2012，45(1): 26-33. DING Qishuo, DING Weimin, PAN Genxing, et al. Soil macro-structure and precision tillage[J]. China Agriculture Science, 2012, 45(1): 26-33. (in Chinese)

35 楊有剛，張宏，馮濤，等. 土壤淺深松聯合松耕機設計和淺松土試驗[J]. 機械工程學報， 2012，48(19): 163-168. YANG Yougang, ZHANG Hong, FENG Tao,et al. Topsoil and subsoil combined cultivator and top-soiling experiment[J]. Journal of Mechanical Engineering,2012,48(19): 163-168.(in Chinese)

36 郭志軍，杜干，周志立，等. 土壤耕作部件宏觀觸土曲面減阻性能研究現狀分析[J]. 農業機械學報，2011，42(6): 47-52. GUO Zhijun, DU Gan, ZHOU Zhili, et al. Actuality analysis of resistance-reducing properties on soil cultivating components with different macroscopic soil-engaging surfaces[J]. Transactions of the Chinese Society for Agricultural Machinery, 2011, 42(6): 47-52. (in Chinese)

37 張強，張璐，于海業，等. 復合形態深松鏟耕作阻力有限元分析與試驗[J/OL]. 農業機械學報， 2012，43(8): 61-65.http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?flag=1&file_no=20120812&journal_id=jcsam.DOI:10.6041/j.issn.1000-1298.2012.08.012. ZHANG Qiang,ZHANG Lu,YU Haiye,et al. Finite element analysis and experiment of soil resistance of multiplex-modality subsoiler[J/OL]. Transactions of the Chinese Society for Agricultural Machinery, 2012,43(8): 61-65.(in Chinese)

38 辛麗麗，李傳琦，梁繼輝，等. 考慮分段土壤作用力的振動減阻分析[J/OL]. 農業機械學報，2014, 45(2): 136-140. http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?flag=1&file_no=20140223&journal_id=jcsam. DOI:10.6041/j.issn.1000-1298.2014.02.023. XIN Lili, LI Chuanqi, LIANG Jihui, et al. Vibrating drag reduction considering acting force of piecewise soil[J/OL]. Transactions of the Chinese Society for Agricultural Machinery, 2014, 45(2): 136-140. (in Chinese)

39 蔣建東，高潔，趙穎娣，等. 土壤旋切振動減阻的有限元分析[J/OL]. 農業機械學報，2012，43(1):58-62.http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?flag=1&file_no=20120112&journal_id=jcsam. DOI:10.6041/j.issn.1000-1298.2012.01.012. JIANG Jiandong, GAO Jie, ZHAO Yingdi, et al. Finite element simulation and analysis on soil rotary tillage with external vibration excitation[J/OL]. Transactions of the Chinese Society for Agricultural Machinery, 2012, 43(1):58-62.(in Chinese)

40 張琰，黃河，任露泉. 挖掘機仿生斗齒土壤切削試驗與減阻機理研究[J/OL]. 農業機械學報，2013, 44(1): 258-261.http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?flag=1&file_no=20130148&journal_id=jcsam. DOI:10.6041/j.issn.1000-1298.2013.01.048. ZHANG Yan,HUANG He,REN Luquan.Experiment and drag reduction mechanism of bionic excavator tooth during soil cutting[J/OL]. Transactions of the Chinese Society for Agricultural Machinery, 2013, 44(1): 258-261.(in Chinese)

41 張金波，佟金，馬云海. 仿生減阻深松鏟設計與試驗[J/OL]. 農業機械學報，2014，45(4): 141-145.http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?flag=1&file_no=20140422&journal_id=jcsam. DOI:10.6041/j.issn.1000-1298.2014.04.022. ZHANG Jinbo, TONG Jin,MA Yunhai. Design and experiment of bionic anti-drag subsoiler [J/OL].Transactions of the Chinese Society for Agricultural Machinery, 2014, 45(4): 141-145.(in Chinese)

DEM Analysis of Subsoiling Process in Wet Clayey Paddy Soil

DING Qishuo1REN Jun1BELAL Eisa Adam1,2ZHAO Jikun1GE Shuangyang1LI Yang1

(1.CollegeofEngineering,NanjingAgriculturalUniversity,Nanjing210031,China2.DepartmentofAgriculturalEngineering,ZalingeiUniversity,Zalingei06,Sudan)

Tool design and working parameters for subsoiling in wet clayey paddy soil are different from those for dry land cropping system. Investigation on the subsoiling process and the governing mechanisms requires detailed description on the underlying influential factors. Field soil parameters were used as a reference to construct a DEM model for subsoiling which was suitable for wet clayey soil analysis and performed in a EDEM software. The proposed DEM model of clayey paddy soil was then implemented to simulate a subsoiling process, in which directly measured draft and macroscopic disturbance of field soil was compared. Microscopic process and mechanisms of soil fragmentation were assessed with simulated stages. Results showed that simulation error with the proposed DEM model was less than 6.63%. Mean error of the resulted micro-relief was 4.39%. Mean error of elevated ridge was 19.22%. The measured results from microscopic approach were re-evaluated and the correctness of related assumptions from microscopic approach was confirmed by DEM results, soil failure boundary evolution and soil failure with particle contacting models were depicted. A soil fragmentation index was proposed for describing soil fragmentation behavior from DEM perspectives. The error of fragmentation index with respect to measured soil fragmentation parameter was less than 3.46%, affirming that the proposed parameter was a suitable tool for microscopic description of soil failure imposed by subsoiling. This research provides a powerful technical basis for the interaction mechanism between soil and soil, and lays a theoretical foundation for the optimization design of the soil contact parts.

subsoiling; soil disturbance; DEM simulation; paddy soil

10.6041/j.issn.1000-1298.2017.03.005

2016-11-14

2016-12-16

國家重點研發計劃項目(2016YFD0300900)、農業部現代農業裝備重點實驗室開放課題項目(201601002)和江蘇省農機基金項目(201-051028)

丁啟朔(1968—)，男，教授，主要從事土壤耕作力學研究，E-mail: qsding@njau.edu.cn

趙吉坤(1977—)，男，副教授，主要從事土壤數值模擬研究，E-mail: jikunzhao_2006@163.com

S222.12+9

A

1000-1298(2017)03-0038-11