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

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

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

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

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

1.1 JKR接觸模型

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

(1)

(2)

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

Δγ=γi+γj-γij

(3)

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

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

1.2 堆積角的測定

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

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

(4)

式中：β為測得的堆積角；k為擬合直線斜率。經(jīng)過5次實(shí)驗(yàn)，求得堆積角的平均值為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%，能夠代表整個(gè)粒子系統(tǒng)參與計(jì)算。將45～255 μm的粒徑劃分為5個(gè)區(qū)間，選取代表粒徑，計(jì)算所占比率，據(jù)此在軟件中創(chuàng)建顆粒，同時(shí)結(jié)合仿真預(yù)實(shí)驗(yàn)確定各顆粒的生成時(shí)間及總數(shù)，如表1所示。

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

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

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

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

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

2 參數(shù)標(biāo)定實(shí)驗(yàn)設(shè)計(jì)及分析

參數(shù)標(biāo)定實(shí)驗(yàn)主要包括3個(gè)部分： 1)Plackett-Burman實(shí)驗(yàn)，用以選擇待標(biāo)定因子中對堆積角影響最顯著的因子； 2)最陡爬坡實(shí)驗(yàn)，用以確定顯著性因子的最優(yōu)區(qū)間； 3)Box-Behnken實(shí)驗(yàn)，用以建立顯著性影響因子與響應(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ù)實(shí)驗(yàn)確定，詳見表2。

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

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

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

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

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

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

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

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

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

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

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

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

2.3 二階回歸方程的建立

選取顯著性影響因素3個(gè)水平進(jìn)行Box-Behnken實(shí)驗(yàn)，設(shè)置5組中心實(shí)驗(yàn)，實(shí)驗(yàn)設(shè)計(jì)及結(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í)驗(yàn)設(shè)計(jì)及仿真實(shí)驗(yàn)結(jié)果Tab.8 Design and simulation results of Box-Behnken test

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

表9 Box-Behnken實(shí)驗(yàn)?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，交互項(xiàng)E×J、G×J以及二次項(xiàng)G2、J2的P值均進(jìn)一步降低，較未優(yōu)化前更加顯著，同時(shí)模型失擬項(xiàng)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)化模型的實(shí)驗(yàn)精密度PA由25.89提高到28.57。綜上，優(yōu)化模型更能真實(shí)反映實(shí)際狀況。

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

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

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

a)交互項(xiàng)E×J的響應(yīng)曲線b)交互項(xiàng)G×J的響應(yīng)曲線圖6 交互項(xiàng)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)。同時(shí)，通過響應(yīng)曲面可以得出，各因子對堆積角影響的顯著性順序?yàn)椋篔>E>G，這與方差分析表顯示的結(jié)果相同。

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

應(yīng)用Design-Expert軟件，以物理實(shí)驗(yàn)實(shí)測所得堆積角度值為目標(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次模擬仿真實(shí)驗(yàn)驗(yàn)證，所得堆積角圖像與實(shí)測堆積角圖像的對比如圖7所示。

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

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

4 結(jié)論

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

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

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