紅外噴動床數值模擬中兩種湍流模型的適用性

侯志昀 段 續,2 任廣躍,2 李琳琳 徐一銘

(1. 河南科技大學食品與生物工程學院，河南 洛陽 471023；2. 糧食儲藏安全河南省協同創新中心，河南 鄭州 450001)

干燥是一個能源密集型過程，消耗了食品加工行業20%～25%的能源。如何節約能耗和保持干燥產品質量是食品干燥過程中最值得關注的兩個問題[1]。當前主流的熱風干燥技術設備操作簡單，但存在能源利用率較低，產品質量較差的問題[2]。真空冷凍干燥技術產品質量高，但能耗大、效率低的問題制約了其推廣應用[3]。因此，食品干燥領域越來越注重利用新型的加熱源與其他傳熱方式耦合的組合干燥技術，如以紅外、微波、射頻等為代表的新型加熱方式的引入，極大地提高了食品干燥領域的能源利用率[4]。紅外輻射在空氣中傳播時損失較小，可將熱直接輻射到被加熱體的表面，由于不存在傳熱界面，且具有穿透性，物料內部形成熱量堆積，同時被干燥的物料水分不斷蒸發帶走熱量，使得物料表面溫度下降，物料內部溫度遠高于外部溫度，所以在干燥過程中物料的濕度梯度和溫度梯度保持一致，極大地提高了干燥速率[5-6]。但紅外輻射干燥極易受到輻射距離、料層厚度、熱源分布等因素影響，從而造成干燥不均勻的問題[7]。噴動床技術被用于小麥(流動性小、顆粒狀)干燥，是由于其干燥過程中物料顆粒與熱空氣的良好接觸和循環運動使得干燥均勻性和傳熱傳質效率極高，非常適合散粒狀熱敏性食品的干燥[8-9]。但噴動床干燥技術最大的問題是需前置空氣加熱器，出風溫度太高，設備的熱效率低。結合紅外干燥技術和噴動床干燥技術的優缺點，將紅外加熱方式用于噴動床，開發新型的紅外噴動床干燥技術，可以利用紅外輻射對噴泉區的物料進行高效直接加熱，省去前置空氣加熱器，從而有效降低熱損失，同時，又消除了單一紅外干燥加熱不均勻的缺點[10-13]。課題組[14]提出了紅外噴動床的設計方法并進行了相關試驗，證實了紅外噴動床具有良好的節能干燥特性。

對紅外噴動床的性能進行準確預測是其設計及應用的必要條件，然而其干燥過程對物料的質熱傳遞過程及其影響因素的研究較復雜。與其他試驗方法相比，使用計算流體力學(CFD)方法為研究和解決結構復雜，操作環境危險等流體力學問題建立了一種新的研究方法[15]。應用CFD分析研究各種干燥機理，有助于提高工藝效率[1]。然而在對紅外噴動床干燥進行數值模擬及性能預測時，湍流模型的選取至關重要。Miltner等[16]為選擇合適湍流模型進行自由射流的模擬，對多種湍流模型 (S-A、Standardk-ε、 RNG、Realizablek-ε、Standardk-ε、SSTk-ε、RSM)進行了試驗驗證及比較； Shives等[17]在對潮流場驅動盤數值模擬中，比較了 Standardk-ε模型與 SSTk-ε的適用性；于曉麗等[18]采用兩方程模型(SSTk-ε與Realizablek-ε)對潮流能水輪機進行了數值模擬適用性研究；李東耀[19]在流化床氣固兩相流模擬中對Standardk-ε、 RNGk-ε、Realizablek-ε3種湍流模型進行了模擬比較研究其適用性。近年來國內外學者用Standardk-ε[20-22]、RNGk-ε[15]湍流模型對噴動床氣固兩相流動進行了數值模擬研究，但紅外噴動床干燥青豆的數值模擬尚未見報道。

研究擬采用不同的湍流模型(Realizablek-ε與SSTk-ε)對紅外聯合噴動床進行數值模擬，分析不同工況下紅外聯合噴動床干燥室流體域內的速度場、溫度場與物料顆粒濃度分布以及氣相和顆粒相速度分布變化情況，并進行紅外聯合噴動床試驗驗證，為選擇合適的湍流模型對紅外聯合噴動床在特定的湍流條件下的三維流場進行描述提供依據。

1 湍流模型

1.1 Realizable k-ε湍流模型

氣相的湍流運輸方程采用Realizablek-ε方程[23]。

k湍流動能方程：

(1)

ε湍流動能耗散率方程：

(2)

式中：

k——氣相的湍動能，J/kg；

ε——氣相湍動能的耗散率，m2/s3；

σk、σε——k和ε相應的普朗特數；

ui——氣相速度矢量；

ρ——氣相的密度，kg/m3；

v——氣相的速度，m/s；

E——平均應變率張量的模量；

μt——氣相湍流對應的黏性系數；

Gk——湍動能增量(由平均速度梯度產生)；

t——時間，s。

Realizablek-ε模型是一個半經驗模型。使用的Realizablek-ε經驗常數[24]為：C1=1.44，C2=1.92，σk=1.0，σε=1.2。

1.2 SST k-ε湍流模型

Menter[25]綜合近壁區Standardk-ε模型的穩定性及邊界層外部Standardk-ε模型獨立性的優點，提出了BSLk-ε湍流模型φ3，可用Standardk-ε模型φ1與Standardk-ε模型φ2的加權相加來表示，在此基礎上又在湍流黏度定義中考慮了湍流剪切應力的傳輸形成了SSTk-ε湍流模型。

φ3=φ1f1+φ2(1-f1)，

(3)

式中：

f1——混合函數。

SSTk-ε湍流模型方程[18,26]為

(4)

(5)

式中：

ω——比耗散率，s-1；

k——湍動能，J/kg；

Gk——湍流動能k的生成項；

Gω——比耗散率ω的生成項。

模型中常數[18,24,26]：β*=0.09，σω=2.0，σω2=0.856。

2 紅外噴動床算例描述

2.1 紅外噴動床模型

采用的紅外輔助熱風噴動試驗設備平臺為研究室自主設計研制(見圖1)。該試驗系統包括風機送風裝置，入口氣體加熱裝置，試驗物料顆粒進、出料設施，試驗數據測量、采集和控制設備等。對試驗設備內所用紅外板在試驗進行前先設定波長(10 μm)和功率(500 W)。出口空氣溫度由位于噴動床蓋子上的溫度感應器測得，風速由風速感應器測得。

1. 噴動床頂蓋 2. 溫度感應器 3. 控制臺 4. 溫度感應器 5. 加熱器 6. 風速感應器 7. 風機 8. 紅外噴動床

結合紅外噴動床具體的設備結構尺寸將其干燥室按照三維幾何形狀模型處理，滿足連續介質模型在計算時耦合多個輸運方程的需求，選擇Solidworks軟件構建三維立體幾何模型，如圖2所示。

圖2 紅外噴動床干燥室幾何模型Figure 2 Infrared spouted bed drying chamber model

網格劃分選擇ANSYS ICEM CFD(2019R1)商用軟件，采用非結構化網格技術干燥室計算流體域整體生成四面體網格，然后將其通過ANSYS Fluent(2019R1)軟件再轉變成多面體網格，經網格無關性驗證最終確定的網格數量是130 861。具體的網格劃分如圖3所示。

圖3 紅外噴動床干燥室網格模型Figure 3 Grid model of infraredspouted bed drying chamber

2.2 邊界條件及初始條件設定

2.2.1 邊界條件 依據紅外噴動床干燥室幾何模型計算流體域特點，邊界條件主要考慮氣流的進入、排出及干燥室固體壁面3個方面條件的設定。

(1) 噴動床入口：采用速度入口，氣相速度為：vx=0；

vy=vg；vz=0。vg為給定值。

(2) 固體壁面：以無滑移考慮干燥室壁面的邊界，不考慮速度因素。

(3) 噴動床出口：依據紅外噴動床干燥設備的干燥室結構情況，將壓力出口確定在干燥室上部蓋子的圓形排氣孔位置，表壓為0。

2.2.2 初始條件 根據紅外噴動干燥設備結構特征，啟動數值計算時必須設置初始條件才能進行，具體設置見表1。紅外噴動床干燥室壁面材料在Fluent中設置為鋁。

表1 初始條件設置

模擬過程中Fluent開啟能量方程，多相流模型選用雙歐拉模型，湍流方程選擇Realizablek-ε或SSTk-ε，輻射方程選擇離散坐標模型(DO)，曳力模型選擇Gidaspow模型，顆粒體積黏度選擇Lun et al.模型，徑向分布函數選擇Lun et al.模型，界面面積選擇Ia-symmetric，顆粒黏度選擇Syamlal-obrien。壓力—速度求解選擇PRESTO算法。

2.3 試驗數據采集

選擇青豆為試驗物料，在風速8 m/s、溫度50 ℃下進行干燥。依據試驗監測所獲得的數據，利用ANSYS Fluent(2019R1)軟件模擬兩種湍流模型(Realizablek-ε模型、SSTk-ε湍流模型)的適用性情況。

3 結果對比及分析

3.1 殘差曲線對比

兩湍流模型在運算時均設置時間步長為0.001 s，時間步數為10 000，每一時間步長最大迭代數為20。由圖4、圖5可知，Realizablek-ε和SSTk-ε最后均達到收斂，說明兩種模型的運算結果均是可信的。運算過程中Realizablek-ε的波動整體上要比SSTk-ε更大一些，而Realizablek-ε的總迭代數比SSTk-ε小，運算時間短。

圖中箭頭處從上到下依次為連續性、u方向空氣相速度、u方向顆粒相速度、v方向空氣相速度、v方向顆粒相速度、w方向空氣相速度、w方向顆粒相速度、空氣相能量方程、顆粒相能量方程、k方程、ε方程、離散坐標模型、顆粒體積分數

圖中箭頭處從上到下依次為連續性、u方向空氣相速度、u方向顆粒相速度、v方向空氣相速度、v方向顆粒相速度、w方向空氣相速度、w方向顆粒相速度、空氣相能量方程、顆粒相能量方程、k方程、ε方程、離散坐標模型、顆粒體積分數

3.2 紅外噴動床噴動過程對比

由圖6、圖7可知，兩模型從初始狀態開始噴動到達穩定噴動狀態運算用時相同。在0.000～0.378 s時兩模型的圖形有明顯差別，SSTk-ε顆粒體積分數變化幅度比Realizablek-ε的更大；但在0.252～0.378 s時能看到噴動區左偏較明顯，可能是不同模型的運算殘差不同。

圖6 Realizable k-ε噴動形成過程顆粒體積分數演變流程圖Figure 6 Flow chart of particle volume fraction evolution during spout formation of Realizable k-ε

圖7 SST k-ε噴動形成過程顆粒體積分數演變流程圖Figure 7 Flow chart of particle volume fraction evolution during spout formation of SST k-ε

3.3 顆粒體積分數徑向及軸向分布對比

由圖8、圖9可知，由室內中心到近壁面處顆粒濃度逐漸增大，室中心位置為噴動區域顆粒濃度小，而近壁面環形空間區域濃度大，噴動和環形空間兩區域的交界面附近形成了一個明顯的過渡區域。在40，70 mm床高處，緊貼壁面處顆粒的體積分數略低，是由于顆粒運動快致使局部區域顆粒濃度降低，相互間空隙率增加。40 mm處，兩模型運算結果差異較大，70，100 mm處兩模型運算結果較接近。

圖8 Realizable k-ε不同床高顆粒體積分數分布圖

圖9 SST k-ε不同床高顆粒體積分數分布圖

3.4 軸向及徑向顆粒速度大小分布對比

圖10、圖11所示分別為同一溫度(50 ℃)同一風速(8 m/s)下，Realizablek-ε和SSTk-ε在不同床高40，70，100 mm處的顆粒速度分布圖。

圖像顯示速度最大的位置是在干燥室的縱軸附近，而速度最小的位于壁面附近。從圖10、圖11可以看出，由干燥室中心處到近壁面處的過程中存在過渡區域，這是噴動區和環隙區之間的交界面。環形空間區域隨著床高高度的上升顆粒速度呈現出由大到小的變化趨勢，歸因于噴動床干燥室由下向上(40 mm→70 mm→100 mm)堆積粒子增多，粒子下降流動受到阻礙更大速度減弱。

從圖10、圖11可以看出，在40 mm處兩模型運算結果差異略大，70，100 mm處兩模型運算結果較接近。

圖10 Realizable k-ε不同床高顆粒速度分布圖

圖11 SST k-ε不同床高顆粒速度分布圖

3.5 不同床層高度顆粒溫度分布對比

由圖12、圖13可知，干燥初期，床高40，70 mm處的顆粒層橫截面的中心噴動區域溫度最低，環隙區域溫度最高。這是由于噴動區的氣流速度高，傳熱傳質速率高，熱傳遞快，而環隙區顆粒密度大、氣流滲透互穿量小，熱傳遞慢，界面處則介于中間；此外，干燥室壁面上有紅外輻射板，所以靠近紅外輻射板的顆粒相溫度和壁面溫度總體比較高。紅外輻射板沿X軸對稱，均勻分布在干燥室兩側，不同床高處溫度最高的區域分布在沿Z軸對稱的兩側，是由于兩個紅外輻射板的熱輻射在此處交匯。不同床高處溫度最低區域在干燥室中心部分且隨床高的增加該區域進一步擴大，這是由于中心部分是噴動區域相對于其周圍的環形空間區域固體物料顆粒分布較少；在床高70 mm處顆粒相中心區域比床高40 mm處小，可能是當前時刻下噴動區裹挾的顆粒剛好到達這一高度的較多。

圖12 Realizable k-ε不同床高顆粒溫度分布圖

圖13 SST k-ε不同床高顆粒溫度分布圖

由圖12、圖13可知，40，70，100 mm處均有一定差異，其中40，70 mm處的差異較明顯。

3.6 兩模型與試驗對比

由圖14可知，Realizablek-ε與SSTk-ε分別與試驗數據的誤差都在合理范圍內，而兩模型的模擬數據略有差異。

圖14 Realizable k-ε與SST k-ε模擬數據和試驗數據對比

4 結論

選取Realizablek-ε與SSTk-ε兩種湍流模型，對紅外噴動床干燥青豆顆粒的干燥室流體域三維湍流情況進行了數值模擬，得到了紅外噴動床干燥青豆顆粒的相應參數及流場(顆粒體積分數分布、顆粒速度分布和顆粒溫度分布)情況，并進行了紅外噴動床干燥青豆試驗。結果表明，在紅外噴動床干燥室進氣速度為8 m/s，溫度為50 ℃的工況下，采用兩種湍流模型得到的干燥室流體域中顆粒體積分數分布、速度分布及溫度分布具有相似性，總體效果較好，均可應用于該設備。考慮試驗所用紅外噴動床干燥設備的實際工況及運算耗時等因素，選擇Realizablek-ε湍流模型進行描述更加合理。后續將進一步選擇其他類型模型如曳力模型、傳熱模型等進行優化對比研究，以獲得更加貼近紅外噴動床實際干燥工況的模型。