屠熹遠(yuǎn) 楊偉東 王媛媛 高翔宇

摘要 在基于微滴噴射的砂型噴墨打印工藝中，砂粒鋪展過程對制件的打印精度具有重要影響，故對砂堆的堆積特性進(jìn)行基于離散元法的仿真研究。首先針對砂粒復(fù)雜形狀和不對稱性，采用了一種近似模擬8種典型砂粒形狀的多球體建模方法，并給出了相應(yīng)的力學(xué)模型；然后基于EDEM軟件進(jìn)行了不同形狀顆粒隨機(jī)混合生成，模擬了不同圓筒和底盤尺寸的休止角試驗(yàn)裝置及顆粒不同形狀對休止角的影響。仿真和實(shí)驗(yàn)結(jié)果表明休止角試驗(yàn)裝置的尺寸選擇和顆粒形狀分布對顆粒堆積特性具有重要影響；當(dāng)圓筒直徑為2.5 mm，底盤直徑為3.5 mm時休止角隨裝置尺寸變化最小；桿狀顆粒比片狀顆粒流動性好，片狀顆粒所占比例增加使砂堆休止角變大。

關(guān) 鍵 詞 噴墨打印；離散元；休止角；尺寸選擇；形狀分布

中圖分類號 TB44? ? ?文獻(xiàn)標(biāo)志碼 A

Accumulation characteristics of non-spherical sand particles based on discrete element method

TU Xiyuan， YANG Weidong， WANG Yuanyuan， GAO Xiangyu

（School of Mechanical Engineering， Hebei University of Technology， Tianjin 300401， China）

Abstract In 3D printing technology based on droplet jet， the spreading process of sand grains has an important influence on the printing accuracy of the workpiece. Firstly， for the complex shape and asymmetry of sand particles， a multi-sphere modeling method is used to approximate eight typical sand particle shapes， and the corresponding mechanical model is given. Then different kinds of particles were randomly mixed and generated based on EDEM software， and the effects of different cylinder and chassis sizes and particle shapes on the repose angle were simulated. The simulation and experimental results show that the size selection of the repose angle device and particle shape distribution have important influence on the particle accumulation characteristics. When the diameter of the cylinder is 2.5 mm and the diameter of the chassis is 3.5 mm， the repose angle of sand piles changes the least with the size of the device. It is also found that the fluidity of rod particles is better than that of flake particles ， and the increase of the proportion of flake particles in the whole ones makes the repose angle of sand pile much larger by the corresponding simulation.

Key words 3D printing; discrete element; repose angle; size selection; shape distribution

0 引言

噴墨打印技術(shù)（3DP）是增材制造技術(shù)中一項(xiàng)高性價比的快速成型工藝，具有非常廣闊的市場前景，針對該技術(shù)的研究工作發(fā)展迅速。其中，鋪粉系統(tǒng)影響粉層的致密度[1]，并最終使制件性能發(fā)生變化。本文以原砂顆粒作為研究對象，目前針對原砂顆粒模擬大多采用球形顆粒作為基本模型[2]，但其實(shí)際顆粒形狀并非球形，顆粒呈現(xiàn)出不規(guī)則形狀，使得模擬仿真的結(jié)果與實(shí)際實(shí)驗(yàn)的偏差較大。大量物理實(shí)驗(yàn)和模擬表明，非球形顆粒的引入對顆粒運(yùn)動和堆積產(chǎn)生重要影響[3-4]，同時顆粒形狀的變化會導(dǎo)致流動性的變化[5]，從而改變了顆粒的分布及平均運(yùn)動速度。非球形顆粒的空間取向分布直接影響砂層的致密度和孔隙率。因此，砂粒形狀的準(zhǔn)確描述和分析砂粒形狀對顆粒堆積特性的影響，對3DP工藝具有重要的指導(dǎo)意義。

近些年來，離散元法（DEM）被廣泛應(yīng)用于顆粒力學(xué)問題，但在使用過程中通常對顆粒形狀和顆粒間接觸模型等方面進(jìn)行簡化[6]。在非球形顆粒的離散模擬過程中，顆粒之間的接觸判斷更為復(fù)雜，接觸力學(xué)模型相對于球形顆粒也要發(fā)生變化。目前描述三維非球形顆粒形狀的方法主要有[3-4]：橢球、超二次、多面體、多球體組合和復(fù)合幾何元素組合。其中多球體組合方法采用球體作為基本單元，原則上能以足夠精度模擬任何復(fù)雜形狀，最大優(yōu)勢是可以直接應(yīng)用為球形DEM開發(fā)的接觸檢測評估程序，是非球形離散元仿真中最通用的方法。Alireza等[7]采用離散元方法模擬漏斗排出顆粒的過程，并用多球體方法模擬真實(shí)顆粒形狀，實(shí)驗(yàn)與仿真的一致性良好。Haeri等[8]在增材制造領(lǐng)域采用離散元方法仿真桿狀顆粒的鋪粉過程，用多個球體組成桿狀顆粒，發(fā)現(xiàn)顆?？v橫比越大或棍子的平移速度越大粉床的表面質(zhì)量越差。Kruggel-Emden等[9]提出了用多球體組合模擬真實(shí)非球體顆粒，并發(fā)現(xiàn)多球體的宏觀碰撞特性與圓球體仿真相比更接近真實(shí)的情況。張程林等[10]認(rèn)為砂土顆粒的形狀對砂土的性質(zhì)具有重要影響，提出根據(jù)顆粒加權(quán)圓度模擬其真實(shí)形狀。

縱觀國內(nèi)外研究發(fā)現(xiàn)，多球體方法適用于增材制造鋪粉過程的顆粒形狀模擬。本文在上述研究基礎(chǔ)上，用多球體的方法解決對實(shí)際顆粒形狀的描述，其中涉及到不同球體（來源于砂粒顯微鏡觀察的結(jié)果）的混合方式（采用等概率方法），研究分析了數(shù)值模擬中試驗(yàn)裝置尺寸及顆粒形狀分布對砂堆流動性的影響，對砂型打印機(jī)的工藝參數(shù)的選擇具有一定的指導(dǎo)意義。

1 力學(xué)模型

本文選用軟球模型作為顆粒系統(tǒng)基本模型，這種方法更適用于碰撞持續(xù)時間與粒子自由飛行的持續(xù)時間相當(dāng)?shù)南到y(tǒng)。顆粒間接觸力學(xué)模型如圖1所示。

Hertz-Mindlin接觸模型[11]認(rèn)為法向力[Fnij]由作用在顆粒i和顆粒[j]上的彈性力和阻尼力組成，分別由彈簧和法向阻尼器產(chǎn)生。對于二維顆粒，彈性力與法向重疊量大小成正比，阻尼力與顆粒相對速度成正比，[Fnij]可表示為

式中：[δ]是法向重疊量；[kn]是顆粒的法向彈性系數(shù)；[cn]是法向阻尼系數(shù)；[vij]是顆粒[i]相對于顆粒[j]的速度；[n]為顆粒[i]的中心到顆粒[j]中心的單位矢量。

切向力[Ftij]可表示為

式中，[kt]、[ct]、[α]、[vct]分別為切向彈性系數(shù)、阻尼系數(shù)、切向位移、接觸點(diǎn)的滑移速度。

仿真中的滾動摩擦[12]可由一作用在接觸面上的力矩表征：

式中：[μ]為滾動摩擦系數(shù)；[wi]是顆粒在接觸點(diǎn)處的單位角速度向量。

Johnson等人充分考慮了表面能對物體接觸的影響，考察了兩彈性球體間粘附效應(yīng)。運(yùn)用能量計(jì)算方法，得到了兩球體的接觸半徑與粘附能、外載荷、球體尺寸之間的關(guān)系，提出了著名的JKR模型[13]，如公式（4）（5）所示：

式中：[a、γ、E*、P]分別是兩球體接觸半徑，表面能密度和外載荷；[E*]和[R]分別為等效彈性模量和等效半徑，由公式（6）和（7）決定：

式中，[Ei，Ej，vi，vj，Ri，Rj]分別為兩球體的彈性模量、泊松比、半徑。JKR模型在DEM的仿真中的可行性已經(jīng)得到了廣泛驗(yàn)證。

本文的休止角數(shù)值模擬使用EDEM軟件，在原砂顆粒間的力學(xué)作用采用上述的Hertz-Mindlin with JKR模型。

2 參數(shù)設(shè)置與方法

2.1 仿真參數(shù)設(shè)置

砂堆休止角是反映原砂材料流動性的數(shù)據(jù)，休止角越大說明該顆粒材料流動性越差。為了研究噴墨打印干砂顆粒的堆積特性，進(jìn)行了球形砂堆和非球形砂堆的仿真模擬，并與原砂砂堆的物理實(shí)驗(yàn)進(jìn)行了對比。仿真參數(shù)的設(shè)定對仿真結(jié)果有重要影響，通過查閱文獻(xiàn)[14-15]、手冊和實(shí)驗(yàn)測定得到所有需要用到的相應(yīng)參數(shù)見表1。仿真網(wǎng)格劃分設(shè)置為[2rmin]，時間步長設(shè)置[16]滿足如下公式：

式中，[rmin]，[mmin]分別為仿真中所用到的最小圓球顆粒的半徑和質(zhì)量。

本文采用英國馬爾文公司型號為 Mastersizer 2000 的激光粒度分析儀，測試未進(jìn)行篩選的焙燒砂的粒度，得到該砂樣的粒徑分布曲線。其平均粒徑為0.173 mm，最大粒徑0.450 mm，粒徑分布如圖2所示，符合正態(tài)分布。

2.2 多球體顆粒建模

多球體方法的一個優(yōu)勢是顆粒間碰撞的建模，法向接觸力不會產(chǎn)生力矩，相對于非球體形狀，球體間的接觸點(diǎn)的判斷較為簡單。但形狀精度與組成單元的數(shù)量息息相關(guān)[17]，隨著球體數(shù)量的增加，模擬得到的非球形顆粒的形狀精度越高。但由于鋪砂過程中砂粒的數(shù)量較多，為了減小計(jì)算工作量和提高效率，模擬計(jì)算中近似逼近砂粒形狀。

為了更好地模擬出砂粒的形狀，本文利用光學(xué)顯微鏡對2 000顆不同形狀的砂粒形貌進(jìn)行了觀測，根據(jù)顆粒形狀的分布選取了8種典型砂粒形狀，如圖3所示。使用Proe軟件按照顆粒形狀，用多個球體組合形成顆粒，近似模擬得到相應(yīng)形狀的外形，如圖4所示。

2.3 多球體顆粒生成混合方法

在EDEM軟件中直接生成仿真所用顆粒的粒子工廠不能實(shí)現(xiàn)多種形狀粒子在同一容器內(nèi)的隨機(jī)混合，通過該軟件可進(jìn)行二次開發(fā)的接口模塊，在立方體容器中設(shè)計(jì)粒子生成算法，實(shí)現(xiàn)過程如圖5所示。首先把立方體內(nèi)部區(qū)域劃分為20 mm×20 mm×30 mm的網(wǎng)絡(luò)單元，每個網(wǎng)格單元為邊長0.3 mm的立方體。然后選取底層左下角落處網(wǎng)格為第1個顆粒的生成區(qū)域，顆粒生成過程完成后循環(huán)到X軸方向的下一個網(wǎng)格生成，直到X軸方向的這一列網(wǎng)格全部生成完成顆粒。之后Y軸坐標(biāo)增加1個邊長長度，生成下一列粒子，直到這一層的所有網(wǎng)格內(nèi)都生成了一類粒子。之后Z軸坐標(biāo)再增加1個邊長長度，生成下一層的粒子，以此類推。直到顆粒的總數(shù)量達(dá)到要求時，顆粒生成算法結(jié)束。

針對非球形顆粒，程序每循環(huán)1次生成1個顆粒，每次循環(huán)生成的顆粒形狀是上述8種類型多球體顆粒的隨機(jī)一種，采用等概率分配方式，概率均為1/8。每個網(wǎng)格內(nèi)的粒子三軸坐標(biāo)有0.1 mm范圍的隨機(jī)浮動。旋轉(zhuǎn)矩陣與歐拉角實(shí)現(xiàn)三維空間坐標(biāo)變換以保證生成顆粒的隨機(jī)擺放角度。使用Box Muller方法令顆粒粒徑分布服從正態(tài)分布。生成的顆粒受重力作用沿著一漏斗流入圓筒中，因EDEM軟件的顆粒的接觸檢測算法在圓筒中形成互不交疊、無不合理孔隙的砂堆，之后去掉立方體容器和過渡漏斗。至此，圓筒仿真試驗(yàn)的多球體砂粒的混合準(zhǔn)備工作完成。

3 結(jié)果與分析

3.1 圓筒試驗(yàn)尺寸選擇

在休止角的測量中一般圓筒的直徑應(yīng)大于顆粒最大直徑尺寸的4～5倍[18]，為了保證足夠的砂粒形成穩(wěn)定的休止角，圓筒高度設(shè)置均大于直徑的3倍。

為了避免尺寸效應(yīng)，取不同尺寸的圓筒和底盤。在仿真過程中，首先采用球形顆粒進(jìn)行仿真，圓筒堆積著所需顆粒，以5 mm/s的速度提升圓筒，最終形成穩(wěn)定的圓錐形砂粒堆。休止角形成過程如圖6所示，球形顆粒砂堆如圖7所示。

本文的休止角測量方法[19]使用Matlab軟件對休止角圖像進(jìn)行圖像處理和曲線擬合，處理過程如圖8所示。首先對剖面圖像進(jìn)行二值化處理，然后使用Zerocross算子對二值圖像進(jìn)行邊緣檢測，并得到所有砂層表面輪廓上點(diǎn)的像素坐標(biāo)，運(yùn)用最小二乘法對這些像素豎直方向的坐標(biāo)進(jìn)行線性擬合，最后得到擬合直線和對應(yīng)的休止角度。不同圓筒直徑和底盤直徑條件下球形顆粒休止角的模擬仿真結(jié)果如表2所示。

由表2可以看出，當(dāng)圓筒的直徑從2 mm增加到4 mm，底盤直徑相應(yīng)從3.5 mm增加到6 mm，此時休止角受裝置尺寸影響產(chǎn)生的變化最小，故選圓筒直徑為2.5 mm，底盤直徑為3.5 mm作為測量砂粒休止角模擬的裝置，采用該大小的裝置進(jìn)行數(shù)值模擬測得的休止角為31.72°。

將圖4中多球體的砂粒通過Proe軟件導(dǎo)入EDEM軟件中，按照相同概率生成8中不同形狀的多球體顆粒，采用圓筒直徑為2.5 mm，底盤直徑為3.5 mm作為測量砂粒休止角模擬的裝置。針對非球形砂粒的仿真結(jié)果如圖9所示，測量得到的多球體顆粒休止角為34.4°。

將數(shù)值模擬得到的裝置尺寸按比例擴(kuò)大搭建試驗(yàn)裝置如圖10所示，底盤直徑為140 mm，圓筒高度為320 mm，直徑為100 mm。使用直流電機(jī)帶動圓筒勻速提升得到砂堆，用同樣的休止角測量方法測得的結(jié)果為33.6°。實(shí)驗(yàn)測量的休止角與球形顆粒仿真相比，誤差為5.6%。而與多球體顆粒休止角相比，誤差為2.4%，其誤差小于球形顆粒時的仿真模擬。

3.2 顆粒形狀對流動性影響

為了進(jìn)一步研究顆粒形狀對顆粒流動性的影響，將圖4中多球體顆粒的1號、2號、4號、5號、7號劃分為桿狀顆粒，3號、6號、8號劃分為片狀顆粒。選取圓球體顆粒、3號片狀顆粒及4號桿狀顆粒用于數(shù)值模擬。

不同形狀顆粒混合后休止角的仿真結(jié)果如表3所示。表3中仿真結(jié)果顯示，不同比例的顆粒形狀組成砂堆對休止角產(chǎn)生重要影響。其中球形顆粒的流動性最好，兩種多球體顆粒相比桿狀顆粒的流動性比片狀顆粒好。當(dāng)采用其中兩種不同形狀的顆?；旌蠒r，圓筒試驗(yàn)形成的砂堆休止角在這兩種單一形狀顆粒形成的休止角之間變化。在顆粒形狀表征精度足夠的前提下，多球體方法會讓數(shù)值模擬的休止角結(jié)果更加趨近于實(shí)驗(yàn)數(shù)據(jù)，流動性更接近于真實(shí)形成的砂堆。

4 結(jié)語

本文基于離散元方法，建立了砂粒間相互作用的力學(xué)模型。通過光學(xué)顯微鏡觀察砂粒，對8種典型砂粒形狀構(gòu)建了多球體顆粒物理模型，得到了如下結(jié)論。

1）基于EDEM軟件采用了一種針對不同形狀顆粒生成混合方法，完成了多球體顆粒組合的混合過程。

2）為避免尺寸效應(yīng)，在休止角模擬仿真實(shí)驗(yàn)中修改圓筒和底盤的尺寸，最終采用圓筒直徑為2.5 mm，底盤直徑為3.5 mm作為測量砂粒休止角模擬的裝置。并按測量休止角模擬裝置比例完成了休止角試驗(yàn)裝置的搭建和實(shí)驗(yàn)。與砂粒休止角實(shí)驗(yàn)數(shù)據(jù)相對比，多球體顆粒仿真模擬的休止角的誤差小于球形顆粒仿真模擬的休止角。砂粒的多球體模型流動性與其球形模型相比更接近實(shí)驗(yàn)數(shù)據(jù)。

3）研究了顆粒形狀變化對休止角產(chǎn)生的影響，發(fā)現(xiàn)球形顆粒的流動性最好。其中2種多球體顆粒相比，桿狀顆粒的流動性比片狀顆粒好。片狀顆粒所占比例增加使砂堆休止角變大。

參考文獻(xiàn)：

[1]? ? ZIEGELMEIER S，CHRISTOU P，W?LLECKE F，et al. An experimental study into the effects of bulk and flow behaviour of laser sintering polymer powders on resulting part properties[J]. Journal of Materials Processing Technology，2015，215：239-250.

[2]? ? 張帥，單忠德，張杰. 基于離散元方法的型砂流動性仿真研究[J]. 鑄造技術(shù)，2016，37（2）：288-291.

[3]? ? LU G，THIRD J R，M?LLER C R. Discrete element models for non-spherical particle systems：from theoretical developments to applications[J]. Chemical Engineering Science，2015，127：425-465.

[4]? ? COETZEE C J. Review：Calibration of the discrete element method[J]. Powder Technology，2017，310：104-142.

[5]? ? BERRETTA S，GHITA O，EVANS K E. Morphology of polymeric powders in laser sintering （LS）：from polyamide to new PEEK powders[J]. European Polymer Journal，2014，59：218-229.

[6]? ? 譚援強(qiáng)，鄭軍輝，張浩，等. 基于離散元法的錐形筒倉中顆粒流體的數(shù)學(xué)模擬[J]. 過程工程學(xué)報，2015，15（6）：916-922.

[7]? ? KHAZENI A，MANSOURPOUR Z. Influence of non-spherical shape approximation on DEM simulation accuracy by multi-sphere method[J]. Powder Technology，2018，332：265-278.

[8]? ? HAERI S，WANG Y，GHITA O，et al. Discrete element simulation and experimental study of powder spreading process in additive manufacturing[J]. Powder Technology，2017，306：45-54.

[9]? ? KRUGGEL-EMDEN H，RICKELT S，WIRTZ S，et al. A study on the validity of the multi-sphere Discrete Element Method[J]. Powder Technology，2008，188（2）：153-165.

[10]? 張程林，周小文. 砂土顆粒三維形狀模擬離散元算法研究[J]. 巖土工程學(xué)報，2015，37（增刊）：115-119.

[11]? 王國強(qiáng)，郝萬軍，王繼新. 離散單元法及其在EDEM上的實(shí)踐[M]. 西安：西北工業(yè)大學(xué)出版社，2010.

[12]? HAN Q Q，GU H，SETCHI R. Discrete element simulation of powder layer thickness in laser additive manufacturing[J]. Powder Technology，2019，352：91-102.

[13]? JOHNSON K，KENDALL K，ROBERTS A D. Surface energy and the contact of elastic solids[J]. Proceedings of the Royal Society of London A Mathematical and Physical Sciences，1971，324：301-313.

[14]? 姚志雄. 三維離散元法解析砂型泥石流細(xì)觀運(yùn)動特征研究[J]. 地下空間與工程學(xué)報，2017，13（5）：1414-1423.

[15]? 頓國強(qiáng)，楊永振，陳海濤，等. 基于EDEM的控位分層施肥開溝器作業(yè)過程的仿真與試驗(yàn)[J]. 湖南農(nóng)業(yè)大學(xué)學(xué)報（自然科學(xué)版），2018，44（1）：95-100.

[16]? MEIER C，WEISSBACH R，WEINBERG J，et al. Modeling and characterization of cohesion in fine metal powders with a focus on additive manufacturing process simulations[J]. Powder Technology，2019，343：855-866.

[17]? FERELLEC J F，MCDOWELL G R. A method to model realistic particle shape and inertia in DEM[J]. Granular Matter，2010，12（5）：459-467.

[18]? 賈富國，韓燕龍，劉揚(yáng)，等. 稻谷顆粒物料堆積角模擬預(yù)測方法[J]. 農(nóng)業(yè)工程學(xué)報，2014，30（11）：254-260.

[19]? 韓燕龍，賈富國，唐玉榮，等. 顆粒滾動摩擦系數(shù)對堆積特性的影響[J]. 物理學(xué)報，2014，63（17）：173-179.