CFD-DEM  

CFD-DEM模擬是一種基于歐拉-拉格朗日參考系的方法，是離散模擬的典型代表(Ge et al., 2017)。該方法的特征是：

1 將流體視為連續介質，對計算空間進行網格劃分，然后求解納維斯托克斯方程得到每個網格上流體的速度、密度、和壓力等；

2 對顆粒相則接計算其受力，然后根據牛頓定律計算其加速度、速度和位置的變化。

流體相的求解一般采用有限體積等傳統計算流體力學方法，相對成熟，并不是制約該類方法的難點。顆粒相的求解，雖然在算法上較為簡單，但巨大的顆粒數目導致計算量巨大，往往限制了該類方法的應用。如何準確計算流體體積分數，如何計算曳力，如何在非規則網格中計算顆粒所在網格等兩相耦合計算也是該方法的難點。

本文簡要介紹CFD-DEM模型方程、算法流程、研究熱點以及典型應用。

1 研究模型

目前在CFD-DEM學術論文中最常見的控制方程都基于A類模型，即顆粒受力直接包含了壓力梯度力，而B類模型只考慮了顆粒排開流體體積所帶來的阿基米德浮力。這導致B類模型無法正確的模擬旋風分離器等旋轉流場(Zhouet al., 2010)。在A類模型中，Ishii形式的動量方程應用較多，其具體形式如下：

• 方程1

• 方程2

其中為流體體積分數，為其密度，u_f 為其速度，為應力，f_drag表示流體對顆粒的曳力。應力的計算方式如下：

• 方程3
  

其中和分別表示動力粘度和第二粘性系數。曳力的計算方式如下：

• 方程4

該公式的含義為，先插值計算顆粒質心處的流體速度，然后計算顆粒受到的曳力，最后根據牛頓第三定律插值計算CFD網格內顆粒對流體的反作用力。

顆粒的運動包括平動和轉動，其平動控制方程如下：

• 方程5
  

公式右側各項分別為重力、壓力梯度力、曳力和碰撞力。

顆粒轉動的控制方程如下：

• 方程6
  

其中，I_p為轉動慣量，其數值為 m_pd_p²/10. W為轉動速度，T為扭矩。

注：碰撞力和扭矩的具體計算方法，可以參考DEM專業書籍或文獻。

2 算法流程

CFD-DEM的計算流程主要包括：初始化、流體相計算和顆粒相計算。

初始化

初始化操作主要包括讀入計算參數，申請內存空間，初始化流場和顆粒，構建顆粒鄰居列表，構建顆粒流體網格映射關系。在初始化顆粒位置時，要保證顆粒-顆粒以及顆粒-壁面無重疊。

流體相計算

流體相的計算采用壓力耦合方程組的半隱式方法（SIMPLE）即根據假定的速度場與壓力場對動量方程離散求解，然后求解壓力修正方程，對壓力和速度進行修正，迭代直到流場收斂。該計算流程和單相流的計算一致，流體網格中的流體體積分數和平均顆粒速度由顆粒位置信息顯式插值計算得到。在進行流體計算時，顆粒的速度、位置等保持不變。

顆粒相計算

當流場收斂后，進行顆粒相的計算，主要包括以下8個子步驟：

1 顆粒-顆粒、顆粒-壁面碰撞力計算；

2 計算流體壓力梯度力和曳力；

3 更新顆粒速度和位置；

4 多進程并行計算時，傳遞進程邊界顆粒信息；

5 更新顆粒鄰居列表；

6 更新顆粒-流體網格映射關系；

7 更新顆粒-流體網格插值系數；

8 更新流體網格中流體體積分數。

一般情況下，由于流體的計算時間步長大于顆粒計算時間步長，以上1-8步需要重復多次，但在每一步的計算中，流體的速度、密度和壓力分布保持不變。當顆粒和流體的時間同步后，停止顆粒的計算轉為流體的計算。重復以上步驟，直到達到預先設定的模擬時間。

3  研究熱點

CFD-DEM研究熱點集中在提高其計算速度和精度。

CFD計算速度的提高主要依賴于高效的線性方程組求解器和合適的流體方程求解算法。

	目前常用的線性方程求解器為代數多重網格和共軛梯度等迭代方法，以減少求解離散方程的時間。
	流體方程求解算法主要包括SIMPLE和PISO。PISO算法一般采用較大的松弛因子或不進行松弛，因此其壓力-速度耦合迭代修正次數小于SIMPLE算法中的迭代次數，該部分的研究一般可直接采用單相流領域的進展，由于CFD計算在CFD-DEM中所占比重較小，多相流研究人員很少直接涉及該領域。

DEM計算速度的提高主要依賴于高效的碰撞搜索算法。

目前應用最廣泛的方法是鄰居列表法：將顆粒映射到搜索網格中，計算顆粒與相鄰網格中其他顆粒的距離并建立鄰居列表。在碰撞檢測時只需檢測其鄰居列表中的顆粒。鄰居顆粒的截斷距離越大，則鄰居顆粒數目越多，可以減少鄰居列表更新的頻率；反之則需要頻繁更新。對于多粒徑系統，采用均勻分布的單一搜索網格時，由于搜索網格需大于系統中最大的顆粒，這導致網格內顆粒數目巨大，從而降低了搜索效率。因此，多采用樹形或多層嵌套搜索網格。除了改進算法外，還可以采用GPU加速(Xu et al., 2011)或者粗粒化方法(Lu et al., 2016)直接提高DEM的計算速度。

CFD-DEM模擬的計算結果受CFD網格質量、離散格式、顆粒碰撞參數、時間步長、曳力模型和耦合方法等多種因素的影響。目前研究熱點集中在曳力模型、耦合方法以及顆粒形狀。

	CFD-DEM模擬可以采用DNS曳力、EMMS曳力或基于試驗數據的曳力。相比于雙流體模擬，CFD-DEM模擬可以獲得每個網格內顆粒的分布信息，因此該信息也可被用來對曳力進行修正。
	耦合方法包括顆粒信息向CFD網格映射和流場信息向顆粒映射。前者用來計算CFD網格內流體體積分數和平均顆粒速度，后者用來計算顆粒質心處的流速和壓力梯度。對于笛卡爾網格，可以直接計算顆粒所在的CFD網格，并采用線性插值的方法計算。對于其他類型的CFD網格，需要根據CFD網格的拓撲結構和顆粒運動速度判定顆粒所在的網格，然后在進行插值計算。此外，還可以在原有CFD網格的基礎上重建一套笛卡爾網格，并在這套網格上進行兩相的耦合。
	大多數CFD-DEM模擬都將顆粒假定為球形，但實際中常常會有圓柱形、膠囊形、方形等規則形狀以及各種非規則形狀。非球形顆粒直接的碰撞、與流體的作用計算都是研究的熱點(Zhong et al., 2016)。

4  典型應用

CFD-DEM常常用來模擬氣固或液固多相流。由于直接跟蹤顆粒的運動，該方法可以方便的用來研究：

顆粒在反應器中的停留時間分布(Lu et al., 2017)；

多粒徑系統的混合和離析(圖1，Lu et al., 2018)；

顆粒之間的靜電力(Korevaar et al., 2014)、粘性力(Zhou et al., 2018)。

由于直接計算顆粒-顆粒，顆粒-壁面的碰撞力，該方法還可以研究顆粒和反應器的磨損(圖2，Christopher et al., 2015)。

采用并行計算以及粗粒化模型提高計算速度后，該方法還可以用來研究煤氣化等(圖3，Hu et al., 2018)。

圖1. 不同粒徑的FCC 顆粒在流化床中的分布，在表觀氣速為0.05m/s的鼓泡床中模擬0.2s后發現粒徑較大的FCC顆粒分布在底部，較小的顆粒分布在床層上部。圖片摘自Elsevier (Lu et al., 2018)

圖2. 由顆粒- 壁面碰撞引起的彎管壁面磨損。圖片摘自Elsevier (Christopher et al., 2015)

圖3. 流化床煤氣化反應器中煤炭顆粒（紅色）和沙子（藍色）的分布。圖片摘自Elsevier(Hu et al., 2018)

本文是面向初步接觸CFD-DEM的讀者所寫的簡介類文章，并不是完整的模型綜述，因此未包含許多該領域的重要文獻，建議對該模型有興趣的讀者自行查閱相關文獻，以便對該方法有更為完整的認識，后續會向大家介紹詳細的CFD-DEM模擬。

作者： 顧龍夏  

來源：多相流在線