CFD-DEM模型并行化及其在流化床氣固流動中的應用

楊春振，陳成敏，劉光霞，段鈺鋒

?

CFD-DEM模型并行化及其在流化床氣固流動中的應用

楊春振1，陳成敏1，劉光霞1，段鈺鋒2

（1山東省科學院能源研究所，山東 濟南 250014；2東南大學能源與環境學院能源熱轉換及其過程測控教育部重點實驗室，江蘇 南京 210096）

摘要:提出廣義區域覆蓋法用于定位計算節點區域，利用Fluent軟件的UDFs實現顆粒信息在計算節點間高效、準確地傳遞，建立了適用于非正交網格的CFD-DEM模型并行求解方法，在探討并行模型計算效率的基礎上對埋管流化床和單孔射流流化床內的氣泡行為及雙支腿流化床半床間的顆粒交換行為進行了數值模擬研究并對比實驗結果。結果表明:CFD-DEM模型并行求解方法能夠快速傳遞及收集計算節點內的顆粒信息；廣義區域覆蓋法降低了顆粒在計算節點內的定位時間，提高了模型的并行計算效率；能夠捕獲與 Rong等一致的埋管周圍氣泡行為（其頻率約為3 Hz），獲得了與Bokkers等結果接近的單孔射流氣泡尺寸及顆粒速度矢量分布，雙支腿流化床半床間顆粒交換行為與實驗結果一致，第一主頻約為1 Hz，第二主頻大于2 Hz。

關鍵詞:并行；DEM模型；兩相流；流化床；數值模擬

引 言

CFD-DEM（computational fluid dynamic and discrete element method）模型是目前研究稠密氣固兩相流流動、傳熱等的重要工具[1-5]。隨著氣固裝置（如流化床）幾何形狀復雜度的提高，諸多學者提出了不同方法實現DEM模型在非正交網格內的使用。Wu等[6-7]提出了精確計算三角形網格和楔形網格內空隙率的方法；趙永志等[8-9]以正交DEM網格作為背景網格完成了DEM模型與Fluent軟件的耦合，并研究了埋管流化床的傳熱行為；Su等[10]將DEM模型與OpenFOAM耦合，研究了多埋管流化床內的顆粒流動行為；楊旭[11]借鑒光滑粒子動力學中插值的思想和權函數的概念實現了 DEM在非規則網格內的使用；本課題組[12-13]提出二維和三維通用的區域覆蓋法（domain coverage method, DCM），實現了顆粒在非正交網格內的高效定位技術，大大提高了 DEM 模型在非規則網格內的應用。

隨著流化床尺寸的增加或顆粒粒徑的減小，CFD-DEM模型需要跟蹤的顆粒量急劇增加，采用串行程序會增加計算量。為了有效解決這一問題，研究者著力于開發 CFD-DEM模型的并行算法。Tsuji等[14]利用16個CPU研究了三維方形流化床內4.5×106個顆粒的氣泡行為；Xu等[15]開發出多尺度CPU-GPU計算模式，采用CPU計算流場、GPU計算顆粒運動的方式模擬了工業流化床內的氣固流動行為；任立波等[16]借助Fluent軟件的MPI技術實現了CFD-DEM模型的并行耦合算法；Wu等[17]借助Fluent軟件在并行計算方面的優勢發展出一套顆粒在多計算節點間信息傳遞的方法，實現了DEM模型的并行化。

本研究在所提出的DCM方法基礎上提出廣義DCM方法，借助Fluent軟件提供的用戶自定義宏函數（user defined function，UDF）實現了CFD-DEM模型的并行化。首先介紹實現顆粒在多計算節點間信息傳遞的方法，借助廣義DCM實現CFD-DEM模型的并行計算；然后數值模擬研究了單埋管流化床和單孔射流流化床的氣泡特性及雙支腿半床間顆粒交換行為，并與相關學者的實驗研究結果進行對比，驗證了所開發的CFD-DEM模型并行計算方法的準確性。本研究所使用的CFD-DEM數學模型與文獻一致，見文獻[12-13, 18-19]。

1　并行CFD-DEM模型

1.1 計算節點間顆粒信息的傳遞

Fluent軟件具有高效的并行計算能力，其提供的 UDF不僅可以實現復雜數學模型的建模，也可通過 Host與 Node間的信息交互 UDF（如PRF_CSEND_###及PRF_CRECV_###等）實現模型的并行化。

CFD-DEM模型并行計算要求Fluent流體網格分布于至少2個Node內，由于每個DEM時間步長內顆粒位置均發生改變，需將各節點內的顆粒信息發送給Host統計并重新分配至其應屬的Node。因此，提高顆粒信息的再分配效率是提高 CFD-DEM模型并行化速度的重點。

圖1為Fluent流體網格在不同Node內的分布。Node 0包含全局網格的1～4列，其中1～3列為內部網格、4列為外部網格（用于 Node間的信息傳遞）；Node 1包含全局網格的3～7列，其中4～7列為內部網格、3列為外部網格。位于Node內部網格的顆粒與分布在全部網格的顆粒必須是時刻相同的，因此為了獲得全部顆粒信息同時保證不重復，可以將Node內部網格內的顆粒信息傳遞給Host統計，并通過host_to_node_###函數實現Host向Node分配顆粒信息的功能；為了保證Node內顆粒碰撞的準確性和完整性，需檢測內部和外部網格內所有顆粒間是否發生碰撞。

圖1　GDCM基本原理Fig. 1 Schematic diagram of GDCM Mesh splitting line; 　 GDCM grid boundary; 　 Fluent grid boundary;particle; particle locate in GDCM grid

為了實現Host將匯總到的顆粒重新合理、高效地分配至Node，本研究基于DCM提出廣義區域覆蓋法（general DCM，GDCM）:首先建立如圖1所示的GDCM網格（GDCC），分析各Fluent流體網格節點坐標獲得各坐標方向上的最大值和最小值，最終確定各節點的計算區域，Host搜集位于該區域內的顆粒并將其信息傳遞給Node；然后在Node內利用DCM方法定位顆粒在Fluent網格內的位置，用于網格空隙率的計算，DCM可見文獻[12]。DCM關鍵在于建立單個Fluent流體網格和多個DEM網格間的對應關系；GDCM借助DCM覆蓋網格的思路，重點在于確定節點內所有Fluent網格的區域，并形成各計算節點區域的覆蓋網格。

圖 2為 CFD-DEM模型并行化流程圖，采用User defined Init、User defined EXECUTE_AT_END 和User defined Adjust這3個UDF實現。User defined Init用于在 Host中初始化顆粒信息，確定各 Node 的GDCC、建立DEM網格和DCC，確定DEM網格與Fluent網格的映射關系；User defined EXECUTE_AT_END可控制流體和顆粒時間步長（Δtg和Δtp），當Δtp<Δtg時判斷DEM模型的執行次數；User defined Adjust實現DEM模型功能，包括顆粒定位、網格空隙率、顆粒曳力、碰撞檢測、計算接觸力及顆粒跟蹤等。

各Node的顆粒信息經過一個時間步長更新后匯總發送至Host。Host將顆粒信息匯總并存儲，借助GDCC重新確定顆粒所屬的Node，通過消息傳遞UDF將顆粒信息回傳至Node。

圖2　CFD-DEM模型并行算法流程Fig. 2 Flowchart of parallel CFD-DEM algorithm

1.2 并行CFD-DEM模型計算效率

為了定量對比并行算法和串行算法的效率，本研究測試了如圖3(a)所示的0.1 m×0.1 m區域（戴爾:i3-2370M）內顆粒定位所耗時間，其中串行算法耗時為全部顆粒定位于各Fluent流體網格內的時間，并行算法耗時包括Host將匯總的顆粒信息傳遞給Node的時間及顆粒定位于各Node的Fluent網格的時間。顆粒總數分別為1萬、4萬和16萬個。按照圖3(a)所示顆粒排列方案，并行模型2個Node的顆粒數分別為0.52萬、2.08萬和8.32萬個。從圖3(b)所示的顆粒定位耗時可以看出定位顆粒所需時間隨顆粒數增加而增加。定位16萬個顆粒時，2個Node的耗時分別為46 ms和51 ms，串行算法需84 ms，可見CFD-DEM模型的并行算法計算效率更高。

圖3　顆粒數對定位顆粒耗時的影響Fig. 3 Effects of particle number on time-consuming

2　二維單埋管流化床氣泡行為

為了驗證并行化的CFD-DEM模型的準確性，采用串行與2個Node的并行CFD-DEM模型分別計算了單埋管流化床氣泡行為，并與Rong等[18]的結果進行對比。圖4(a)給出了二維單埋管流化床幾何結構，寬和高分別為180mm和420mm，埋管直徑為30 mm，采用非規則網格劃分，網格數為8262個。本研究模擬20000個顆粒，采用重力沉降法獲得顆粒的初始分布，如圖4(c)所示。流化風通過在布風板上設置的20個Φ3 mm孔進入流化床，采用速度入口邊界條件(u=1.1 m·s-1)、充分發展出口邊界條件。

串行模型的DEM網格數為180×210。由于并行計算區域邊緣呈鋸齒狀，并行模型2個Node分別設置90×210個DEM網格。模擬物理時間為5 s。表1列出了模擬參數，DEM模型見文獻[18]。

表1　模擬參數Table 1 Simulation parameters

圖4　二維單埋管流化床Fig. 4 2D fluidized bed with single immersed tube

圖 5給出了采用串行 CFD-DEM模型和并行CFD-DEM模型獲得的埋管右側水平位置為基準上下2°范圍內空隙率隨時間的變化。從圖中可以看出串行模型和并行模型均出現空隙率為1的時間段，這表明此時大氣泡穿過埋管，而且完全包裹住埋管，阻礙顆粒與埋管發生碰撞。在空隙率不為1的時間段（由于采用質心定位法及擬三維法求解空隙率，圖中呈現出空隙率小于 0.25的情形），呈現小氣泡流動、顆粒包裹埋管及兩種并存的3種階段。

為了定量分析氣泡生成、破碎的周期行為，對圖5中的空隙率變化進行功率譜密度（PSD）分析，如圖6所示。串行算法和并行算法獲得的氣泡頻率分別為3.084 Hz和3.203 Hz，與Rong等[18]（u=1.0 m·s-1）計算的2～3 Hz是相符的。

圖5　埋管右側上下2°范圍內空隙率隨時間的變化Fig. 5 Variations of void fraction at ±2° of right-side of tube simulated by serial and parallel CFD-DEM models

3　三維單孔射流流化床氣泡行為

本研究模擬了三維矩形單孔射流流化床內的氣泡行為，并與 Bokkers等[19]的實驗和數值模擬結果進行對比。流化床寬深高分別為0.15 m、0.015 m、0.45 m；采用正交結構化網格劃分，網格數為4050個；噴口寬10 mm，入口速度為20 m·s-1，流化速度為1.2 m·s-1。模擬3萬個粒徑為2.5 mm的顆粒。為了精確計算顆粒碰撞過程，顆粒運動的時間步長設為2×10-6s。采用2個Node。表2列出了模擬參數。本研究采用的CFD-DEM模型與文獻[19]一致。

圖6　對隨時間變化的空隙率的功率譜密度分析Fig. 6 PSD analysis of variations of void fraction

表3　模擬參數Table 2 Simulation parameters

圖7對比了t=150 ms時刻本研究并行算法獲得的氣泡與Bokkers等串行算法的模擬及實驗結果。對比圖7(a)和圖 7(c)可以看出，本研究模擬得到的氣泡較 Bokkers等的模擬結果小（本研究采用 De Felice曳力模型；Bokkers等采用Ergun/Wen-Yu曳力模型）。然而對比圖 7(a)與圖 7(b)所示的實驗結果，本研究模擬結果與實驗氣泡相似。

為了定量描述射流氣泡行為，圖8給出了并行模擬和Bokkers等實驗獲得的氣泡尺寸及其周圍顆粒的速度分布。并行模擬的最大處的氣泡尺寸約為80 mm×150 mm，實驗得到的氣泡尺寸約為 77 mm×120 mm，兩者尺寸在高度方向存在較大的差別，但在寬度方向基本相同。此外，結合顆粒的速度矢量分布可以看出本研究并行算法能夠在細節上較為準確地捕獲氣泡周圍顆粒的運動規律，如氣泡上部周圍顆粒及其速度矢量分布一致性較好。但并行模擬結果在氣泡頂端及氣泡尾渦區域與實驗結果存在一定的差異，這與顆粒的碰撞參數及曳力模型的選擇有關。

圖7　單孔射流氣泡 (t=150 ms)Fig. 7 Bubble in single jet fluidized bed at t=150 ms

4　二維雙支腿流化床顆粒交換行為

雙支腿結構是大型循環流化床鍋爐的一種重要設計方案，然而如果操作不當半床間易發生不穩定的顆粒交換并造成“翻床”現象。本課題組采用Yule-Walker方法對雙支腿流化床（dual leg fluidized bed，DL-FB）半床間的壓力脈動信號進行功率譜分析，研究顆粒交換的主頻及其對應的物理特性[20]，并采用CFD-DEM模型串行求解方法模擬了半床間的顆粒交換行為[13]。為進一步驗證本研究并行求解方法的準確性，采用4個Node模擬顆粒交換行為。

圖9給出了DL-FB的幾何尺寸、網格分布（網格共10240個）以及重力沉降法獲得的顆粒初始分布。雙支腿入口氣體流速為4.16 m·s-1，頂部出口氣體充分發展，床料高80 mm，模擬時間為60 s。表3列出了模擬參數。

圖9　雙支腿流化床Fig. 9 Dual leg fluidized bed

表3　模擬參數Table 3 Simulation parameters

圖8　氣泡尺寸及氣泡周圍顆粒的速度矢量分布Fig.8 Bubble size and velocity vectors of solids

圖10(a)為60 s內DL-FB半床內顆粒數（Np）的變化，可以看出顆粒在半床間發生劇烈的交換，而且呈現出兩種明顯的交換行為:顆粒數變化尖峰（顆粒以團簇的形式交換）、較緩和區域（顆粒以單顆粒的形式交換）并存。這與本課題組[20]的實驗結果一致。為了定量分析兩種行為的頻率及強度，圖10(b)給出了顆粒數變化的PSD圖，可以看出存在兩個主頻，分別為1.152 Hz和2.324 Hz，這與實驗得到的0.938 Hz和2.792 Hz[20]基本一致。

圖10　DL-FB半床間顆粒交換行為及功率譜分析Fig. 10 Solid exchange behavior and PSD diagram in DL-FB

5　結 論

本研究提出廣義區域覆蓋法，借助Fluent軟件的用戶自定義函數實現顆粒信息在多計算節點間的快速、準確傳遞，建立了可適用于非規則網格的CFD-DEM模型并行求解方法。在對比了并行算法與串行算法的計算效率基礎上，為了驗證并行化CFD-DEM模型的計算精度，對比分析了單埋管流化床和單孔射流流化床內的氣泡行為以及雙支腿流化床內的顆粒交換行為，并采用功率譜分析方法定量獲得了信號的脈動頻率。結論如下。

（1）廣義區域覆蓋法能夠快速將顆粒分配至多計算節點內；并行化CFD-DEM模型能夠高效地保證顆粒信息在Host及Node間傳遞。

（2）模型并行算法具有較高的計算效率，雙Node并行計算效率約為串行模型的1.6～1.8倍。

（3）并行化CFD-DEM模型能夠較好地捕獲流化床埋管周圍的氣泡行為，埋管周圍氣泡頻率約為3 Hz；獲得了與Bokkers等實驗結果一致的單孔射流氣泡尺寸及顆粒速度矢量；雙支腿流化床半床間顆粒交換行為存在兩個主頻，數值模擬得到的第一、二頻率分布與實驗結果一致。

符 號 說 明

d——直徑，m

e——顆粒碰撞恢復系數

F——頻率，Hz

k——彈性系數，N·m-1

PN——顆粒數，個

t——時間，s

Δt——時間步長，s

u——速度，m·s-1

α——空隙率

μg, μpp, pw——分別為氣體黏度和顆粒摩擦系數

ρ ——密度，kg·m-3

下角標

g——氣體

n——法向

p——顆粒

t——切向

References

[1] 彭麗, 吳迎亞, 李佳瑤, 等. 基于 DEM 模擬氣固鼓泡床中顆粒碰撞參數對流場間歇性的影響[J]. 化工學報, 2015, 66(6): 2041-2048.

PENG L, WU Y Y, LI J Y, et al. Effect of granular collision parameters on DEM simulation of flow field intermittency in gas-solids bubbling fluidized bed [J]. CIESC Journal, 2015, 66(6): 2041-2048.

[2] 王猛, 朱衛兵, 孫巧群, 等. 提升管內氣固流動特性的離散元模擬[J]. 化工學報, 2013, 64(7): 2436-2445.

WANG M, ZHU W B, SUN Q Q, et al.. Discrete element simulation of gas-solids flow behavior in riser[J]. CIESC Journal, 2013, 64(7): 2436-2445.

[3] 周池樓, 趙永志. 離散單元法及其在流態化領域的應用[J]. 化工學報, 2014, 65(7): 2520-2534.

ZHOU C L, ZHAO Y Z. Discrete element method and its applications in fluidization[J]. CIESC Journal, 2014, 65(7): 2520-2534.

[4] 卜昌盛, 陳曉平, 劉道銀, 等. 基于顆粒尺度的離散顆粒傳熱模型[J]. 化工學報, 2012, 63(3): 698-704.

BU C S, CHEN X P, LIU D Y, et al. Heat transfer model for particles with dicrete element method[J]. CIESC Journal, 2012, 63(3): 698-704.

[5] 吳迎亞, 藍興英, 高金森. 基于DEM模擬的氣固鼓泡床內流場間歇性及顆粒相干結構的分析[J]. 化工學報, 2014, 65(7): 2724-2732.

WU Y Y, LAN X Y, GAO J S. Analysis of flow field intermittencyand coherent structure of particles based on DEM simulation of gas-solids bubbling bed[J]. CIESC Journal, 2014, 65(7): 2724-2732.

[6] WU C L, BERROUK A S, NABDAJYMAR K. Three-dimensional discrete particle model for gas-solid fluidized beds on unstructured mesh[J]. Chemical Engineering Journal, 2009, 152(2/3): 514-529.

[7] WU C L, ZHAN J M, LI Y S, et al. Accurate void fraction calculation for three-dimensional discrete particle model on unstructured mesh[J]. Chemical Engineering Science, 2009, 64(6): 1260-1266.

[8] ZHAO Y Z, JIANG M Q, CHENG Y. Particle-scale simulation of the flow and heat transfer behaviors in fluidized bed with immersed tube[J]. Frontiers of Chemical Engineering in China, 2008, 2(3): 341-345.

[9] 趙永志, 程易. 沉浸管式流化床的顆粒尺度模擬[J]. 化學工程, 2007, 35(11): 21-24.

ZHAO Y Z, CHENG Y. Particle-scale simulation of fluidized bed with immersed tubes[J]. Chemical Engineering (China), 2007, 35(11): 21-24.

[10] SU J W, GU Z L, XU X Y. Discrete element simulation of particle flow in arbitrarily complex geometries[J]. Chemical Engineering Science, 2011, 66(23): 6069-6088.

[11] 楊旭. 微型流化床的數值模擬及混合特性研究[D]. 長春: 吉林大學, 2014.

YANG X. Numerical simulation on the mixing characteristics in micro fluidized bed[D]. Changchun: Jilin University, 2014.

[12] 楊春振, 段鈺鋒, 孫榮峰, 等. 埋管流化床顆粒流動行為的數值模擬[J]. 化工學報, 2013, 64(8): 2788-2793.

YANG C Z, DUAN Y F, SUN R F, et al. Numerical study of solid behavior in a fluidized bed with multi-immersed tubes[J]. CIESC Journal, 2013, 64(8): 2788-2793.

[13] YANG C Z, DUAN Y F. CFD-DEM model for simulating solid exchange in a dual-leg fluidized bed[J]. Chemical Engineering & Technology, 2013, 36(11): 1907–1914.

[14] TSUJI T, YABUMOTO K, TANAKA T. Spontaneous structures in three-dimensional bubbling gas-fluidized bed by parallel DEM–CFD coupling simulation[J]. Powder Technology, 2008, 184(2): 132-140.

[15] XU M, CHEN F G, LIU X H, et al. Discrete particle simulation of gas-solid two-phase flows with multi-scale CPU-GPU hybrid computation[J]. Chemical Engineering Journal, 2012, 207/208: 746-757.

[16] 任立波, 韓吉田. 基于 CFD-DEM 耦合并行算法的錐形噴動床內離散顆粒數值模擬[J]. 東南大學學報(自然科學版), 2014, 44(5): 993-998.

REN L B, HAN J T. Numerical simulation of discrete particles in conical-base spouted bed based on parallel coupled CFD-DEM model[J]. Journal of Southeast University (Natural Science Edition), 2014, 44(5): 993-998.

[17] WU C L, AYENI O, BERROUK A S, et al. Parallel algorithms for CFD-DEM modeling of dense particulate flows[J]. Chemical Engineering Science, 2014, 118: 221-244.

[18] RONG D G, TAKAFUMI M, MASAYUKI H. Particle and bubble movements around tubes immersed in fluidized beds — a numerical study[J]. Chemical Engineering Science, 1999, 54(23): 5737-5754.

[19] BOKKERS G A, VAN SINT ANNALAND M, KUIPERS J A M. Mixing and segregation in a bidisperse gas-solid fluidised bed: a numerical and experimental study[J]. Powder Technology, 2004, 140(3): 176-186.

[20] YANG C Z, DUAN Y F, HU H T, et al. Pressure fluctuation analysis of solid exchange in a dual-leg fluidized bed[J]. Powder Technology, 2012, 224(224): 69-75.

2016-01-07收到初稿，2016-03-14收到修改稿。

聯系人:陳成敏。第一作者:楊春振（1984—），男，博士。

Received date: 2016-01-07.

中圖分類號:TK 1

文獻標志碼:A

文章編號:0438—1157（2016）07—2748—08

DOI:10.11949/j.issn.0438-1157.20160026

基金項目:山東省自主創新及成果轉化專項（2014ZZCX04215）；山東省科學院青年基金項目（2014QN016，2013QN016）。

Corresponding author:CHEN Chengmin, chencm@sderi.cn supported by the Project of Independent Innovation and Achievement Transformation of Shandong Province (2014ZZCX04215) and the Youth Fund of Shandong Academy of Sciences (2014QN016, 2013QN016).

Parallel algorithm of CFD-DEM model and applications on gas-solids flow in fluidized beds

YANG Chunzhen1, CHEN Chengmin1, LIU Guangxia1, DUAN Yufeng2
(1Energy Research Institute of Shandong Academy of Sciences, Jinan 250014, Shandong, China;2Key Laboratory of Energy Thermal Conversion and Control, Ministry of Education, School of Energy and Environment, Southeast University, Nanjing 210096, Jiangsu, China)

Abstract:General domain coverage method (GDCM) was proposed to determine the region for each computing node efficiently by comparing all cell nodes, where these regions were helpful for Host to identify which computing node solid particle belongs to. Particle information was passed among host and nodes with the help of Fluent UDFs. The parallel algorithm for CFD-DEM model was developed successfully to simulate gas-solids flows in the fluidized beds (FBs) which were meshed in arbitrary shapes. Based on the comparison of calculation efficiency between parallel algorithm and serial algorithm, the bubble behaviors were studied in a FB with an immersed tube in a single jet FB, and the solids exchange behaviors between two half beds were simulated and compared with experimental results in a dual leg FB. The simulation results showed that the particle information was passed between host and nodes efficiently by parallel algorithm. It was efficient for GDCM to determine which computing node particles are in. The simulated bubbling frequency (about 3 Hz) was similar as Rong’s simulation result. Compared with Bokkers’ simulated and experimental jetting bubble size, the more similar bubble size and solid velocity vectors around the bubble were obtained by the parallel algorithm. The solidexchange between two half beds in both forms of the first (～1 Hz) and second (>2 Hz) domain frequency was simulated, which was similar as the experimental result.

Key words:parallel; DEM model; two-phase flow; fluidized bed; numerical simulation