異質(zhì)柱形顆粒與球形顆?；旌狭鲃犹匦缘腃FD-DEM數(shù)值模擬分析

王恒，仲兆平，王佳，王澤宇，王肖祎，朱玲莉

?

異質(zhì)柱形顆粒與球形顆?；旌狭鲃犹匦缘腃FD-DEM數(shù)值模擬分析

王恒1，仲兆平1，王佳1，王澤宇1，王肖祎2，朱玲莉1

(1. 東南大學(xué)能源與環(huán)境學(xué)院，能源熱轉(zhuǎn)換及其過程測控教育部重點實驗室，江蘇南京，210096；2. 東南大學(xué)建筑設(shè)計院有限公司，江蘇南京，210096)

針對成型生物質(zhì)顆粒與流化介質(zhì)在氣固流化床中的混合流動過程，采用數(shù)值模擬對該過程進行描述和分析。應(yīng)用Fortran語言編程，構(gòu)建復(fù)雜組分的三維氣固兩相流模型，在介觀尺度下對流化床中異質(zhì)異形顆粒的混合流動模擬問題進行研究。在歐拉和拉格朗日框架下，分別采用不同方法構(gòu)造柱形顆粒：在CFD模型中，采用虛擬小球法描述柱形顆粒，引入體積分數(shù)改進曳力系數(shù)公式；在DEM部分，采用球元疊加法構(gòu)造柱形顆粒，引入軟球模型計算顆粒間的碰撞力?；谏鲜龇椒ǎM表觀氣速分別為1.0 m /s和2.0 m/s時的柱形生物質(zhì)與石英砂混合流動的物理過程，并將模擬結(jié)果與試驗結(jié)果進行對比。研究結(jié)果表明：該方法能較好地模擬在鼓泡床狀態(tài)下異質(zhì)柱形顆粒與球形顆?；旌狭鲃拥倪^程。

CFD-DEM模擬；柱形顆粒；氣固流化床

氣固流化床中的顆粒流動模擬是多相流研究中一個重要的領(lǐng)域。氣固流化床中的顆粒流動是氣流與顆粒之間、顆粒與顆粒之間以及顆粒與壁面之間相互作用的過程，由于固相物料通常不是由單一組分組成，該過程為復(fù)雜的異形異質(zhì)顆粒混合流動過程。針對顆粒在氣固流化床中的流動過程，國內(nèi)外學(xué)者已經(jīng)進行了一些實驗和模擬研究，但研究對象多是單組分顆粒。TSUJI等[1?3]將DEM與CFD耦合的方法引入氣固流化床領(lǐng)域，成功模擬了單組分球形顆粒的流動過程。桂南[4]應(yīng)用DEM-LES方法對兩相流中的顆粒碰撞問題進行研究，研究對象為單組分的球形顆粒。對于單組分非球形顆粒流動的研究也有相關(guān)報道。盧洲等[5]模擬了柱狀顆粒在彎管輸送過程中的碰撞特性；ZHONG等[6]采用球元疊加的方法構(gòu)造柱形顆粒并模擬其在噴動床中的流動；REN等[7?9]采用同樣的構(gòu)造方法分別研究了圓柱形顆粒、玉米型顆粒在噴動床中的流動過程。目前針對氣固流化床中2種不同顆?；旌狭鲃訂栴}的模擬研究較少。LIU等[10]采用實驗方法研究了柱形顆粒與球形顆?；旌狭鲃犹匦?。張勇[11]針對異質(zhì)球形顆粒和異徑球形顆粒在流化床中的混合特性分別進行了探索，YU等[12]采用微動力建模法研究了異徑球形顆粒的混合與分離特性。任冰[13]模擬了等密度的非球形顆粒與球形顆?；旌狭鲃忧闆r。上述研究或針對同質(zhì)異形顆?；蜥槍Ξ愘|(zhì)同形顆粒，而氣固流化床中的物料如固體垃圾、生物質(zhì)燃料等，其形狀多為非球形，通常會與煤、河沙等摻混，從而呈現(xiàn)相對復(fù)雜的流動特性，因此開展異質(zhì)非球形顆粒與球形顆粒的混合流動模擬研究工作具有十分重要的意義。

1 柱形顆粒模型

1.1 虛擬小球法

在CFD部分計算氣固作用力時，需要考慮柱形顆粒與球形顆粒的差異，采用虛擬小球法模擬CFD部分的柱形顆粒。將柱形顆粒視為多個虛擬小球堆積組成，如圖1所示。

圖1 虛擬小球構(gòu)造柱形顆粒示意圖

虛擬小球的直徑定義為f，虛擬小球相對于柱形顆粒的體積分數(shù)為f，當f足夠小，f足夠大時，柱形顆粒的體積近似為兩者的乘積。此時空隙率為單位網(wǎng)格中除去球形顆粒及虛擬小球體積的部分：

在計算氣固作用力時，固相平均速度s計算公式如下：

(2)

式中：c和s分別為柱形顆粒和球形顆粒的體積分數(shù)；c和s分別為柱形顆粒和球形顆粒在單位網(wǎng)格內(nèi)的平均速度。

式中：f和s分別為構(gòu)成柱形顆粒虛擬小球的數(shù)量和球形顆粒的數(shù)量；f和s分別為構(gòu)成柱形顆粒虛擬小球的直徑和球形顆粒的直徑。

1.2 球元疊加法

在DEM部分的計算中，對柱形顆粒采用球元疊加法進行構(gòu)造。球元疊加法[14]即采用多個小尺寸球形顆粒順序排列重疊構(gòu)造成非球形顆粒的方法。圖2所示為球元疊加法的示意圖，球元的數(shù)量和尺寸決定了柱形顆粒的準確程度和計算所用的時間。

(a) 2球元；(b) 32球元

2 數(shù)值模擬方法

在歐拉?拉格朗日框架下，基于CFD-DEM耦合算法，分別建立氣相、固相物理模型。由于固相顆粒由異質(zhì)柱形顆粒和球形顆粒組成，在計算氣固作用力以及顆粒碰撞時，不能按照單組分顆粒模型進行處理，本文對有關(guān)參數(shù)進行了改進。

2.1 氣相物理模型

氣相流體采用連續(xù)性方程與納維?斯托克斯方程進行描述：

(5)

式中：為單位網(wǎng)格內(nèi)的空隙率；g為氣相密度；g為氣相流速；為氣相壓力；為重力加速度；為應(yīng)力張量。g為單位網(wǎng)格的氣固相互作用轉(zhuǎn)換項：

g=d+m+s(6)

式中：d為曳力；m為Magnus升力；s為Saffman升力，與曳力相比，后兩者相對較小，可忽略。

d=(s?g) (7)

式中：為氣固曳力系數(shù)，可根據(jù)空隙率，由Ergun公式[15]或Wen&Yu關(guān)系式[16]進行求解，具體表達式如下：

式中：g為氣相黏度；為顆粒的平均直徑；D為單個球形顆粒曳力系數(shù)。

(9)

顆粒雷諾數(shù)定義如下：

2.2 固相物理模型

固相顆粒的運動在模擬條件下遵循牛頓第二定律，在空間中進行平動和轉(zhuǎn)動，計算公式如下：

(12)

式中：s為顆粒轉(zhuǎn)動慣量；s為顆粒力矩；s為顆粒質(zhì)量；s為重力；c顆粒所受的碰撞力；其他作用力在本文計算中不考慮。異質(zhì)柱形顆粒與球形顆粒的混合運動較傳統(tǒng)的DEM模型更為復(fù)雜，下面對顆粒碰撞力c和氣固曳力d進行詳細描述。

2.3 顆粒碰撞模型

顆粒碰撞采用軟球碰撞模型，軟球模型先計算顆粒下一時間步長的位移，對于球形顆粒，通過比較顆粒質(zhì)心與兩球半徑之和(壁面可以看作靜止顆粒)，判定是否發(fā)生碰撞。對于發(fā)生碰撞的顆粒，進一步計算其與周圍顆?；蛘弑诿骈g的作用力。相對于硬球模型，軟球模型允許顆粒同時與多個顆粒及壁面發(fā)生碰撞。

將球形顆粒的碰撞受力c分為切向力ct和法向力cn：

(14)

式中：t和n為切向和法向彈性系數(shù)；t和n表示切向和法向的相對位移；t和n分別為切向和法向的阻尼系數(shù)；s為顆粒相對角速度；r為顆粒相對線速度。若切向力和法向力之間存在ct＞|s?g|，則發(fā)生相對滑移，此時，切向力的計算公式為：

式中：為顆粒碰撞后的摩擦因數(shù)。

柱形顆粒由球元構(gòu)成，同一個柱形顆粒的球元位置相對固定，不存在碰撞。柱形顆粒的碰撞判定思路是判定表面球元的球心之間的距離與半徑之和的關(guān)系(假設(shè)壁面為靜止顆粒)。柱形顆粒的碰撞可分為以下6種情況：1) 單個球元與單個球元的碰撞；2) 單個球元與多個球元的碰撞；3) 單個球元與壁面的碰撞； 4) 多個球元與多個球元的平行碰撞；5) 多個球元與多個球元的交叉碰撞；6) 多個球元與壁面的碰撞。

柱形顆粒的碰撞方式雖然不同，但各種碰撞受力分析的基本思路是一致的，即柱形顆粒碰撞受力為其所有球元受力的疊加。

2.4 氣固曳力

DEM方法對流場中每一個顆粒的受力都進行了計算，對于球形顆粒，其所受的氣固曳力為

式中：s為球形顆粒的體積；si為球形顆粒的速度。

柱形顆粒所受的氣固曳力，為球元受力的疊加，計算公式為

式中：c為柱形顆粒的體積；fi為柱形顆粒的速度。

3 異質(zhì)柱形顆粒與球形顆?；旌狭鲃幽Ｐ偷慕?/h2>

基于上述CFD-DEM模擬方法，本課題組以成型生物質(zhì)顆粒與石英砂混合流動為研究對象，建立了異質(zhì)柱形顆粒與球形顆粒混合流動的并行計算模型，建立截面長×寬為150 mm×20 mm，高度為1 000 mm的流化床模型，分別計算了表觀氣速為1.0 m/s和2.0 m/s的混合流動模型，實際計算中，顆粒運動集中在流化床下部，為節(jié)省時間，在高度方向上僅計算500 mm以下的顆粒流化過程。

網(wǎng)格尺寸決定了計算域內(nèi)離散節(jié)點的位置，若網(wǎng)格尺寸過大，則計算精度偏低；若網(wǎng)格尺寸過小，則可能使網(wǎng)格尺寸與顆粒粒徑的關(guān)系將不能滿足曳力公式的要求。柱形顆粒由小球或球元構(gòu)成，小球或球元個數(shù)越多，則小球或球元的粒徑就越小。因此，本文選擇如下計算條件：CFD計算部分的柱形顆粒由2 072個直徑為0.8 mm的虛擬小球組成，計算網(wǎng)格尺寸為2 mm，DEM計算部分的柱形顆粒由768個球元疊加而成，球元直徑為3 mm，，計算網(wǎng)格尺寸為3 mm。

壁面邊界假設(shè)為無滑移邊界條件；入口采用速度入口邊界條件，出口采用壓力出口邊界條件。布風方式為均勻布風。相關(guān)模擬參數(shù)見表1。

表1 模擬參數(shù)

4 模擬結(jié)果分析

4.1 模擬結(jié)果與試驗瞬態(tài)對比

表觀氣速1.0 m/s的工況模擬結(jié)果如圖3所示。流化開始時，柱形顆粒與球形顆粒按照坐標進行排列布置，試驗測得該工況下的最小流化速度為0.68 m/s。在表觀氣速1.0 m/s的工況中，固相顆粒有明顯的床層界面，顆粒揚起的最大高度約為靜止床高的2倍。圖3中，顆粒從0 s開始至0.48 s為鼓泡初期，氣體進入床層并形成氣泡，攜帶球形顆粒上升，并將頂層柱形顆粒與球形顆粒一并揚起。顆粒隨著氣泡到達床層表面，氣泡破裂后，顆粒被拋向床層上部，達到最高點后回落并再次與上升的氣泡相遇，從0.6 s至1.34 s，柱形顆粒與球形顆?；旌狭己?。圖4所示為成型生物質(zhì)顆粒與石英砂混合流動特性試驗瞬態(tài)圖。圖4(a)所示為試驗中流化前的顆粒在床中的布置，與模擬對應(yīng)，所用顆粒沒有將柱形顆粒與球形顆?；旌?，圖4(b)~(e)所示為表觀氣速1 m/s的成型生物質(zhì)顆粒與石英砂混合流動瞬態(tài)圖。圖3中0.6 s后的顆?；旌闲袨榕c圖4中試驗拍攝的顆粒流動行為十分相似，顆粒混合效果良好。

4.2 壓力波動特性分析

壓力波動分析方法是對流化床中的壓差進行實時測量，并通過分析其時域和頻域變化特性，對氣固混合流動特性進行定性分析的一種經(jīng)典的研究方法。實驗中采樣頻率為100 Hz。將模擬得到的瞬時壓力與試驗測得的壓力脈動信號進行處理，得到柱形顆粒與球形顆?；旌狭鲃釉囼灪湍M的功率譜密度(PSD)，如圖5所示。壓力脈動信號的功率譜表征的是壓力信號的能量在頻域上的分布，能量最大的點對應(yīng)的頻率為主頻。在1.0 m/s的表觀氣速下，試驗與模擬得到了相似的結(jié)果，主頻均在2 Hz左右，表明兩者具有相似的流化特性。模擬和實驗得到的功率譜均集中在0~5 Hz，表明氣泡行為對流化的影響占主導(dǎo)地位。

4.3 柱形顆粒與球形顆?；旌咸匦苑治?/p>

通過分別統(tǒng)計模擬結(jié)果中柱形顆粒和球形顆粒的平均高度，可分析異質(zhì)柱形顆粒與球形顆粒隨時間的混合情況。不同氣速下柱形顆粒與球形顆粒平均高度隨時間的變化如圖6所示。在0.6 s后，不同氣速下的柱形顆粒和球形顆粒呈現(xiàn)了相似的分布特性，即柱形顆粒與球形顆粒的高度趨于一致，這表明2種顆粒在0.6 s后混合較好，沒有發(fā)生分層現(xiàn)象。不同表觀氣速下，顆粒平均高度在起始流化過程中差異較大：在表觀氣速較大的工況下柱形顆粒與球形顆粒所達到的最大高度也較高，且達到最大的高度需要的時間也更長，該現(xiàn)象符合本文對顆粒運動的描述。

DEM模擬結(jié)果保存了每一個顆粒的運動信息，按照不同的床層高度，統(tǒng)計顆粒沿方向和沿方向的時均速度，可以宏觀的分析顆粒的運動狀態(tài)。圖7所示為混合顆粒沿床高方向上的時均速度分布，為了更準確地描述混合顆粒的特性，數(shù)據(jù)未包括0.6 s前起始流化狀態(tài)的部分。床寬為150 mm，因此，75 mm處為流化床的中心位置，從75 mm處向兩側(cè)觀察可以發(fā)現(xiàn)，混合顆粒的氣速在方向大致對稱分布，靠近兩側(cè)壁面的顆粒以靜止或下落為主，中心區(qū)域附近的顆粒則以向上運動為主。不同床層高度的顆粒具有不同的運動規(guī)律：低床層(=5~60 mm)的顆粒距離布風板較近，顆粒速度受氣相射流的影響較大，在方向中間位置呈現(xiàn)規(guī)律性波動；=0 mm和=150 mm附近，顆粒速度為0 m/s或略小于0 m/s，其原因是布風板開口離兩側(cè)壁面距離略遠，氣相對顆粒作用力較小。位于中床層的顆粒(=80~120 mm)在方向中間位置的速度分布與低床層的顆粒類似；但是=0 mm和=150 mm附近，顆粒以一定的速度下落，因為處于中間位置的顆粒向上運動到最大高度時開始從向上作用力較小的兩側(cè)位置回落。高床層顆粒的時均速度較小，表明該位置的顆?？赡苁且堰\動到最大高度，處于逐漸停止向上運動和開始向下運動的狀態(tài)。

圖3 柱形顆粒與球形顆?；旌狭鲃幽M瞬態(tài)圖

(a) 流化前；(b)~(g) 氣泡運動過程

(a) 模擬結(jié)果；(b)實驗結(jié)果

表觀氣速/(m?s?1)：(a) 1.0；(b) 2.0

H/mm：1—5；2—20；3—40；4—60；5—80；6—120；7—140。

圖8所示為混合顆粒在方向上的時均速度分布。顆粒在低床層時沿方向向床層中心位置(坐標75 mm處為流化床中心位置)移動，即在中心顆粒向上運動后，周圍的顆粒向中心位置補充，顆粒在中床層高度(=80~120 mm)沿方向背離中心位置運動，即向兩側(cè)壁面處移動，對照圖7可知該高度的顆粒呈中間上升兩側(cè)下落的狀態(tài)，圖8表明該下落位置的顆粒同時向兩側(cè)的壁面位置移動。高床層位置的顆粒在方向沒有比較明顯的運動。

H/mm：1—5；2—20；3—40；4—60；5—80；6—120；7—140。

結(jié)合圖7和圖8可知，方向中間區(qū)域(=75 mm)的顆粒從低床層位置向上運動，兩側(cè)顆粒向中間區(qū)域補充，中間床層位置的顆粒在方向中間區(qū)域向上運動同時向遠離中間區(qū)域的壁面處運動，在兩側(cè)靠近壁面的區(qū)域向下運動，高床層顆粒呈現(xiàn)到達高點開始回落。床內(nèi)顆粒整體呈現(xiàn)從中部向上到達高點后兩側(cè)回落的運動狀態(tài)。

圖9所示為不同床層高度處的時均空隙率，整體上呈現(xiàn)床高越高，空隙率越大的趨勢，顆粒集中在中低床層。從圖9可見：低床層的空隙率沿方向波動較大，因為低床層顆?？拷鼩饬魅肟冢軞怏w射流的影響顯著。

H/mm：1—5；2—20；3—40；4—60；5—80；6—120；7—140。

5 結(jié)論

1) 通過構(gòu)建CFD-DEM模型，采用球元疊加法和虛擬小球法描述柱形顆粒，引入體積分數(shù)改進曳力系數(shù)公式，實現(xiàn)了柱形顆粒與球形顆粒的混合流動。

2) 模擬結(jié)果與試驗結(jié)果接近，說明該方法能較好地模擬在鼓泡床狀態(tài)下異質(zhì)柱形顆粒與球形顆?；旌狭鲃拥倪^程。

3) 隨著表觀氣速的增大，顆粒上升高度增大，異質(zhì)柱形顆粒與球形顆粒的混合流動整體呈現(xiàn)從中部向上到達高點后兩側(cè)回落的周期運動狀態(tài)，符合實際流化過程中顆粒運動規(guī)律。

[1] TSUJI Y, KAWAGUCHI T, TANAKA T. Discrete particle simulation of two-dimensional fluidized bed[J]. Powder Technology, 1993, 77: 79?87.

[2] KAWAGUCHI T, TANAKA T, TSUJI Y. Numerical simulation of two-dimensional fluidized beds using the discrete element method (comparison between the two- and three-dimensional models)[J]. Powder Technology, 1998, 96: 129?138.

[3] GERA D, GAUTAM M, TSUJI Y, et al. Computer simulation of bubbles in large-particle fluidized beds[J]. Powder Technology, 1998, 98: 38?47.

[4] 桂南. 復(fù)雜兩相流動中顆粒碰撞的DEM-LES/DNS耦合模擬研究[D]. 杭州: 浙江大學(xué)熱能工程研究所, 2010: 87?107. GUI Nan. Simulation study on DEM-LES/DNS coupling of particle collisions in complex two-phase flows[D]. Hangzhou: Zhejiang University. Institute of Thermal Energy Engineering, 2010: 87?107.

[5] 盧洲, 劉雪東, 潘兵. 基于CFD-DEM方法的柱狀顆粒在彎管中輸送過程的數(shù)值模擬[J]. 中國粉體技術(shù), 2011, 17(5): 65?69. LU Zhou, LIU Xuedong, PAN Bing. Numerical simulation of the process of cylindrical particles in elbow pipe based on CFD-DEM method[J]. China Powder Technology, 2011, 17(5): 65?69.

[6] ZHONG Weiqi, ZHANG Yong, JIN Baosheng, et al. Discrete element method simulation of cylinder-shaped particle flow in a gas-solid fluidized bed[J]. Chem Eng Technol, 2009, 32(3): 386?391.

[7] REN Bing, ZHONG Wenqi, CHEN Yu, et al. CFD-DEM simulation of spouting of corn-shaped particles[J]. Particuology, 2012, 10(5): 562?572.

[8] REN Bing, ZHONG Wenqi, JIANG Xiaofeng, et al. Numerical simulation of spouting of cylindroid particles in a spouted bed[J]. Canadian Journal of Chemical Engineering, 2014, 92(5): 928?934.

[9] REN Bing, ZHONG Wenqi, JIN Baosheng, et al. Numerical simulation on the mixing behavior of corn-shaped particles in a spouted bed[J]. Powder Technology, 2013, 234(1): 58?66.

[10] LIU Xuejiao, ZHONG Wenqi, JIANG Xiaofeng, et al. Spouting behaviors of binary mixtures of cylindroid and spherical particles[J]. Aiche J, 2015, 61(1): 58?67.

[11] 張勇. 氣固流化床非球異質(zhì)顆粒介觀混合特性的實驗研究與三維DEM直接數(shù)值模擬[D]. 南京: 東南大學(xué)能源與環(huán)境學(xué)院, 2010: 15?34. ZHANG Yong. Experimental study of mesoscopic mixing characteristics of non spherical particles in gas-solid fluidized bed and direct numerical simulation of three-dimensional DEM[D]. Nanjing: Southeast University. College of Energy and Environment, 2010: 15?34.

[12] YU A B, FENG Y Q. Microdynamic modelling and analysis of the mixing and segregation of binary mixtures of particles in gas fluidization[J]. Chem Eng Sci, 2007, 62(1/2): 256?268.

[13] 任冰. 稠密氣固系統(tǒng)非球形顆粒運動機制的試驗和CFD-DEM耦合模擬研究[D]. 南京: 東南大學(xué)能源與環(huán)境學(xué)院, 2013: 91?121. REN Bing. Experimental study on the mechanism of non spherical particle motion in dense gas solid system and CFD-DEM simulation[D]. Nanjing: Southeast University. College of Energy and Environment, 2013: 91?121.

[14] KODAM M, BHARADWAJ R, CURTIS J, et al. Cylindrical object contact detection for use in discrete element method simulations. Part I: contact detection algorithms[J]. Chem Eng Sci, 2010, 65(22): 5852?5862.

[15] ERGUN S. Fluid flow through packed columns[J]. Chem Eng Prog, 1952, 48(2): 89?94.

[16] WEN C Y, YU Y H. Mechanics of fluidization[J]. Chemical Engineering Progress Symposium Series, 1955, 62: 100?111.

(編輯 趙俊)

Simulation study on characteristics of cylinder-shaped particles and sphere-shaped particles flow

WANG Heng1, ZHONG Zhaoping1, WANG Jia1, WANG Zeyu1, WANG Xiaoyi2, ZHU Lingli1

(1. Key Laboratory of Energy Thermal Conversion and Control of Ministry of Education,School of Energy and Environment, Southeast University, Nanjing 210096, China;2. Engineering Design & Research Institute, Southeast University, Nanjing 210096, China)

A numerical simulation method was used to describe and analyze the mixing flow of biomass particle and quartz sands in gas-solid fluidized bed. Three-dimensional gas-solid flow models were constructed to simulate the mixed flow of heterogeneous particles in mesoscopic scale. The integrated model was built in an Eulerian-Lagrangian approach and the constructed methods of cylinder-shaped particles were different when it came to different numerical methods. Each cylinder-shaped particle was constructed as an agglomerate of fictitious small particles in CFD part, while in DEM method, cylinder-shaped particles were built by multi-sphere method, in which small sphere element merged with each other. Soft sphere model was used to get the connect force between particles. The total connect force of cylinder-shaped particle was calculated as the sum of the small sphere particles’ forces. Two models with different superficial gas velocities (1.0 m/s and 2.0 m/s) were built and the results of simulation were compared with the experimental results. The results show that the present work provides an effective approach to simulation the flow of two component particles.

CFD-DEM simulation; cylinder-shaped particle; gas-solid fluidized bed

10.11817/j.issn.1672?7207.2017.06.034

TQ051

A

1672?7207(2017)06?1667?07

2016?06?01；

2016?09?28

國家自然科學(xué)基金資助項目(51276040，U1361115)(Projects(51276040, U1361115) supported by the National Natural Science Foundation of China)

仲兆平，博士，教授，從事能源與環(huán)境工程研究；E-mail：zzhong@seu.edu.cn