基于Plackett-Burman試驗設計與響應面法優(yōu)化玉米秸稈離散元模型

朱惠斌 錢 誠 白麗珍* 李 慧 牟丹磊 李駿杰

(1.昆明理工大學 農(nóng)業(yè)與食品學院，昆明 650500； 2.山東省農(nóng)業(yè)機械科學研究院，濟南 250010)

我國的玉米種植長期采用傳統(tǒng)耕作模式，存在土壤有機碳含量降低、耕層厚度減小、土壤潛在生產(chǎn)力降低等問題[1-3]。保護性耕作具有改善土壤理化性質(zhì)、增強土壤蓄水保墑能力和提高土壤微生物群落數(shù)量等顯著效益[4-6]。免耕播種機是在玉米種植區(qū)進行保護性耕作的重要農(nóng)機具。防堵裝置是免耕播種機的核心部件，其主要作用是切斷或撥離秸稈從而保證免耕播種施肥作業(yè)的順利進行。為提高免耕播種機防堵裝置的作業(yè)性能，應對其進行數(shù)值模擬，從而對防堵裝置的結構進行優(yōu)化改進。防堵裝置數(shù)值模擬的關鍵點是建立可破碎的玉米秸稈模型，從而對切割或撥離秸稈的過程進行數(shù)值模擬。

目前，機械部件的數(shù)值模擬方法主要為有限元法[7-8]和離散元法[9-10]，尤其在農(nóng)業(yè)機械設計與理論分析方面，如對冰草種子[11]和花生[12]進行建模，對深松鏟[13]、開溝器[14]、旋耕刀[15]和排種器[16]等進行優(yōu)化改進。離散元法可用于模擬免耕播種機防堵裝置的作業(yè)過程，為防堵裝置優(yōu)化提供理論依據(jù)，如已有研究建立玉米秸稈層離散元模型[17]和小麥秸稈離散元模型[18-19]對防堵裝置進行數(shù)值模擬。但上述基于Hertz-Mindlin (no slip)模型建立的秸稈離散元模型存在局限性，這些秸稈模型不會因外力而發(fā)生斷裂或破壞，無法對防堵裝置切割秸稈的過程進行數(shù)值模擬。為對防堵裝置切割或撥離秸稈過程進行數(shù)值模擬，以Hertz-Mindlin with bonding模型為接觸模型建立可破碎的玉米秸稈離散元模型。

采用離散元法對免耕播種機防堵裝置進行數(shù)值模擬時，需對秸稈離散元模型進行參數(shù)標定從而減小數(shù)值模擬的誤差。已有研究主要采用堆積角試驗對各類離散元模型進行標定，如苜蓿秸稈[20]、油菜莖稈[21]、水稻種子[22]、胡麻籽粒[23]和三七種子[24]。堆積角試驗標定法不能準確反映出離散元模型的力學特性，因此采用單軸壓縮試驗法對玉米秸稈離散元模型進行標定，以減小免耕播種機防堵裝置數(shù)值模擬的誤差。

本研究擬對玉米秸稈進行物理單軸壓縮試驗，基于Hertz-Mindlin with bonding接觸模型建立玉米秸稈離散元模型；以物理試驗與仿真試驗的臨界載荷相對誤差為試驗指標，采用單因素試驗法、Plackett-Burman試驗設計和Box-Behnken響應面法確定離散元模型的最優(yōu)仿真參數(shù)組合；以期為免耕播種機防堵裝置的數(shù)值模擬提供仿真參數(shù)依據(jù)。

1 材料與方法

1.1 玉米秸稈物理單軸壓縮試驗

試驗材料選用昆明理工大學保護性耕作試驗田種植的“紅單6號”玉米，在玉米收獲2周之后采集秸稈。選取粗壯筆直的節(jié)間玉米秸稈為試樣，標距為60 mm，直徑為25.67 mm，含水率約為53.6%。試驗設備為CSS-44100萬能材料試驗機，最大載荷100 kN。以10 mm/min加載速度和10 kN加載力進行單軸壓縮試驗，試驗重復9次。根據(jù)試驗結果，玉米秸稈單軸壓縮試驗的臨界載荷為935.4 N。

1.2 玉米秸稈仿真單軸壓縮模型

1.2.1The Hertz-Mindlin with bonding接觸模型

The Hertz-Mindlin with bonding接觸模型是基于Potyondya提出的顆粒粘結模型(Bonded particle model，BPM)，主要用于物料破碎、磨損的模擬(圖1)[25]。該接觸模型采用有限尺寸的粘結鍵將顆粒粘結，粘結鍵可以抵抗切向和法向位移，直至達到最大法向和切向剪應力使粘結鍵斷裂。此后，粒子以硬球的形式相互作用。在粘結鍵產(chǎn)生前，通過Hertz-Mindlin接觸模型計算顆粒間的相互作用。

在某一時刻粘結鍵產(chǎn)生后，將粘結力和力矩設置為零，求出在每個時間步所受粘結力和力矩的疊加增量。表達式如下：

Fn=-vnSnAδt

(1)

Ft=-vtStAδt

(2)

Mn=-ωnStJδt

(3)

(4)

式中：Fn為法向粘結力，N；Ft為切向粘結力，N；δt為時間步長，s；Mn為法向力矩，N·m；Mt為切向力矩，N·m；Sn為法向剛度，N/m3；St為切向剛度，N/m3；vn為法向速度，m/s；vt為切向速度，m/s；ωn為法向角速度，rad/s；ωt為切向角速度，rad/s；A為

ti表示切線方向；ni表示法線方向。 Mt，切向力矩；Mn，法向力矩；Fi，粘結力；Rb，粘結鍵半徑；Lb，粘結鍵長度。 ti, tangent direction; ni, normal direction. Mt, tangential moment; Mn, normal moment; Fi, bonding force; Rb, radius of bond; Lb, length of bond.圖1 The Hertz-Mindlin with bonding接觸模型Fig.1 The Hertz-Mindlin with bonding contact model

當模型中粘結鍵的法向應力σmax和切向應力τmax到達臨界值時，粘結鍵產(chǎn)生破壞，表達式如下：

(5)

(6)

1.2.2玉米秸稈離散元模型

玉米秸稈可分為內(nèi)瓤和外表皮，屬于各向異性材料。本研究的玉米秸稈離散元模型主要用于免耕播種機防堵裝置的數(shù)值模擬，只需保證玉米秸稈離散元模型能體現(xiàn)秸稈的力學特性。在實際建模時，對玉米秸稈進行簡化，將秸稈等效為各向同性材料?；陔x散元法建立玉米秸稈離散元模型，以顆粒之間產(chǎn)生粘結鍵來體現(xiàn)玉米秸稈的力學性能。采用顆粒快速填充法建立玉米秸稈離散元模型，建模流程見圖2。

基于Gambit軟件對玉米秸稈進行幾何建模，設置其半徑為13 mm，高為60 mm。對幾何模型進行網(wǎng)格劃分，設置spacing為3、網(wǎng)格類型為四面體網(wǎng)格，完成網(wǎng)格劃分。將已劃分網(wǎng)格的幾何建模導入刀Fluent軟件中，以用戶自定義函數(shù)下Compiled UDFs命令導入編寫好的C文件，對C文件進行編譯，得到包含顆粒坐標的文件。C文件可根據(jù)用戶定義的最大填充系數(shù)和基礎半徑，求出幾何體所有網(wǎng)格的中心坐標，最后生成有關顆粒坐標、速度和粒徑倍數(shù)等信息的文檔。設置最大填充系數(shù)為0.61，基礎半徑為1 mm。

圖2 離散元模型建立流程Fig.2 Process of building discrete element model

建立EDEM文件，對顆粒與幾何模型的本征參數(shù)和接觸參數(shù)進行設置，如表1所示[26]。設置顆粒的半徑為1 mm，接觸半徑根據(jù)試驗方案進行設置。把顆粒坐標數(shù)據(jù)復制到顆粒工廠數(shù)據(jù)文件“Block_Factory_Data.txt”中，將顆粒工廠API文件“myFactory.dll”導入EDEM。設置顆粒之間的接觸模型為“The Hertz-Mindlin with bonding”，相關參數(shù)根據(jù)試驗方案進行設置。將劃分好網(wǎng)格的玉米秸稈模型導入，設置時間步長為Rayleigh的時間步長2%，完成玉米秸稈離散元模型的建立(圖3)。

基于后處理功能將粘結完成的秸稈離散元模型輸出為“.edm”文件，將壓板與支撐底座的幾何模型導入輸出文件，為壓板添加Force Controller，設置壓縮力為10 kN，速度為10 mm/min，玉米秸稈仿真單軸壓縮試驗模型見圖4。

表1 顆粒與幾何模型的仿真參數(shù)Table 1 Simulation parameters of particle and geometric model

圖3 玉米秸稈離散元模型Fig.3 Discrete element model of corn stalk

1.壓板；2.玉米秸稈；3.支撐底座 1.Compression; 2.Corn stalk; 3.Support base圖4 仿真壓縮試驗模型Fig.4 Model of simulated compression test

1.3 試驗因素與試驗指標

選取接觸半徑和Hertz-Mindlin with bonding接觸模型的相關參數(shù)為試驗因素，根據(jù)前期試驗結果，選取因素的試驗水平見表2。在仿真單軸壓縮試驗完成之后，導出壓板在每個時間步所受到力的數(shù)據(jù)文件，根據(jù)牛頓第三定律，該力等效為玉米秸稈在壓縮過程中所受到的力，以其峰值為玉米秸稈所受到的臨界載荷。試驗以仿真試驗與物理試驗中玉米秸稈所受到的臨界載荷的相對誤差為試驗指標，根據(jù)式(7)進行計算：

(7)

表2 Hertz-Mindlin with bonding接觸模型參數(shù)Table 2 Parameters contact model of Hertz-Mindlin with bonding

式中：ER為物理試驗與仿真試驗的相對誤差，%；Fsmax為在仿真單軸壓縮試驗中玉米秸稈所受到的臨界載荷，N；Fpmax為在物理單軸壓縮試驗中玉米秸稈所受到的臨界載荷，N。

1.4 試驗設計

1.4.1單因素試驗設計

以仿真試驗與物理試驗中玉米秸稈所受到的臨界載荷的相對誤差為試驗指標，對6個試驗因素進行單因素試驗，分析各個因素對相對誤差的變化規(guī)律。單因素試驗采用的固定參數(shù)分別為：接觸半徑1.2 mm、單位面積法向剛度1×108N/m3、單位面積切向剛度1×108N/m3、臨界法向應力1×108Pa、臨界切向應力1×108Pa和粘結半徑2 mm。

1.4.2Plackett-Burman試驗設計

Plackett-Burman試驗設計法在可在試驗因素較多的條件下，選取試驗因素的兩個高低水平進行試驗設計來篩選出對試驗結果有顯著影響的因素。在單因素試驗的基礎上，對接觸半徑、單位面積法向剛度、單位面積切向剛度、臨界法向應力、臨界切向應力和粘結半徑共6個試驗因素進行篩選，得到對相對誤差影響顯著的因素，設計了試驗次數(shù)N=12的Plackett-Burman試驗。

1.4.3最陡爬坡試驗

根據(jù)Plackett-Burman試驗結果，對具有顯著性的因素(接觸半徑、單位面積法向剛度和單位面積切向剛度)進行最陡爬坡試驗，使試驗因素水平快速逼近最優(yōu)范圍，以進行Box-Behnken響應面法試驗。

1.4.4Box-Behnken響應面法試驗設計

根據(jù)最陡爬坡試驗試驗的結果，以接觸半徑、單位面積法向剛度和單位面積切向剛度為自變量，相對誤差為響應值，進行Box-Behnken響應面試驗。設計3因素3水平共17次的響應面試驗，探究變量對響應值的影響，因素試驗因素和水平見表3。

表3 Box-Behnken試驗因素與水平Table 3 Experimental factors and levels of Box-Behnken

2 結果與分析

2.1 單因素試驗

1)隨著接觸半徑增大，相對誤差先減小然后增大，這是由于接觸半徑增大使粘結鍵的數(shù)量增加，導致臨界載荷增大(圖5(a))。接觸半徑為1.2 mm時，相對誤差最小。

2)隨著單位面積法向剛度增加，相對誤差先減小然后增大(圖5(b))。這是由于單位面積法向剛度增大使顆粒之間粘結力增加，導致臨界載荷增大。單位面積法向剛度為1×108N/m3時，相對誤差最小。

3)隨著單位面積切向剛度增加，相對誤差先減小然后增大(圖5(c))。這是因為單位面積切向剛度增大使顆粒之間粘結力增加，從而提高玉米秸稈離散元模型抵抗變形和破壞的能力，最終導致臨界載荷增大。單位面積切向剛度為1×108N/m3時，相對誤差最小。

4)隨著臨界法向應力增加，相對誤差的變化不大(圖5(d))。臨界法向應力為1×108Pa時，相對誤差最小。

5)隨著臨界切向應力增加，相對誤差先減小后增大，但變化幅度不大(圖5(e))。臨界切向應力為1×108Pa時，相對誤差最小。

6)隨著粘結鍵半徑增加，相對誤差先減小后增大(圖5(f))。這是因為粘結鍵半徑增大使顆粒粘結力增大，導致臨界載荷變大。粘結鍵半徑為2 mm 時，相對誤差最小。

X1，X2，…，X6的含義見表2，下表、圖同。 The meanings of X1， X2， …， X6 are shown in Table 2, and the following Tables and Figures are the same.圖5 各因素對物理試驗與仿真試驗相對誤差的影響Fig.5 Effect of factors on relative error of physical test and simulation test

2.2 Plackett-Burman試驗

Plackett-Burman試驗結果見表4，對試驗結果進行方差分析，結果見表5?？芍簺Q定系數(shù)R2為0.97，說明模型適用于97.44%的試驗數(shù)據(jù)；X1和X3

表4 Plackett-Burman試驗結果Table 4 Test result of Plackett-Burman

表5 Plackett-Burman試驗方差分析Table 5 ANOVA for the Plackett-Burman design

的P<0.01，表明接觸半徑和單位面積切向剛度對相對誤差的影響極顯著；X2的P<0.05，即單位面積法向剛度對相對誤差的影響顯著；X4、X5和X6的P>0.05，說明粘臨界法向應力、臨界切向應力和粘結鍵半徑對相對誤差的影響不顯著。因此，選取因素X1、X2和X3進行最陡爬坡試驗和中心組合響應面試驗。在進行最陡爬坡試驗和中心組合響應面試驗時，根據(jù)單因素試驗結果取不顯著因素的水平為：臨界切向應力1×108Pa、臨界切向應力1×108Pa和粘結鍵半徑2 mm。

2.3 最陡爬坡試驗

最陡爬坡試驗的結果見表6，可知：試驗5與試驗6的臨界載荷分別為775.73和986.85 N。玉米秸稈物理單軸壓縮試驗的臨界載荷為935.4 N，在試驗5與試驗6的結果區(qū)間之內(nèi)，而且試驗5與試驗6的相對誤差較小，因此選取試驗5和試驗6的因素水平分別為Box-Behnken試驗低水平和高水平。

表6 最陡爬坡試驗結果Table 6 The path of steepest ascent design and results

2.4 Box-Behnken試驗

為得到離散元模型參數(shù)的最優(yōu)組合，進行Box-Behnken試驗，結果見表7。對試驗結果進行回歸擬合后，得到以相對誤差為響應值，X1(接觸半徑)、X2(單位面積法向剛度)和X3(單位面積切向剛度)為變量的回歸方程：

(8)

根據(jù)回歸方程得到各因素交互作用對物理試驗與仿真試驗相對誤差影響的響應面圖(圖6)。

表7 Box-Behnken試驗結果Table 7 Experimental results of Box-Behnken

表8 Box-Behnken回歸模型的方差分析Table 8 ANOVA for Box-Behnken quadratic model

Y為玉米秸稈臨界載荷試驗值與仿真值之間的相對誤差。 Y was relative error between the critical load test value and simulation value of corn stalk.圖6 交互作用對物理試驗與仿真試驗相對誤差的影響Fig.6 Effect of interaction factors on relative error of physical test and simulation test

當單位面積法向剛度為定值時，相對誤差隨著接觸半徑的增大，呈減小趨勢，變化幅度明顯，接觸半徑對相對誤差的影響高于單位面積法向剛度的影響(圖6(a))。

接觸半徑與單位面積切向剛度的交互作用對相對誤差的影響明顯(圖6(b))。當接觸半徑為定值時，隨著單位面積切向剛度的增大，相對誤差呈先減小后增大的趨勢；當單位面積切向剛度的增大為定值時，隨著接觸半徑的增大，相對誤差呈減小趨勢，且變化幅度明顯。

由圖6(c)可知，當接觸半徑為固定水平值時，單位面積法向剛度與單位面積切向剛度的交互作用對相對誤差的影響不顯著。當單位面積法向剛度為定值時，相對誤差隨著單位面積切向剛度的增大，呈減小趨勢，單位面積切向剛度對相對誤差的影響高于單位面積法向剛度的影響。

2.5 最優(yōu)參數(shù)與模型驗證

為了使仿真單軸壓縮試驗中玉米秸稈所受到的臨界載荷接近物理試驗的935.4 N，采用Design-Expert軟件對回歸模型進行優(yōu)化分析，得到最優(yōu)仿真參數(shù)為：接觸半徑1.2 mm、單位面積法向剛度9.361×107N/m3和單位面積切向剛度9.845×107N/m3，預測理論相對誤差為2.94%。根據(jù)最優(yōu)仿真參數(shù)進行驗證試驗，玉米秸稈離散元模型在驗證試驗中臨界載荷為950.2 N，與物理試驗的相對誤差為1.58%，試驗結果與預測值基本一致，進一步證明了模型的準確性與可靠性。

3 結 論

本研究基于Plackett-Burman試驗設計和Box-Behnken響應面法，以單軸壓縮試驗法對玉米秸稈離散元模型進行標定，探究了玉米秸稈離散元模型的最優(yōu)參數(shù)組合，主要結論如下：

1)進行物理試驗測定得到玉米秸稈單軸壓縮試驗的臨界載荷為935.4 N。開展Plackett-Burman試驗，篩選出對臨界載荷試驗值與仿真值之間的相對誤差有顯著影響的因素為：接觸半徑、單位面積法向剛度和單位面積切向剛度。

2)進行最陡爬坡試驗縮小顯著性因素的水平范圍，再進行Box-Behnken響應曲面試驗，對試驗結果進行方差分析，得到二次多項式回歸方程。方差分析表明，回歸方程極顯著，且擬合度高。對回歸方程進行優(yōu)化求解，得到最優(yōu)參數(shù)為：接觸半徑1.2 mm、單位面積法向剛度9.361×107N/m3和單位面積切向剛度9.845×107N/m3。進行驗證試驗，驗證試驗的臨界載荷為950.2 N，與物理試驗935.4 N的相對誤差為1.58%，驗證了玉米秸稈離散元模型的仿真參數(shù)的可靠性。