基于JKR模型的辣椒籽離散元參數標定*

徐莊威，王士林，易中懿，潘健，呂曉蘭

(1. 江蘇省農業科學院農業設施與裝備研究所,南京市,210014; 2. 南京林業大學機械電子工程學院,南京市,210037; 3. 農業農村部長江中下游設施農業工程重點試驗室,南京市,210014)

0 引言

機械化播種是農業智能化發展的一種趨勢,其對種子的質量與形狀有一定的要求,因此小粒不規則種子的丸化加工便成為機械化播種的前提[1]。辣椒籽扁平、輕質,在造粒過程中,辣椒籽之間以及辣椒籽與釜壁之間的相互作用非常復雜,對于其涉及的一些接觸參數難以定量,從而導致數值模擬時模型不精確的問題。

針對以上問題,目前學者主要以離散元法(DEM)研究散體顆粒的力學性質[2-3],以及對顆粒參數的標定[4-5]。石林榕等[6]以安息角為響應值,調節胡麻籽模擬模型時的滾動摩擦系數,最終得到胡麻籽顆粒真實滾動的摩擦系數。邢潔潔等[7]結合EDEM軟件中Hertz-Mindlin with JKR模型研究海南磚紅土壤的仿真參數,同樣以休止角作為優化目標,最終獲得仿真參數的最佳組合。侯占峰等[8]結合物理試驗和無滑移模型標定了冰草籽顆粒表面參數。Coetzee等[9]為了精確顆粒在DEM仿真時參數的精確度,利用堆積角試驗標定顆粒的摩擦系數。然而,關于辣椒籽的參數報告較少,且由于辣椒籽表面有粘結力的作用,簡單的無滑移模型難以揭示丸化顆粒的運動機理。

因此,本文提出一種JKR粘力模型,引入表面能參數用于表征辣椒籽表面的粘結力,通過一系列物理試驗和圖像處理技術得到辣椒籽的基本參數,然后運用優化設計試驗Plackett-Burman對所有影響因子進行顯著性篩選,建立3個顯著影響參數與響應值休止角的多元回歸模型,以最小相對誤差對模型進行優化,得到最佳參數組合。為辣椒籽丸?；峁﹨狄罁?提高離散元數值模擬的準確度和可靠性。

1 材料與方法

1.1 試驗材料和物性參數測定

1.1.1 基本參數測定

隨機選擇1 000粒辣椒籽,使用精密電子天平測量其質量,將辣椒籽等分成5組,計算其質量均值。辣椒籽的規格尺寸、大小和密度由數字游標卡尺(精度為0.01 mm)和量筒(精度為1 mL)測量,含水率使用干燥法測得。共進行10次試驗,結果如表1所示。

表1 辣椒籽基本物性參數Tab. 1 Basic physical parameters of pepper seeds

1.1.2 泊松比

辣椒籽的形狀不同,在長度和厚度方向上也有許多差異,泊松比很難通過傳統檢測方法確定。使用專業質構儀(TMS-Touch),在20 ℃的試驗溫度下對辣椒籽進行壓縮試驗,并使用靜態測量法獲得辣椒籽彈性模量。

為避免重復,此處僅以編號為1的辣椒籽彈性模量的測試過程為例。將1號辣椒籽水平放置在試驗的測試平臺上,并選擇直徑20 mm的探針。使用Texture Lab Pro軟件將傳感器范圍設置為0～250 N,將探頭下降速度設置15 mm/min、運行時間設置5 s,選擇90%形變量,開始后探頭緩慢下降,得到時間—載荷變化曲線,試驗結果如圖1所示。

圖1 時間—載荷散點圖Fig. 1 Time-load scatter diagram

由圖1可知,辣椒籽的時間—載荷曲線分為初始載荷段AB、線彈性段BC、峰值段CD和下降段DE。當質構儀探頭剛開始擠壓辣椒籽時,由于辣椒籽顆粒表面彎曲,豎直方向形變受壓力的影響較大;當處于線彈性變形階段時辣椒籽應力與應變比值為常數,如果解除應力,形變則會消失。當在BC段時,根據胡克定律,在線彈性變形下,辣椒籽的應變與應力成正比,且應力與應變的比例常數即為辣椒籽的彈性模量。

選取10次試驗數據,分別對AB線彈性變形階段進行分析。橫坐標為時間變量,縱坐標為壓力指數,截取3～4 s的數據擬合,獲得高度為0.6 mm、橫截面積為12.25 mm2辣椒籽壓力隨時間影響的折線圖,如圖2所示。

選取10組時間—載荷曲線的斜率的平均值代入式(1)計算可得到該類辣椒籽的彈性模量E為20.98 MPa。

(1)

式中: ΔF——壓力增量,N;

A——厚度方向橫截面的面積,mm;

Δt——時間增加量,s;

v——探頭的移動速度,mm/s;

v·Δt——位移增加量,mm;

h——原厚度,mm;

k——線彈性變形階段時辣椒籽的壓力增加量與時間增加量的比值,即壓縮時間—載荷曲線圖上線彈性形變階段的斜率。

泊松比在種子物性參數中極其重要,能夠體現種子一定的物理性能。通過泊松比原理對其求解,過程如下。

(2)

式中:μ——泊松比;

εx——橫向應變;

εy——縱向應變;

ΔL——橫向絕對變形量,mm;

L——橫向原長度,mm;

ΔH——厚度方向絕對變形量,mm;

H——厚度方向原長度,mm。

采用專業質構儀(TMS-Touch)作為加載裝置進行泊松比試驗。用游標卡尺量取單粒辣椒籽軸向和縱向的原尺寸。試驗開始時,將探頭下降速度設置成30 mm/min、運行時間設置成3 s,質構儀探頭開始緩慢下降,給辣椒籽施加載荷直至停止。記錄10次辣椒籽的縱向變形量并求均值。根據方程式(2)計算可得辣椒籽的泊松比為0.351±0.139,由式(3)可得辣椒籽的剪切模量為7.765 MPa。

(3)

式中:G——剪切模量,Pa;

E——彈性模量,Pa。

1.2 辣椒籽接觸參數測定

當辣椒籽進行一些模擬試驗時所涉及的參數有:辣椒籽間、辣椒籽與鋼板之間的恢復系數、靜摩擦系數和滾動摩擦系數。為了精確化所涉及的參數,以物理試驗為參照,后用模擬試驗標定。

1.2.1 碰撞恢復系數的測定試驗

辣椒種子丸?；鳂I時,顆粒間以及顆粒與轉釜之間會碰撞接觸,因此恢復系數的大小對丸?；鳂I影響較大,能夠直接影響辣椒種群的運動規律。根據牛頓碰撞定律,如果兩個物體材質不變,碰撞后與碰撞前相對速度的比值便是恢復系數[10]。

(4)

式中:v1′、v2′——碰撞后物體1和2的速度,m/s;

v1、v2——碰撞前物體1和2的速度,m/s。

圖3 恢復系數測試原理示意圖Fig. 3 Schematic diagram of restoration coefficient test principle1.落種孔 2.刻度尺 3.承載架 4.底座

(5)

式中:g——重力加速度,m/s2。

H′——種子下落高度;

h′——種子反彈高度。

根據以上原理進行試驗,利用UX100型高速攝像儀捕捉辣椒籽自由落體的過程,調節拍攝距離與試驗臺相距100 cm,整個系統主要有高速攝像儀、計算機和輔助光源。在本文中,辣椒籽均勻平鋪粘合在鋼板(轉釜材料)上,進行辣椒籽碰撞試驗,以確定辣椒籽與鋼板之間、辣椒籽間的恢復系數。辣椒籽從20 cm高的承載架上自由落下,與被測物體碰撞后反彈,并使用高速攝像儀拍攝整個辣椒籽碰撞和移動過程。

利用Matlab軟件將籽與鋼板碰撞過程的視頻進行分析,進行視頻關鍵幀提取和三維灰度處理,通過旋轉三維灰度曲面可以得到辣椒籽所在的區域以及環境的灰度值大概范圍,通過閾值分割和降噪處理得到辣椒籽的區域圖,接著運用region-props函數得到圖像中辣椒籽區域的面積、質心,最后使用循環語句分析所有幀得到辣椒籽的運動軌跡和速度如圖4所示。

(a) 關鍵幀提取

(b) 三維灰度處理

(c) 曲面旋轉尋找閾值

(d) 閾值分割

(e) 質心標定

(f) 種子軌跡曲線圖4 辣椒籽與鋼板碰撞視頻分析過程Fig. 4 Video analysis process of collision between seed and steel plate

圖4(f)為辣椒籽顆粒和鋼板接觸碰撞后的路徑軌跡,h′表示辣椒籽與鋼板碰撞時反彈的高度,自由下落高度為20 cm,由式(5)可得到辣椒籽的恢復系數。

按照以上過程,進行10次試驗,統計辣椒種間、辣椒籽與鋼板間的恢復系數,計算其平均值,可知辣椒籽之間恢復系數0.377±0.039,辣椒籽與鋼板之間恢復系數0.597±0.086。

1.2.2 靜摩擦系數測定試驗

靜摩擦系數對辣椒籽丸?；茸鳂I影響顯著。此處選用表面附有鋼板的自制測試儀,如圖5所示。

圖5 摩擦角測量儀Fig. 5 Friction angle measuring instrument

將重量為G1的辣椒種群穩定放置在斜面儀上(此時測試儀為水平靜止),搖動斜面儀的手柄,使種群所在平面的傾角α變大,當物體靜止于斜面上時,可以得到

G1sinα=γG1cosα

(6)

當傾角增大,使物體表面正好有滑動的趨勢時,即

tanα=γ

(7)

式中:γ——靜摩擦系數。

此刻,傾角α即為辣椒籽的摩擦角。通過摩擦角與摩擦系數之間的關系便可得到辣椒籽在剛體上的靜摩擦系數[11]。

辣椒籽間靜摩擦系數的測量,先使用粘性材料將辣椒籽平鋪固定在平面上,然后將辣椒籽群均勻置于該平面上,緩慢旋轉抬起測試儀,當有少數辣椒籽在平面上滑動時,記錄此刻測試儀角度,根據對應關系計算出辣椒籽間的靜摩擦系數。辣椒籽與鋼板間靜摩擦系數原理相同,將鋼板置于測量儀面上,操作一致。共進行10次試驗,并統計記錄摩擦系數均值。辣椒籽間的靜摩擦系數均值為0.564±0.045,辣椒籽與鋼板間的靜摩擦系數均值為0.325±0.016。

1.2.3 滾動摩擦系數測定試驗

當多數辣椒籽在傾斜提升時開始滾動,測量儀角度的正切值則是被測對象之間的滾動摩擦系數,原理與靜摩擦系數測定相同。但是,由于辣椒籽較小,形狀扁平,籽面粗糙,試驗過程中滾動系數變化顯著。經過10次試驗得辣椒籽間的平均滾動摩擦系數為0.734±0.036,辣椒籽與鋼板之間的平均滾動摩擦系數為0.513±0.036。但是由于物理試驗與仿真試驗之間存在或多或少的誤差,物理試驗測定的結果在仿真試驗中難以具備說服力。為此,使用休止角進行模擬和標定滾動摩擦系數[12]。

1.3 物理試驗休止角測定

由于辣椒籽規格較小,形狀扁平輕質,在形成自然的堆垛時,其斜邊是折線,休止角難以計算。因此,在累積時,使用無底圓筒填充材料,此為注射法[13]。試驗材料包括一個圓柱體和一個載荷板,圓筒由亞克力板制成,根據辣椒籽大小定義尺寸?？招膱A筒內徑至少是辣椒籽尺寸的4～5倍,其高度應是直徑得3倍[13]。為保證種子積累穩定,亞克力空心圓筒尺寸為:d=50 mm,H=150 mm。承載板由45鋼制成,規格為400 mm×400 mm×3 mm。

水平拍攝試驗堆積形成的自然休止角圖像,截取種子堆的單側圖像通過MATLAB圖像處理,經二值化、閾值分割以及膨脹腐蝕等處理,用邊界函數提取邊界信息,用最小二乘法調整邊界[14],得到方程的斜率kθ,處理過程如圖6所示,最終休止角θ表達式如式(8)所示。

(a) 原始圖像

(b) 圖像灰度化

(d) 圖像腐蝕膨脹

(e) 邊緣檢測

(f) 線性擬合圖6 休止角處理過程Fig. 6 Repose angle processing

θ=arctan|kθ|

(8)

利用無底圓筒法進行堆積試驗,除了辣椒籽的大小、形狀、自然參數和圓筒尺寸外,圓筒的上升速度也影響堆積角,為了防止堆積體過度散開并影響后續角度的準確性,圓通上升速度應小于0.07 m/s[15],上升速度越低,堆積體越穩定,但速度過低會延長仿真時間。鑒于上述情況,對圓筒不同上升速度下的辣椒籽開展了研究。如圖7所示,模擬與試驗中的休止角變化趨勢一致。當上升速度低于0.02 m/s時,休止角與上升速度成反比;當上升速度由0.02 m/s增加到0.03 m/s時,休止角急速減小;當上升速度在0.03 m/s后,休止角幾乎不變,趨于穩定。據上述分析,為提高模擬效率,在研究注入法模型在標定辣椒籽參數方面的有效性時,圓筒上升速度應為0.03 m/s。

圖7 圓筒上升速度對休止角的影響Fig. 7 Influence of cylinder rising speed on angle of repose

試驗時取1 000粒;將辣椒籽裝入圓筒內,然后放入水平裝料板上,以0.03 m/s的速度提高圓筒,在支撐板上形成穩定的堆垛,重復5次試驗,求出休止角的平均值。對得到的堆垛進行圖像處理,擬合結果,計算求得方程斜率kθ,得到自然休止角為28.88°、27.05°、26.77°、28.71°、28.59°,計算得均值為28.00°,標準偏差為1.01°。

2 基于JKR模型辣椒籽離散元參數標定

2.1 EDEM休止角仿真試驗

通過初步物理試驗測定了辣椒籽的主要物理參數,通過EDEM軟件進行模擬,利用多球體重疊球模型[16-18],重疊球模型能有效減緩“自動阻塞”的發生,更好的擬合邊界[19-20]。利用三維Catia建模軟件,根據辣椒籽的實際全局尺寸建立幾何模型(粒子模型)。然后將籽模型轉換為IGES中性格式并導入EDEM作為顆粒模板。在EDEM軟件中,辣椒籽顆粒模型填充獨特的球形顆粒,建立辣椒籽顆粒的模擬模型,如圖8所示。

圖8 辣椒籽球填充模型Fig. 8 Filling model of pepper seed ball

在EDEM中選擇接觸模型“Hertz-Mindlin with JKR”進行模擬試驗。模擬中參考物理試驗圓筒的內徑和高度建模,如圖9所示。

(a) 種子生成

(c) 種子堆形成圖9 辣椒籽堆積仿真過程Fig. 9 Chili seed accumulation simulation process

在圓筒頂端定義虛擬Polygon顆粒生成平面,選擇dynamic生成模式,效率為12 000個/s,共10 000個。根據Rayleigh時間步長規則,固定時間步長在5%～30%為佳,為考慮效率選擇30%固定步長,采取Euler積分,0.02 s記錄一次數據。模擬共進行10 s,顆粒生成總時間為1.5 s,網格尺寸為最小顆粒直徑的1.5倍[21]。辣椒籽穩定后,圓筒以0.03 m/s升起,待支撐板上辣椒籽堆穩定后,采用相同的圖像處理方法,獲得邊界斜率kδ,并使用式(16)計算模擬辣椒籽休止角δ。

δ=arctan|kδ|

(9)

2.2 優化試驗設計及分析

2.2.1 Plackett-Burman試驗設計及顯著參數篩選

試驗設計模擬中使用的參數對丸?；Ч胁煌挠绊?并非所有使用的參數都對丸?；囼炗酗@著影響,具有顯著影響的參數可用于模擬分析。由于影響參數太多,且與響應變量相關的參數的顯著性尚未確定,因此本研究使用軟件Design-Expert中Plackett-Burman Design設計模塊過濾選擇對辣椒籽運動軌跡具有較大影響的參數[22]。為避免試驗分析過于繁瑣,不計辣椒籽能直接精確測量的參數外,選取了以下7個參數:辣椒籽的表面能(A′)、種子—種子恢復系數(B′)、種子—種子靜摩擦系數(C′)、種子—種子滾動摩擦系數(D′)、種子—鋼板恢復系數(E′)、種子—鋼板靜摩擦系數(F′)、種子—鋼板滾動摩擦系數(G′)等作為變量,考慮誤差因素,在試驗中添加了四個虛擬參數H′、J′、K′和L′。因此選用N=16的Plackett-Burman設計表,并以辣椒籽的休止角作為響應值。表2中參數設置高低2個水平,即+1和-1,共在試驗的四個中心點進行了16次測試,試驗項目和結果如表3所示。

表2 Plackett-Burman Design因素水平Tab. 2 Factors and levels of Plackett-Burman design

表3 Plackett-Burman方案及結果Tab. 3 Scheme and results of Plackett-Burman design

對結果作方差分析,獲得所選參數的顯著性,如表4所示。辣椒籽—辣椒籽滾動摩擦系數D′和辣椒籽表面能A′的P<0.01,表明有極大影響;辣椒籽—鋼板滾動摩擦系數G′的P<0.05,表明有影響;其余仿真試驗參數的P>0.05,表明影響忽略不計。

表4 Plackett-Burman試驗參數顯著性分析Tab. 4 Significance analysis of Plackett-Burman test parameters

Plackett-Burman試驗參數顯著性分析表明辣椒籽的表面能A′、辣椒籽—辣椒籽滾動摩擦系數D′、辣椒籽—鋼板滾動摩擦系數G′對休止角的影響在0.05水平達顯著,P值越小影響越大,對休止角的影響順序是D′>A′>G′。

由圖10可知,D′、A′和G′因素對休止角影響在界值2.365以上,影響顯著。

圖10 PB Design帕累托圖Fig. 10 Pareto chart of PB Design

2.2.2 最陡爬坡試驗設計

為更有效地進行響應面試驗,找到最大梯度下降的中心點,以誤差為指標,對范圍進行尋優[19],爬坡試驗設計及結果如表5所示。由表5可知,相對誤差呈現先減小后增大的趨勢,在水平3時,準確度最高,可知水平3左右為理想范圍。因此,水平3被視為中心點,水平2、4分別為響應面試驗的低、高水平。其中相對誤差φ的表達式如式(10)所示。

表5 最陡爬坡試驗設計及結果Tab. 5 Design and results of steepest climbing test

(10)

式中:δ′——休止角測量值;

θ′——休止角目標值。

2.2.3 Box-Behnken試驗設計及回歸模型分析

1) Box-Behnken響應曲面試驗設計。在Box-Behnken試驗中,由爬坡試驗的結果可知,爬坡試驗中水平3為中心水平,水平2、水平4被視為設計試驗的-1和+1水平。根據BBD設計原則,由于有12個因素點以及5個中心點故選擇3因子3水平響應面,如表6所示,進行17次試驗,并通過分析獲得二次回歸方程以及參數的最優組合。

表6 Box-Behnken試驗設計及結果Tab. 6 Box-Behnken test design and results

2) 休止角回歸模型分析。對表6的試驗結果再次分析,建立完全二次模型,表7為模型方差分析表。

表7 完全二次模型分析結果Tab. 7 Complete quadratic model analysis results

模型各項方差分析表明,辣椒籽的表面能A′、辣椒籽—辣椒籽滾動摩擦系數D′、辣椒籽—鋼板滾動摩擦系數G′對休止角影響極其顯著,同時因素A′的曲線效應也顯著,交互項D′G′的影響較為顯著,其余項均不顯著。

仿真試驗休止角二階回歸方程為

δ′=28.098+3.153A′+2.263D′+1.955G′-

0.103A′D′+0.688A′G′-1.283D′G′-

2.293A′2-0.653D′2-0.628G′2

(11)

3) 因素效應分析。如圖11所示,隨著表面能A′水平增加休止角先快速增加后緩慢下降,辣椒籽間滾動摩擦系數、籽與鋼板間摩擦系數水平與休止角正相關。

(b) 辣椒籽—辣椒籽滾動摩擦系數

(c) 辣椒籽—鋼板滾動摩擦系數圖11 單因素影響趨勢圖Fig. 11 Effect-tendency chart of single factor

4) 不同因素間交互作用分析。響應面分析圖是由休止角和3個顯著因素構成的三維曲面圖(圖12),顯示了當任一因素(A′、D′和G′)處于零水平時,其余兩個因素對休止角的影響,傾角越大,交互作用越強[23]。等高線圖是曲面上相同的因素值在底面上形成的曲線(圖12),等高線輪廓圓扁程度體現了兩因素之間的交互作用的強弱,等高線越密集,交互作用越強[23]。

圖12 響應面交互作用分析圖Fig. 12 Interaction analysis diagram of response surface

由圖12可知,辣椒籽—辣椒籽滾動摩擦系數D′和辣椒籽—鋼板滾動摩擦系數G′都與辣椒籽的表面能A′有交互作用。在表面能A′低水平時,交互作用較強;在表面能A′高水平時,交互作用較弱。

5) 參數尋優。本研究對休止角建立二次多項式模型,利用Design-Expert根據實際休止角28.00°為指標對回歸模型進行優化,優化方案唯一,如表8所示。

表8 相對誤差最優解Tab. 8 Optimal solution of relative error

由表8可知,在該參數組合下得到休止角為27.3°,與試驗的相對誤差為2.5%,可取性(Desirability)為1,說明模型可靠。因此辣椒種子丸化加工時的表面能A′為0.31 J/m2,辣椒籽間滾動摩擦系數為0.75,辣椒籽與鋼板的滾動摩擦系數為0.60。

2.2.4 仿真模擬驗證

為了驗證標定參數的準確度,利用辣椒籽標定后的參數進行模擬,并用辣椒籽仿真的結果圖像與試驗圖像驗證,如圖13所示。反復模擬三次得到辣椒籽自然休止角為25.90°、26.98°、27.25°,平均值26.71°,與模擬休止角誤差為2.16 %,與試驗休止角誤差為4.61%,說明辣椒籽所標定的物理參數具有可靠性及參考意義。

(a) 試驗積累

(b) 模擬積累圖13 辣椒籽仿真模擬與試驗堆積體圖像對比Fig. 13 Comparison between simulation and experimental accumulation of chili seed

3 結論

1) 首先通過試驗確定辣椒籽便于測量的基本參數(總粒徑、千粒重、密度、含水率、彈性模量、剪切模量和泊松比);利用UX100型高速攝像儀和自制的斜面儀測試辣椒籽種間碰撞恢復系數、靜摩擦系數、滾動摩擦系數的均值分別為0.377±0.039、0.564±0.045、0.734±0.036;辣椒籽與鋼板間碰撞恢復系數、靜摩擦系數、滾動摩擦系數的均值分別為0.597±0.086、0.325±0.016、0.513±0.036。

2) 以辣椒籽為研究對象,基于離散元仿真軟件EDEM中的“Hertz-Mindlin with JKR”粘力模型,選擇注入法,進行休止角累積試驗,同時對辣椒籽在接觸過程中的物理性能參數進行標定。

3) 利用Plackett-Burman試驗過濾出對休止角影響較大的參數(辣椒籽種間的滾動摩擦系數、靜摩擦系數、JKR表面能),而后通過最陡爬坡試驗尋找顯著項的最優區間,并通過響應面Box-Behnken試驗以最佳休止角目標值(28.00°)進行顯著項交互作用分析。最后通過模型優化得到最優擬合度(辣椒籽種間滾動摩擦為0.75、辣椒籽表面能為0.31 J/m2、辣椒籽—鋼板滾動摩擦系數為0.60)及擬合方程,用于預測休止角。利用試驗結果得出的參數組合進行驗證試驗,與試驗休止角相對誤差為4.61%,說明了模型的可靠性。