半夏塊莖物理特性研究及離散元仿真參數標定

楊貴川，張富貴，鄭樂，王震，孔曼曼，章鑫鵬

（貴州大學機械工程學院，貴陽 550025）

半夏（Pinellia ternata）具有燥濕化痰、降逆止嘔、消痞散結的功效，全國第三次中藥資源普查顯示，在588 個抽樣處方中，半夏在處方中出現的頻率位于第22 位；國家210 個法定漢方中含半夏的有46 個，其應用量約占第21 位，臨床應用廣泛，具有十分重要的藥用價值[1-3]。半夏在貴州廣泛種植，種植面積、產量均位居全國前列，還擁有“大方圓珠半夏”和“赫章半夏”2 個國家地理標志保護產品。半夏塊莖具有極高的經濟價值，目前成品半夏塊莖市場價格在80～100 元·kg-1，極大地助力了當地脫貧攻堅工作[4-5]。半夏的種植方式主要為條播和壟上撒播，為了保證半夏品質，貴州地區采用壟上撒播種植模式（在半夏播種初期施用有機肥和復合肥，后期不再對其進行除草、殺蟲作業，通過高密度播種的方式抑制雜草生長，保證產量）模擬野生生長環境，撒播作業通過人工進行，費時、費力且種植成本居高不下。因此為了省工降本，進一步擴大生產規模，需盡快實現機械化作業。

部分半夏產區采用鏈勺式排種器、外槽輪式排種器等設備進行條播作業，存在易傷種、作業幅寬過小、播種均勻性差、播種密度無法達到作業要求等問題。離散元仿真廣泛應用于分析顆粒物料與設備交互作用下的運動、受力等，可以縮短設備的研發周期，降低研發成本。為了設計適應貴州地區半夏種植特點的半夏撒播排種器，有必要對半夏塊莖參數進行仿真標定。文恩楊等[6-8]設計了大蒜播種機，并通過離散元仿真優化了關鍵部分的結構參數；劉文政等[9-11]對微型馬鈴薯進行了離散元仿真標定，并確定了基于振動排序的馬鈴薯微型種薯播種機的最佳振動頻率范圍；溫翔宇等[12]通過離散元仿真試驗優化了氣力變量配比施肥裝置；劉彩鈴等[13-14]基于EDEM 離散元仿真對離心甩盤撒肥器的性能進行了分析與試驗，確定了最佳甩盤轉速、喂入量、喂入角和喂入位置角；吳孟宸等[15]對花生顆粒進行了離散元仿真標定，并驗證了花生排種器的排種性能；賴慶輝等[16-19]基于三七種苗的離散元仿真參數設計了滾筒板齒式三七種苗分離裝置、超窄行氣吸式三七精密排種器及凸包異形孔窩眼輪式人參精密排種器，并進行了優化設計。

本研究采用實際接觸試驗確定仿真參數試驗范圍，結合最陡爬坡試驗和Box-Behnken 試驗對半夏的離散元仿真參數進行標定；通過壓縮試驗對半夏的本征物理特性進行測定分析，建立了半夏離散元仿真模型，以期為半夏撒播排種器提供離散元仿真參數。

1 材料與方法

1.1 試驗材料

本研究以貴州省赫章縣種植基地的半夏塊莖為試驗材料，塊莖近似橢球形，隨機選取100 顆，利用游標卡尺對其三軸尺寸進行測量，根據下式求出半夏的等效直徑。

式中，d為等效直徑，L為長度，W為寬度，T為厚度。

采用軟件Origin 軟件擬合得出半夏等效直徑的分布規律，半夏等效直徑為正態分布（圖1），峰值為16.8 mm。根據所測結果，對比實際半夏塊莖在EDEM 中組合創建半夏的仿真模型（圖2），模型的物理特性參數根據測定參數設定。

圖1 半夏等效直徑分布Fig.1 Equivalent diameter distribution of P.ternata stems

圖2 半夏塊莖及仿真模型Fig.2 Tuber and simulated particle model

1.2 半夏基本物理特性參數測定

用電子秤稱取播種期半夏質量（上1 年度10月份收獲的半夏，次年3月播種），將半夏放入DHG-9741A 型電熱恒溫干燥箱（上海精宏實驗設備有限公司），105 ℃烘烤至恒重，計算半夏含水率。在500 mL 的量筒中注入300 mL 清水，將鮮半夏放入量筒中使其自然沉入水底，讀取量筒內水的體積變化量，計算半夏密度。

根據泊松比的定義，對半夏進行壓縮試驗。由于半夏形狀不規則，因此將半夏用內徑10 mm、高40 mm 的環刀和美工刀切割成直徑10 mm、高10 mm 的規整試樣。將試樣用QJ211S 型農作物專用試驗機（上海傾技儀器儀表科技有限公司）對試樣進行壓縮試驗，并測量試樣受壓變形后的橫向及縱向應變量[15,20-21]。分別以1、3、5 mm·min-1的速率壓縮試樣，每個速率重復9 次，試樣產生塑性變形前停止壓縮，用電子游標卡尺測量其橫向與縱向變形尺寸，計算泊松比。

式中，μ1為泊松比；εl為橫向應變；ε為縱向應變；L0為試樣原始長度，mm；L1為壓縮變形后試樣長度，mm；d0為試樣初始直徑，mm；d1為試樣變形后直徑，mm。彈性模量E計算公式如下。

式中，E為彈性模量，MPa；ΔF為壓力變化量，N；A0為試樣橫截面積，mm2。

1.3 碰撞恢復系數標定

分別用尼龍繩綁住球形半夏胚芽，尼龍繩另一端固定在同一點，固定點背景為A3 坐標紙；為便于觀察碰撞過程，應盡量選擇大的擺動半徑，取擺動半徑為360 mm（由于半夏中心點難以確定，因此試驗過程中，測量拉繩與半夏的捆綁點高度變化，代替中心點高度變化）。如圖3 所示，將半夏a 沿圓周方向拉起一定高度并釋放，使其自由擺動最終與半夏b 發生碰撞，a、b 沿周向運動，最終動能全部轉換為重力勢能（試驗過程由高速攝像機記錄），通過背景坐標紙測量a、b 到零點（a、b 自然懸垂點）的垂直距離，根據式（4）計算碰撞恢復系數。

圖3 半夏-半夏碰撞仿真試驗Fig.3 Simulation test of pinellia stalk and pinellia stalk collision

式中，e1為碰撞恢復系數；vb為半夏b 的瞬時分離速度，mm·s-1；va為半夏a 的瞬時分離速度，mm·s-1；v0為半夏a的瞬時接近速度，mm·s-1；H0為半夏a初始釋放高度，mm；Ha為半夏a極限爬升高度，mm；Hb為半夏b極限爬升高度，mm。

試驗過程中發現，將半夏從初始位置拉升抬高100 mm處時將超出測量范圍，為確保高速相機能夠清晰地記錄碰撞全過程及坐標紙標尺，取H0分別為30、50 mm，重復30 次碰撞試驗，測得Ha分別為（3.2±0.5）和（3.75±1）mm，Hb分別為（16.56±3）和（25.70±6）mm，碰撞恢復系數e＇1分別為（0.469 0±0.05）和（0.474 0±0.05）mm，實際試驗測得的e＇1變異系數為0.75%，表明H0對碰撞恢復系數并無影響，因此取50 mm 作為仿真試驗初始釋放高度（H0）。

在EDEM 仿真試驗（圖3）中，采用中心半徑為360 mm 的半圓形空心圓筒以控制半夏的擺動軌跡（代替拉繩）。H0為初始釋放高度，ya為顆粒a 的極限爬升高度，yb為顆粒b 的極限爬升高度，圓筒內徑為15 mm，分別在圓筒的軌道鉛錘方向管道中心點和距離該點50 mm 處各生成1 粒球形半夏顆粒模型a、b,為避免摩擦力對仿真結果產生影響，將顆粒模型與管道模型間的摩擦系數設置為0。根據仿真結果分別建立回歸模型，結合實際試驗結果標定出半夏-不銹鋼間的碰撞恢復系數e1。

同理，分別在40、80、120 mm 處釋放半夏，測得半夏-鋼板的反彈平均高度H1分別為(16.47±3.52)(29.33±2.43)(45.67±2.21) mm，根據式（5）計算出半夏-不銹鋼實碰撞恢復系數e＇2分別為（0.63±0.05）（0.60±0.03）（0.61±0.03），變異系數為2.49%，變異系數較小，表明半夏的釋放高度對半夏-不銹鋼板間的碰撞恢復系數影響越小。

式中，H＇0為半夏跌落高度，mm；H1為半夏極限回彈高度，mm。

選擇120 mm 作為仿真試驗釋放高度開展標定，通過實際試驗表明半夏與不銹鋼板間的碰撞恢復系數在0.58～0.70 之間。在EDEM 仿真過程中，不銹鋼的密度8 000 kg·m-3，泊松比為0.29、彈性模量7 500 MPa[22]，為避免摩擦系數的影響，將摩擦系數設定為0，將半夏-不銹鋼板的碰撞恢復系數設定在0.58～0.70，步長為0.02。根據仿真結果分別建立回歸模型，結合實際試驗結果標定出半夏-不銹鋼間的碰撞恢復系數e2。

1.4 摩擦因素標定

1.4.1 半夏-不銹鋼靜摩擦因素標定 將4 個體積相當的半夏用透明膠帶粘貼在一起（防止半夏發生滾動），通過斜面滑動試驗測定半夏-不銹鋼的靜摩擦系數，利用數顯角度尺（精度0.05°，量程0°～360°）測量鋼板與水平面間的夾角，靜摩擦系數與抬升平面傾角間的關系如下。

式中，μ2為半夏-不銹鋼板靜摩擦系數；f為摩擦力，N；F為顆粒平行于斜面的分力，N；β為斜面與水平面的夾角，測量5次取均值，β為31.55°。

由于在EDEM 軟件仿真過程中，其他仿真參數對靜摩擦因素并無影響，因此引用文獻[11]對離散元仿真顆粒模型靜摩擦系數標定的回歸方程。

式中，y2為傾角，x2為塊莖與鋼板間的靜摩擦系數，取y2=31.55°代入式（7），可得半夏-不銹鋼靜摩擦因素μ2。

1.4.2 滾動摩擦因素及半夏間靜摩擦系數標定由于滾動摩擦系數、顆粒間的靜摩擦系數難以直接進行測定，因此采用堆積角試驗進行間接測量，顆粒堆積角是物料流動性、物料摩擦系數等本征物理特性的宏觀體現，物料的靜摩擦系數、滾動摩擦系數等物理特性參數對堆積角的角度大小影響顯著[23-26]，因此通過無底圓筒（內徑120 mm，高度360 mm）提升試驗測定半夏在不銹鋼板上的堆積角[27]。圓筒提升試驗獲取堆積角，重復5 次取平均值，測得堆積角為（37.55°±1°）。

用EDEM 軟件對上述試驗作重復仿真試驗，設計最陡爬坡試驗確定半夏-半夏摩擦系數、半夏-不銹鋼滾動摩擦系數取值范圍，初始值為0，以0.03為步長逐漸增加，進行爬坡試驗（仿真過程中其他參數采用文中已標定的數值）。用Box-Behnken 試驗確定顯著性影響因數，并根據試驗結果建立回歸模型，對回歸模型進行二次優化，最終結合實際試驗標定摩擦系數。

2 結果與分析

2.1 半夏本征物理特性參數分析

通過烘干試驗測得半夏含水率為62.23%，通過浸液試驗測得半夏的體積，并測量半夏的質量，最終求得半夏的密度為1 210 kg·m-3。經壓縮試驗得到不同壓縮速率下的泊松比及彈性模量，結果如表1 所示。通過對比分析得知，不同壓縮速率下的半夏泊松比μ1變異系數為1.85%、彈性模量E變異系數為9.01%，泊松比及彈性模量變異系數較小，因此取3 種壓縮速率下的平均值，半夏的實際泊松比μ1=0.373 1、彈性模量E=5.757 1 MPa。

表1 半夏壓縮試驗結果Table 1 Result of P.ternata tuber compression

2.2 半夏-半夏碰撞恢復系數仿真參數分析

為保證仿真試驗不受摩擦力的影響，將圓筒軌道-半夏模型、半夏模型間的摩擦系數、碰撞恢復系數設為0；用EDEM 后處理功能測量半夏a、b碰撞后的最大提升高度ya、yb，結果如表2 所示。以半夏碰撞恢復系數x為自變量，半夏a、半夏b碰撞后的最大爬升高度ya、yb為因變量。用Origin對仿真結果進行二次多項式擬合，得到式（8）。

表2 半夏-半夏碰撞仿真試驗結果Table 2 Simulation test results of P.ternata-P.ternata collision

擬合方程決定系數分別為=0.957 1、=0.992 2，均接近1，表明擬合方程準確可靠。將Ha=3.75 mm、Hb=25.70 mm 代入式(8)，聯立可得x=0.472 0。將0.472 0 作為半夏間恢復系數，在H0=50 mm 為初始條件下各做3 次重復仿真試驗，得半夏a 碰撞后上升的最大高度平均值ya=4.83 mm，與實際試驗平均值的相對誤差為22.36%；半夏b 碰撞后上升的最大高度平均值yb=26.27 mm，與實際試驗測量值相對誤差為2.17%；仿真碰撞恢復系數x=0.472 0 與真實試驗測量中得碰撞恢復系數e＇1=0.474 0，二者相對誤差僅為0.42%；表明仿真參數與實際試驗結果差異較小，因此取半夏間仿真碰撞恢復系數e1=0.472 0。

2.3 半夏-不銹鋼的碰撞恢復系數仿真參數分析

半夏-不銹鋼板的碰撞仿真試驗結果如表3所示。以恢復系數x1為試驗因數，回彈高度y1為評價指標，對結果進行曲線擬合，建立擬合方程。

表3 半夏-不銹鋼板碰撞仿真試驗結果Table 3 Simulation test result of P.ternata-stainless steel plate collision

決定系數R2=0.999 9，接近于1，表明擬合方程可信度高。將y1=45.67 帶入方程（9）得到x1=0.635 8。以0.635 8 為半夏-不銹鋼板碰撞恢復系數，分別在H0為40、80、120 mm 釋放高度下各做3 組重復仿真試驗，得塊莖平均仿真回彈高度H1分別為16.77、31.24、45.82 mm，與真實試驗條件下得到的反彈高度相對誤差分別為1.82%、6.51%、0.33%，x1與e2＇間的相對誤差為3.89%誤差較小，故取半夏-不銹鋼板間的仿真碰撞恢復系數為e2=0.635 8。

2.4 半夏-不銹鋼板靜摩擦系數仿真參數分析

將y2=31.55°帶入式（7）中，得x2=0.615 4，將0.615 4 設置為半夏-不銹鋼板間的靜摩擦系數，其他摩擦系數設置為0，進行3次重復斜面滑動仿真試驗，得到不銹鋼板傾角分別為31.58°、31.57°、31.57°，取平均值為31.573 0°，與實際試驗值31.550 0°的相對誤差為0.643%。表明標定后的仿真結果與實際試驗結果差異較小，因此取半夏-不銹鋼板將的靜摩擦系數μ2=0.615 4。

2.5 半夏-半夏靜摩擦系數對滾動摩擦素的影響

2.5.1 半夏-半夏靜摩擦系數對堆積角的影響爬坡試驗結果如表4所示，隨著A、B、C增加，堆積角逐漸增大，真實試驗與仿真試驗得到的半夏堆積角的相對誤差先減小后增加。在6 號水平時，該相對誤差最小，由此可知最優值區間在6 號水平附近，因此選取6 號水平為中心點，5 號、7 號水平為低、高水平進行后續的響應面設計。

表4 最陡爬坡試驗Table 4 Steepest climbing test

利用Design-Expert 軟件設計Box-Behnken 試驗，根據仿真試驗結果（表5）建立回歸模型（式10）。

表5 堆積角Box-Behnken試驗結果Table 5 Result of stacking angle Box-Behnken test

方差分析（表6）表明，B的二次項以及BC交互項對堆積角影響極為顯著；擬合方程的決定系數R2=0.995 6，校正決定系數R2adj=0.987 8，二者均接近1，表明擬合方程可靠度高；精密度為37.422 2，表明該模型具有良好的精確度。

表6 堆積角試驗方差分析Table 6 Analysis of variance in the stacking angle test

剔除表6 中對堆積角仿真試驗影響不顯著的項，對二階回歸模型進行優化，結果如表7 所示，得到優化后的回歸模型（式11）。

表7 二次優化后回歸模型Table 7 Regression model after second optimization

優化后的回歸模型變異系數為1.390 0%，決定系數R2=0.979 4，校正決定系數R2adj=0.967 9，二者都接近1，模型精密度增大到56.052 9，表明擬合方程可靠度高，因此可利用該回歸方程可以用來預測半夏的堆積角。

2.5.2 回歸模型尋優求解及仿真驗證 通過表7可知，交互項BC對半夏的堆積角影響顯著。當A=0.15時，應用Design-Expert軟件繪制BC交互項的響應曲面，如圖5 所示，當BC同時增大或減小時半夏堆積角變化顯著。通過Design-Expert軟件對優化后的回歸模型進行目標尋優求解，以實際試驗時測得的半夏平均堆積角37.55°為目標，由于A不顯著，取中間水平0.15，通過Design-Expert求解可得：B為0.554，C為0.157。為驗證標定參數的準確性，采用上文中標定的各參數數值進行無底圓筒提升仿真試驗，獲得半夏堆積角，如圖6所示，得到堆積角分別為37.80°、37.29°、37.43°，取平均值為37.51°，與實際堆積角的相對誤差為0.107%，表明仿真結果與實際試驗無顯著差異。

圖5 BC交互項響應曲面Fig.5 BC interaction term response surface

圖6 半夏堆積角Fig.6 Pile of stems of P.ternata

3 討論

本文以貴州赫章半夏為研究對象，建立了半夏的EDEM 仿真模型，并測定了半夏的含水率為62.23%；采用浸液試驗測定了半夏的密度為1 210 kg·m-3，半夏的泊松比μ=0.373 1、彈性模量E=5.757 1 MPa。當壓縮應力超過標定參數時，半夏將發生塑性形變，因此在設計半夏排種器時，應考慮半夏受力狀態，避免破壞半夏導致發芽率降低。

利用實際物理試驗結合EDEM 仿真試驗的方法測定半夏間的碰撞恢復系數、半夏-不銹鋼板間的碰撞恢復系數、半夏-不銹鋼板間的靜摩擦系數，對試驗結果進行回歸分析，并建立回歸方程，分析結果表明：半夏-半夏間碰撞恢復系數e1=0.472 0，半夏-不銹鋼板間碰撞恢復系數e2=0.635 8，半夏-不銹鋼板間的靜摩擦系數μ2=0.615 4。

由于A、B、C難以直接進行測量，采用實際圓筒提升試驗結合最陡爬坡仿真試驗，篩選出中心水平，并用Design-Expert軟件設計Box-Behnken試驗，根據試驗結果建立回歸方程，對回歸方程進行方差分析，并對回歸模型進行二次優化，最終獲得到A=0.150、B=0.554、C=0.157。

顆粒物料的結拱與物料的直徑、摩擦因、含水率等因數有關[28-29]，目前對物料結拱條件的理論性研究較少，離散元仿真可以用于研究半夏結拱形成條件，對半夏進行離散元仿真參數標定及物理特性研究，可以為半夏排種器結構尺寸確定提供依據，同時也為半夏排種器離散元仿真分析提供了參數依據。