基于動態雙重網格下噴動床滯止區流動特性CFD-DEM數值模擬

王洪遠，紀律，孟繁旭，李斌，楊建蒙，陳海生

（1 華北電力大學能源動力與機械工程學院，河北保定071003； 2 中國科學院工程熱物理研究所，北京100190）

引 言

噴動床作為流化床的一種，其應用領域十分廣泛，如農業顆粒物料（稻谷、油菜籽、小麥等）的干燥和去皮、去殼，農作物種子以及制藥工業中的包衣和化學肥料的造粒，以及固態物料的混合等[1]。噴動床在工作時床層受到通過小孔或噴嘴射入氣流的影響，中央射流內的固體顆粒被氣流抬升到一定高度，此時顆粒的重力大于射流作用在顆粒上的力，顆粒開始下降，之后又被氣流重新朝上帶起而呈噴泉狀，從而導致顆粒在噴動床內往復運動。噴動床的工作效果主要受床體內的氣固兩相流運動影響，而床體內的氣固兩相流內部的流動規律又是一個極為復雜的過程組合，關于內在規律的研究還在進行。常規的實驗觀察對流化床內氣固兩相流動特性的研究難以獲得顆粒詳細的運動信息，而計算機模擬技術的發展，有助于獲得微觀層次上顆粒的運動信息[2-7]。

離散單元法(discrete element method,DEM)是目前較為先進的研究方法，憑借能追蹤每個顆粒的位置、受力等運動信息的優點，其被廣泛應用到氣固兩相流研究中[8-9]。國內外學者分析流化床內氣固兩相流動特性時側重于CFD-DEM 和LBM-DEM 方法，但由于LBM-DEM 方法較于CFD-DEM 方法只能對床層尺度較小、顆粒直徑較小的工況進行模擬，所以對于尺寸較大的噴動床常采用CFD-DEM方法開展。王敬哲[10]建立了氣流-噴動床的三維數值模型，對床內的氣固兩相流動特性進行了模擬，較為完整地揭示了噴動的形成過程及氣固兩相在床內的速度分布情況。彭麗等[11]利用計算流體力學和離散元(CFD-DEM)方法模擬了擬二維矩形鼓泡床內的氣固流動特性，分析了固相彈性、恢復系數對流場間歇性的影響。Olaofe 等[12]則模擬了顆粒在矩形噴動流化床內的混合過程。Hilton等[13]在CFDDEM 模型下，選取改進的壓力梯度力模型與非球型顆粒（橢圓形、多邊形、多球形）對氣固噴動床的動力特性進行研究。楊建蒙等[14]基于MFIX 軟件，利用MFIX 開源程序中的MFIX-DEM 模型，比較了Gidaspow 和Syamlal-O’Brien 2 種曳力模型下對壓力源項處理的2種方式Model A和Model B的三維噴動流化床氣固兩相流動特性，并與同工況下的實驗結果進行了對比。

諸多研究者在對氣固兩相流領域的研究中各有建樹，他們的研究推動了氣固兩相流技術的進步[15-19]，但在現有的研究中均未對噴動床內滯止區現象進行深入研究，因此本文首先建立了基于動態雙重網格技術的CFD-DEM 方法，然后進行了噴動床徑向混合實驗與模擬研究，同時與單網格模型進行對比，驗證了動態雙網格方法結果的準確性，然后利用動態雙重網格方法建立了不同進口氣速和不同初始堆積高度下的計算模型，同時對滯止區內顆粒進行示蹤標記，分析示蹤區域的顆粒流動情況，得出顆粒在不同模型的滯止區流動混合的規律。

1 數學模型

1.1 氣固兩相的運動模型

氣固兩相運動過程在流化床內非常復雜，而CFD-DEM 數值模擬技術已成為解決眾多領域顆粒材料問題的有力工具。

顆粒之間不斷碰撞并受氣體夾帶作用運動，同時近壁面顆粒還與床體碰撞，另外重力、電場力和磁場力等外界物理場產生的作用力也會對顆粒造成影響。因此本研究顆粒相的碰撞模型選擇基于時間驅動原理的軟球模型，其運動方程見文獻[20]。氣相運動模型則采用氣固兩相耦合的Navier-Stocks方程，湍流模型采用k-ε模型，詳細的控制方程參見文獻[20]。

1.2 氣固兩相間耦合

1.2.1 曳力模型 氣固相間曳力表達式通常表示為曳力系數乘以速度差的形式：

根據Gidaspow 的曳力模型，氣相對顆粒的曳力系數[21]如下：

1.2.2 氣固兩相間動量耦合 氣固兩相間的作用力應該遵循牛頓第三定律，顆粒受到了流體對它的曳力，相應的流體也會受到顆粒對它的反作用力，這個力應該與曳力大小相等，方向相反。所以氣固兩相間的動量耦合公式為：

1.2.3 動態雙重網格模型 圖1為三種網格劃分示意圖，為了便于分清氣相網格和顆粒相網格，用紅色的網格代表顆粒相粗網格，用黑色的網格代表氣相細網格。由圖可見，在單網格模型下氣相和固相的計算都在同一網格內進行，即氣相和顆粒相的計算共用一套網格[22-27]。在固定雙重網格模型下氣相和顆粒相分別有一套網格，氣相的細網格被包含在顆粒相粗網格中，所以如顆粒的速度、溫度、床體空隙率等在顆粒相的粗網格中計算；而氣相的流場和溫度場等在細網格內計算，在計算過程中通過氣固兩相間的耦合使兩套網格之間即時映射，在同一時刻下通過網格間數據的傳遞來交換相間能量和動量等信息。

圖1 網格模型劃分示意Fig.1 Schematic diagram of grid model division

動態雙重網格模型與前兩種網格劃分情況不同之處在于，動態雙重網格方法主要在顆粒相的網格的設置有極大改變，同樣氣相采用與固定雙重網格相同的細網格，用來在計算過程中映射氣相的計算信息，同樣根據計算對象的不同可以改變氣相網格的疏密程度；而顆粒相網格劃分則不以傳統網格劃分方法為依托，其以顆粒的球形中心為網格中心，以圖中顆粒i為例，在某時刻下以i顆粒的中心作為i顆粒網格的中心，同樣其網格尺寸可以根據i顆粒的尺寸來改變，其網格尺寸設置值的范圍為顆粒直徑的3～4 倍[28]，同時判斷統計位于i顆粒網格內的其他顆粒個數，以及氣相網格編號和個數，在該時刻其他顆粒的網格也可以通過此方法劃分，所以稱隨顆粒運動的顆粒相網格為動態粗網格，而氣相網格稱為動態細網格。在計算過程中通過氣固兩相耦合進行數據映射和相間的能量交換[29]。此外如何判定顆粒是否在所計算網格之內，是根據該顆粒的質心是否在計算的網格之內。只要該顆粒的質心在所計算單元之內，則可以認為該顆粒整體都屬于此計算單元。

1.2.4 網格間數據傳遞 在固定細網格中計算的氣體速度等物理量可通過求解氣相方程來得到[30]。而在動態粗網格中計算的顆粒相的速度等物理量則是由其包含的所有細網格中對應物理量的加權平均得到。

動態粗網格中氣體速度計算公式如下：

粗網格映射給其包含的所有細網格空隙率為：

2 模型驗證

為了探究噴動床混合規律和驗證動態雙重網格的準確性，實驗過程選取直徑4 mm的紅黃兩種顏色的亞克力球形顆粒作為實驗物料，在噴動床流化過程中來觀察顆粒混合特性。

圖2 為進口表觀氣速1.5 m/s、初始堆積高度0.16 m的流化床內顆粒實驗和模擬的徑向混合序列圖，圖3 為數值計算的顆粒徑向混合過程垂直速度分布云圖。從圖中可知，在t=0 s 時，兩種顏色顆粒沿床高方向彼此分離，隨后進口氣流帶動底部中間區域顆粒向上運動，而靠近壁面兩側區域的顆粒則向下運動，在顆粒往復運動過程中氣流沖擊床體內部而形成氣泡，隨后氣泡不斷擴張帶動床體顆粒向上膨脹，紅色顆粒開始向上混合形成一個“樹”狀圖形，而黃色顆粒則向下運動來填補缺失的空間。當顆粒在擴散進程中自身重力大于氣流的曳力時，顆粒便開始下行，同時與其他顆粒不斷碰撞混合，緩慢下移至床層底部后，又被卷吸至噴口附近重新進入下一流化循環。顆粒在此循環運動中和其周圍顆粒不斷進行著混合與擴散，其混合程度也不斷加深，最后隨著時間的推移，會達到顆粒間混合的動態平衡。

圖2 顆粒徑向混合過程實驗與模擬對比Fig.2 Comparison between experiment and simulation of radial mixing process of particles

圖3 顆粒徑向混合過程垂直速度云圖Fig.3 Vertical velocity nephogram of particle radial mixing process

此外觀察0～2 s內實驗與模擬圖片可以發現，在床層底部壁面與左右側壁面由于空間位置的限制，存在一定的滯止區域，該區域內顆粒由于遠離噴口位置，氣流不能及時擴散，造成流動緩慢，在流化混合過程與壁面形成“三角”區域。該區域在流化時間內的面積隨時間推移不斷縮小，同時觀察利用動態雙重網格計算的模擬區域也可發現此規律，因此通過徑向混合實驗和模擬對比也論證了動態雙重網格的計算結果與實驗混合流化過程結果高度一致。在此模型基礎上，本次研究還利用數值模擬中顆粒追蹤技術開展滯止區顆粒流動特性的研究。

3 數值模擬及結果分析

模擬對象為單孔射流準三維矩形噴動床，其長寬高分別為150 mm、4 mm、900 mm。氣流入口設在床層底部中心，寬度為0.01 m，噴動床模型如圖4 所示。固相采用了1904 個球形顆粒，顆粒直徑4 mm、密度2700 kg/m3。氣相選擇空氣，密度和黏度分別取值1.205 kg/m3和1.8×10-5kg/(m·s)。具體模擬計算參數見表1。將流體在壁面處設置為無滑移邊界條件，顆粒在壁面處設置為滑移邊界條件。

圖4 噴動床幾何模型Fig.4 Geometric model of spouted bed

表1 數值模擬計算參數Table 1 Calculation parameter in numerical simulation

3.1 動態雙網格模型與單網格模型對比分析

3.1.1 滯止區混合特性對比 為了體現動態雙網格技術的計算準確性，在計算過程中將床體內顆粒分為四個區域，利用黑色區域作為滯止區參考，其距離左側和右側壁面為0.008 m，紅色區域作為噴動區參考，亮青色區域用來代表床層上半部分區域，綠色區域用來代表與滯止區相鄰的區域。

圖5和圖6給出了初始堆積高度為0.18 m，在噴口氣流速度為27 m/s時的單網格模型和動態雙網格模型在2 s 內床體顆粒的混合流動情況。從圖中可以看出，氣流由底部進入床體，導致床層出現松動，在初始時刻由于空間狹小，導致噴口附近的局部壓力瞬時增大，由此產生一條裂縫，推動顆粒往上運動。而整個床體中心處的顆粒速度較大，并且由于顆粒之間的動能傳遞，導致床體上部的部分顆粒躍起，并隨著時間的推移，氣流往上繼續流動。在0.2～0.3 s 時刻附近，床層內部產生一個較大氣泡，隨后噴動區底部的顆粒在顆粒之間產生氣泡的帶動下一起朝上運動；到達床層表面附近即噴泉區后顆粒向床體壁面移動分散；由于位于環隙區上部的兩側壁面附近的氣體瞬時速度較小，導致顆粒受到的曳力減小，在重力作用下沿著床體壁面向下部運動，并在這個過程中參與了環隙區內顆粒的混合，顆粒受氣流卷吸影響，下降到床層底部后會重新進入噴動區，從而實現了顆粒在流化床內循環。

圖5 單網格模型2 s內混合序列圖Fig.5 Mixed sequence diagram within 2 s of single grid model

圖6 動態雙網格模型2 s內混合序列圖Fig.6 Mixed sequence diagram within 2 s of dynamic double grid model

在這個循環過程中，由于壁面附近的顆粒遠離噴口，其氣流速度較小，致使床底與壁面附近的顆粒被壓在床的底部角落處，形成流動停滯區，也稱“死區”。由圖5 可知，利用單網格方法在初始時刻，滯止區內黑色顆粒受到床層的膨脹后回落沖擊，其高度隨著時間推移逐漸變低，t=0.2 s 時,滯止區高度變為原來的一半；t=0.4 s 時,滯止區內顆粒完全隨著噴口氣流進入噴動區；t=0.6 s 時，觀察可以發現在床層底部黑色顆粒所占比例大大降低；在隨后的時間內整個床層內顆粒分布均勻。觀察圖6，則發現在0～2 s 內床層底部的滯止區示蹤區一直存在，黑色顆粒在滯止區的含量隨著時間的推移逐漸減少，t=0.5 s 時其滯止區黑色顆粒的高度變為原來的一半，t=2 s 時滯止區內仍然存在黑色顆粒，此時刻的滯止區內被其他顆粒所占據，重復此循環。

3.1.2 滯止區示蹤顆粒流動高度對比 為了對滯止區流動特性有更為清晰的描述，在研究中對黑色標記的顆粒進行統計追蹤，取每個時刻的黑色示蹤顆粒的平均高度作為參考，得到0～2 s內示蹤顆粒的流動高度如圖7所示。

圖7 兩種網格方法下滯止區示蹤顆粒平均高度Fig.7 Average height of tracer particles in stagnation zone under two grid methods

圖7 中初始床層高度為0.18 m,進口速度為27 m/s。可以看出，隨著時間的推移，兩種方法下的滯止區示蹤顆粒平均高度不斷下降，利用單網格方法比利用動態雙網格方法滯止區示蹤顆粒平均高度要低，并且利用單網格方法滯止區示蹤顆粒平均床高變化較快，在t=0.5 s 滯止區內示蹤顆粒平均床高趨于0，這說明在此后滯止區示蹤顆粒含量很少，而利用動態雙網格模型滯止區示蹤顆粒平均床高變化較為緩慢，滯止區現象更為明顯。此外，對動態雙網格滯止區示蹤顆粒平均高度在t=0.5 s 和0.9 s 出現波動情況進行分析，可能是由于床層波動造成擠壓，致使示蹤顆粒又往滯止區運動，造成滯止區示蹤顆粒數量增多，進而使平均床高增加。

3.1.3 滯止區示蹤顆粒流動寬度對比 為了探究滯止區黑色示蹤顆粒受噴口擾動的影響，在研究中取每個時刻所有黑色示蹤顆粒與壁面的距離作為參考，得到0～2 s 內示蹤顆粒的平均寬度如圖8所示。

圖8 兩種網格方法下示蹤顆粒平均寬度Fig.8 Average width of tracer particles under two grid methods

圖8 中初始床層高度為0.18 m,進口速度為27 m/s。從圖中可以看出，隨著示蹤顆粒的下滑，造成示蹤顆粒在床層寬度方向逐漸延伸，從而逐漸接近噴口附近最后隨同噴口氣流往上運動。對比動態雙網格方法與單網格方法可以發現，單網格方法示蹤顆粒往噴口方向延伸速度更快，在0.8 s就延伸至噴口附近，而此時刻利用動態雙網格方法其示蹤顆粒才延伸至距左側壁面0.04 m 的位置。而在2 s 時，兩種方法的示蹤顆粒平均寬度又接近一致，可能是由于隨著時間的推移，兩種網格方法下的示蹤顆粒向噴口擴散的數量接近一致造成的。

通過對比單網格模型和動態雙網格模型內關于示蹤顆粒的運動可以明顯發現，利用動態雙網格模型的計算結果中滯止區區域更為明顯，而在采用單網格模型模擬大顆粒稠密氣固兩相流時，由于氣相網格不夠精確，使得氣流速度分布均勻，阻礙了該現象的形成。

3.2 不同速度下滯止區示蹤顆粒流動高度分析

圖9給出了基于動態雙重網格方法下初始床層高度為0.18 m,進口速度為27、30、32 和35 m/s 時的床層滯止區示蹤顆粒在2 s 時刻內的平均高度變化情況。從圖中可以看出，隨著時間的推移，滯止區示蹤顆粒平均高度不斷下降，四種速度下的床層滯止區示蹤顆粒平均高度下降趨勢基本保持一致,在0～0.5 s 內，可以看出隨著進口速度的增加，其滯止區示蹤顆粒平均高度下降速度也隨之增快。此外，動態雙網格滯止區示蹤顆粒平均高度在t=0.5 s 和0.9 s 出現波動，可能是由于床層波動造成顆粒間擠壓，致使示蹤顆粒又往滯止區運動，造成滯止區示蹤顆粒數量增多，進而使平均床高增加。整體來看，在計算時刻內，四種速度下床層內部都有良好的滯止區現象表現，速度的增加對滯止區顆粒高度下降影響不明顯。

圖9 不同速度下滯止區示蹤顆粒平均高度Fig.9 Average height of tracer particles in stagnation zone at different speeds

3.3 不同速度下滯止區示蹤顆粒流動寬度分析

圖10 給出了基于動態雙重網格方法下初始床層高度為0.18 m,進口速度為27、30、32 和35 m/s 時的床層滯止區示蹤顆粒在2 s 時刻內的示蹤顆粒平均寬度變化情況。從圖中可以看出，隨著示蹤顆粒的下滑，示蹤顆粒在床層寬度方向逐漸延伸，從而逐漸接近噴口附近最后隨同噴口氣流往上運動。從圖中可以明顯看出，在0～1 s 內，32 m/s 和35 m/s工況下，其延伸速度一致，較27 m/s和30 m/s的工況更快。32 m/s和35 m/s工況下在0.8 s就延伸至噴口附近，而此時27 m/s 和30 m/s 的工況才延伸至距左側壁面0.04 m的位置。

圖10 不同速度下示蹤顆粒平均寬度Fig.10 Average width of tracer particles at different speeds

3.4 不同初始堆積高度下滯止區示蹤顆粒流動高度分析

圖11 給出了基于動態雙重網格方法下初始床層高度為0.13、0.18、0.23 和0.28 m,進口速度為27 m/s時的床層滯止區示蹤顆粒在2 s時刻內的平均高度變化情況。從圖中可明顯看出，隨著床層初始堆積高度的增加，滯止區示蹤顆粒平均高度下降速度也隨之增加，但是初始堆積高度為0.28 m 的情況下在0～0.3 s下降速度最快，在0.3～2 s下降速度趨于平緩，這可能是由于初始堆積床層較大，而氣體速度更接近最小流化速度造成的。同時對四種初始堆積床高進行對比，發現在0～0.3 s 快速下降后，在之后的下降速度都較為平緩，床層高度較大的床層，在0.3 s后存在較為明顯的波動現象。

圖11 不同初始高度下示蹤顆粒平均高度Fig.11 Average height of tracer particles at different heights

3.5 不同初始堆積高度下滯止區示蹤顆粒流動寬度分析

圖12 給出了基于動態雙重網格方法下初始床層高度為0.13、0.18、0.23 和0.28 m,進口速度為27 m/s 時的床層滯止區示蹤顆粒在2 s 時刻內的示蹤顆粒平均寬度變化情況。可以看出0～0.4 s，床層高度為0.18 m 的示蹤顆粒寬度增長速度最大，床層高度為0.28 m 的變化則最慢。0.4～1.0 s，初始堆積高度為0.13 m 的寬度增長速度最快，床層初始堆積高度為0.28 m 的變化最慢。從四種初始堆積高度下其滯止區寬度變化情況來看，當初始堆積越高其寬度變化越慢，可能是由于床層高度增加導致滯止區顆粒下降速度增快，使其更接近噴口速度導致的。

圖12 不同初始高度下示蹤顆粒平均寬度Fig.12 Average width of tracer particles at different heights

4 結 論

本研究基于動態雙網格方法，使用足夠精密的連續相網格來獲得獨立于網格的解，得到如下主要結論。

（1）徑向混合實驗結果與數值模擬結果有很好一致性。

（2）動態雙重網格CFD-DEM 方法結果表明，由于壁面附近的顆粒遠離噴口，其氣流速度較小，致使床底與壁面附近的顆粒被壓在床的底部角落處，形成流動停滯區，滯止區內的顆粒流動性較差。

（3）當床層初始堆積高度一定時，隨著速度的增加，滯止區顆粒高度下降速率和向噴口延伸速度無明顯變化。

（4）當進口速度不變，隨著初始堆積高度的增加，滯止區顆粒下降速度也隨之增加，但其向噴口延伸速度逐漸變慢。

符 號 說 明

CD,Wen-Yu——有效曳力系數

dp——顆粒直徑，m

Fp——顆粒對氣體的單位體積作用力，N

Fy,i——氣體對顆粒i的曳力，N

Kc——網格內顆粒總數目

kc——計算粗網格內顆粒的數目

n——粗網格包含的細網格數目

tc——粗網格內氣體溫度，K

tf,j——細網格內氣體溫度，K

uc——粗網格內氣體速度，m/s

uf,j——細網格內氣體速度，m/s

uG,uS——分別為氣體和顆粒的速度，m/s

Vc——計算粗網格的體積，m3

Vp——單顆粒體積，m3

ΔV——網格體積，m3

β——曳力系數

εG,εS——分別為氣體和顆粒的空隙率

μG——氣體剪切黏度，m2/s

ρG——氣體密度，kg/m3