風(fēng)速對毛豆顆粒流態(tài)化行為影響的仿真分析

楊曉燕，錢 丹，樓曉華，李林建，許 彤，劉冬瑞

（1.四方科技集團(tuán)股份有限公司，江蘇南通 226371；2.南通四方罐式儲運(yùn)設(shè)備制造有限公司，江蘇南通 226371）

0 引言

流態(tài)化技術(shù)具有混合效率高、傳熱傳質(zhì)快等優(yōu)點，在生物質(zhì)能源、礦業(yè)和食品等行業(yè)應(yīng)用廣泛［1-2］。目前在食品行業(yè)，流態(tài)化技術(shù)主要用于蔬果類食品的快速干燥、凍結(jié)等。毛豆速凍常采用流態(tài)式速凍設(shè)備，將毛豆顆粒置于水平網(wǎng)篩上，在高速低溫氣流作用下，毛豆顆粒層產(chǎn)生“懸浮”，從而實現(xiàn)快速凍結(jié)效果［3］。

毛豆顆粒較大，屬于Geldart D顆粒范疇，在食品技術(shù)領(lǐng)域，對于毛豆流態(tài)化的研究報道不多見。在生物質(zhì)、能源領(lǐng)域中，有許多關(guān)于Geldart D顆粒流態(tài)化行為的研究［4-5］。在食品顆粒流態(tài)化速凍領(lǐng)域，常見有通過布風(fēng)板的結(jié)構(gòu)［6-7］、不同的風(fēng)機(jī)頻率和控溫方式［8］、脈動氣流［9］等方式改善食品顆粒流態(tài)化品質(zhì)的舉措。對于C類（粒徑小于50 μm）和D類（粒徑大于2 mm）顆粒，流化速度高，混合質(zhì)量受到形狀影響較大［10］。

在毛豆顆粒速凍過程中，為解決顆粒邊緣部分難以流化、易粘連等問題，本文利用CFD-DEM氣固耦合方法，探索毛豆顆粒流態(tài)化的行為特征，分析在不同進(jìn)口風(fēng)速下，毛豆顆粒的床層膨脹比和混合效率，為其速凍選擇合理的進(jìn)口風(fēng)速，保障其能夠充分混合均勻，并進(jìn)行相關(guān)試驗驗證。為食品行業(yè)流態(tài)化設(shè)備改進(jìn)和參數(shù)優(yōu)化提供一定的技術(shù)支持。

1 數(shù)學(xué)模型

在CFD-DEM求解過程中，氣相使用Navier-Stokes方程進(jìn)行求解［11］；顆粒作為離散相，使用力學(xué)平衡方程求解［12］。對于1組具有相同大小的顆粒，顆粒的平動和轉(zhuǎn)動約束方程如下：

式中 mp——顆粒質(zhì)量；

vp——顆粒平動速度；

Fc—— 顆粒之間以及顆粒和邊界間的接觸力；

Ff→p——流體對顆粒的作用力；

g——重力加速度；

Jp——顆粒轉(zhuǎn)動慣量；

ωp——顆粒轉(zhuǎn)動速度；

Mc——顆粒旋轉(zhuǎn)產(chǎn)生的扭矩；

Mf→p—— 流體相的速度梯度對顆粒產(chǎn)生的額外扭矩。

相互作用力方程：

式中 FD——阻力項；

FN→D—— 非阻力項，包括壓力梯度力F?P、升力FL、附加的虛質(zhì)量力FVM以及其他力Fothers。

式中VP——顆粒體積；

?P——局部壓力梯度。

使用阻力系數(shù)計算流體對顆粒的作用力［13］：

式中 CD——阻力系數(shù)；

ρf——?dú)怏w密度；

A'——顆粒在流體流動方向上的投影面積；

u-vp——顆粒和流體之間的相對速度。

根據(jù)不同的顆粒形狀以及顆粒相的密度，阻力系數(shù)有多種計算方法。針對流化床稠密相顆粒，使用Huilin & Gidaspow阻力系數(shù)模型，計算式如下：

其中，ψ為混合參數(shù)，是與流體體積分?jǐn)?shù)αf相關(guān)的函數(shù)，具體表示為：

CDErgun和CDWen&Yu為另外2種計算阻尼系數(shù)的方法，分別適用于較高粒子濃度（αs≥0.2）和較低粒子濃度（αs＜0.2），且CDErgun考慮了壓力損失的情況，其具體公式為［14］：

式中φ——顆粒的球形度。

2 建模和方法

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

流態(tài)化范圍設(shè)置為500 mm×200 mm×1 000 mm（長×寬×高），毛豆初始堆積高度為100 mm，顆粒填充率為20%，見圖1。

圖1 毛豆顆粒初始堆積Fig.1 Initial stacking status of soybean particles

劃分流體網(wǎng)格，網(wǎng)格大小為20 mm，氣體密度為1.225 kg/m3，黏度為1.789 4×10-5kg/（m·s），環(huán)境溫度為293.15 K；風(fēng)力進(jìn)口采用速度進(jìn)口，方向垂直于截面；出口與大氣相通，采用壓力出口；CFD時間步長為0.002 s。

給定毛豆的最小過篩尺寸為10 mm，為能使顆粒通過的最小方形孔，孔尺寸縱向長寬比y/z=1.8，橫向長寬比z/x=0.85，顆?？倲?shù)為7 175。在顆粒仿真軟件中，設(shè)置顆粒與顆粒之間的靜摩擦系數(shù)為0.2，動摩擦系數(shù)為0.3，恢復(fù)系數(shù)為0.9；顆粒與邊界之間的靜摩擦系數(shù)為0.3，動摩擦系數(shù)為0.3，恢復(fù)系數(shù)為0.3；仿真總時長為10 s，時間步長為0.02 s。不考慮顆粒的破碎、侵蝕等。

在食品流態(tài)凍結(jié)過程中，食品顆粒被放置在網(wǎng)孔均勻分布的網(wǎng)帶上，如圖2所示。

圖2 網(wǎng)帶進(jìn)風(fēng)示意圖Fig.2 Schematic diagram of net belt air inlet

在仿真的過程中，如考慮網(wǎng)孔的情況，則網(wǎng)格的數(shù)量將呈指數(shù)級的增加，大大延長計算時間。FALAH［15］將9個內(nèi)徑為2 mm的圓孔等效為1個相同面積的矩形孔，進(jìn)行相應(yīng)的耦合計算，得到較為準(zhǔn)確的仿真結(jié)果。因此，將有網(wǎng)孔的情況簡化為無網(wǎng)孔進(jìn)行討論。

根據(jù)最小流態(tài)化速度vmf預(yù)測公式［16］：

式中 μf——?dú)怏w黏度；

Ar——阿基米德數(shù)，

根據(jù)WEN等［17］的研究，C1=33.67；C2=0.036 5，計算得到毛豆顆粒vmf=1.584 m/s。

考慮壁面、顆粒間摩擦以及預(yù)測公式偏差等因素，結(jié)合相關(guān)食品流態(tài)化測試數(shù)據(jù)，設(shè)定無網(wǎng)帶工況下進(jìn)風(fēng)速度為2.4 m/s。參考FALAH［15］簡化網(wǎng)孔的做法，已知網(wǎng)帶孔隙率為17.3%，則在有網(wǎng)帶情況下，為保證流量相等，網(wǎng)孔的進(jìn)口風(fēng)速應(yīng)為2.4/17.3%=13.8 m/s。其余條件保持相同，進(jìn)行耦合仿真。

床層高度和壓降是評價流態(tài)化過程的重要參數(shù)。圖3為有、無網(wǎng)孔時壓降、床層高度隨時間變化曲線。在10 s的流態(tài)化過程中，床層高度呈周期性的波動。在仿真中，無網(wǎng)孔的平均壓降峰值和吹起高度峰值均要大于有網(wǎng)孔的情況，但是整體偏差較小。因此，可以使用無網(wǎng)孔的工況代替有網(wǎng)孔的工況進(jìn)行相似的仿真計算。

圖3 壓降、床層高度隨時間變化曲線Fig.3 Curve of pressure drop and bed height over time

2.2 仿真結(jié)果

2.2.1 顆粒流動分布特征

0～2 s內(nèi)毛豆顆粒的分布和流速分布（+z方向）如圖4所示。該時間段剛好為床層高度波動的第1個周期。由于初始堆積較密集，毛豆顆粒在0.1 s呈現(xiàn)出整體向上懸浮的狀態(tài)。隨后在0.2～0.3 s出現(xiàn)顆粒松散分層的現(xiàn)象。隨著床層的繼續(xù)膨脹，顆粒內(nèi)部中心區(qū)域出現(xiàn)橢圓形的氣泡，0.4～0.7 s氣泡逐漸增大。伴隨著氣泡內(nèi)部翻騰的顆粒個數(shù)減少，和氣泡邊緣的顆粒不斷碰撞，橢圓形氣泡逐漸縮減，顆粒也逐漸從氣泡邊緣回到床層底部。在1～2 s期間，毛豆顆粒出現(xiàn)分層涌動的現(xiàn)象，在兩側(cè)出現(xiàn)小型氣泡，氣泡的邊緣顆粒數(shù)較多，氣泡較小，出現(xiàn)時間短，床層的上表面不斷波動。

圖4 毛豆顆粒分布（+z方向）Fig.4 Soybean particle distribution （+z direction）

2.2.2 最大床層高度和壓降

設(shè)置氣體流速區(qū)間為0.2～3.2 m/s，速度增量為0.2 m/s，進(jìn)行16組不同進(jìn)口風(fēng)速下氣固耦合仿真，得到壓降均值和最大床層高度變化曲線，如圖5所示。在固定床階段，壓降值與進(jìn)口風(fēng)速呈正相關(guān)；在流化床階段，壓降值增長幅度較小。擬合曲線的轉(zhuǎn)折點在2.0 m/s，即毛豆顆粒的最小流態(tài)化速度為2.0 m/s，略大于預(yù)測公式計算值（1.584 m/s）。

圖5 毛豆顆粒在不同風(fēng)速下的最大床層高度和壓降變化Fig.5 Maximum bed height and pressure drop changes of soybean particles under different wind speeds

在小于2.0 m/s風(fēng)速工況下，床層最大高度在100 mm左右，近似于顆粒的初始堆積高度，此時未發(fā)生流化。從2.0 m/s開始，床層高度隨著進(jìn)口風(fēng)速的增大而增大，進(jìn)入流化狀態(tài)。

2.2.3 膨脹比和混合率

在流態(tài)化的過程中，當(dāng)氣體不斷注入床層，顆粒間的間隙增大，出現(xiàn)床層膨脹的現(xiàn)象。為了尋找毛豆顆粒的最佳流態(tài)化速度，引入膨脹比與混合效率。采用CAXA軟件對圖像進(jìn)行處理，計算得到毛豆在流態(tài)化過程中的膨脹比：

以+z方向為正方向選取圖像，在CAXA軟件中使用樣條曲線的方法，對顆粒床層外緣輪廓進(jìn)行描繪，并測量計算出堆積面積。

在風(fēng)速2.0 m/s之前，床層高度未發(fā)生明顯變化，因此以此速度值為起始，計算不同風(fēng)速下的膨脹比，如圖6所示。隨著進(jìn)口風(fēng)速增大，毛豆顆粒的膨脹比增大，即床層高度增高，與最大床層高度的變化情況相似。盡管風(fēng)速越大，床層膨脹越明顯，然而進(jìn)口風(fēng)速越大，對風(fēng)機(jī)要求越高，能耗越高。

圖6 不同風(fēng)速下膨脹比Fig.6 Expansion ratio under different wind speeds

將初始時刻的毛豆顆粒均分為上、下2層，在最后時刻記錄上、下2層顆粒在測量區(qū)域的占比，如圖7（a）所示。從圖7（b）毛豆顆粒的混合情況看，在2.0 m/s時，上層與下層顆?；緵]有混合；在2.6 m/s時，上層與下層顆粒的混合率均在50%左右，說明在床層高度方向上，顆?；旌线_(dá)到比較均勻的狀態(tài)。隨著風(fēng)速的增大，混合率基本保持不變。因此，綜合考慮膨脹比、風(fēng)機(jī)能耗以及混合效果，在初始堆積高度為100 mm的工況下，2.6 m/s為毛豆顆粒的最佳流態(tài)化速度。

圖7 不同風(fēng)速下混合效果Fig.7 Mixing effect under different wind speeds

3 流態(tài)化試驗設(shè)置

流化床實驗平臺如圖8所示。風(fēng)機(jī)將蒸發(fā)器提供的冷量通過空氣輸送至風(fēng)道中，經(jīng)過網(wǎng)帶傳遞到食品顆粒上，使食品顆粒產(chǎn)生翻滾、沸騰甚至懸浮的狀態(tài)，實現(xiàn)食品顆粒的冷凍，最后熱空氣從網(wǎng)帶上方流回蒸發(fā)器進(jìn)行熱交換。蒸發(fā)器外部連接冷凝器、壓縮機(jī)等設(shè)備，完成制冷循環(huán)。在網(wǎng)帶上用塑料板隔出500 mm×200 mm×400 mm的范圍進(jìn)行試驗，將毛豆平鋪在網(wǎng)帶上。該網(wǎng)帶的孔隙率為17.3%，毛豆的堆積厚度為100 mm。通過工業(yè)相機(jī)進(jìn)行拍攝，以50幀/s的幀數(shù)拍攝10 s。

圖8 流化床實驗平臺Fig.8 Fluidized bed experimental platform

圖9是毛豆的流態(tài)仿真和試驗工況的對比結(jié)果，截取5個時刻的圖像。試驗拍攝圖片顯示，毛豆床層表面出現(xiàn)連續(xù)的波動起伏，在風(fēng)力作用下局部區(qū)域突起，隨后伴隨著少量毛豆顆粒的噴射行為。這與仿真現(xiàn)象較為一致，證明仿真結(jié)果較為準(zhǔn)確。

圖9 仿真和試驗結(jié)果對比Fig.9 Comparison of simulation and experiment results

4 結(jié)語

本文以毛豆顆粒為研究對象，基于歐拉-拉格朗日方法對毛豆流態(tài)化行為進(jìn)行描述，通過CFD-DEM仿真手段探究不同風(fēng)速下毛豆的最大床層高度、壓降變化，得到毛豆的初始流態(tài)化速度為2.0 m/s。從膨脹比和混合率指標(biāo)出發(fā)，確定在初始堆積高度為100 mm時，毛豆顆粒的最佳流態(tài)化速度為2.6 m/s。研究對食品流化床設(shè)備的優(yōu)化設(shè)計，以及進(jìn)一步提高食品速凍過程中的混合傳熱效率具有一定的參考價值。