噴動床粉體混合機混合性能的數(shù)值模擬

王一澤，熊桂龍，蘇文康，劉振峰

( 1.南昌大學a.資源環(huán)境與化工學院；b.鄱陽湖環(huán)境與資源利用教育部重點實驗室，江西 南昌 330031；2.宜春萬申制藥機械有限公司，江西 宜春 336000)

離散物料的混合設備主要有兩大類：機械式混合設備與氣力式混合設備[1]。與傳統(tǒng)的機械式混合設備相比，氣力式粉體混合機具有混合體積大、均勻性好、能耗低、密閉無泄漏、清潔環(huán)保等優(yōu)點，在能源、化工、制藥，食品等行業(yè)領域應用廣泛[2-3]。氣力式混合設備主要有噴動床、流化床混合設備，其工作原理是通過壓縮空氣使離散物料流態(tài)化并向上膨脹，壓縮空氣膨脹到一定程度后物料由于重力作用向下沉降，離散物料在氣流作用下產(chǎn)生上下翻滾和激烈的沸騰攪拌，在相對較短的時間內達到混合均勻[4-5]。因此，對噴動床粉體混合機的粉體混合過程及運動行為規(guī)律進行分析具有重要意義。

對于噴動床粉體混合的研究常采用實驗方法與數(shù)值模擬的手段，但由于受到實驗條件與設備的限制，通過實驗難以獲得顆粒微觀尺度上的全面運動信息，例如顆粒速度、顆粒混合質量、基于顆粒尺度的碰撞信息等。隨著計算機技術快速發(fā)展，數(shù)值模擬在研究氣固兩相流動問題的優(yōu)勢日益凸顯，可以獲得更為全面的顆粒運動信息[6-7]。李斌等[8]采用CFD-DEM耦合的數(shù)值模擬方法對循環(huán)流化床內顆粒的混合行為進行了研究，認為增大流化風速有助于加速顆粒混合。張俊強等[9-10]采用數(shù)值模擬的方法分別研究了單孔與雙孔射流流化床顆粒的混合效果，結果表明，顆粒的軸向混合速率大于徑向混合速率。陳程等[11]基于CFD-DEM的方法對旋轉流化床粉體混合機進行了研究，結果表明，進氣管傾斜一定角度布置，可全面提高混合強度。Sharma等[12]對鼓泡流化床內不同顆粒的混合與分離過程進行了數(shù)值模擬。Olaofe等[13]采用數(shù)值模擬的方法研究了流化床內氣固兩相運動特性，分析了顆粒在氣流作用下的混合過程。Ren等[14]采用數(shù)值模擬的方法對噴動床中不同密度顆粒的混合行為進行了研究，研究結果表明，隨著顆粒密度差的增大，粉體混合均勻性指數(shù)增大。已有的研究大多以實驗室小規(guī)模裝置為研究對象，所用的顆粒大多為毫米級，然而，在實際工業(yè)過程中設備的尺寸相對較大、顆粒的粒度較小，顆粒粒徑通常在微米或亞微米的數(shù)量級，顆粒數(shù)量較大，對實際工程尺度的噴動床進行數(shù)值模擬研究是噴動床粉體混合研究的熱點與難點。

本文采用數(shù)值模擬的方法對工業(yè)級噴動粉體混合機中直徑為165 μm的同種淀粉顆?；旌线^程進行了研究，根據(jù)相似原理與流動相似理論，對顆粒直徑進行了工程縮放，以減少顆粒數(shù)量，降低計算資源與計算量。建立了噴動床中粉體混合的模型，分析了噴動床中粉體微觀混合過程，噴動床中氣體與顆粒沿床寬和床高方向的速度分布，顆粒循環(huán)運動周期，評價粉體混合均勻性的混合指數(shù)，為工業(yè)級噴動床粉體混合機的研發(fā)與應用提供理論基礎與參考。

1 數(shù)值模擬

1.1 模型建立

噴動床粉體混合機的幾何模型和流體網(wǎng)格模型如圖1(a)～(c)所示，采用底部軸向進氣，頂部為出氣口，坐標原點設置在底部的圓心處，初始粉體堆積狀態(tài)如圖1(a)所示，床層表面顆粒1與床層底部顆粒2的比例為1:10；噴動床幾何尺寸參數(shù)如圖1(b)與表1所示；噴動床中流體運動區(qū)域的網(wǎng)格劃分采用結構化網(wǎng)格，流體網(wǎng)格模型如圖1(c)所示，經(jīng)網(wǎng)格無關性驗證后，最終Fluent網(wǎng)格劃分的節(jié)點數(shù)為45 903，網(wǎng)格數(shù)37 660。

(a) 噴動床結構圖 (b) 幾何尺寸

表1 噴動床幾何尺寸參數(shù)Tab.1 Geometry parameters of spouted bed

1.2 數(shù)學模型

1.2.1 流體相控制方程

對于不可壓縮流體，利用N-S方程來描述氣相的流動，連續(xù)性方程和動量方程分別為：

(1)

(2)

式中：ρf為流體密度，ρf=1.205 kg·m-3；αf為空隙率；p為壓力；uf為流體速度；τ為應力張量；Fpf為顆粒相與流體相的動量交換。

1.2.2 顆粒相控制方程

暮色四合，潼關泣血，天策府的女將軍曹雪陽綽長槍，領著將士們百戰(zhàn)斷頭，蹚在血海里，奮力殺賊，饒是如此，也難挽哥舒翰力戰(zhàn)敗北，歸降安祿山。大戰(zhàn)之余，朗朗晴天風云變幻，電閃雷鳴中，暴雨如注，蕩滌谷中血肉，將數(shù)十萬腥臭的尸骨沖進黃河。大唐經(jīng)此一役，由盛轉衰，奄忽百年遂亡，這是后話不提。

采用離散單元法對顆粒相進行模擬計算，該方法可描述系統(tǒng)中的每個單個粒子i通過接觸力Fc,ij與其他粒子或壁面j的相互作用。接觸力Fc,ij可以被分解為法向分力Fn,ij和切向分力Ft,ij。顆粒運動過程中所受的力主要包括：重力、顆粒之間、顆粒與壁面間的接觸力、摩擦力、曳力、浮力。顆粒的運動包括平動和轉動，由牛頓第二定律確定[15]：

(3)

(4)

式中：mp,i、Ip,i分別為顆粒i的質量、轉動慣量；vp,i、ωp,i分別為顆粒的平動速度、轉動速度；Fg,i和Ff,i分別表示顆粒重力和氣相對顆粒的作用力；Tt,ij和Tr,ij分別是由切向力和滾動摩擦產(chǎn)生的扭矩。顆粒間接觸采用Hertz-Mindlin模型。

1.2.3 曳力模型

氣體對顆粒的曳力Fd,i采用Gidaspow曳力模型，在空隙率αf>0.8時，顆粒的曳力計算方程為：

(5)

在空隙率αf≤0.8時，顆粒的曳力計算方程為：

(6)

Cd為單個顆粒的曳力系數(shù)，其取值與顆粒雷諾數(shù)Rep有關，Rep≤1 000時，Cd計算方程為：

(7)

當Rep=1 000時，

Cd=0.4

(8)

粉體物料混合后，其混合的均勻性可采用目前國內外使用較多的Lacey指數(shù)[11]來表征與評價，當粉體處于完全分離的狀態(tài)時，Lacey指數(shù)為0，粉體在理想混合均勻狀態(tài)時，Lacey指數(shù)為1，在實際混合過程中，Lacey指數(shù)在0～1之間變化，因此，可以利用Lacey指數(shù)的變化來分析粉體在整個噴動過程中的混合情況。表征粉體混合均勻性的Lacey指數(shù)為：

(9)

(10)

(11)

(12)

式中：P為一種顆粒所占比例；(1-P)為另一種顆粒所占比例；n為每一個樣本內平均顆粒數(shù)量；NS為樣本總數(shù)；ai為兩類顆粒中任意一類在樣本i中所占的比例；a為相應的顆粒在噴動床中所占比例；ki為樣本i的權重，表示為：

(13)

(14)

式中:Ni為樣本i內的顆粒數(shù)；Nt為顆粒數(shù)總和。

1.3 模型參數(shù)

顆粒系統(tǒng)中有兩種顆粒，初始堆積狀態(tài)下床層表面淀粉顆粒1的質量5 kg，床層底部淀粉顆粒2的質量為50 kg，兩種顆?？倲?shù)量為146億個，顆粒直徑為165 μm，顆粒密度為1 600 kg·m-3。氣流進口邊界條件設為速度進口，并采用UDF函數(shù)控制入口處氣流速度，前5 s為勻加速進氣狀態(tài)，初速度為0 m·s-1，加速度為1 m·s-2，在第5秒時入口處進氣速度達到并穩(wěn)定在5 m·s-1，直到第30秒時停止進氣；氣流出口邊界條件設置為一個標準大氣壓。為了減小計算量，縮短計算時間，基于相似原理與流動相似理論進行工程縮放，在工程縮放過程中，保證顆粒總體積、顆粒最小流化速度Umf、顆粒雷諾數(shù)Rep、阿基米德數(shù)Ar為常數(shù)不變[16-18]，可將計算的粒子總數(shù)Np由146億個減少到128 227個。

為了保證縮放前后顆粒床層體積的相似性，縮放后的顆粒直徑可由下式得到：

(15)

式中:mbed、N、ρp、dp分別為顆粒床層總質量、顆粒數(shù)目、顆粒密度和顆粒直徑。下角標1和2分別為初始顆粒系統(tǒng)與縮放后的顆粒系統(tǒng)。

縮放因子kc可由下式得到：

(16)

對于縮放前后的系統(tǒng)，氣體密度ρg和重力加速度g保持不變。

在數(shù)值模擬計算中，DEM時間步長為10 μs，CFD時間步長為500 μs，顆粒與顆粒間發(fā)生碰撞時恢復系數(shù)、靜摩擦系數(shù)、滾動摩擦系數(shù)分別為0.75、0.5、0.06，顆粒與壁面間發(fā)生碰撞時恢復系數(shù)、靜摩擦系數(shù)、滾動摩擦系數(shù)分別為0.8、0.3、0.05?？s放前后的顆粒與氣體的相關參數(shù)如表2所示。計算中所需的其他參數(shù)設置如表3所示。

表2 縮放前后顆粒和氣體的參數(shù)Tab.2 Parameters of particles and gases before and after scaling

表3 模擬參數(shù)設置Tab.3 Simulation parameter setting

采用EDEM2017與Fluent17.0耦合的方法，通過UDF(用戶自定義函數(shù))實現(xiàn)耦合，模擬計算粒子運動與氣相流場，離散相可以與流體相交換動量、質量和能量。

2 結果與討論

2.1 粉體的混合過程

圖2(a)～(f)是在不同時刻的噴動床混合機中粉體的混合過程，模擬時長共31 s。初始狀態(tài)下，床層底部淀粉顆粒2與床層表面淀粉顆粒1的體積比為10:1。

(a) t=0.32 s (b) t=5 s (c) t=9 s (d) t=13 s (e) t=19 s (f) t=31 s

由圖2(a)～(f)的混合過程可看出，根據(jù)噴動床中粉體的運動特征，噴動床可劃分為噴射區(qū)、噴泉區(qū)、環(huán)隙區(qū)。從圖2(a)可看出，在第0.32秒時刻，粉體顆粒依然處于初始堆積狀態(tài)，兩種粉體顆粒之間有明顯的分層；隨著床層底部進口氣流速度的增大，噴動床底部氣體壓強增大，床層開始向上膨脹，如圖2(b)所示；當進口氣體流速在第5秒時刻達到5 m·s-1并穩(wěn)定之后，噴射區(qū)逐漸形成，噴動氣體開始穿透粉體床層，噴動床層底部的粉體開始在氣流曳力作用下向上運動，粉體在向上運動的過程中會與周圍粉體發(fā)生相互碰撞；從圖2(c)可看出，當向上運動粉體在到達噴泉區(qū)時，在氣流和其他粉體的碰撞作用下，噴泉區(qū)的粉體開始向四周擴散并向下運動；從圖2(d)可知，粉體在噴泉區(qū)會發(fā)生強烈的混合，環(huán)隙區(qū)的粉體顆粒運動較緩慢，混合程度較弱，由圖2(e)可知，當向下運動的粉體沉降到床層表面時，由于粉體所受重力大于流體曳力，粉體沿壁面向下運動，之后粉體從環(huán)隙區(qū)進入噴射區(qū)，在氣體曳力作用下又重新向上運動，形成一個完整的循環(huán)運動過程。在循環(huán)運動的過程中最終達到粉體混合的動態(tài)平衡，如圖2(f)所示。

2.2 粉體與氣流速度分布

圖3為噴動床達到穩(wěn)定噴動時，在YOZ平面床層高度分別為0.2、0.4、0.6、0.8、1.0 m的粉體軸向速度和徑向速度沿床層寬度的分布，軸向速度正方向和Z軸正方向一致，徑向速度正方向和Y軸正方向一致。

床寬/m(a) 粉體軸向速度

由圖3(a)粉體軸向速度分布可知，噴射區(qū)不同高度處粉體的軸向速度為正值，該區(qū)域粉體向上運動，并且粉體軸向速度隨著床層高度的增加先增加后減小，因為在噴射區(qū)床層高度為0.2 m時，粉體受到氣體曳力作用加速運動，此時粉體軸向速度小于流體軸向速度，曳力為動力；當床層高度為1.0 m時，粉體軸向速度大于氣體軸向速度，此時曳力為阻力，粉體減速運動。在環(huán)隙區(qū)的粉體軸向速度為負值，粉體向下運動，壁面附近粉體的流動性差，其軸向速度越小。由圖3(b)粉體徑向速度分布可知，床層高度為0.2～0.4 m時，噴射區(qū)左側粉體徑向速度為正值，粉體向右運動進入噴射區(qū)；右側粉體徑向速度為負值，粉體向左運動。床層高度為0.6～1.0 m時，噴射區(qū)左側粉體徑向速度為負值，粉體向左運動，離開噴泉區(qū)后向四周擴散，且粉體徑向速度隨高度的降低而增大；噴射區(qū)右側粉體徑向速度為正值，粉體向右運動。由此可知，在噴射區(qū)粉體受到氣體的卷吸作用不斷向中間聚集，之后在氣流作用下向上運動，在噴泉區(qū)粉體不斷混合，之后向四周擴散，在重力作用下向下運動進入環(huán)隙區(qū)，與環(huán)隙區(qū)的粉體重新混合后向噴射區(qū)運動，如此循環(huán)反復，粉體混合程度不斷加深。

圖4為噴動床形成穩(wěn)定噴動后，在YOZ平面床層高度分別為0.2、0.4、0.6、0.8、1.0 m的氣流運動的軸向速度和徑向速度沿床層寬度的分布。

床寬/m(a) 氣流軸向速度

由圖4可知，在軸向上，噴射區(qū)氣流向上運動，環(huán)隙區(qū)的氣流向下流動；隨著床層高度的增加，在噴射區(qū)氣流軸向速度先增大后減小。沿徑向床層底部氣流向噴射區(qū)流動而噴泉區(qū)氣體向四周擴散。由圖3和圖4對比可知，粉體與氣流的分布規(guī)律相似，這也證明氣流曳力作用下的粉體運動是影響噴動床粉體混合的主要原因，此外，氣流軸向速度大于粉體軸向速度，流體動能轉化為粉體動能和流體穿透床層的壓力能。

2.3 粉體的循環(huán)運動周期

粉體在噴動床中依次經(jīng)過噴射區(qū)、噴泉區(qū)和環(huán)隙區(qū)，做周期性的循環(huán)運動，其循環(huán)運動周期可定義為顆粒完成從床層表面開始向上運動，經(jīng)過噴泉區(qū)，然后向下經(jīng)過環(huán)隙區(qū)，然后進入噴射區(qū)并回到床層表面所需要的時間。為了分析粉體在每個循環(huán)周期內的運動過程，分別提取了床層表面顆粒1與床層底部顆粒2的空間位置隨時間的變化，如圖5所示。

圖5 顆粒的三維空間坐標隨時間的變化Fig.5 Changes of particle's three-dimensional space coordinates over time

由圖5可知，床層底部顆粒2在第5.2秒時從床層表面開始向上運動，其在Z向的坐標值快速增大，運動至噴泉區(qū)，在第6秒左右的時刻，回落床層表面，進入環(huán)隙區(qū)，在第14.6秒左右的時刻，進入噴射區(qū)重新上升至床層表面，經(jīng)歷了一個完整的循環(huán)運動，其循環(huán)運動的周期約為9.4 s。其中，粉體在環(huán)隙區(qū)運動的時間為8.6 s，占整個循環(huán)運動周期的90%以上。顆粒1在第12.4秒時，完成第1次循環(huán)運動，顆粒2在第14.6秒時，完成第1次循環(huán)運動，床層表面和床層底部的粉體循環(huán)周期相差約2 s左右。顆粒1與顆粒2在第5.2秒時，進入噴泉區(qū)向四周擴散，其對應的X向與Y向坐標值突然增大；在從環(huán)隙區(qū)運動至噴射區(qū)的過程中，其對應的X與Y向坐標值會突然減小。在整個模擬計算時長內，粉體顆粒大約經(jīng)歷了3個循環(huán)運動周期：第1個運動周期在第14秒左右的時刻結束；第2個運動周期在第23秒左右的時刻結束；第3個運動周期在第30秒左右的時刻結束。由上述分析可知，顆粒在噴動床中循環(huán)運動周期約為9.4 s，在循環(huán)運動的過程中粉體不斷混合并最終達到動態(tài)平衡。

2.4 粉體混合的均勻性

在粉體顆?；旌喜僮鬟^程完成后，可利用大小一致的立方體虛擬網(wǎng)格對混合后的粉體區(qū)域進行劃分。網(wǎng)格大小對Lacey指數(shù)有較大影響，因此，需選擇合適的網(wǎng)格大小。如圖6(a)所示，分別采用邊長為32、40、50 mm的立方體虛擬網(wǎng)格進行劃分，并對相應的混合均勻性指數(shù)進行計算，提取混合均勻性指數(shù)隨時間的變化規(guī)律。

t/s(a) 樣本網(wǎng)格大小對混合指數(shù)的影響

由圖6(a)可知，在0～5 s時，選用不同大小的取樣網(wǎng)格，其對應的混合均勻性指數(shù)差別明顯，這是因為粉體在初始堆積狀態(tài)下，在兩種粉體顆粒的交界面取樣時，會認為兩種粉體有一定的混合度，因此，在初始階段，混合均勻性指數(shù)越小，越能反映粉體物料混合的真實情況，因此，可以認為，在0～5 s內，選用邊長為50 mm的立方體取樣更為合理；在5～15 s時，不同時刻混合指數(shù)的取值較為接近，在15 s之后，不同時刻混合指數(shù)的取值極為接近，綜合考慮認為邊長50 mm的立方體網(wǎng)格更為合理。由圖6(b)可知在0～5 s內，混合均勻性指數(shù)較小且會發(fā)生波動，這是因為在起始階段，床層膨脹，整體粉體床層整體向上運動，隨著兩種粉體分界面向上提升，在虛擬取樣網(wǎng)格內會認為粉體進行了一定程度的混合。在5～15 s內，由于噴射區(qū)的形成，兩種粉體在噴泉區(qū)發(fā)生劇烈混合，此階段的混合均勻性指數(shù)快速增大。在15 s之后，隨著混合過程的進行，兩種粉體的混合已經(jīng)到了較高的水平，混合指數(shù)隨時間的變化開始變得緩慢。當混合均勻性指數(shù)達到0.95時，可以認為兩種粉體的混合已經(jīng)基本均勻，所需要的時間約為20 s，之后隨著混合過程的進行，混合指數(shù)略微增加，最終可達到0.99以上，認為兩種粉體混合均勻。

3 結論

1) 噴動床中粉體在噴射區(qū)、噴泉區(qū)、環(huán)隙區(qū)循環(huán)運動，最終達到粉體混合的動態(tài)平衡，粉體混合主要發(fā)生在噴泉區(qū)。

2) 粉體運動的軸向速度大于徑向速度，氣流運動的軸向速度大于粉體運動的軸向速度，氣流曳力作用下的粉體運動是影響噴動床粉體混合的主要原因，噴射區(qū)對粉體有卷吸作用，噴泉區(qū)促進了粉體擴散運動。

3) 不同的粉體在噴動床中循環(huán)運動規(guī)律一致，循環(huán)周期約為9.4 s，在一個循環(huán)周期中，粉體在環(huán)隙區(qū)運動的時間占整個循環(huán)運動周期的90%以上，在噴射區(qū)和噴泉區(qū)的時間極短，整個混合過程經(jīng)歷了3個循環(huán)運動周期。

4) 粉體混合初期，混合均勻性指數(shù)快速增加，在第20秒時快速達到0.95，之后混合過程減緩，在第31秒時，混合均勻性指數(shù)增加至0.99以上，兩種粉體混合均勻。