基于DEM的玉米秸稈離散元模型參數標定*

童世合，邵明璽，曹猛，趙軍杰

(1. 青海省特種設備檢驗檢測院，西寧市，810005； 2. 青海大學，西寧市，810016)

0 引言

玉米秸稈是畜牧飼料的主要來源之一，也是一種重要的生物質能源[1]，玉米秸稈的深度加工對于飼料增收和生物質能源發展具有深遠意義[2-3]。離散元法(DEM)用于研究離散物料在外部載荷下的運動學特性[4]，具有設計效率高、仿真精確度高和實際誤差低等優點，尤其在農業領域得到廣泛應用。

秸稈加工機械關鍵部件在設計時存在模具制作成本高、設計周期長、操作不靈活等問題，秸稈與機構之間的相互作用特性和規律難以通過常規技術手段進行觀測，機構設計缺少可靠依據，導致已有的結構可靠性差、作業效果不理想，借助DEM法建立秸稈精細化仿真模型，開展機構的參數優化設計，有效提高秸稈生產效率。然而，秸稈離散元仿真模型在外形和材質等方面仍與實際具有一定差異，通過力學試驗獲得彈性模量、摩擦系數等關鍵參數通常不能夠直接應用，需要對相關參數重新進行標定。

在農業物料仿真參數標定方面，國內外研究人員基于EDEM軟件開展了大量研究，EDEM是一款基于離散元法的顆粒物料通用仿真軟件。Adajar等[5]開展了油菜、玉米、小麥等5種作物的秸稈剪切特性標定，與實際值相對誤差小于2%；Schramm等[6]利用EDEM開展了小麥秸稈彈性模量和泊松比標定，與物理試驗結果的相對誤差僅為3.11%；Fang等[7]利用Plackett-Burman試驗對玉米碎秸稈的堆積角進行了仿真參數的標定，結果表明仿真試驗和物理試驗的平均相對誤差僅為0.29%，說明利用EDEM進行參數標定獲得的結果準確性好，具有較高的精度；張濤等[8]為了提高玉米秸稈揉碎過程的仿真精度，開展了秸稈—揉碎機、秸稈—秸稈接觸參數的標定研究。

上述研究驗證了離散元仿真方法的可靠性和準確性，但仿真過程中為了提高運算速度，不同程度地簡化了秸稈模型，導致模型精細化程度不足。為此，本文針對秸稈外表皮質地堅硬、纖維排列規律且均勻、內瓤木質部質地偏軟、組織緊密不規則的顯著差異，建立了反映真實秸稈有莖節式樣和無莖節式樣的外表皮和內瓤性狀的離散元仿真模型，結合單軸壓縮試驗確定仿真參數，并通過剪切試驗進行驗證，為后續秸稈加工機構的研制提供數據參考。

1 材料與方法

1.1 Hertz接觸理論

Hertz接觸理論是EDEM軟件計算顆粒—顆粒和顆粒—幾何體接觸特性的理論基礎，通過Hertz-Mindlin無滑移模型、Hertz-Mindlinwithbonding模型和Hertz-Mindlin JKR模型等基礎模型建立真實情況下兩物體之間的接觸規則[9-11]。考慮到秸稈內部組織間的連接方式和表現出的各向異性特性，采用Hertz-Mindlinwithbonding將秸稈模型中的顆粒進行粘結，顆粒與顆粒之間接觸作用如圖1所示，法向、切向應力應變方程如式(1)所示。

(1)

式中：A——顆粒之間的接觸面積，mm2，A=πRB2；

J——顆粒的轉動慣量，kg·m2，J=πRB4/4；

RB——顆粒的接觸半徑，mm；

vn——顆粒法向運動速度，m/s；

vt——顆粒切向運動速度，m/s；

Fn——顆粒受到的法向作用力，N；

Ft——顆粒受到的切向作用力，N；

Mn——顆粒法向力矩，N·m；

Mt——顆粒切向力矩，N·m；

Sn——法向剛度系數；

St——切向剛度系數。

圖1 球型顆粒接觸示意圖

式(1)反映了顆粒所受作用力與顆粒自身屬性的關系，法向剛度系數、顆粒的彈性模量和半徑相關，切向剛度系數、顆粒的剪切模量和半徑相關，綜上所述，顆粒粘結力的大小主要取決于法向剛度系數、切向剛度系數、臨界法向應力和臨界切向應力。

1.2 離散元仿真模型的建立

本文基于EDEM提出了一種快速構建精細化散粒物料顆粒的方法，具體構建方法如下。

1) 規定顆粒生成規則。考慮到秸稈外表皮和內瓤在形態和結構上存在較大差異，用于單軸壓縮試驗和剪切試驗的秸稈直徑為30 mm、長度為30 mm，外表皮模型采用直徑為3 mm球型顆粒建立，內瓤模型采用0.5～1.0 mm球型顆粒建立。外表皮模型采用“分層疊加”的方法進行建立，秸稈直徑根據測量結果設置為30 mm，因此每個底層顆粒的質心位置均勻分布在直徑為30 mm的圓上，在空間直角坐標系中捕獲顆粒的質心在X和Y軸上的坐標，Z軸坐標則以間隔4 mm 遞增，制定顆粒坐標規則(部分)如圖2(a)所示，將表中所有數據復制到Meta-particle面板中，建立外表皮模型如圖2(b)所示。內瓤模型采用顆粒工廠模板內靜態生成的方法進行建立，顆粒直徑遵循正態分布規則，能夠保證區域內顆粒填充密度與實際相符，建立內瓤模型如圖2(c)所示。

(a) 外表皮模型生成規則表(部分)

(b) 秸稈外表皮仿真模型(左圖為無莖節，右圖為有莖節)

2) 設置粘結參數。在EDEM中，合理的粘結參數是保證仿真模型成功建立的重要前提，粘結參數數值過大，會導致模型發生“爆炸”，粘結參數數值過小，模型出現坍縮、變形和塌陷等不穩定現象，通過反復調試，觀測生成的仿真模型的穩定性，最終確定了外表皮—外表皮、內瓤—內瓤和外表皮—內瓤的粘結參數，如表1所示。

表1 顆粒粘結參數設置Tab. 1 Particle bonding parameter settings

3) 固定顆粒運動狀態為靜止，在計算域內添加一個平面，設置為Virtual并創建顆粒工廠，生成顆粒數量為1，生成顆粒類型為Meta-particle，生成方向沿Z軸，生成外表皮模型后，在中心位置生成圓柱形顆粒工廠，生成方式為填充整個區域，構建秸稈式樣完整模型如圖2(d)所示。

1.3 仿真環境參數設置

為了提高仿真精確性，防止漏算，建立加載件模具的三維模型，劃分網格后共計1 080個網格單元，導入EDEM軟件中，加載件壓縮速度為2 mm/s，與實際壓縮試驗條件一致。合適的時間步長是成功建立顆粒之間粘結鍵的關鍵參數，時間步長過大，顆粒運動誤差計算結果精度差，導致粘結鍵生成失敗，顆粒將脫粒粘結鍵向四周散開，從而影響模型精度；時間步長過小會增大運算總時長，仿真效率低，綜合考慮計算精度和仿真效率，時間步長設置為瑞麗步長的5%。網格大小為顆粒半徑的2.5倍，共計8 489 664個網格。國內外研究人員針對玉米秸稈與加載件的接觸參數進行了大量工作，為本研究提供了充足的數據參考，各項材料屬性參數[9, 12-14]設置如表2所示，玉米無莖節秸稈單軸壓縮仿真試驗和萬能實驗機試驗如圖3所示。

表2 材料屬性參數設置Tab. 2 Material property parameter settings

(a) 仿真式樣

(b) 壓縮式樣

1.4 Plackett-Burman仿真試驗

Plackett-Burman試驗能夠在多個因素同時影響的條件下，選取每個因素的高低極限水平進行單因素試驗，從而得到影響最為顯著的因素[15-17]。本研究針對玉米秸稈外表皮和內瓤的仿真參數開展標定試驗。需要標定的因素包括外表皮—內瓤碰撞恢復系數A、外表皮—內瓤靜摩擦系數B、外表皮—內瓤滑動摩擦系數C、內瓤—內瓤碰撞恢復系數D、內瓤—內瓤靜摩擦系數E和內瓤—內瓤滑動摩擦系數F共6個仿真參數進行篩選，試驗水平如表3所示。仿真結束后和移除載荷后的玉米秸稈式樣如圖4所示。

表3 Plackett-Burman試驗水平表Tab. 3 Plackett-Burman test level table

(a) 仿真式樣

(b) 壓縮式樣

1.5 最陡爬坡仿真試驗

由Plackett-Burman仿真試驗結果可知，對玉米秸稈抗壓力表現出顯著影響的因素有外表皮—內瓤靜摩擦系數、內瓤—內瓤碰撞恢復系數和內瓤—內瓤靜摩擦系數。為了使上述3個試驗因素水平快速逼近最優取值范圍，開展最陡爬坡試驗。

1.6 Box-Behnken正交試驗

根據最陡爬坡仿真試驗結果，以外表皮—內瓤靜摩擦系數B、內瓤—內瓤碰撞恢復系數D和內瓤—內瓤靜摩擦系數E3個參數為試驗因素，以顆粒所受抗破壞力N為響應指標開展3因素3水平正交試驗，共17組，探明上述3個試驗因素對響應指標的影響規律，并進行最優參數優化[18-20]。

1.7 最優參數組合驗證試驗

根據正交試驗結果，利用Design-Expert軟件進行最優參數優化設計，設計規則如式(2)所示。為了驗證仿真結果的準確性，將最優參數組合作為新的玉米秸稈仿真參數輸入，分別進行剪切仿真試驗和物理試驗，對比相同參數組合下有莖節玉米秸稈式樣的抗剪切力的大小，如圖5所示。

(2)

圖5 玉米秸稈剪切驗證試驗

2 結果與討論

2.1 Plackett-Burman仿真試驗結果

Plackett-Burman仿真試驗結果如表4所示，方差分析結果如表5所示。該模型的P值為0.022 5<0.05，表明模型對顆粒抗破壞力的影響為極顯著，因素B、D和E的P值<0.05，說明外表皮—內瓤靜摩擦系數、內瓤—內瓤碰撞恢復系數和內瓤—內瓤靜摩擦系數對顆粒抗破壞力的影響為極顯著，而其他因素的影響均為不顯著。

因此，將因素B、D和E作為最陡爬坡試驗和Box-Behnken正交試驗的變量，同時取外表皮—內瓤靜摩擦系數A為0.5、外表皮—內瓤滑動摩擦系數C為0.055、內瓤—內瓤滑動摩擦系數F為0.055，均為其區間的中位數。

表5 Plackett-Burman試驗分析結果Tab. 5 Plackett-Burman test analysis results

2.2 最陡爬坡仿真試驗結果

最陡爬坡仿真試驗結果如表6所示，第2組試驗和第3組試驗的顆粒抗破壞力分別為69.55 N和80.34 N，秸稈單軸壓縮試驗的抗破壞力為75.54 N，在第2組和第3組試驗結果的區間范圍內，因此將這兩組試驗因素的水平作為正交試驗自變量因素的+1和-1水平。

表6 最陡爬坡試驗結果Tab. 6 Steepest climbing test results

2.3 Box-Behnken正交試驗結果

為了獲得顯著影響因素(外表皮-內瓤靜摩擦系數、內瓤-內瓤碰撞恢復系數和內瓤-內瓤靜摩擦系數)的最優參數組合，正交試驗安排表與顆粒抗破壞力結果如表7所示，方差分析結果如表8所示，得到變量的回歸方程如式(3)所示。

F=-281.83+446.57B+895.38D+24.34E-

1 326.31BD+7.29BE+2.19DE+

259.73B2-716.43D2-1.44E2

(3)

由表8可知，因素擬合度為0.99，說明該模型具有較高的可信度，整個模型的P值<0.000 1，說明該模型的影響為極顯著，能夠可靠反映自變量與因變量之間的關系，模型失擬項的P值=0.475 8>0.10，說明該模型中不存在無關因素的影響。另外，自變量B、D和E的P值均小于0.05，說明選取的3個因素對顆粒抗破壞力的影響均為極顯著，交互因素BD、二次項因素D2和E2的影響也均為極顯著，其他項均為不顯著。

表7 正交試驗安排與結果Tab. 7 Orthogonal test arrangement and results

表8 正交試驗方差分析結果Tab. 8 Results of ANOVA

剔除了影響不顯著的因素后，各項顯著性因素的影響程度次序為D、C2、BD、E、E2、B，根據模型的回歸方程得到各自變量的交互作用對顆粒抗破壞力的響應面分析如圖6所示。

(a) F=(B,D,0)

(b) F=(B, 0, E)

(c) F=(0, D, E)

由圖6(a)可知，當外表皮—內瓤靜摩擦系數為中間水平時，顆粒的抗破壞力隨著內瓤—內瓤碰撞恢復系數的增大而平緩增大；當其處在低水平時，抗破壞力則隨著內瓤—內瓤碰撞恢復系數的增大而迅速增大；當其處于高水平時，抗破壞力隨著內瓤—內瓤碰撞恢復系數的增大先增大，后減小，變化幅度較為明顯；由圖6(b)可知，顆粒的抗破壞力隨著外表皮—內瓤靜摩擦系數和內瓤—內瓤靜摩擦系數的增大而增大，變化幅度明顯；由圖6(c)可知，顆粒的抗破壞力隨內瓤—內瓤碰撞恢復系數和內瓤—內瓤靜摩擦系數的增大而增大，變化趨勢明顯，2個因素對抗破壞力的影響程度基本相同。

由式(3)得到實現抗破壞力最小值的參數組合為外表皮—內瓤靜摩擦系數0.28、內瓤—內瓤碰撞恢復系數0.29、內瓤—內瓤靜摩擦系數0.23，此時加載件受力為76 N，與目標值誤差為0.46 N。

2.4 驗證試驗

根據最優仿真參數組合進行仿真試驗，并隨機取玉米有莖節秸稈在萬能實驗機進行剪切試驗的5組抗破壞力曲線與仿真試驗中顆粒抗破壞力曲線進行對比，結果如圖7所示。玉米有莖節秸稈離散元模型在驗證試驗中的抗破壞力為132.29 N，在物理試驗中抗破壞力平均值為131.36 N，相對誤差僅為0.93 N，說明仿真試驗結果與物理試驗結果基本一致，驗證了玉米秸稈精細化離散元仿真模型的可靠性與準確性。

圖7 玉米有莖節秸稈剪切驗證試驗對比結果

3 結論

針對玉米秸稈外表皮與內瓤存在較大差異，導致離散元仿真模型精確性不足的問題，提出了一種玉米無莖節和有莖節秸稈的精細化離散元仿真模型的建立方法，探究了該模型下的仿真參數最優組合。

1) 萬能實驗機物理試驗測得玉米無莖節秸稈的單軸壓縮抗破壞力為75.54 N，為了獲得與實際值相近的仿真參數，利用Plackett-Burman試驗和Box-Behnken回歸正交試驗開展6個參數的仿真試驗，結果表明對抗破壞力具有顯著影響的因素包括外表皮—內瓤滑動摩擦系數、內瓤—內瓤碰撞恢復系數和內瓤—內瓤靜摩擦系數。

2) 根據方差分析結果，進行最優參數優化，得到滿足要求的參數組合為外表皮—內瓤靜摩擦系數0.28、內瓤—內瓤碰撞恢復系數0.29、內瓤—內瓤靜摩擦系數0.23，此時顆粒抗破壞力的預測值為76 N，與實際值75.54 N的誤差為0.46 N。

3) 建立玉米有莖節秸稈離散元精細化仿真模型，利用最優參數組合開展剪切驗證試驗，仿真試驗與物理試驗的抗剪切力分別為132.29 N和131.36 N，誤差僅為0.93 N，表明本文提出的離散元仿真模型建立方法和仿真結果可靠有效，能夠為其他作物離散元仿真研究提供參考。