謝方平，吳正陽，王修善，劉大為，鄔 備，張正中

基于無側(cè)限抗壓強(qiáng)度試驗(yàn)的土壤離散元參數(shù)標(biāo)定

謝方平1,2，吳正陽1，王修善1,2，劉大為1,2，鄔 備1,2，張正中1

（1.湖南農(nóng)業(yè)大學(xué)機(jī)電工程學(xué)院，長(zhǎng)沙 410128；2.智能農(nóng)機(jī)裝備湖南省重點(diǎn)實(shí)驗(yàn)室，長(zhǎng)沙 410128）

為標(biāo)定基于Edinburgh Elasto-Plastic Adhesion (EEPA)模型的土壤離散元參數(shù)，該研究通過單軸密閉壓縮試驗(yàn)和無側(cè)限抗壓強(qiáng)度試驗(yàn)，以無側(cè)限抗壓強(qiáng)度σ為黏性指標(biāo)，軸向應(yīng)變ε為塑性指標(biāo)，軸向壓力-軸向應(yīng)變曲線特征參數(shù)（，）為彈性指標(biāo)，基于響應(yīng)面標(biāo)定了土壤離散元仿真參數(shù)。結(jié)合文獻(xiàn)與實(shí)際情況確定參數(shù)試驗(yàn)范圍，應(yīng)用Plackett-Burman試驗(yàn)對(duì)7個(gè)初始參數(shù)進(jìn)行篩選，發(fā)現(xiàn)塑性變行比λ和表面能Δ對(duì)黏塑性指標(biāo)影響顯著。根據(jù)試驗(yàn)結(jié)果，進(jìn)行2次Central Composite試驗(yàn)，建立σ、ε與黏塑性參數(shù)的二次回歸模型和σ、ε、、與彈性參數(shù)二次回歸模型，并以實(shí)測(cè)值為目標(biāo)優(yōu)化求解，獲得的優(yōu)化參數(shù)組合為：塑性變形比為0.36、表面能為15.6 J/m2、恢復(fù)系數(shù)為0.37、靜摩擦系數(shù)為0.6、滾動(dòng)摩擦系數(shù)為0.26、黏性分支指數(shù)為4.24、切向剛度因子為0.52。最后將該參數(shù)組合下的仿真值與實(shí)測(cè)值進(jìn)行對(duì)比驗(yàn)證，發(fā)現(xiàn)仿真ε和與實(shí)測(cè)值無顯著差異，σ和與實(shí)測(cè)值存在較大差異。結(jié)果表明基于響應(yīng)面法標(biāo)定的EEPA模型參數(shù)可表征試驗(yàn)土壤的軸向塑性變形和3%～45%軸向應(yīng)變內(nèi)的應(yīng)力-應(yīng)變行為。

0 引 言

離散元法（Discrete Element Method，DEM）是一種將介質(zhì)整體視為若干顆粒單元集合的數(shù)值模擬方法，利用牛頓第二定律和接觸模型，以單個(gè)顆粒的運(yùn)動(dòng)學(xué)和動(dòng)力學(xué)行為描述整體特征，被廣泛用于谷物脫粒[1-2]與清選[3-5]、土壤與機(jī)具相互作用[6-10]等領(lǐng)域。研究對(duì)象外形各異、尺寸不一、粒度分布復(fù)雜、孔隙度不均勻，DEM模擬這些復(fù)雜的物理特性存在一定的誤差，該誤差會(huì)導(dǎo)致仿真對(duì)象的物理特性與實(shí)際情況存在較大差異[11-14]，為在DEM中再現(xiàn)研究對(duì)象的實(shí)際物理特性，需要確定合理的DEM參數(shù)以減小差異。因此，DEM參數(shù)的準(zhǔn)確標(biāo)定是成功仿真的關(guān)鍵，也是大部分離散元仿真的重要組成部分。近年來，DEM相關(guān)文獻(xiàn)出版數(shù)量激增[15]，詳細(xì)的DEM參數(shù)標(biāo)定方法研究不斷增多[16]。

在黏性介質(zhì)的DEM參數(shù)標(biāo)定研究中，一定固結(jié)應(yīng)力條件下，無側(cè)限抗壓強(qiáng)度能直觀地反映固結(jié)介質(zhì)間的黏聚力水平[17-18]。土壤的無側(cè)限抗壓強(qiáng)度是指固結(jié)土壤在無側(cè)限條件下，抵抗軸向壓力的極限強(qiáng)度，是影響土壤抗破壞能力的重要指標(biāo)[17]。土壤的無側(cè)限抗壓強(qiáng)度試驗(yàn)一般從單軸密閉壓縮試驗(yàn)中獲取固結(jié)試樣。

Karkala等[17]在無側(cè)限抗壓強(qiáng)度仿真試驗(yàn)中，以100 mm/s的恒定速率（應(yīng)變率=4 s-1）加載，直至試樣失效；試驗(yàn)中記錄壓盤的壓力和試樣端面的變形，分析發(fā)現(xiàn)JKR表面能是無側(cè)限抗壓強(qiáng)度的極敏感參數(shù)，在JKR(Johnson-Kendall-Roberts, JKR）表面能參數(shù)的可行域內(nèi)等分5組并逐一嘗試，選擇一組與實(shí)測(cè)值誤差最小的JKR表面能參數(shù)作為標(biāo)定結(jié)果。類似地，王憲良等[19]發(fā)現(xiàn)休止角、黏聚力和內(nèi)摩擦角對(duì)顆粒半徑、顆粒間靜摩擦系數(shù)和滾動(dòng)摩擦系數(shù)的敏感度都比其他參數(shù)高出一個(gè)數(shù)量級(jí)，利用響應(yīng)面法得到這3項(xiàng)參考指標(biāo)與3個(gè)主要參數(shù)的回歸模型并優(yōu)化求解得到最優(yōu)的參數(shù)組合。Xia等[20]為模擬煤礦機(jī)械與煤礦之間的相互作用，根據(jù)黏性物料特性使用Johnson-Kendall-Roberts（JKR）模型，通過Plackett-Burman（PB）試驗(yàn)獲得對(duì)休止角敏感的DEM參數(shù)，并將基于響應(yīng)面法原理的Box-Behnken設(shè)計(jì)應(yīng)用于標(biāo)定。從已有文獻(xiàn)發(fā)現(xiàn)，在多數(shù)離散元仿真研究中，僅對(duì)幾個(gè)敏感參數(shù)進(jìn)行標(biāo)定是DEM參數(shù)標(biāo)定常見的方法。其余未標(biāo)定參數(shù)，如摩擦系數(shù)[21]、恢復(fù)系數(shù)[22]等，通常采用直接測(cè)量法確定，但JKR表面能、塑性變形比、粘結(jié)強(qiáng)度等參數(shù)難以通過直接測(cè)量法確定。

Coetzee[23]基于Hertze-Mindlin（no slip）模型以休止角為響應(yīng)，對(duì)顆粒間摩擦系數(shù)進(jìn)行標(biāo)定。Ucgul等[24]將線性黏附內(nèi)聚模型整合到線性Hysteretic Spring模型，并通過改變土壤靜摩擦系數(shù)、滾動(dòng)摩擦系數(shù)和恢復(fù)系數(shù)，獲得與實(shí)際休止角和實(shí)測(cè)累計(jì)能量-貫入深度曲線相近的DEM參數(shù)。李俊偉等[25]以不同含水率的黑黏土為研究對(duì)象，基于JKR模型，通過休止角試驗(yàn)對(duì)顆粒間JKR表面能、靜摩擦系數(shù)、滾動(dòng)摩擦系數(shù)和恢復(fù)系數(shù)進(jìn)行標(biāo)定。Thakur等[26]首次提出Edinburgh Elasto-Plastic Adhesion (EEPA)接觸模型，EEPA模型在Hertz接觸理論的基礎(chǔ)上，包含了顆粒塑性和黏性。Janda等[18]進(jìn)一步證明EEPA模型適宜模擬黏重、可塑性強(qiáng)的農(nóng)田土壤?；贘KR模型的黏彈性物料[17,25]、基于Hertze-Mindlin（no slip）模型的無黏性物料[14-16]和基于線性Hysteretic Spring模型的塑性變形物料[7,11,24]的DEM參數(shù)標(biāo)定已有許多研究，而基于EEPA模型，對(duì)具備黏彈性的非線性塑性物料DEM參數(shù)標(biāo)定鮮有報(bào)道。

Roessler等[16]認(rèn)為，對(duì)于不同的應(yīng)用，需要不同的標(biāo)定試驗(yàn)，單獨(dú)的宏觀響應(yīng)（例如堆積角）可能導(dǎo)致多個(gè)可行的參數(shù)組合，而多個(gè)參數(shù)組合可能無法用于單個(gè)響應(yīng)以外的其他應(yīng)用。本文以具有黏性且易產(chǎn)生塑性變形的土壤為研究對(duì)象，以軸向應(yīng)變、無側(cè)限抗壓強(qiáng)度和壓縮土壤時(shí)的非線性應(yīng)力-應(yīng)變特征為響應(yīng)，通過Plackett-Burman和Central Composite試驗(yàn)進(jìn)行EEPA模型的全部接觸參數(shù)標(biāo)定，并通過仿真與實(shí)測(cè)值進(jìn)行對(duì)比驗(yàn)證，以期為基于EEPA模型的離散元仿真研究提供參數(shù)依據(jù)和可行的參數(shù)標(biāo)定方法

1 材料與方法

本研究以單軸密閉壓縮試驗(yàn)獲得的固結(jié)土壤試樣為材料，通過無側(cè)限抗壓強(qiáng)度試驗(yàn)，對(duì)EEPA模型參數(shù)進(jìn)行標(biāo)定。參數(shù)標(biāo)定時(shí)，首先以恢復(fù)系數(shù)、靜摩擦系數(shù)μ、滾動(dòng)摩擦系數(shù)μ、表面能Δ、塑性變形比λ、黏結(jié)分支指數(shù)和切向剛度因子k為影響因素進(jìn)行Plackett-Burman仿真試驗(yàn)，得到對(duì)軸向應(yīng)變ε和無側(cè)限抗壓強(qiáng)度σ影響顯著的因素，并以其為變量進(jìn)行二次正交旋轉(zhuǎn)組合仿真試驗(yàn)，得到ε和σ關(guān)于影響因素的回歸方程，并以實(shí)測(cè)值為目標(biāo)進(jìn)行求解。在此基礎(chǔ)上，以、μ、、k為影響因素進(jìn)行4因素Central Composite仿真試驗(yàn)，得到4個(gè)影響因素對(duì)ε、σ和單軸密閉壓縮階段軸向壓力-軸向應(yīng)變曲線特征參數(shù)的方程組。在Design Expert 8.0.5軟件中以實(shí)測(cè)的ε、σ和曲線特征參數(shù)為目標(biāo)，對(duì)方程組求解，選擇Desirability值最接近1(可靠度最高)的一組解。

1.1 土壤基本參數(shù)

利用圓孔篩篩分干土，取直徑為2～3 mm的土壤顆粒。取樣地點(diǎn)為湖南農(nóng)業(yè)大學(xué)工程實(shí)訓(xùn)中心試驗(yàn)土槽，所用土壤來自湖南農(nóng)業(yè)大學(xué)試驗(yàn)田，土質(zhì)黏重、孔隙度大、可塑性強(qiáng)。為模擬土壤潤(rùn)濕過程，將土樣放置于土槽內(nèi)，噴灑適量水并覆蓋地膜防止土樣失水太快。放置12 h后進(jìn)行含水率測(cè)定和后續(xù)試驗(yàn)，測(cè)定的試樣含水率為26.25%。通過比重瓶試驗(yàn)測(cè)量土壤密度為2 670 kg/m3。由于可塑狀態(tài)黏土的泊松比一般大于粉質(zhì)黏土和沙土[8-10]，本文選取土壤泊松比為0.38。應(yīng)力應(yīng)變測(cè)試三軸剪切儀測(cè)量土壤剪切模量為1 MPa。

1.2 土壤物理特性參數(shù)

1.2.1 塑性和彈性指標(biāo)

如圖1所示，單軸密閉壓縮試驗(yàn)時(shí)，內(nèi)徑40 mm、外徑60 mm、高度150 mm的有機(jī)玻璃圓筒內(nèi)壁涂抹適量潤(rùn)滑油，以減小限制壁與試樣摩擦的影響。為便于后續(xù)取出試樣進(jìn)行無側(cè)限抗壓強(qiáng)度試驗(yàn)，將圓筒沿軸線均勻分割成2塊并用抱箍固定。將240 g土樣注入到圓筒中進(jìn)行試驗(yàn)，萬能力學(xué)試驗(yàn)裝置以8 mm/s的恒定速率推動(dòng)底部直徑39 mm的壓盤向下加載。加載至軸向壓力為376.8 N（固結(jié)應(yīng)力約為300 kPa）時(shí)以相同速率回車，直至完全卸載。使用線性EEPA模型達(dá)到300 kPa固結(jié)應(yīng)力所用時(shí)間約為500 kPa所用時(shí)間的三分之一[18]，而非線性模型理論上將耗時(shí)更少，固結(jié)應(yīng)力設(shè)置為300 kPa有利于節(jié)約仿真時(shí)長(zhǎng)。

1.萬能力學(xué)試驗(yàn)裝置 2.壓盤 3.抱箍 4.有機(jī)玻璃圓筒

利用試驗(yàn)獲得土壤的軸向壓力-軸向應(yīng)變關(guān)系曲線（圖2）表征土壤彈性。軸向壓力-軸向應(yīng)變擬合函數(shù)如下：

式中F是軸向壓力，N；是試樣的軸向應(yīng)變，%；和是單軸密閉壓縮階段軸向壓力-軸向應(yīng)變曲線特征參數(shù)。

顯然式（1）是經(jīng)過原點(diǎn)的曲線，而萬能力學(xué)試驗(yàn)裝置以一定的入口力為起點(diǎn)采集數(shù)據(jù)，且加載初期，殘余應(yīng)變影響[27]導(dǎo)致式（1）與測(cè)量值的擬合程度不佳。因此，單軸密閉壓縮期間試樣的軸向壓力-軸向應(yīng)變曲線以固結(jié)應(yīng)力為1 kPa即軸向壓力約為1.256 N時(shí)的數(shù)據(jù)點(diǎn)為起點(diǎn)，將擬合曲線向左平移，使擬合曲線通過（0，1.256），為此，將式（1）調(diào)整為：

由于實(shí)際軸向壓力-軸向應(yīng)變曲線在45%的應(yīng)變以后與仿真曲線有較大差別，殘余應(yīng)變約為軸向應(yīng)變的0%～3%，所以式（2）僅擬合實(shí)測(cè)試驗(yàn)中加載階段3.0%～45%的軸向應(yīng)變。

利用試驗(yàn)后固結(jié)試樣發(fā)生的軸向應(yīng)變ε表征土壤塑性。ε定義為

式中0是單軸密閉壓縮試驗(yàn)中試樣的初始高度，mm；1是到達(dá)預(yù)設(shè)固結(jié)應(yīng)力并完全卸載后試樣的高度，mm。

測(cè)量塑性指標(biāo)時(shí)，將單軸密閉壓縮期間軸向壓力為1.2 N（固結(jié)應(yīng)力約為1 kPa）時(shí)的試樣高度記為0，完全卸載后的試樣高度記為1，依據(jù)式（3）計(jì)算ε。測(cè)量彈性指標(biāo)時(shí)，依據(jù)式（2）擬合3.0%～45%軸向應(yīng)變的軸向壓力-軸向應(yīng)變數(shù)據(jù)，并根據(jù)式（2）求解曲線特征參數(shù)和。

試驗(yàn)重復(fù)3次，結(jié)果為：為0.038、0.039、0.046；為3.667、3.782、3.702；ε為36.19%、37.22%、34.17%。取平均值有：=0.041，=3.717，ε=35.86%。

圖2 單軸密閉壓縮壓力曲線

1.2.2 土壤黏性指標(biāo)

利用無側(cè)限抗壓強(qiáng)度試驗(yàn)測(cè)量土壤無側(cè)限抗壓強(qiáng)度σ，試驗(yàn)如圖3所示。

1.萬能力學(xué)試驗(yàn)裝置 2.壓盤 3.失效試樣

σ定義為

試驗(yàn)時(shí)，利用切土刀將單軸密閉壓縮試驗(yàn)結(jié)束后的試樣高度切削至81 mm。修整試樣端面至平整，最后試樣高度約為80 mm，進(jìn)行無側(cè)限抗壓強(qiáng)度試驗(yàn)。萬能力學(xué)試驗(yàn)機(jī)以3.2 mm/s的恒定速率推動(dòng)底部直徑為100 mm的壓盤向下加載，直至?xí)r間-壓力曲線出現(xiàn)峰值。設(shè)置軸向壓力衰減至95%的峰值壓力時(shí)停止并回車。記錄峰值壓力和回車后的試樣端面校正半徑。試驗(yàn)重復(fù)3次，依據(jù)式（4）計(jì)算得到σ分別為30.32，35.46和42.61 kPa，結(jié)果取平均值為36.13 kPa。

1.3 仿真參數(shù)設(shè)置

1.3.1 單軸密閉壓縮仿真試驗(yàn)

單軸密閉壓縮仿真試驗(yàn)中，選擇顆粒-顆粒接觸模型為EEPA模型，使用高300 mm（以保證以圓筒為顆粒工廠生成足夠多的顆粒。實(shí)際試驗(yàn)中圓筒的高度為150 mm）、底部直徑40 mm的圓柱模擬有機(jī)玻璃圓筒。仿真參數(shù)設(shè)置如表1所示。試驗(yàn)中土壤顆粒與限制壁之間黏性相比于土壤顆粒之間黏性較小，故不考慮限制壁與土壤顆粒黏性的影響，顆粒-接觸部件模型選擇Hertze-Mindlin (no slip) 模型。

表1 單軸密閉壓縮仿真參數(shù)

Cleary[28]研究表明，無論顆粒間摩擦的角度如何，球形顆粒都不能代表“真實(shí)”固體，而非球形顆粒造成的顆粒群排列方式不同，可能導(dǎo)致試驗(yàn)數(shù)據(jù)的散射性較大。因此，單軸密閉壓縮期間以圓筒為顆粒工廠，生成25 960個(gè)半徑為1 mm的球形顆粒，以匹配與實(shí)際試驗(yàn)相同的初始填充質(zhì)量。

仿真中對(duì)同一初始填充試樣分別使用8（實(shí)際加載速率）、32.4和50 mm/s的加載速率（萬能力學(xué)試驗(yàn)機(jī)能達(dá)到的加載速率范圍為8～50 mm/s），發(fā)現(xiàn)3種加載速率下的試樣應(yīng)力-應(yīng)變行為具有良好的一致性，如圖4所示。此外，在50 mm/s的加載速率條件下，保持相同的仿真時(shí)間步長(zhǎng)，較大k和組合下的單軸密閉壓縮仿真試驗(yàn)容易發(fā)生穿模和崩潰。因此，為縮短仿真時(shí)間，單軸密閉壓縮仿真試驗(yàn)的加載速率均為32.4 mm/s。

圖4 不同加載速率下試樣的軸向壓力-軸向應(yīng)變曲線

1.3.2 無側(cè)限抗壓強(qiáng)度仿真試驗(yàn)

單軸密閉壓縮仿真試驗(yàn)結(jié)束后，去除限制壁，試樣松弛0.2 s，以消除由于移除限制壁而增加的動(dòng)能[18]。側(cè)限移除將導(dǎo)致試樣產(chǎn)生一定的徑向膨脹加速度，為平衡徑向加速度，防止無側(cè)限抗壓強(qiáng)度試驗(yàn)期間試樣滑移出計(jì)算域，在去除限制壁之前，適當(dāng)增加顆粒-幾何體靜摩擦系數(shù)至0.3。調(diào)整計(jì)算域高度以保留80 mm高度的試樣，試樣繼續(xù)松弛0.2 s以消除由于部分顆粒移除而導(dǎo)致的動(dòng)能變化。壓盤以3.2 mm/s的恒定速率向下移動(dòng)直至壓力衰減至95%峰值壓力，記錄峰值壓力和峰值壓力時(shí)的校正半徑，并根據(jù)式（4）計(jì)算σ。所有仿真試驗(yàn)的時(shí)間步長(zhǎng)均為1×10-5s，Rayleigh時(shí)間步長(zhǎng)約為6%。

2 參數(shù)標(biāo)定試驗(yàn)與結(jié)果分析

根據(jù)EEPA模型可同時(shí)表征塑性和黏性的特點(diǎn)，利用塑性和黏性指標(biāo)為響應(yīng)值，依據(jù)仿真設(shè)置，通過Plackett-Burman試驗(yàn)篩選出與黏塑性顯著相關(guān)的模型參數(shù)以減少參數(shù)標(biāo)定試驗(yàn)因素。通過二次正交旋轉(zhuǎn)組合試驗(yàn)，得到黏塑性指標(biāo)與顯著因素的回歸模型，并求解與實(shí)測(cè)黏塑性指標(biāo)對(duì)應(yīng)的黏塑性參數(shù)。在標(biāo)定黏塑性參數(shù)的基礎(chǔ)上，通過二次正交旋轉(zhuǎn)組合試驗(yàn)，得到彈性指標(biāo)與待標(biāo)定參數(shù)的回歸模型，并求解與實(shí)測(cè)黏塑性指標(biāo)和彈性指標(biāo)對(duì)應(yīng)的模型參數(shù)。

2.1 Plackett-Burman 試驗(yàn)

為從待標(biāo)定參數(shù)中獲取對(duì)ε和σ影響顯著的參數(shù)，依據(jù)參數(shù)范圍，以1個(gè)中心點(diǎn)估計(jì)試驗(yàn)隨機(jī)誤差和7個(gè)待標(biāo)定參數(shù)為影響因素，以ε和σ為指標(biāo)，Plackett-Burman仿真試驗(yàn)方案及結(jié)果如表2所示。

表2 Plackett-Burman試驗(yàn)方案與結(jié)果

對(duì)Plackett-Burman試驗(yàn)得到的ε和σ進(jìn)行方差分析，以評(píng)價(jià)各參數(shù)對(duì)黏塑性指標(biāo)的顯著性。方差分析結(jié)果如表3所示。以置信度0.05與各因素值比較可知，λ對(duì)σ影響顯著，Δ對(duì)σ影響極顯著，λ對(duì)ε影響極顯著，所以Δ和λ是影響土壤黏塑性的顯著因素。

從表3可以看出，靜摩擦系數(shù)μ的對(duì)σ和ε的值分別為0.695 7（>0.05）和0.804 7（>0.05），且當(dāng)顆粒形狀為球形時(shí)，靜摩擦系數(shù)對(duì)單軸密閉壓縮期間試樣的應(yīng)力-應(yīng)變行為無顯著影響[26]，說明μ在其試驗(yàn)范圍（0.2～1）內(nèi)對(duì)軸向壓力-軸向應(yīng)變曲線特征參數(shù)（，）、ε、σ均無顯著影響。因此后續(xù)試驗(yàn)將μ取中間值0.6。

2.2 黏塑性參數(shù)標(biāo)定結(jié)果與分析

為標(biāo)定與實(shí)測(cè)黏塑性指標(biāo)（σ和ε）對(duì)應(yīng)的Δ和λ，基于Plackett-Burman試驗(yàn)因素范圍，依據(jù)二次正交旋轉(zhuǎn)組合試驗(yàn)原理，以3個(gè)中心點(diǎn)、2個(gè)因素（Δ和λ）、2個(gè)指標(biāo)（σ和ε）進(jìn)行二次正交旋轉(zhuǎn)組合試驗(yàn)，試驗(yàn)因素水平如表4所示，試驗(yàn)方案和結(jié)果如表5所示。仿真試驗(yàn)中對(duì)黏塑性影響不顯著的其他參數(shù)均取Plackett-Burman試驗(yàn)因素范圍中間值，具體為：=0.5、μ=0.6、μ=0.3、=2.25，k=0.75。

表3 εn和σu的方差分析

注：<0.01表示極顯著，0.01≤<0.05表示顯著，≥0.05表示不顯著。下同。

Note:<0.01 means extremely significant, 0.01≤<0.05 means significant,≥0.05 means insignificant. The same as below.

表4 Δγ和λp二次正交旋轉(zhuǎn)組合試驗(yàn)因素水平編碼表

根據(jù)表5，利用Design Expert 8.0.5得到以下2個(gè)回歸方程：

表5 Δγ和λp的二次正交旋轉(zhuǎn)組合試驗(yàn)方案及結(jié)果

為驗(yàn)證所得ε和σ與Δ和λ關(guān)系模型的可靠性，進(jìn)行3組驗(yàn)證試驗(yàn)。結(jié)果表明，ε的預(yù)測(cè)值與實(shí)測(cè)值的平均誤差約為0.77%，標(biāo)準(zhǔn)差約為1.27%；σ的預(yù)測(cè)值與實(shí)測(cè)值的平均誤差約為1.20%，標(biāo)準(zhǔn)差約為10.59%。所得模型可準(zhǔn)確且穩(wěn)定預(yù)測(cè)試樣軸向應(yīng)變ε，可準(zhǔn)確預(yù)測(cè)σ但穩(wěn)定性略差。

將式（5）和式（6）聯(lián)立求解，得到與σ和ε對(duì)應(yīng)的Δ和λ分別為15.6 J/m2和0.36。

2.3 彈性參數(shù)標(biāo)定結(jié)果與分析

為減少與實(shí)測(cè)彈性指標(biāo)（，）對(duì)應(yīng)的彈性參數(shù)、μ、、k的解的數(shù)量，同時(shí)將黏彈塑性指標(biāo)作為評(píng)價(jià)指標(biāo)，以3個(gè)中心點(diǎn)、4個(gè)因素（、μ、和k）設(shè)計(jì)二次正交旋轉(zhuǎn)組合試驗(yàn)，星號(hào)臂長(zhǎng)1.546 71，因素水平如表6所示。將加載階段試樣的軸向壓力-軸向應(yīng)變數(shù)據(jù)與式（2）擬合，試驗(yàn)結(jié)果如表7所示。

表6 彈性參數(shù)的二次正交旋轉(zhuǎn)組合試驗(yàn)因素水平編碼表

根據(jù)表7，利用Design Expert 8.0.5軟件，保留不顯著項(xiàng)的條件下得到以下4個(gè)回歸方程：

表7 e,μr, X和ktm的二次正交旋轉(zhuǎn)組合試驗(yàn)方案及結(jié)果

對(duì)得到的ε、σ、和進(jìn)行方差分析，結(jié)果如表8所示。根據(jù)表8，ε、σ、、的模型值分別為0.001 1（<0.01）、0.005 9（<0.01）、<0.000 1（<0.01）和<0.000 1（<0.01），說明標(biāo)定彈性參數(shù)的各指標(biāo)擬合模型均極顯著，擬合模型能夠準(zhǔn)確描述因變量與自變量之間的關(guān)系。切向剛度因子k對(duì)ε、σ、、的值均為<0.000 1（<0.01），且k二次項(xiàng)對(duì)ε、σ影響極顯著。

以實(shí)測(cè)的ε、σ、和為目標(biāo)對(duì)式（7）～（10）方程組求解，獲得27組解。這些解幾乎無差異，以Desirability值為參考，選擇一組Desirability值最接近1（可靠度最高）的參數(shù)=0.37、μ=0.26、=4.24和k=0.52。

表8 黏彈塑性指標(biāo)方差分析

續(xù)表

3 標(biāo)定結(jié)果驗(yàn)證

根據(jù)黏塑性參數(shù)和彈性參數(shù)標(biāo)定結(jié)果，利用EDEM_v2018軟件建立離散元仿真模型，進(jìn)行單軸密閉壓縮仿真試驗(yàn)和無側(cè)限抗壓強(qiáng)度仿真試驗(yàn)，以驗(yàn)證標(biāo)定結(jié)果的可信度。試驗(yàn)重復(fù)3次，結(jié)果如表9所示。

表9 單軸密閉壓縮和無側(cè)限抗壓強(qiáng)度驗(yàn)證試驗(yàn)結(jié)果

將表9中第一組單軸密閉壓縮仿真試驗(yàn)的軸向壓力-軸向應(yīng)變曲線與實(shí)測(cè)的土壤軸向壓力-軸向應(yīng)變曲線對(duì)比，如圖5所示。利用Orign Pro8.5軟件，將仿真值與實(shí)測(cè)值對(duì)軸向應(yīng)變積分求差，取其相對(duì)軸向應(yīng)變的平均值，得到仿真軸向壓力與實(shí)測(cè)值平均誤差約為1.95 N。積分差值相對(duì)實(shí)測(cè)值積分誤差約為3.3%。故可以認(rèn)為，在3%～45%的軸向應(yīng)變內(nèi)仿真試樣的應(yīng)力-應(yīng)變行為與實(shí)際情況基本一致。

圖5 標(biāo)定區(qū)間的應(yīng)力-應(yīng)變曲線對(duì)比

4 結(jié) 論

本研究通過Plackett-Burman試驗(yàn)和Central Composite試驗(yàn)標(biāo)定了EEPA模型參數(shù)。結(jié)合標(biāo)定試驗(yàn)與驗(yàn)證試驗(yàn)結(jié)果與分析，主要結(jié)論如下：

1）Plackett-Burman試驗(yàn)結(jié)果表明，模型參數(shù)中的塑性變形比對(duì)軸向應(yīng)變影響極顯著，對(duì)無側(cè)限抗壓強(qiáng)度影響顯著；表面能對(duì)無側(cè)限抗壓強(qiáng)度影響極顯著。

2）參數(shù)標(biāo)定試驗(yàn)結(jié)果表明，軸向應(yīng)變、無側(cè)限抗壓強(qiáng)度關(guān)于黏塑性參數(shù)的擬合模型均顯著。軸向應(yīng)變、無側(cè)限抗壓強(qiáng)度、曲線特征參數(shù)關(guān)于彈性參數(shù)的擬合模型均極顯著。對(duì)回歸模型求解，得到與實(shí)測(cè)值對(duì)應(yīng)的表面能、塑性變形比、恢復(fù)系數(shù)、滾動(dòng)摩擦系數(shù)、黏結(jié)分支指數(shù)、切向剛度因子分別為15.6 J/m2、0.36、0.37、0.26、4.24、0.52。

3）標(biāo)定結(jié)果驗(yàn)證試驗(yàn)結(jié)果表明，軸向應(yīng)變、曲線特征參數(shù)的仿真值與實(shí)測(cè)值的相對(duì)誤差分別為-3.40%和0.54%，軸向壓力仿真值與實(shí)測(cè)值誤差為1.95 N，與實(shí)測(cè)值積分差值相對(duì)誤差為3.3%。無側(cè)限抗壓強(qiáng)度、曲線特征參數(shù)與實(shí)測(cè)值的差異略大。

對(duì)土壤黏塑性指標(biāo)影響顯著的因素不只是表面能和塑性變形比，切向剛度因子對(duì)其影響也較大，因而造成標(biāo)定結(jié)果的驗(yàn)證試驗(yàn)中，無側(cè)限抗壓強(qiáng)度的仿真值與實(shí)測(cè)值存在較大差異，后續(xù)研究應(yīng)在Plackett-Burman試驗(yàn)前通過單因素試驗(yàn)確定切向剛度因子的取值范圍。

[1] Coetzee C J, Lombard S G. The destemming of grapes: Experiments and discrete element modelling[J]. Biosystems Engineering, 2013, 114(3): 232-248.

[2] Yu Y, Fu H, Yu J. DEM-based simulation of the corn threshing process[J]. Advanced Powder Technology, 2015, 26(5): 1400-1409.

[3] 王立軍，馮鑫，鄭招輝，等. 玉米清選組合孔篩體設(shè)計(jì)與試驗(yàn)[J]. 農(nóng)業(yè)機(jī)械學(xué)報(bào)，2019，50(5)：104-113.

Wang Lijun, Feng Xin, Zheng Zhaohui, et al. Design and test of combined sieve of maize screening[J]. Transactions of the Chinese Society for Agricultural Machinery, 2019, 50(5): 104-113. (in Chinese with English abstract)

[4] 王立軍，武振超，馮鑫，等. 玉米收獲機(jī)清選曲面篩設(shè)計(jì)與試驗(yàn)[J]. 農(nóng)業(yè)機(jī)械學(xué)報(bào)，2019，50(2)：90-101.

Wang Lijun, Wu Zhenchao, Feng Xin, et al. Design and experiment of curved screen for maize grain harvester[J]. Transactions of the Chinese Society for Agricultural Machinery, 2019, 50(2): 90-101. (in Chinese with English abstract)

[5] 王立軍，張傳根，丁振軍. 玉米收獲機(jī)清選篩體結(jié)構(gòu)優(yōu)化[J]. 農(nóng)業(yè)機(jī)械學(xué)報(bào)，2016，47(9)：108-114.

Wang Lijun, Zhang Chuangen, Ding Zhenjun. Structure optimization of cleaning screen for maize harvester[J]. Transactions of the Chinese Society for Agricultural Machinery, 2016, 47(9): 108-114. (in Chinese with English abstract)

[6] 趙淑紅，劉宏俊，譚賀文，等. 仿旗魚頭部曲線型開溝器設(shè)計(jì)與性能試驗(yàn)[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2017，33(5)：32-39.

Zhao Shuhong, Liu Hongjun, Tan Hewen, et al. Design and performance experiment of opener based on bionic sailfish head curve[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2017, 33(5): 32-39. (in Chinese with English abstract)

[7] 鄭侃，何進(jìn)，李洪文，等. 基于離散元深松土壤模型的折線破土刃深松鏟研究[J]. 農(nóng)業(yè)機(jī)械學(xué)報(bào)，2016，47(9)：62-72.

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]. Transactions of the Chinese Society for Agricultural Machinery, 2016, 47(9): 62-72. (in Chinese with English abstract)

[8] 張銳，李建橋，周長(zhǎng)海，等. 推土板表面形態(tài)對(duì)土壤動(dòng)態(tài)行為影響的離散元模擬[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2007，23(9)：13-19.

Zhang Rui, Li Jianqiao, Zhou Changhai, et al. Simulation of dynamic behavior of soil ahead of the bulldozing plates with different surface configurations by discrete element method[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2007, 23(9): 13-19. (in Chinese with English abstract)

[9] 于建群，錢立彬，于文靜，等. 開溝器工作阻力的離散元法仿真分析[J]. 農(nóng)業(yè)機(jī)械學(xué)報(bào)，2009，40(6)：53-57.

Yu Jianqun, Qian Libin, Yu Wenjing, et al. DEM Analysis of the resistances applied on furrow openers[J]. Transactions of the Chinese Society for Agricultural Machinery, 2009, 40(6): 53-57. (in Chinese with English abstract)

[10] 馬躍進(jìn)，王安，趙建國(guó)，等. 基于離散元法的凸圓刃式深松鏟減阻效果仿真分析與試驗(yàn)[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2019，35(3)：16-23.

Ma Yuejin, Wang An, Zhao Jianguo, et al. Simulation analysis and experiment of drag reduction effect of convex blade subsoiler based on discrete element method[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2019, 35(3): 16-23. (in Chinese with English abstract)

[11] Ucgul M, Fielke J M, Saunders C. 3D DEM tillage simulation: Validation of a hysteretic spring (plastic) contact model for a sweep tool operating in a cohesionless soil[J]. Soil and Tillage Research, 2014, 144: 220-227.

[12] Bravo E L, Tijskens E, Suarez M H, et al. Prediction model for non-inversion soil tillage implemented on discrete element method[J]. Comput. Electron. Agric, 2014, 106: 120-127.

[13] Chen Y, Munkholm L J, Nyord T. A discrete element model for soil-sweep interaction in three different soils[J]. Soil and Tillage Research, 2013, 126: 34-41.

[14] Li B, Chen Y, Chen J. Modeling of soil–claw interaction using the discrete element method (DEM)[J]. Soil and Tillage Research, 2016, 158: 177-185.

[15] Katterfeld A, Wensrich C. Understanding granular media: from fundamentals and simulations to industrial application[J]. Granular Matter, 2017, 19(4).

[16] Roessler T, Richter C, Katterfeld A, et al. Development of a standard calibration procedure for the DEM parameters of cohesionless bulk materials-part I: Solving the problem of ambiguous parameter combinations[J]. Powder Technology, 2019, 343: 803-812.

[17] Karkala S, Davis N, Wassgren C, et al. Calibration of discrete- element-method parameters for cohesive materials using dynamic-yield-strength and shear-cell experiments[J]. Processes, 2019, 7(5): 278.

[18] Janda A, Ooi J Y. DEM modeling of cone penetration and unconfined compression in cohesive solids[J]. Powder Technology, 2016, 293: 60-68.

[19] 王憲良，胡紅，王慶杰，等. 基于離散元的土壤模型參數(shù)標(biāo)定方法[J]. 農(nóng)業(yè)機(jī)械學(xué)報(bào)，2017，48(12)：78-85.

Wang Xianliang, Hu Hong, Wang Qingjie, et al. Calibration method of soil contact characteristic parameters based on DEM theory[J]. Transactions of the Chinese Society for Agricultural Machinery, 2017, 48(12): 78-85. (in Chinese with English abstract)

[20] Xia R, Li B, Wang X, et al. Measurement and calibration of the discrete element parameters of wet bulk coal[J]. Measurement, 2019, 142: 94-95.

[21] Ucgul M, Fielke J M, Saunders C. Three-dimensional discrete element modelling of tillage: Determination of a suitable contact model and parameters for a cohesionless soil[J]. Biosystems Engineering, 2014, 121:105-117.

[22] Barrios G K P, Carvalho R M, Kwade A, et al. Contact parameter estimation for DEM simulation of iron ore pellet handling[J]. Powder Technology, 2013, 248:84-93.

[23] Coetzee C. Calibration of the discrete element method and the effect of particle shape[J]. Powder Technol, 2016, 297: 50-70.

[24] Ucgul M, Fielke J M, Saunders C. Three-dimensional discrete element modelling (DEM) of tillage: Accounting for soil cohesion and adhesion[J]. Biosystems Engineering, 2015, 129: 298-306.

[25] 李俊偉，佟金，胡斌，等. 不同含水率黏重黑土與觸土部件互作的離散元仿真參數(shù)標(biāo)定[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2019，35(6)：130-140.

Li Junwei, Tong Jin, Hu Bin, et al. Calibration of parameters of interaction between clayey black soil with different moisture content and soil-engaging component in northeast China[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2019, 35(6): 130-140.

[26] Thakur S C, Morrissey J P, Sun J, et al. Micromechanical analysis of cohesive granular materials using the discrete element method with an adhesive elasto-plastic contact model[J]. Granular Matter, 2014, 16(3): 383-400.

[27] Tomas J. Adhesion of ultrafine particles: A micromechanical approach[J]. Chemical Engineering Science, 2007, 62(7): 1997-2010.

[28] Cleary P W. DEM prediction of industrial and geophysical particle flows[J]. Particuology, 2010, 8(2): 106-118.

Calibration of discrete element parameters of soils based on unconfined compressive strength test

Xie Fangping1,2, Wu Zhengyang1, Wang Xiushan1,2, Liu Dawei1,2, Wu Bei1,2, Zhang Zhengzhong1

(1.,410128,;2.410128,)

In order to calibrate the discrete element parameters of viscoplastic soil based on the Edinburgh Elasto-Plastic Adhesion (EEPA) model, two constants to describe the stress-strain behavior of the sample soils were defined in this study. The axial strain of the soil under a load of 300 kPa was used to characterize the plasticity of the soil in the uniaxial closed compression test. First, two EEPA model parameters based on the physical properties of the soil and other scholars’ research were determined, i.e. the constant pull-off force and the load branch index. Then, another two factors, surface energy and plastic deformation ratio, that had significant effects on axial strain and unconfined compressive strength based on the Plackett-Burman test results were described. Next, the central composite test based on the response surface method was designed, and the two factors that matched the actual measured axial strain and unconfined compressive strength based on the test results were determined.According to the test results. the four discrete element parameters corresponding to the two measured constants, unconfined compressive strength and axial strain were solved. Finally, the soil discrete element parameters were calibrated based on the EEPA model,and those were that plastic deformation ratio of 0.36, surface energy of 15.6 J/m2, static friction coefficient of 0.6, rolling friction coefficient of 0.26, recovery coefficient of 0.37, adhesion branch index of 4.24 and tangential stiffness factor of 0.52.Verification test results showed that the EEPA model parameters calibrated based on the response surface method could simulate the plastic deformation of the sample soil under a load of 300 kPa and the stress-strain behavior within 3%-45% of the axial strain. In addition, the results of Quadratic orthogonal rotation combined test were analyzedand it showed that the tangential stiffness factor was one of the key parameters affecting axial strain and unconfined compressive strength, and the random error of the unconfined compressive strength simulation test was also one of the reasons. Moreover, it was found that the limit of the value range of the tangential stiffness factor was the cause of the huge difference between the simulation and the measured values of the unconfined compressive strength.

soils;stress; strain;discrete element method; calibration; plastic deformation; unconfined compressive strength

謝方平，吳正陽，王修善，等. 基于無側(cè)限抗壓強(qiáng)度試驗(yàn)的土壤離散元參數(shù)標(biāo)定[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2020，36(13)：39-47.doi：10.11975/j.issn.1002-6819.2020.13.005 http://www.tcsae.org

Xie Fangping, Wu Zhengyang, Wang Xiushan, et al. Calibration of discrete element parameters of soils based on unconfined compressive strength test[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2020, 36(13): 39-47. (in Chinese with English abstract) doi：10.11975/j.issn.1002-6819.2020.13.005 http://www.tcsae.org

2019-12-18

2020-06-07

湖南省農(nóng)機(jī)裝備創(chuàng)新發(fā)展項(xiàng)目（湘財(cái)農(nóng)指（2018）175號(hào)）

謝方平，教授，主要從事農(nóng)業(yè)機(jī)械性能試驗(yàn)與創(chuàng)新設(shè)計(jì)。Email：hunanxie2002@163.com

10.11975/j.issn.1002-6819.2020.13.005

O347.7

A

1002-6819(2020)-13-0039-09