離散元聚苯乙烯顆粒細觀參數辨識

龐紹銳，楊麗紅，王雙園

(上海理工大學機械工程學院，上海 200000)

0 引 言

真空擔架由聚氯乙烯（PVC）軟膜熱熔后經密封形成腔體，內部填充變形量大，密度小的聚苯乙烯(EPS)顆粒。腔體被抽真空時，腔體內的EPS 顆粒會發生流動并轉變成“類固體”狀態[1]，隨之產生可變的剛度[2]。由于人體各部分最大承受壓力不同，通過研究顆粒與軟膜抽真空過程的力學特性，可以研制不同結構的腔體，以便擔架抽真空后可以更適合包裹和固定人體。離散元法（discete element method）是一種基于分子動力學的顆粒物料分析方法[3]，通過使用離散元軟件EDEM 與多體動力學仿真軟件Recurdyn 耦合，可以分析腔體壓縮顆粒后的結構變化特征，為擔架設計提供參考依據。由于顆粒的細觀參數直接影響著宏觀的力學特性變化，準確的離散元參數是仿真模擬的基礎，仿真之前需辨識出EPS 顆粒與PVC 軟膜的細觀參數。

一般離散元理論主要研究對象為巖石、沙土、混泥土等[4]，對于顆粒的的細觀參數測定，國內外學者基于離散元方法對土壤、糧食之類顆粒物質進行了大量的研究。GONVIA[5]采取休止角為模擬指標發現，采用具有摩擦等特性的離散元能很好地增加顆粒模擬效果，顆粒物質的細觀參數對其宏觀運動有著直接影響。以宏觀休止角為響應值，韓燕龍[6]通過對休止角的仿真，對稻谷之間的滾動摩擦系數進行了仿真標定。王美美[7]等通過使用響應面優化方法對玉米籽粒進行離散元參數標定，提出多因數對響應休止角的最優參數標定。針對于顆粒的參數組合對休止角的影響，張銳等[8]使用EDEM 軟件對顆粒之間相互作用參數進行標定，提出一種系統標定沙土顆粒相互作用參數的方法，顆粒的細觀參數組合直接影響著堆積的休止角值。賈富國[9]等人利用離散元方法對稻谷顆粒堆積的休止角進行預測稻谷的細觀參數。目前很多類型的顆粒測定采用經驗值，這種方法往往導致實驗數據不準確，研究結果的誤差較大。EPS 顆粒相比于常見的沙土、巖石、糧食等物質特性更加特殊，需尋找一種不同的測定方法，且在研究顆粒散體動態加載的過程時，理想化使用有限元方法研究EPS 顆粒會忽略顆粒的不連續性[10]，造成散體運動研究的不準確性。

結合前文所述的相關研究，本文以EPS 顆粒與PVC 軟膜為對象，設計一套方案來測定辨識EPS顆粒之間，EPS 顆粒與PVC 膜間的細觀接觸參數。研究基于Hertz-Mindlin 接觸模型，以多種試驗和仿真模擬結合的方法對顆粒進行離散元參數的測定，旨在為仿真研究顆粒在真空壓縮過程產生的力學機理提供準確的接觸參數。

1 離散元模型

1.1 顆粒模型

真空擔架填充的EPS 顆粒，顆粒外觀為球形，通過對顆粒進行網格篩選可知：粒徑在3～3.5 mm占總量的20.5%，3.5～4 mm 占總量的21.7%，4～4.5 mm 為35.6%，4.5～5 mm 為22.2%。顆粒模型選擇EDEM 軟件內置的球形顆粒模型如圖1 所示。在仿真時，將整個粒徑分布分為4 個粒徑段，可認為每個粒徑段均勻分布，在仿真時以粒徑3.25 mm、3.75 mm、4.25 mm、4.75 mm 作為顆粒的直徑，分別按照所屬粒徑段的比例進行顆粒工廠設置。

圖1 EPS 顆粒模型

1.2 顆粒接觸模型

EPS 顆粒之間沒有粘結，且顆粒接觸面很小，在接觸面上僅發生彈性變形。EPS 顆粒的運動是法向力與切向力耦合的一個狀態[11]，很難得到兩個接觸顆粒的法向力和切向力分析解。因此，EPS 顆粒仿真分析時采用EDEM 內置的Hertz-Mindlin(no slip)接觸模型來描述。EPS 顆粒接觸模型如圖2 所示，該模型將顆粒之間的法向和切向接觸力簡化成為彈簧和阻尼器的并聯。法向力的計算采用Hertz 接觸理論計算，切向力與切向位移關系采用Mindlin 接觸理論[11]。

圖2 接觸力學模型

在EPS 顆粒發生運動時，顆粒之間碰撞擠壓產生的法向力N為

式中：σ——法向重疊量；

R?——有效顆粒半徑；

E?——有效彈性模量。

其中R?、E?計算公式為

式中：E1、E2——顆粒1 和顆粒2 的彈性模量；

ν1、ν2——顆粒2 和顆粒2 泊松比。

當兩顆粒間重疊量為 ?σ 時，法向力增量 ?N由式(1)可得

兩接觸面表面切向位移δ 時，切向力為T

其中Sτ為切向剛度。

在接觸表面產生 ?δ的切向位移增量，兩者之間切向力增量 ?T為

式中：μ——顆粒表面靜摩擦因數；

G?——有效剪切模量；

a——接觸半徑；

k=0，1，2——對應切向力的加載，卸載和卸載后重新加載，當

顆粒和顆粒，顆粒和PVC 軟膜之間法向阻尼力Fnd和 切向阻尼力Fτd滿足

式中：m?——等效質量；

vτrel——切向相對速度；

vnrel——法向相對速度；

Sn——法向剛度，N/mm；

ε——恢復系數；

β——阻尼比。

顆粒之間存在的滾動運動，該運動主要通過接觸表面滾動摩擦力矩Ti和 切向力矩Tτ決定。

式中：fr——顆粒之間的滾動摩擦因數；

Ri——接觸點到質心的距離；

ω——顆粒在接觸點處的角速度。

2 接觸參數的測定

離散元的部分接觸參數在誤差允許范圍內可通過物理實驗測量，為提高參數的準確性，利用離散元軟件EDEM 進行仿真模擬和實驗相結合來測定參數。所需的EPS 顆粒PVC 軟膜本征參數[12-13]如表1 所示。

表1 本征參數

2.1 EPS 顆粒與PVC 軟膜碰撞恢復參數測定

碰撞恢復系數是表現一個物體在受到碰撞時變形恢復的能力，在速度和沖量定義方法中，恢復系數定義為兩物體碰撞前后在接觸點相對分離速度與相對接近速度之比，理論計算公式為

式中：H1——反彈之后的最大高度；

H2——初始顆粒降落高度。

相關研究表明，顆粒在下落實驗測恢復系數時，下落高度影響著碰撞能量損失，恢復系數測定偏小[14]。為減小測定 ε1時能量損失帶來的誤差，如圖3（a）選擇將EPS 顆粒從高度H2=10 cm 處進行自由落體，與實驗臺上平鋪的PVC 軟膜進行碰撞。通過高速攝像機記錄顆粒最高的彈起高度H1，實驗進行六次，計算平均值，得到反彈起的高度實驗均值為H1=2.55 cm。

圖3 碰撞恢復系數

以實驗值為依據進行多次如圖3（b）仿真實驗,設置顆粒不同的碰撞恢復系數，得到顆粒反彈之后的軌跡，實驗結果如圖4 所示，可以看出隨著EPS顆粒與PVC 軟膜之間恢復系數 ε1降低，EPS 顆粒的初次反彈高度H1逐漸減小。

圖4 顆粒降落位置曲線

根據實驗的數據，設計仿真彈跳實驗，以彈性恢復系數 ε1為 變量，ε1范圍選擇0.44～0.52，步長為0.02。每組進行三次仿真，以反彈高度H1平均值為結果，仿真實驗的結果擬合方程為

決定系數r2=0.998，該值接近1 表明擬合可靠度高。將實驗測量值2.55 cm 代入式（14）中，得到彈性恢復系數 ε1=0.506，由公式（14）驗證可得顆粒反彈的預測高度為2.56 cm，與實驗的誤差為0.38%，誤差極小，則可確定仿真參數 ε1=0.506。

2.2 EPS 顆粒與PVC 軟膜靜摩擦因數測定

使用斜面滑測定EPS 顆粒與PVC 軟膜之間的靜摩擦系數 μ1時，顆粒的重力在斜面上分解為平行于斜面的力F1和 垂直于斜板的力F2，顆粒與軟膜之間的靜止摩擦力為f。當斜面的傾角 α大于靜摩擦臨界角的時候，F1>f時發生滑移，此時通過公式(15)可計算μ1

在實驗與仿真的時候采取多顆粒粘結成塊的方式，并將其放置在軟膜表面如圖5（a）。實驗與仿真時，斜面旋轉速度一致保持為0 .05 rad/s，仿真中斜面傾角通過時間乘以轉速來計算。設置仿真時間步長為Rayleigh（TR）時間的8%，設置采集點間隔為0.001 s，可以精確采集每個點顆粒板的運動狀態并計算斜面角度α。

圖5 EPS 顆粒和PVC 表面靜摩擦實驗

進行斜面實驗時，當顆粒板塊發生滑動，記錄此時斜面傾斜角，重復測量5 次之后，取平均值為42.15°。在EDEM 仿真軟件Creator 界面下創建New Section (0.06 m×0.06 m)來模擬角測量儀斜面上的PVC 軟膜，仿真過程中改變顆粒與PVC 軟膜接觸的靜摩擦因數，進行6 組仿真實驗，記錄每組仿真實驗滑動斜角，實驗結果如表2 所示。

表2 μ 1選取與傾角結果

對仿真實驗中的數據進行曲線擬合，得到 μ1與y的二次多項式如圖6 所示，曲線方程為：

圖6 靜摩擦系數與傾斜角擬合曲線

該方程式的決定系數為r2=0.99891，接近1，表明此擬合多項式極其符合EPS 顆粒斜面實驗的系數與摩擦因數之間的關系。將室內斜面實驗所得到的平均值42.15°代入方程可得 μ1=0.885。由公式(15)可得靜摩擦角的預測值為0.905，仿真與理論誤差2.2%，表明仿真結果與試驗結果一致，則可確定仿真參數 μ1=0.885。

2.3 EPS 顆粒之間滾動摩擦系數測定

離散元的出現源于分子動力學，其將整個介質看作是一系列離散的獨立運動的單元組成，單元本身就具有一定的幾何形狀、物理、化學性質。基于分子動力學仿真模擬出EPS 顆粒單元的滾動摩擦系數f1，模擬過程如圖7 所示，兩個相同的聚苯乙烯球形顆粒單元作為對象，對顆粒進行點接觸滾動運動模擬。聚苯乙烯是一種只含有碳氫元素的聚合物(C8H8)，使用Material Studio 軟件建立顆粒運動模型，構建兩個材料相同，半徑20 nm，密度均為20 kg/ m3的球形單元。在分子動力學模擬中，將10個聚苯乙烯分子首尾相連組成分子鏈，球形單元由聚苯乙烯分子鏈填充，分子鏈填充量使球形單元密度達20 kg/ m3為 止。對顆粒1 施加載荷Fz=10 nN，方向為z軸負方向，賦予顆粒1 沿V方向速度為0.1 m/s 的運動，模擬運動過程中小球1 和小球2 在點接觸時發生滾動運動。

圖7 顆粒單元滾動碰撞工況

本實驗中的滾動摩擦系數K計算公式為

式中：T——顆粒單元接觸時的切向力；

N——接觸時候的法向力；

R——顆粒1 的半徑。

模擬溫度通過Nose-Hoover 熱浴法設定在300 K，耦合時間為0.1 ps，壓強保持為1 atm，分子動力學模擬通過Lammps 軟件完成。輸出仿真過程中顆粒之間的滾動摩擦系數值f1，除去突變的數值，剩余滾動摩擦系數一直維持在0.02，則在接下來的離散元仿真模擬時f1為0.02。

2.4 堆積實驗分析和測定參數

前文已經對 μ1、ε1、f1完成測定。下文對EPS 顆粒與PVC 滾動摩擦因數f2，EPS 顆粒之間碰撞恢復系數 ε2、靜摩擦因數 μ2進行分析和測定。為避免物理測定帶來的較大誤差，通過顆粒堆積實驗與仿真相結合，以休止角為響應，根據最陡爬坡實驗確定參數 μ2、ε2、f2最佳范圍，結合Box-Behnken 實驗建立休止角與參數的二階回歸模型，對參數進行求解。

2.4.1 堆積實驗休止角測量

堆積實驗如圖8(a)所示，材料選擇EPS 顆粒，底板選用PVC 軟膜，圓柱管為內徑41 mm，長度為190 mm 玻璃管。通過電機伸縮桿末端固定玻璃管并勻速提升，EPS 顆粒自然滑落并堆積如圖8（b）。

圖8 顆粒堆積實驗

堆積實驗采集到的圖像如圖9(a)所示，為保證測量休止角的精準性，使用Matlab 讀取實驗圖像，依次對圖像進行去噪處理，灰度處理，二值化處理，最后提取輪廓邊界點，利用最小二乘法對邊界點進行直線擬合，擬合的直線斜率即為所要測得的休止角正切值。利用同種測量方法實驗五次，求斜率均值為0.25594，則EPS 顆粒在PVC 軟膜上的休止角θ為14.356°。

圖9 休止角測量

2.4.2 最陡爬坡實驗

使用離散元仿真軟件EDEM 進行如圖10 顆粒堆積仿真實驗，通過改變 ε2、f2、μ2三個參數，對仿真實驗的休止角θ 進行測定，顆粒的其余參數均采用已測定數值。

圖10 顆粒堆積仿真實驗

在一定范圍內f2、μ2增大會使休止角增加[15]，ε2的增大會導致休止角減小，結合細觀參數對休止角的影響研究結論[16]，確定爬坡試驗的方向和步長。通過前期的預仿真，ε2、f2、μ2的參數選擇區間以及實驗結果如表3 所示?？芍抡嫘葜菇桥c試驗測量休止角最小誤差在第c、d、e 組區域，由d 組實驗可以看出當 μ2為 0 .4，f2為 0.44，ε2為0.35 時，仿真休止角值與實驗值比較接近。

表3 最陡爬坡實驗設計及其結果

2.4.3 響應面分析

結合最陡爬坡實驗選定d 組參數為中心點，確定三個參數合理范圍，根據Box-Behnken 中心組合設計原理[17]，以–1,0,1 為低、中、高水平進行三因素三水平的響應面實驗，響應值為休止角θ，如表4 所示。

表4 仿真實驗因素與水平

響應面實驗設計及結果如表5 所示，用軟件Design-expert 對實驗結果進行回歸分析得到回歸方程為：

表5 Box-Behnken 實驗設計及結果

對休止角的模型進行方差分析如表6 所示，方程模型極顯著（P<0.0001），失擬項不顯著（P=0.7612>0.05），模型穩定。相關系數r2值為0.9853，校正決定系數=0.9665，預測決定系數=0.9282，三者都大于0.8，信噪比為26.142>4 表明模型精確度很高。

表6 實驗模型方差分析

根據表7 結果可知，在保證模型的準確性前提下，剔除對休止角影響不顯著的項（P>0.05），優化后得到新的模型為

表7 優化模型方差分析

r2=0.9804，，信噪比41.86，模型較之前精確度更高。

2.4.4 最優參數的選取

利用Design-Expert 三維響應表面，設置優化目標為休止角接近14.356°，誤差值優化目標為最小，對模型進行求解，選取最優參數值,顆粒與顆粒靜摩擦因數 μ2，顆粒與PVC 滾動摩擦f2，顆粒與顆粒之間恢復系數 ε2的約束函數是：

得到最優的參數組合為：顆粒與顆粒之間靜止摩擦因數 μ2=0.40，恢復系數 ε2=0.43，顆粒與PVC軟膜之間的滾動摩擦因數f2=0.45。用最佳的參數組合進行三次休止角仿真實驗，得到的仿真休止角度數均值14.352°。由于堆積試驗休止角的角度為14.356°，仿真與實測的休止角誤差為0.124%，表明測定的參數組合與試驗參數相符。

3 試驗驗證

為驗證前文測定的顆粒與顆粒之間靜止摩擦因數、滾動摩擦因數、碰撞恢復系數，顆粒與PVC 軟膜之間的靜止摩擦因數、滾動摩擦因數、碰撞恢復系數的準確性，將六個參數值：0.40、0.02、0.43、0.885、0.45、0.506 代入EDEM 仿真軟件，進行側壁坍塌仿真實驗[18]。

如圖11(a)所示，實驗器材由箱體（長90 mm,寬90 mm，高90 mm）和玻璃擋板組成，為控制實驗對象的接觸參數數量，箱體墻壁以PVC 軟膜覆蓋。實驗時以較慢速度向上抬升玻璃擋板，使EPS 顆粒在PVC 軟膜上充分擴散，當顆粒散體穩定之后進行水平拍照，將采集到的圖片進行2.4 節相同步驟圖像處理，重復實驗5 次并測量休止角，取均值為18.08°。將上述測定的參數帶入側壁坍塌仿真實驗如圖11(b)，測量仿真實驗的坍塌休止角，實驗進行5 次。測量休止角均值為18.43°，與試驗所測量的休止角相對誤差為1.9%，表明仿真實驗中的結果和實際情況吻合，參數驗證合理。由此可得上述參數f1、μ1、ε1、f2、μ2、ε2適用于EPS 顆粒和PVC 軟膜的離散元分析，在類似的彈性顆粒和柔性體的接觸該測定方法同樣適用。

圖11 側壁崩塌實驗

4 結束語

1）以EPS 顆粒和PVC 軟膜為研究對象，基于Hertz-Mindlin 接觸模型，利用實驗和仿真進行顆粒與軟膜之間的參數測定，得出EPS 顆粒與PVC 軟膜的靜摩擦因數為0.885，恢復系數為0.506。通過分子動力學的仿真模擬方法，標定出顆粒與顆粒之間的滾動摩擦因數為0.02。

2）對顆粒與顆粒的碰撞恢復系數，顆粒與PVC軟膜的滾動摩擦系數，顆粒與顆粒的靜止摩擦系數進行參數測定：通過設計最陡爬坡法逼近最大響應區域，結合Box-Behnken 響應面實驗，以顆粒堆積實驗休止角為響應值，選取與實驗測量休止角誤差值最小的參數組合為最優解，得三者參數分別為0.43、0.45、0.40。

3）通過側壁崩塌實驗驗證參數的可靠性，發現崩塌實驗的仿真值和實驗值休止角結果接近，誤差為1.9%，說明接觸參數的測定準確，也為研究柔性體和顆粒細觀接觸參數提供了具體的測定方法。