基于JKR接觸模型的微米級顆粒離散元參數(shù)標(biāo)定

韓 偉，王紹宗，張 倩，田宇航

(機(jī)械科學(xué)研究總院 先進(jìn)成形技術(shù)與裝備國家重點(diǎn)實驗室，北京 100044)

微米級粉體顆粒物質(zhì)作為一種常見的工業(yè)原料，在生產(chǎn)應(yīng)用中不可避免地涉及其倉儲、 輸送、 稱量、 包裝等過程[1-3]，但是，由于部分微米級顆粒物質(zhì)流動性較差，在進(jìn)行倉儲、 輸送等過程中極易出現(xiàn)架橋、 搭拱、 堵結(jié)輸送管路等情況[4]，不利于生產(chǎn)活動順利進(jìn)行，因此，有必要探究此類微米級粉體顆粒物質(zhì)的物性參數(shù)，準(zhǔn)確把握其流動規(guī)律，以便于有針對性地開發(fā)相應(yīng)的機(jī)械裝置來解決此類問題。

近年來，在顆粒物質(zhì)的物性參數(shù)、流動規(guī)律研究方面，離散單元法應(yīng)用較為廣泛[5]。離散單元法基于牛頓第二定律描述顆粒系統(tǒng)中每個粒子單元的運(yùn)動狀態(tài)，進(jìn)而推導(dǎo)整個系統(tǒng)的流動規(guī)律，可為顆粒物質(zhì)倉儲、 輸送、 稱量等裝置的設(shè)計提供理論參考。對于顆粒物質(zhì)而言，開展離散元模擬的首要問題是獲取物理實驗難以測得的顆粒本征參數(shù)、 基本接觸參數(shù)以及接觸模型參數(shù)等離散元參數(shù)。目前，比較通用的方法是通過“參數(shù)標(biāo)定”實驗確定，其實質(zhì)是建立物理實驗與仿真實驗之間的二階回歸數(shù)學(xué)模型，并基于某一響應(yīng)值求取該模型最優(yōu)解的過程，其中最優(yōu)解即可視為待標(biāo)定的離散元參數(shù)。針對沙土顆粒、 小麥顆粒、 蚯蚓糞基質(zhì)這類粒徑較大的顆粒物質(zhì)，張銳等[6]、 劉凡一等[7]、 羅帥等[8]分別以堆積角為響應(yīng)值展開了參數(shù)標(biāo)定實驗：而對于諸如面粉顆粒等粒徑較小的顆粒物質(zhì)，李永祥等[9]基于Feng等[10]提出的精確縮尺理論，對顆粒進(jìn)行粒徑放大并開展了參數(shù)標(biāo)定實驗；針對不同含水率的顆粒物質(zhì)，Grim等[11]以堆積角為響應(yīng)值，分別對干、濕顆粒的滾動摩擦系數(shù)進(jìn)行了參數(shù)標(biāo)定；此外，Li等[12]還研究了顆粒摩擦系數(shù)的滑板滑動實驗測定方法，并進(jìn)行離散元模擬驗證，取得了較為理想的結(jié)果。

目前，應(yīng)用離散單元法進(jìn)行的參數(shù)標(biāo)定實驗，其對象更傾向于大顆粒物質(zhì)，粒徑多為毫米級及以上，即使是微米級顆粒也多進(jìn)行精確縮放等操作，然而嚴(yán)格的精確縮放過程并不會提升計算機(jī)仿真效率，反而增加了理論計算的工作量。為了探究微米級顆粒原始尺寸下的離散元參數(shù)標(biāo)定問題，本文中選取微米級粉體顆粒物質(zhì)，印染行業(yè)中粒徑分布在10～400 μm的一種粉體活性染料3BF Hong開展參數(shù)標(biāo)定實驗。由于物料堆積角能夠反映物料流動特性，因此參數(shù)標(biāo)定實驗選用物料堆積角作為響應(yīng)值。此外，實驗中接觸模型的選擇是一個重要問題。粉體活性染料在流動過程中，由于顆粒間的粘結(jié)作用，物料流動不暢，適用于JKR模型加以描述。JKR模型即“Hertz-Mindlin with JKR”接觸模型，作為一種凝聚力接觸模型，JKR模型在計算中引入了顆粒間表面能的概念，可以最大程度地表征濕顆粒之間的接觸特性，多用于表示濕式和細(xì)小顆粒間的粘結(jié)作用[13]。本文中應(yīng)用離散元分析軟件EDEM，選取JKR模型作為接觸模型，對粉體活性染料3BF Hong進(jìn)行參數(shù)標(biāo)定實驗，以實測得到的粉體活性染料堆積角角度為響應(yīng)值，通過Plackett-Burman實驗、 2次最陡爬坡實驗以及Box-Behnken實驗，推導(dǎo)顯著影響仿真實驗中染料堆積角大小的不同離散元參數(shù)與堆積角的二階回歸模型，求取回歸模型最優(yōu)解即為染料顆粒的離散元參數(shù)，所得數(shù)據(jù)可為后續(xù)研究顆粒流動規(guī)律，設(shè)計顆粒物質(zhì)輸送、稱量裝置提供數(shù)據(jù)參考。

1 參數(shù)測定與模型建立

選取粉體活性染料3BF Hong為研究對象，依據(jù)參數(shù)標(biāo)定實驗要求，測定其堆積角作為參數(shù)標(biāo)定實驗的響應(yīng)值；并對該種粉體活性染料進(jìn)行理化分析，根據(jù)其密度、粒徑分布以及表面微觀特征等信息建立離散元模型。

1.1 JKR接觸模型

JKR模型的切向彈性接觸力、切向耗散力和法向耗散力的計算基于“Hertz-Mindlin (No Slip)”接觸模型理論，而法向彈性接觸力基于Johnson-Kendall-Roberts理論，用法向重疊量、相互作用參數(shù)以及顆粒表面能計算法向彈性接觸力為

(1)

(2)

式中：FJ為法向彈性接觸力，N；E*為當(dāng)量楊氏模量，Pa；R*為當(dāng)量半徑，m；α為接觸半徑，m；Δγ為接觸顆粒表面能，J/m2；δ為法向重疊量，m。其中

Δγ=γi+γj-γij

(3)

式中：γi、γj分別為2個顆粒表面的單位面積粘附表面能，γij為界面表面能，當(dāng)相接觸的顆粒材料相同時，γi=γj且γij= 0，所以有Δγ=2γ。

JKR模型對于摩擦力的計算取決于 JKR 法向力的正向排斥部分，這與“Hertz-Mindlin”接觸模型不同，JKR模型在接觸力的凝聚力分量更大時，提供一個更大的摩擦力，Gilabert等[14]對強(qiáng)凝聚力情況下摩擦力模型修正的重要性和優(yōu)勢已做過相應(yīng)探討。JKR接觸理論中顆粒間粘結(jié)示意圖如圖1所示。

1.2 堆積角的測定

參考目前廣泛使用的堆積角測定方法，基于點(diǎn)源法設(shè)計了如圖2所示堆積角測定裝置，主要由鐵架臺、防塵罩、漏斗及堆積底座組成。對所得粉體活性染料料堆進(jìn)行數(shù)字圖像處理，讀取料堆中部比較平直的邊界，通過最小二乘法等數(shù)學(xué)手段，對邊界進(jìn)行擬合計算得出料堆的堆積角[15]，處理過程如圖3所示。

擬合直線函數(shù)斜率與堆積角角度值的轉(zhuǎn)換為

(4)

式中：β為測得的堆積角；k為擬合直線斜率。經(jīng)過5次實驗，求得堆積角的平均值為38.57°。

Ri、 Rj—顆粒半徑; α—接觸半徑; δ—法向重疊量。圖1 JKR接觸理論中顆粒間粘結(jié)示意圖Fig.1 Schematic diagram of adhesion between particles圖2 堆積角測定儀Fig.2 Accumulation angle tester

a)原始圖像b)圖像灰度化c)圖像二值化d)圖像濾波e)邊緣檢測f)線性擬合圖3 堆積角圖像處理Fig.3 Image processing process of accumulation angle

1.3 染料顆粒參數(shù)測定及離散元模型建立

經(jīng)北京市理化分析測試中心測定，粉體活性染料3BF Hong的平均密度為0.570 g/cm3，顆粒粒徑分布如圖4所示。其中D10=42.36 μm、D50=123.09 μm、D90=229.80 μm，平均粒徑為123.09 μm，體積中值直徑為130.42 μm。粒徑分布在45～255 μm區(qū)間的顆粒數(shù)量分?jǐn)?shù)為84.66%，能夠代表整個粒子系統(tǒng)參與計算。將45～255 μm的粒徑劃分為5個區(qū)間，選取代表粒徑，計算所占比率，據(jù)此在軟件中創(chuàng)建顆粒，同時結(jié)合仿真預(yù)實驗確定各顆粒的生成時間及總數(shù)，如表1所示。

圖4 染料顆粒粒徑分布圖Fig.4 Dye particle size distribution chart

表1 顆粒參數(shù)信息表Tab.1 Parameters of created particles

對染料顆粒進(jìn)行顯微鏡拍照，得到如圖5所示染料顯微鏡圖像。由圖可知，其主要由大小不均、較為完整的球狀顆粒組成，偶有夾雜不規(guī)則的顆粒碎片，因此，在建立染料顆粒離散元模型時，可以使用單一球狀顆粒代替染料顆粒進(jìn)行仿真計算。

圖5 染料顆粒顯微鏡圖像Fig.5 Micrograph of dye particles

由于粉體物料堆積角的大小與測定儀器尺寸及待測堆積物料質(zhì)量無關(guān)[6]，為了簡化計算，按照1∶5的比例縮放堆積角測定儀，且只在仿真軟件中導(dǎo)入測定儀的漏斗和底座模型。按照表1給出的具體數(shù)據(jù)設(shè)置顆粒工廠，顆粒生成方式為“Dynamic”。由于不同的參數(shù)組合會對應(yīng)力波在顆粒中的傳播產(chǎn)生影響，導(dǎo)致各仿真中瑞利時間步不同[7]，因此在所有仿真中，時間步統(tǒng)一取瑞利時間步的30%，網(wǎng)格尺寸選取為最小球形單元尺寸的2倍，為保證仿真結(jié)束后顆粒充分靜止，仿真時間設(shè)定為0.8 s。待所有模型建立完成后進(jìn)行后續(xù)實驗。

2 參數(shù)標(biāo)定實驗設(shè)計及分析

參數(shù)標(biāo)定實驗主要包括3個部分： 1)Plackett-Burman實驗，用以選擇待標(biāo)定因子中對堆積角影響最顯著的因子； 2)最陡爬坡實驗，用以確定顯著性因子的最優(yōu)區(qū)間； 3)Box-Behnken實驗，用以建立顯著性影響因子與響應(yīng)值的二階回歸方程，從而求解出顯著性因子的具體參數(shù)值。

2.1 待標(biāo)定參數(shù)設(shè)置及顯著性參數(shù)判定

目前，通過已有的研究無法直接獲得染料顆粒各參數(shù)的具體范圍，本文中根據(jù)EDEM軟件自帶的“GMEE”數(shù)據(jù)庫，參考國內(nèi)外文獻(xiàn)中粉體顆粒、不銹鋼材料的離散元參數(shù)設(shè)置，初步確定了待標(biāo)定參數(shù)表及部分參數(shù)的取值范圍，對于難以查詢且受顆粒粒徑、含水率影響較大的JKR 表面能數(shù)據(jù)，其取值范圍由仿真預(yù)實驗確定，詳見表2。

表2 離散元模擬參數(shù)表Tab.2 Parameters required in discrete element methodsimulation

Plackett-Burman實驗可以通過考察目標(biāo)響應(yīng)與各因子間的關(guān)系，比較各個因子兩水平間的差異，從多個因子中選取對實驗指標(biāo)有顯著影響的因素。故此，以染料顆粒的堆積角為響應(yīng)值，設(shè)計Plackett-Burman實驗，按照待標(biāo)定參數(shù)顯著性高低進(jìn)行篩選，實驗中插入2個虛擬參數(shù)K、L用于誤差估計，根據(jù)表2各參數(shù)范圍進(jìn)行參數(shù)高低水平設(shè)置，擬定各參數(shù)的參數(shù)編碼，如表3所示。

表3 Plackett-Burman實驗各參數(shù)水平表Tab.3 Parameters of Plackett-Burman test

將表3各參數(shù)水平導(dǎo)入到Design-Expert 軟件中，設(shè)置4個中心點(diǎn)，進(jìn)行16組實驗，所得Plackett-Burman實驗設(shè)計及對應(yīng)的仿真結(jié)果如表4所示。

表4 Plackett-Burman實驗設(shè)計及仿真結(jié)果Tab.4 Design and results of Plackett-Burman test

對表4所示Plackett-Burman實驗結(jié)果進(jìn)行方差分析，得到各參數(shù)對堆積角的影響顯著性如表5所示。由表可以看出，JKR表面能、 染料-染料滾動摩擦系數(shù)的P<0.01，對顆粒堆積角的影響極其顯著；染料-不銹鋼靜摩擦系數(shù)的P<0.05，對顆粒堆積角的影響顯著；而其余參數(shù)P值均大于0.05，對顆粒堆積角影響極小。在最陡爬坡以及Box-Behnken實驗中只考慮這3個影響顯著的參數(shù)。

2.2 最優(yōu)區(qū)間判定

表5為各參數(shù)對堆積角影響顯著性的方差分析表。由表可知，3個顯著性因素對堆積角的影響效應(yīng)均是正效應(yīng)，因此，在爬坡實驗中設(shè)定它們的爬坡步長逐步增加，其余因素選擇表3所列高、低水平的平均值。設(shè)計實驗方案從Plackett-Burman實驗中心點(diǎn)開始，為了較快逼近最優(yōu)值，爬坡步長取較大值，此處取步長ΔG為0.2，剩余2個參數(shù)根據(jù)Plackett-Burman 實驗確定的回歸模型的效應(yīng)值計算得出。綜合考慮參數(shù)的取值區(qū)間，經(jīng)過近似優(yōu)化后，各爬坡步長可確定為ΔE=0.05、 ΔJ=0.001，實驗方案設(shè)計及仿真實驗結(jié)果分析如表6所示。

表5 Plackett-Burman實驗參數(shù)顯著性分析表Tab.5 Analysis of significance of parameters in Plackett-Burman test

表6 最陡爬坡實驗方案設(shè)計及仿真實驗結(jié)果分析Tab.6 Design and analysis of simulation resultsof steepest ascent test

根據(jù)最陡爬坡實驗，可以發(fā)現(xiàn)在2號水平處堆積角的相對誤差較小，由1號水平到3號水平其相對誤差由大變小再變大，由此可知，最優(yōu)值在1號和3號水平代表的參數(shù)值之間。響應(yīng)面擬合方程只在考察的緊接鄰域中才充分近似真實情形，只有逼近最大響應(yīng)區(qū)域后，才能建立有效的響應(yīng)面擬合方程。在第1次爬坡實驗后，1號和3號水平的相對誤差值較大，因此，有必要進(jìn)行二次爬坡實驗，以便進(jìn)一步縮小最優(yōu)解范圍。二次爬坡實驗的爬坡步長設(shè)置與第1次類似，實驗設(shè)計方案及結(jié)果如表7所示。由表7可以看出，2號水平至4號水平相對誤差由大變小再變大，故此，可以選取3號水平為中心點(diǎn)，分別選擇2、 4號水平為低、高水平進(jìn)行后續(xù)Box-Behnken實驗設(shè)計。

表7 二次最陡爬坡實驗方案設(shè)計及仿真實驗結(jié)果分析Tab.7 Design and analysis of simulation resultsof thesecond steepest ascent test

2.3 二階回歸方程的建立

選取顯著性影響因素3個水平進(jìn)行Box-Behnken實驗，設(shè)置5組中心實驗，實驗設(shè)計及結(jié)果如表8所示。顯著性參數(shù)與堆積角的二階回歸方程為

β=50.43+447.91E+121.42G+32 545.50J-155.00EG-85 200.00EJ-

20 850.00GJ-340.80E2-70.55G2-2 812 000J2。

(5)

表8 Box-Behnken實驗設(shè)計及仿真實驗結(jié)果Tab.8 Design and simulation results of Box-Behnken test

方差分析結(jié)果如表9所示。由表可知，交互項E×G與二次項E2的P>0.05，影響不顯著。

表9 Box-Behnken實驗?zāi)Ｐ头讲罘治鯰ab.9 ANOVA for model of Box-Behnken test

為進(jìn)一步優(yōu)化模型，除去表9中影響不顯著的因素，得到顯著性參數(shù)與堆積角的二階回歸方程為

β=-41.27+317.75E+106.81G+32 724.87J-85 200.00EJ-

20 850.00GJ-71.67G2-2 856 840J2。

(6)

參數(shù)優(yōu)化后模型的方差分析結(jié)果如表10所示。由表可知，模型的P<0.000 1，優(yōu)化模型極其顯著。參數(shù)E、G、J，交互項E×J、G×J以及二次項G2、J2的P值均進(jìn)一步降低，較未優(yōu)化前更加顯著，同時模型失擬項PL=0.128 1，擬合性較優(yōu)化前有提高；變異系數(shù)CV為1.46%，模型決定系數(shù)R2為0.981，校正決定系數(shù)為0.966，預(yù)測決定系數(shù)Rp為0.903，均在合理范圍之內(nèi)；優(yōu)化模型的實驗精密度PA由25.89提高到28.57。綜上，優(yōu)化模型更能真實反映實際狀況。

表10 Box-Behnken實驗優(yōu)化模型方差分析Tab.10 ANOVA for modified model of Box-Behnken test

2.4 回歸模型交互效應(yīng)分析

根據(jù)優(yōu)化模型方差分析結(jié)果，可知交互項E×J、G×J的P值均小于0.01，影響極其顯著，為探究各交互項中兩參數(shù)影響預(yù)測響應(yīng)值的交互效應(yīng)，分別選定G為0.4、E為0.1，繪制交互項E×J、G×J的響應(yīng)曲面，如圖6所示。

a)交互項E×J的響應(yīng)曲線b)交互項G×J的響應(yīng)曲線圖6 交互項E×J、 G×J的響應(yīng)曲面Fig.6 Interaction effect diagram of E×J、 G×J

由圖可以直觀得出，當(dāng)固定參數(shù)G的值不變后，在單位范圍內(nèi)，相對參數(shù)E，參數(shù)J的效應(yīng)曲線較陡，對堆積角影響較顯著，且參數(shù)J和參數(shù)E對堆積角的影響均為正效應(yīng)；當(dāng)固定參數(shù)E的值不變后，參數(shù)G和參數(shù)J對堆積角的影響均相對平緩，單位范圍內(nèi)參數(shù)J的效應(yīng)曲線相對參數(shù)G的曲線較陡峭，影響更顯著，兩參數(shù)堆積角的影響均為正效應(yīng)。同時，通過響應(yīng)曲面可以得出，各因子對堆積角影響的顯著性順序為：J>E>G，這與方差分析表顯示的結(jié)果相同。

3 最優(yōu)參數(shù)組合確定及實驗驗證

應(yīng)用Design-Expert軟件，以物理實驗實測所得堆積角度值為目標(biāo)值，對優(yōu)化后所得二階回歸方程求最優(yōu)解，所得參數(shù)值分別為染料-染料顆粒間滾動摩擦系數(shù)E=0.075 2、 染料-不銹鋼靜摩擦系數(shù)G=0.493以及JKR表面能J=0.002 01，將所得參數(shù)代入EDEM軟件中進(jìn)行3次模擬仿真實驗驗證，所得堆積角圖像與實測堆積角圖像的對比如圖7所示。

a)實驗實測結(jié)果b)仿真結(jié)果圖7 實驗實測結(jié)果與仿真結(jié)果的對比Fig.7 Comparison of test result and simulation result

測量仿真實驗所得堆積角分別為38.41°、 37.86°、 38.06°，標(biāo)準(zhǔn)差為0.278°，說明該組合參數(shù)仿真結(jié)果穩(wěn)定；對仿真結(jié)果和實驗實測結(jié)果進(jìn)行秩和檢驗分析，表明仿真結(jié)果與實驗實測結(jié)果無顯著性差異；仿真結(jié)果平均值為38.11°，與實驗實測結(jié)果平均值38.57°的相對誤差為1.19%，相對誤差較小，結(jié)果可靠。

4 結(jié)論

1)選擇粉體活性染料3BF Hong作為研究對象進(jìn)行離散元參數(shù)標(biāo)定實驗，通過Plackett-Burman實驗得出，JKR接觸模型下，影響粉體活性染料堆積角的3個顯著性參數(shù)為染料-染料顆粒間滾動摩擦系數(shù)、 染料-不銹鋼靜摩擦系數(shù)以及JKR表面能，Box-Behnken實驗得到的二次回歸模型經(jīng)過優(yōu)化后，其方差分析結(jié)果表明該模型擬合度良好，方程擬合精確度高，可以用來預(yù)測堆積角。

2)開展最陡爬坡實驗判定所考察因素的最大響應(yīng)區(qū)域，若考察因素在高低水平處誤差較大，為保證后續(xù)響應(yīng)面擬合方程的精度，可以進(jìn)行二次最陡爬坡實驗降低相對誤差值，以進(jìn)一步提高擬合方程的精確度。

3)對優(yōu)化后的二次回歸模型求最優(yōu)解，可得染料-染料顆粒間滾動摩擦系數(shù)、 染料-不銹鋼靜摩擦系數(shù)、 JKR表面能參數(shù)的數(shù)值分別為0.075 2、 0.493、 0.002 01，通過仿真驗證，所得堆積角度與實測值38.57°之間無顯著性差異，其相對誤差為1.19%，模擬值與實測值之間誤差較小，因此判定此實驗標(biāo)定的粉體活性染料顆粒離散元參數(shù)是可靠的，應(yīng)用離散單元法進(jìn)行微米級顆粒的離散元參數(shù)標(biāo)定是可行的。