丘陵地區剛性鎮壓輪性能仿真與試驗

劉宏俊 韓濟遠 陳佳奇 呂金慶 趙淑紅

(東北農業大學工程學院， 哈爾濱 150030)

0 引言

在丘陵地區鎮壓作業時，由于地塊小，坡度大，在保證播種機平穩性的前提下[1-2]，對鎮壓要求較高。丘陵地區的土壤條件差異較大，如地勢高的土壤含水率與地勢低的土壤含水率差異明顯；陽面坡地與陰面坡地的土壤含水率差異較大[3]。當土壤含水率變化時，需要控制載荷和加載時間(前進速度)，達到較優的剛性鎮壓輪壓實效果。載荷、前進速度和土壤含水率共同影響剛性鎮壓輪壓實效果。由上述可知，在研究丘陵地區剛性鎮壓輪與土壤相互作用時，不能像平原地區采用準靜態原則研究鎮壓與土壤相互作用機理[4-5]，而應該加入速度變量，研究動態作用機理。本文對丘陵地區鎮壓對應的動態下陷量進行機理研究，完善適應于丘陵地區的鎮壓技術理論，逐步完善提升地區配套農業機械技術及裝備，縮短與平原地區配套農業機械技術水平的差距[6-7]。

農業機械關鍵部件的研究多采用傳統試驗法(田間試驗方法和土槽試驗法)[8-10]。隨著土壤本構模型及光滑粒子流體動力學(SPH)算法的日益完善[11-14]，基于SPH算法的仿真試驗法在農業機械中的觸土部件與土壤的相互作用研究中越來越受到重視[15-18]。上述研究表明，SPH算法應用于剛性鎮壓輪與土壤的相互作用的研究具有可行性。

本文通過機理分析確定丘陵地區的動態土壤下陷量與土壤條件、工作參數的關系，并在此基礎上建立剛性鎮壓輪與土壤相互作用的三維動態有限元模型，采用SPH算法研究剛性鎮壓輪與土壤的動態過程。以土壤含水率、載荷和前進速度為試驗因素，選擇下陷量和作業阻力為試驗指標，旨在揭示丘陵地區剛性鎮壓輪-土壤的工作機理，以期為丘陵地區播種機及配套技術的研究提供理論及依據。

1 丘陵地區鎮壓作業機理分析

1.1 動態下陷量機理分析

鎮壓過程中土壤的力學性質和變形過程較為復雜，目前對某些參數(流動性、時間效應及應力應變等)的綜合分析及相互作用關系的研究還不完善，而且丘陵山區條件下的相互作用研究鮮有報道。由于丘陵地區的特殊地形，鎮壓作業時土壤的物理性質、載荷、速度相互作用是個動態的過程，動態過程中剛性鎮壓輪下陷量的機理模型還未有人提出，而目前研究該過程多是在一定的假設條件下進行分析推算或由試驗方法測得，往往是一種準靜態預測。本文在基于準靜態原則的鎮壓與土壤的機理模型中加入速度變量，建立動態鎮壓與土壤的機理模型。

圖1 剛性鎮壓輪受力分析Fig.1 Force analysis of rigid press wheel

dQ=pBds

(1)

式中B——剛性鎮壓輪寬度

p——土壤抗壓強度

剛性鎮壓輪受到的載荷與土壤的總抗壓反作用力大小相等，則

(2)

(3)

由Bekker經驗公式確定土壤抗壓強度p與下陷量z0關系方程為

(4)

式中k——土壤變形模量n——土壤變形指數

土壤變形模量k與速度v2z的關系方程為[19]

(5)

式中k0——土壤變形靜態模量

m——速度變化指數

由基點法及三角形相似定理得

(6)

聯立式(3)～(6)可得

(7)

式(7)中有兩個變量：x和z0。式(7)積分前，需將x替換z0。由圖1的幾何關系找出x與z0的關系式為

x2=r2-[r-(z1-z0)]2

(8)

簡化式(8)可得

x2=2r(z1-z0)-(z1-z0)2

(9)

下陷量差值(z1-z0)相對剛性鎮壓輪直徑2r較小，一般忽略不計[4]，簡化式(9)可得

x2=2r(z1-z0)

(10)

對式(10)進行變換得

(11)

將式(11)代入式(7)中，可得

(12)

(13)

對式(13)進行積分可得

(14)

由式(11)可知，當z0=0時，x1與z1的關系式為

(15)

將式(15)代入式(14)得

(16)

由式(16)可知，剛性鎮壓輪的結構參數、土壤條件、前進速度和載荷共同影響土壤下陷量。實際鎮壓狀態時，不同下陷情況下，剛性鎮壓輪與土壤的相互作用及運動趨勢均不確定，而土壤作用模型可以獲得之間的關系，前進速度越大，下陷越小，反之相反。載荷越大，則下陷量越大，反之相反。土壤相關參數對下陷量有影響，而具體的影響規律有必要通過試驗進一步確定。

1.2 作業阻力機理分析

鎮壓作業時，由剛性鎮壓輪與土壤在接觸處變形和摩擦而產生的力為作業阻力[20-21]，其由3部分組成：剛性鎮壓輪壓實土壤過程中前進阻力；推移土壤，形成波浪狀凸起所需的推土阻力；土壤粘附在剛性鎮壓輪表面上所需的粘附阻力。由于形成土壤粘附的因素較多且復雜，無法用合適的數學公式表達，同時仿真模擬時無法模擬土壤粘附，在研究過程中重點分析前進阻力和推土阻力。由文獻[22]和圖1可確定前進阻力F1和推土阻力F2為

(17)

(18)

其中

Kc1=(Nc1-tanφ)cos2φ

式中Nc1、Nγ——太沙基承接能力系數

φ——土壤內摩擦角

c1——土壤內聚力

γ——土壤密度

剛性鎮壓輪鎮壓過程中的作業阻力F3為

F3=F1+F2

(19)

由式(16)、(17)可知，前進阻力F1與結構參數、土壤條件、前進速度和載荷有關；由式(18)可知，推土阻力F2主要與土壤條件和結構參數有關。故確定前進阻力F1和推土阻力F2的合力(作業阻力F)與結構參數、土壤條件、前進速度和載荷等均有關系，而具體的影響規律有必要通過試驗進一步研究。

2 剛性鎮壓輪動力學模型建立

剛性鎮壓輪在丘陵山區的市場占有率很高，同時由于本文的丘陵地區鎮壓作業機理分析是基于剛性鎮壓輪建立的，為了保證研究的連貫性且具有代表性，本文選擇剛性鎮壓輪作為研究對象。

2.1 土壤本構模型

針對東北丘陵地區黑壤土的特性，仿真采用LS-DYNA中的MAT147(MAT_FHWA_SOIL)土壤材料模型。該土壤材料模型采用修正的Mohr Coulomb屈服準則，其應力不變量等式表示為

(20)

式中P——靜水壓力

J2——應力偏張量的第二不變量

K(θ)——應力羅德角函數

α——修正后屈服面和標準Mohr-Coulomb屈服面之間貼近程度的參數

c——黏聚力

將土壤視為各向同性材料，及各部分的含水率、堅實度、密度等物理參數一致。為了解土壤條件對丘陵地區鎮壓輪的影響規律，選擇3種典型東北地區土壤。東北丘陵地區黑壤土的主要參數如表1所示，其余參數采用軟件默認值。

表1 土壤材料參數Tab.1 Parameters of soil

2.2 剛性鎮壓輪-土壤相互作用有限元模型

2.2.1剛性鎮壓輪有限元模型

為了控制仿真時間等問題，忽略影響仿真精度較弱的部位,簡化剛性鎮壓輪物理模型為中空圓柱體(直徑為300 mm，寬度為210 mm)。在LS-PrePost軟件直接進行三維建模，并采用sweep方式對其進行網格劃分，其網格單元類型為Solid164，如圖2所示。剛性鎮壓輪材料參數：密度為7 800 kg/m3；彈性模量為2.1×1011Pa；泊松比為0.3。

圖2 剛性鎮壓輪-土壤有限元模型Fig.2 Finite element model of wheel and soil

2.2.2土壤有限元模型

在LS-PrePost軟件中直接按土壤模型長度為1 200 mm、寬度為400 mm、高度為400 mm的尺寸生成SPH土壤粒子，臨近粒子之間的距離為10 mm，如圖2所示。

2.3 邊界條件

邊界約束：約束模型中土壤底部的全部自由度；約束剛性鎮壓輪除繞軸轉動以外的所有旋轉自由度，同時僅約束剛性鎮壓輪橫向移動自由度；整個模型加入重力，重力加速度為9.8 m/s2。接觸算法為contact-automatic-nodes-to-surface。載荷加載：按照仿真要求，對鎮壓輪施加載荷及前進速度。LS-PrePost軟件生成相應的K文件，然后遞交給MPP LS-DYNA R8.1.0處理器求解，生成求解報告。

2.4 鎮壓作業過程分析

鎮壓輪在不同工況下，雖然指標數值不同，但是狀態及規律相類似[23-24]，故本文以前進速度為3 km/h、載荷為600 N、土壤含水率為12%仿真模型為例，重點分析剛性鎮壓輪與土壤相互作用的動態過程。

圖3為不同時刻剛性鎮壓輪與土壤相互作用狀態。由圖3中的整個運動過程可知，不同深度處土壤豎向應力均有所增加，但不同土壤深度處，土壤豎向應力相差較大，且振幅的變化對表層土壤的豎向應力影響較大，而對深層土壤的豎向應力影響較小。

圖3 剛性鎮壓輪鎮壓土壤動態過程Fig.3 Soil dynamic process of wheel

圖4為下陷量的時間過程曲線。由圖4可知，當載荷一定時，只要時間足夠長，下陷量就能達到最大值，同時相比運動狀態時，下陷量更大。由圖3b到圖3d過程中，剛性鎮壓輪的下陷量降低。由圖3的定性分析和圖4的定量分析確定前進速度影響下陷量。圖4中包括3階段：①a1階段為靜壓階段，土壤被向下和向兩側擠壓，發生土壤變形，此階段下陷量先陡增，后期平穩，這與平板試驗的結果相符合，該階段可證明仿真試驗結果的可靠性。②b1階段屬于靜壓穩定狀態，土壤阻力和載荷及重力幾乎相等。③c1階段為剛性鎮壓輪運動階段，在c1階段的初始階段，由于前兩階段的下陷，造成輪子下方的土壤變硬。由b1階段向c1階段過渡時，在給定剛性鎮壓輪前進速度后，剛性鎮壓輪運動被向上抬起，故剛性鎮壓輪向上運動更大，進而土壤下陷量降低。在c1階段的中后期，剛性鎮壓輪運動趨于相對平穩，但相比于b1階段，下陷量圍繞某一值上下波動，出現這種現象的主要原因是鎮壓作業時，土壤阻力發生變化，造成剛性鎮壓輪激振，這與實際情況相符合。在0.55～1.2 s，剛性鎮壓輪向前運動，輪子下方的土壤發生擠壓變形，后方的土壤發生塑性變形，對距離土壤表面越遠的土壤作用越輕。土壤向斜向下方擠壓土壤，擠壓的土壤發生變形，變形量較大。

圖4 下陷量與仿真時間關系曲線Fig.4 Relationship curve between sinkage and simulation time

3 多目標優化試驗

3.1 中心復合試驗設計

借助數值模擬法研究試驗因素(前進速度、載荷和土壤含水率)對試驗指標(土壤下陷量和作業阻力)的影響規律。為了保證試驗具有很好的序貫性、避免水平增加帶來的復雜性和減少試驗誤差，采用中心復合表面試驗設計(Central composite face-centered design，CCF)方法優化設計試驗方案。基于CCF方案的編碼值如表2所示,每組試驗重復3次取平均值，在P=0.05 水平進行F檢驗。

表2 因素編碼Tab.2 Experimental factors and codes

3.2 中心復合試驗結果與分析

按照試驗方案安排仿真試驗，仿真試驗結果如表3所示。其中x1、x2、x3分別表示載荷W、前進速度v和土壤含水率σ編碼值。

3.2.1回歸方程

采用Design-Expert 8.0.5b軟件對試驗數據進行分析，在置信度α=0.05下進行F檢驗，將不顯著的因素和交互作用的偏平方和及自由度并入剩余平方和后，重新進行分析，保證所有因素都達到顯著或極顯著水平，得到優化的回歸方程

z=14.46+5.75x1+2.69x2+3.81x3-2.24x1x3+ (21)

3.2.2回歸方程方差分析

表4 土壤下陷量回歸模型方差分析Tab.4 Variance analysis of regression equation for sinkage

注：*表示較顯著(0.05

表5 作業阻力回歸模型方差分析Tab.5 Variance analysis of regression equation for working resistance

圖5 因素之間的交互作用影響試驗指標的響應曲面Fig.5 Response surfaces between interactions of factors and experimental indicators

3.2.3各因素對試驗指標影響規律

因素交互與試驗指標之間的關系規律如圖5所示。由圖5a可知，在載荷一定時，前進速度與下陷量的關系趨于線性，同時隨著前進速度的增加，下陷量增加；在前進速度一定時，下陷量隨著載荷的增加，先降低后增加。當載荷較大和前進速度較大時，下陷量較大。當載荷較小和前進速度較小時，下陷量較小。由圖5b可知，在載荷一定時，土壤含水率與下陷量的關系趨于線性；在土壤含水率一定時，載荷與下陷量的關系也趨于線性。下陷量均隨著因素增加而增加。由圖5c可知，在載荷一定時，當前進速度較小時，前進速度增加，作業阻力增加不明顯，而當前進速度超過某一前進速度后，作業阻力增加迅速；在前進速度一定時，載荷對作業阻力影響較明顯，隨著載荷遞增，工作阻力增加。由圖5d可知，在載荷一定時，當土壤含水率較大或較小，作業阻力較小，而當土壤含水率處于中間水平附近，作業阻力較大。在一定土壤條件下，以一定的載荷壓實土壤時，由于土壤中含水率由少增多，土壤下陷量逐漸增加。當土壤中含水率達到某一值時，下陷量最大。

3.3 基于R語言蟻群算法的多目標優化

蟻群算法已經成功應用于若干領域，其中最成功的應用是組合優化問題[25-26]。作為一種全局搜索算法，蟻群算法能夠有效地避免局部極優，但對于大空間的多點全局搜索，卻不可避免地增加搜索時間。為了使算法能夠更好更快的找到問題獲得最優解，在其搜索過程中加入針對具體問題的局部搜索算法。

基于R語言的蟻群算法對多目標進行優化。采用加權距離法將兩個目標方程合并成一個全局方程，具體轉換方程y1為

(23)

式中t1x、t2x——第1、2項指標最大值與最小值的差值

t1d、t2d——第1、2項指標的平均值

μ1、μ2——第1、2項指標權重

下陷量直接影響出苗率與幼苗生長情況，并且最終影響作物產量；作業阻力直接影響功耗，即作業成本。由此可知，在鎮壓作業過程中，下陷量是主要考慮指標，而作業阻力是次要考慮指標。根據各指標的重要性及專家經驗法，確定下陷量、作業阻力的權重μ1、μ2分別為0.75、0.25。

目標函數

對目標函數進行蟻群算法多目標優化后，獲得目標函數之間的權衡曲線，一組非支配解(Pareto最優解)如圖6所示。

圖6 螞蟻算法優化后的最優解前沿分布Fig.6 Distribution of Pareto front from ant colony algorithm

丘陵地區實際作業時，土壤含水率隨著地形及地勢改變而改變，屬于不可控因素。當土壤含水率改變時，若不改變相應的載荷和前進速度，鎮壓效果較差。同時土壤含水率的實際變化是一個連續變化的過程，若不是智能鎮壓系統，則無法實時改變載荷和前進速度，來保證鎮壓效果。由文獻[27]可知，丘陵地區的土壤含水率一般在12%～20%波動。本文將土壤含水率范圍分成5份，由優化結果得到每個土壤含水率對應的載荷和前進速度，進而保證丘陵地區鎮壓作業時鎮壓效果一直處于較優水平。由該原則從Pareto最優解集中篩選出土壤含水率為(12±0.1)%、(14±0.1)%、(16±0.1)%、(18±0.1)%、(20±0.1)%的5組最優解，并將求解時的水平與實際值進行編碼轉換，得到最優水平組合，如表6所示。

表6 最優水平合集Tab.6 Optimal level set

3.4 試驗驗證

驗證試驗在黑龍江省農業機械工程科學研究院的室內土槽中進行，如圖7所示。試驗設備：雙向仿形鎮壓裝置[19]；剛性鎮壓輪；TCC-3型土槽試驗車；土槽；六分力測試裝置；SC-900型土壤堅實度儀。試驗條件：土壤為典型東北黑土；土壤容重為1 311 kg/m3；土壤內摩擦角為18.1°。

圖7 土槽試驗Fig.7 Soil bin tests1.雙向仿形鎮壓裝置 2.剛性鎮壓輪

試驗前，按照含水率由小到大依次排列對土槽分段，共5段，每段長度為5 m，如圖7b所示。在土槽內每段土壤含水率均達到對應表6中土壤含水率后，鎮壓輪按照圖8的箭頭方向進行運動。同時按照表6中土壤含水率對應的工作參數進行試驗，每組驗證試驗重復4次。試驗驗證結果與仿真優化結果如圖9所示。

圖8 試驗方案分布及實施圖Fig.8 Distribution and implementation diagram of test plans

圖9 試驗驗證圖Fig.9 Chart of verification test result

由圖9可知，仿真結果與實測結果相差較小，出現這種現象可能是土槽中的土壤情況更復雜，數值模擬無法完全模擬出來，但是模擬結果的誤差均小于12%，由此表明試驗優化結果的可靠性，同時也驗證了仿真的可行性。

4 結論

(1)針對丘陵地區鎮壓作業特點，在傳統剛性鎮壓輪與土壤作用模型的基礎上加入前進速度變量，建立剛性鎮壓輪作業時剛性鎮壓輪與土壤的動態相互作用模型，為研究丘陵地區鎮壓動態過程中

土壤阻力和變形過程提供理論支撐。

(2)采用SPH算法建立剛性鎮壓輪與土壤相互作用的三維動態有限元模型，重點分析剛性鎮壓輪與土壤的相互作用，預測不同的作業參數對剛性鎮壓輪作業效果的影響。

(3)選擇土壤含水率、載荷、前進速度為試驗因素，選擇下陷量和作業阻力為試驗指標，采用中心復合表面試驗方案獲得回歸模型，然后采用蟻群算法對回歸模型進行多目標優化，獲得不同含水率下的最優水平組合，然后進行土槽試驗對最優水平組合加以驗證，對比結果表明試驗結果合理。