任德志，張露籍，宮元娟，孟 軍，赫天一

（沈陽農(nóng)業(yè)大學(xué)a.農(nóng)學(xué)院/國家生物炭研究院/農(nóng)業(yè)農(nóng)村部生物炭與土壤改良重點實驗室，b.工程學(xué)院，沈陽 110161）

生物質(zhì)炭化還田是秸稈等農(nóng)林廢棄物資源化利用的可行途徑，在土壤改良和作物栽培領(lǐng)域備受關(guān)注[1-2]。目前生產(chǎn)中生物炭主要用于加工生物炭基肥料，施用生物炭基肥料可在不增加額外作業(yè)流程的情況下，實現(xiàn)生物炭的還田，并發(fā)揮其改良土壤、促進作物生產(chǎn)的積極作用[3]。當(dāng)前關(guān)于生物炭還田的研究多為探究炭土混合后對土壤及各類農(nóng)作物生長的影響[4-5]，而對于以何種方式高效還田及如何提高生物炭還田的機械化程度鮮有研究。離散元法（Discrete Element Method）是一種用于預(yù)測散裝固體行為的數(shù)值計算模擬技術(shù)，離散元模型已被廣泛用于農(nóng)業(yè)領(lǐng)域，它可以幫助縮短和改善農(nóng)業(yè)機械及其部件的設(shè)計過程[6-9]。利用離散元法進行仿真試驗，是反映生物炭炭基肥料還田過程的可行途徑，而進行離散元仿真試驗前，需對生物炭基肥料固體顆粒進行參數(shù)標(biāo)定。

國內(nèi)外已有相關(guān)生物質(zhì)材料方面的參數(shù)標(biāo)定理論與試驗方面的研究，侯俊銘等[10]以蓖麻休止角為響應(yīng)值，設(shè)計中心組合試驗得到了蓖麻接觸參數(shù)的最優(yōu)組合；PACHóN-MORALES 等[11]利用DEM 校準(zhǔn)框架，以木制纖維素生物質(zhì)的休止角、堆積密度和附著率為響應(yīng)值，得到了生物質(zhì)顆粒流動和與設(shè)備交互作用的離散元仿真模型；王韋韋等[12]利用EDEM 軟件得到了玉米秸稈粉料致密成型離散元仿真模型最優(yōu)化參數(shù)。然而值得注意的是，生物炭基肥料作為一種備受關(guān)注新興農(nóng)業(yè)投入品，尚未有離散元法參數(shù)標(biāo)定和實驗方面的研究。

本研究在前人研究的基礎(chǔ)上，提出一種真實試驗與仿真試驗結(jié)合的方法測定炭基肥顆粒的物理特性，以物理試驗為參照建立生物炭炭基肥顆粒離散元模型并標(biāo)定仿真參數(shù)，旨在為實現(xiàn)生物炭科學(xué)高效還田、仿真生物炭基肥料還田過程和設(shè)計生物炭基肥料專用還田機械提供新的理論依據(jù)。

1 材料與方法

1.1 試驗材料與物料特性

試驗材料為沈陽農(nóng)業(yè)大學(xué)國家生物炭研究院監(jiān)制的生物炭基肥料，該顆粒材料平均含水率為5%，百粒質(zhì)量5.26 g，球形度92%。從炭顆粒試驗樣品中隨機抓取百粒，測定記錄顆粒粒徑，重復(fù)8 次，炭基肥顆粒堆俯視圖和粒徑分布如圖1，2 mm粒徑占比35%；3 mm顆粒占比42%；4 mm顆粒占比23%。

圖1 生物炭基肥顆粒及粒徑分布Figure 1 Biochar-based fertilizer particles

1.2 堆積角的測量

堆積角是表征顆粒物料流動、摩擦等特性的宏觀參數(shù)，是仿真試驗中的關(guān)鍵響應(yīng)指標(biāo)。參考球形顆粒肥料的參數(shù)標(biāo)定研究[13-16]，選擇運用圓筒提升法測量顆粒堆積角（圖2）。

試驗所用裝置和材料包括:物理特性萬能試驗機、內(nèi)徑40 mm 套筒、游標(biāo)卡尺和載物平板。試驗過程為：圓筒置于平板中心位置，測量物理特性的886 粒炭基肥充滿于圓筒內(nèi)；試驗機夾具固定圓筒后，以0.1 m·s-1的速度將圓筒提起，使生物炭形成顆粒堆，靜止5 min，剔除滾出顆粒堆外的顆粒，利用游標(biāo)卡尺測量顆粒堆底部最大直徑d和顆粒堆最大高度h，生物炭炭基肥顆粒堆積角為：

試驗重復(fù)5次，結(jié)果如表1，最后計算得到堆積角平均值為23.65°。

1.3 接觸模型

離散元軟件ESSS Rocky DEM 提供了專用于預(yù)測農(nóng)業(yè)設(shè)備中顆粒行為所需的模型，仿真模型的選取直接影響仿真結(jié)果，模型包括重力模型、接觸模型以及熱模型。炭基肥顆粒之間的作用存在粘結(jié)現(xiàn)象，合理的選擇接觸模型可以得到與實際試驗相近的結(jié)果。本研究設(shè)置接觸參數(shù)如表2。

表2 接觸參數(shù)模型設(shè)置Table 2 Contact parameter model setup

1.3.1 法向力模型 DEM模型中的法向力模型需要滿足兩個條件：該力是一個反作用力；法向力模型滿足高能量耗散。本研究中法向力模型選擇Hertzian Spring Dashpot 模型，該模型與Linear Spring-dashpot模型區(qū)別在于Hertzian模型中法向力的彈性和阻尼分量是重疊的非線性函數(shù)[17]，其計算公式為：

式中：Fn為法向力值；為剛度系數(shù)；為阻尼系數(shù)；sn為接觸法線重迭量；為觸法線重迭量時間導(dǎo)數(shù)。

根據(jù)TSUJI 等[18]研究，Hertzian Spring Dashpot 模型中剛度系數(shù)阻尼系數(shù)的定義與Linear Spring-dashpot 模型中的定義類似，法向阻尼系數(shù)的值可通過黏性能量耗散與非彈性碰撞的能量耗散相匹配來確定，而非彈性碰撞的能量耗散則由恢復(fù)系數(shù)的值來確定，則剛度系數(shù)和阻尼系數(shù)為：

式中:E*為降低后的楊氏模量；R*為有效或等效半徑；η為阻尼比，其值與恢復(fù)系數(shù)有關(guān)；m*為接觸的有效質(zhì)量或等效質(zhì)量。其中計算公式為：

式中：m1和m2分別為相互接觸粒子的質(zhì)量，而m是與邊界接觸的粒子的質(zhì)量。

Hertzian Spring Dashpot 模型適用于球體顆粒，并且該模型適用于所有切向力模型，與3種附著力模型也有較好的適配性，綜上考慮選擇Hertzian Spring Dashpot 模型。

1.3.2 切向力模型 切向力模型用于計算接觸力的切向分量，包括Linear Spring Coulomb Limit 模型、Coulomb Limit 模型和Mindlin-Deresiewicz 模型。本研究選擇Linear Spring Coulomb Limit 模型，該模型中的切向力是彈性摩擦力，若切向力視為純彈性，其與時間t的關(guān)系式為：

式中：為模型切向力為前一瞬間切向力的值；Δsτ為粒子的切向相對位移；Kτ為切向剛度。

式（6）中的切向剛度與粒子材料相關(guān)。在此模型中，切向力不能超過庫倫極限，因此切向力的完整表達式為：

式中：為隨時間變化的接觸法向力；μ為摩擦系數(shù)，若在接觸處無滑動選擇靜摩擦系數(shù)μs，存在滑動選擇動摩擦系數(shù)μd。

在模擬計算炭基肥顆粒中引用Coulomb Limit 模型不適用于顆粒少精度高的仿真試驗；選擇Mindlin-Deresiewicz 模型則會禁用滾動阻力模型。綜上考慮選擇適用于所有模型的Linear Spring Coulomb Limit 模型。

1.3.3 附著力模型 上述模型中描述的正常接觸模型適合于模擬非粘性和干燥的顆粒材料，實際上，絕大多數(shù)顆粒狀材料都會含有不同的水分，所以會產(chǎn)生粒子與粒子、粒子與其他表明黏附的現(xiàn)象。為了捕捉這種行為，必須補充附著力模型，以準(zhǔn)確預(yù)測其流動特性。需要注意的是，顆粒之間的粘附力是應(yīng)力的一個函數(shù)，Rocky中的附著力模型通過將黏附力和接觸力縮放來捕捉這一現(xiàn)象。

本研究中選擇的附著力模型為JOHNSON等[19]提出的JKR模型，在JKR模型中，粒子間的接觸面積比Hertzian理論中所所預(yù)測的略大，這是表面能加入到系統(tǒng)總能量中的結(jié)果。選擇JKR模型后，法向力的附著力公式為：

式中：Fn,adh為法向附著力；Γ為模型參數(shù)，該系數(shù)的值與材料有關(guān)；E*為有效的楊氏模量；a為粒子與粒子或與表明之間的接觸半徑。

JKR 模型會影響法向和切向力模型，Hertzian Spring Dashpot 模型中的彈性部分需要修正，并以接觸半徑a的函數(shù)表達，即：

式中:R*為有效粒子半徑。

JKR 模型擁有夯實的理論基礎(chǔ)，被廣泛應(yīng)用于附著力模型中，由于表面能可以通過實驗測量，該模型不需要校準(zhǔn)并且可以用于模擬球形顆粒。在Rocky軟件中JKR模型需要與Hertzian Spring-dashpot 模型一起使用。

1.3.4 滾動阻力模型 滾動阻力是在建模中引入與粒子滾動運動相反的力矩時使用的通用名稱，通常將該力矩作為一種實用的方法來表示滾動球體上的非球形效應(yīng)或其他類型粒子上的表面不規(guī)則性效應(yīng)。本研究選擇的滾動阻力模型為Type C:Linear Spring Rolling Limit[20]。

Linear Spring Rolling Limit 是一個彈塑性模型，該模型將黏性阻尼項替代為滾動剛度值以更好的模擬滾動阻力行為，其中滾動剛度為：

式中：Kr為滾動剛度；Kτ為切向剛度；Rr為滾動半徑，定義為在給定時間內(nèi)連接粒子中心點和接觸點的向量，若滾動阻力為純彈性力，則滾動阻力為：

Rocky軟件中提供兩種滾動阻力模型，另一模型Type A 只適用于光滑滾動，由于生物炭顆粒之間存在黏連作用，所以滾動阻力模型選擇可配合附著力模型的Type C:Linear Spring Rolling Limit模型。

1.4 離散元仿真模型

Rocky軟件材料本構(gòu)提供了用于描述農(nóng)業(yè)顆粒行為所需的數(shù)學(xué)模型，通過Rocky建立生物炭粒子?；谝淹瓿傻念w粒物理特性試驗得到生物炭顆粒球形度為92%，根據(jù)離散元法理論[21]，當(dāng)顆粒非球形度小于10%時，設(shè)置球體顆粒對仿真分析結(jié)果影響極小，所以建立Biochar粒子，形狀類型設(shè)為球體顆粒，粒徑分布參考物理特性試驗設(shè)置為:4 mm 直徑的粒子數(shù)目占總顆粒數(shù)25%、3 mm 直徑的粒子數(shù)目占總顆粒數(shù)40%、2 mm 直徑的粒子數(shù)目占總顆粒數(shù)35%。在SolidWorks 中建立圓筒和鋼板的三維模型，導(dǎo)出為STL 格式并導(dǎo)入Rocky 軟件中（圖3）。

圖3 仿真模型Figure 3 Simulation model

為更好地對照真實試驗，參考設(shè)置圓筒上升速度為0.1 m·s-1，粒子生成數(shù)目886個，生成時間2 s，總仿真時長8 s。為記錄設(shè)置不同的接觸參數(shù)組合下所產(chǎn)生的炭顆粒堆積角，在每次仿真結(jié)束后創(chuàng)建在t=8 s 處粒子的交叉散點圖，它可以觀測到在某一瞬間粒子堆積在X-Y平面上的投影（圖4），根據(jù)交叉散點圖可確定本次接觸參數(shù)組合所得到的堆積角，計算方法同式（1）。

圖4 交叉散點圖Figure 4 Cross-scatter plot

2 試驗設(shè)計與結(jié)果分析

2.1 Placket-Burman試驗

由于物料的接觸參數(shù)中并非所有參數(shù)都對堆積角有顯著影響，所以，本研究設(shè)計Placket-Burman 試驗，篩選出關(guān)鍵仿真接觸參數(shù)，P-B 試驗因素水平如表3。

表3 接觸參數(shù)模型設(shè)置Table 3 Contact parameter model setup

P-B 試驗方案及結(jié)果如表4，其中H、J、K、L 為虛擬參數(shù)，即空白列。對試驗結(jié)果進行方差分析如表5，結(jié)果表明炭-炭滾動摩擦系數(shù)（C）、炭-鋼滾動摩擦系數(shù)（F）和炭顆粒表面能（G）對生物炭顆粒的堆積角影響顯著。

表4 Plackett-Burman試驗方案及結(jié)果Table 4 Plackett-Burman test protocol and results

表5 Plackett-Burman試驗方差分析Table 5 Variance analysis of Plackett-Burman test

2.2 最陡爬坡試驗

對選擇的顯著影響參數(shù)進行最陡爬坡試驗以確定顯著性參數(shù)最優(yōu)值區(qū)間，試驗方案及結(jié)果如表6，結(jié)果表明:隨著炭-炭滾動摩擦系數(shù)（C）、炭-鋼滾動摩擦系數(shù)（F）和炭顆粒表面能（G）的值逐漸增大，仿真堆積角逐漸增大，相對誤差先減小再增大，其中第3組試驗的相對誤差最小，則選擇第3組顯著性參數(shù)值為中心水平，建立回歸模型求解參數(shù)組合最優(yōu)值。

表6 最陡爬坡試驗方案及結(jié)果Table 6 Steepest climb test protocol and results

2.3 Box-Behnken響應(yīng)面試驗

根據(jù)最陡爬坡試驗結(jié)果，選擇第3 組試驗因素值為中心水平值，第2 組、第4 組試驗因素值為-1 和1 水平值，設(shè)計Box-Behnken響應(yīng)面試驗，試驗方案和結(jié)果如表7。

表7 Box-Behnken試驗方案及結(jié)果Table 7 Box-Behnken test protocol and results

通過Design Expert試驗分析軟件對試驗結(jié)果進行顯著性檢驗，得到關(guān)于堆積角的方差分析如表8。通過檢驗結(jié)果可發(fā)現(xiàn)，模型R1的顯著值p＜0.001，并且失擬項p＞0.05，說明本次多因素試驗的回歸方程檢驗達到高度顯著，并且擬合度較好，因此可以對本次試驗數(shù)據(jù)結(jié)果進行進一步的分析和優(yōu)化。去除不顯著項后得到的多項式回歸方程為：

表8 Box-Behnken試驗方差分析Table 8 Variance analysis of Box-Behnken test

式中:A為炭-炭滾動摩擦系數(shù)（C）編碼值；B為炭顆粒表面能（G）編碼值；C為炭-鋼滾動摩擦系數(shù)（F）編碼值。

將Design Expert 中得到的多項式回歸方程導(dǎo)入繪圖軟件Origin 的矩陣方程中，通過3D 映射曲面圖繪制出三因素對堆積角的響應(yīng)面圖如圖5。通過方差分析和響應(yīng)面圖可知，炭-炭滾動摩擦系數(shù)（C）和炭顆粒表面能（G）單獨對堆積角影響顯著；炭顆粒表面能（G）和炭-鋼滾動摩擦系數(shù)（F）的交互影響對堆積角影響顯著，并且當(dāng)炭顆粒表面能（G）增大時，堆積角越趨近于真實試驗值23.65°。綜上可知炭顆粒表面能（G）對堆積角的影響最為關(guān)鍵。

圖5 多因素交互作用的響應(yīng)面圖Figure 5 Multi-factor response surface diagram

通過Design-Expert 軟件中試驗數(shù)據(jù)參數(shù)優(yōu)化模塊，設(shè)置邊界條件為堆積角趨近于目標(biāo)值23.65°，得到最優(yōu)參數(shù)組合為:炭-炭滾動摩擦系數(shù)0.36，炭顆粒表面能0.30，炭-鋼滾動摩擦系數(shù)0.20。

2.4 驗證試驗

為驗證所得最優(yōu)化接觸參數(shù)值的準(zhǔn)確性，以炭-炭滾動摩擦系數(shù)0.36、炭顆粒表面能0.30、炭-鋼滾動摩擦系數(shù)0.20、其他接觸參數(shù)取中間水平值（炭-炭恢復(fù)系數(shù)0.5、炭-炭靜摩擦系數(shù)0.65、炭-鋼恢復(fù)系數(shù)0.6、炭-鋼靜摩擦系數(shù)0.7，其余設(shè)置不變進行離散元仿真試驗，測得堆積角為23.48°，與實際試驗測定堆積角23.65°的相對誤差為0.72%。仿真試驗圖與實際試驗圖對比如圖6，兩者無明顯差異。

圖6 試驗對比圖Figure 6 Test comparison figure

3 討論與結(jié)論

國內(nèi)關(guān)于生物炭的研究基本來源于由陳溫福院士帶領(lǐng)成立的生物炭工程技術(shù)研究中心，其研制的Ⅰ級生物炭適用于東北地區(qū)土壤，但該團隊未對生物炭如何還田以及如何更好的還田做出詳細的研究。關(guān)于離散元參數(shù)標(biāo)定的研究已經(jīng)趨于成熟，與其他研究相比，生物炭還田的仿真過程類似于肥料顆粒入土，但兩者物理特性的差異使得需要對生物炭顆粒進行參數(shù)標(biāo)定，以為后續(xù)的生物炭還田技術(shù)研究提供理論參考，其結(jié)果的準(zhǔn)確性直接影響后續(xù)對于生物炭還田工作過程仿真的真實性。

本研究結(jié)果表明，基于實際物理實驗堆積角測定值23.65°，標(biāo)定生物炭顆粒接觸參數(shù)。設(shè)定離散元接觸模型，通過Placket-Burman 試驗篩選得到對生物碳顆粒堆積角有顯著影響的參數(shù):炭-炭滾動摩擦系數(shù)、炭顆粒表面能和炭-鋼滾動摩擦系數(shù)。根據(jù)最陡爬坡試驗結(jié)果確定參數(shù)的最優(yōu)質(zhì)區(qū)間，設(shè)計Box-Behnken響應(yīng)面試驗建立堆積角與3個顯著參數(shù)的回歸模型。顯著性檢驗和響應(yīng)面分析表明:炭-炭滾動摩擦系數(shù)和炭顆粒表面能對堆積角均有顯著影響，并且炭顆粒表面能和炭-鋼滾動摩擦系數(shù)之間存在明顯交互作用。根據(jù)最優(yōu)化參數(shù)組合:炭-炭滾動摩擦系數(shù)0.36、炭顆粒表面能0.30、炭-鋼滾動摩擦系數(shù)0.20，重新設(shè)置接觸參數(shù)并進行驗證試驗。試驗表明，在此接觸參數(shù)組合下仿真得到的堆積角23.48°，與實際測定值相對誤差為0.72%，所標(biāo)定的生物炭顆粒離散元模型參數(shù)真實可靠。