濃相氣固流動模擬方法 

目前，濃相氣固流動模擬方法根據對顆粒離散相的描述方式主要分為兩類：

• 歐拉-歐拉（Eulerian-Eulerian）方法；

• 歐拉-拉格朗日（Eulerian-Lagrangian）方法。

歐拉-歐拉方法主要以采用顆粒動力學模型（KineticTheory of Granular Flows-KTGF）作為固相封閉模型的雙流體模型（TFM）為代表，其將氣相和固相作為相互滲透的連續相，因而在計算規模和計算精度上都得到一定程度的保證。

歐拉-拉格朗日方法將顆粒作為離散相描述，由于描述尺度的不同以及顆粒作用間的模型不同，歐拉-拉格朗日細分出眾多方法。其中最具代表性的是計算流體力學耦合離散元模型的方法CFD-DEM，將氣相用連續相描述，而對每個固相顆粒進行追蹤并詳細計算顆粒間的碰撞過程。雙流體模型可以模擬較大規模的裝置，但不能對顆粒尺度的流動進行解析。離散元模型能準確模擬顆粒間作用力，適用于返料系統中各種復雜流型，但計算非常耗資源，能模擬的顆粒數較少，無法做大尺度的氣固流動模擬，限制了其在工業領域的應用。近年來，兼顧了雙流體模型和離散元模型優點的MP-PIC方法得到越來越多的應用。

1 MP-PIC

多相質點網格模型MP-PIC（Multiphase Particle-in-cell）作為混合型歐拉-拉格朗日方法，氣相的流動采用體積平均N-S控制方程組描述，顆粒相用數值粒子進行離散，對每個粒子都進行追蹤從而獲得顆粒相在不同時刻的空間分布，在需要知道顆粒相整體特性時則將所有顆粒屬性都插值到流體網格進行求解，從而使顆粒相同時具有離散相和連續相的性質。

MP-PIC特點

1	不同于CFD-DEM，MP-PIC使用數值粒子（Parcel）概念表示顆粒相，并使用連續的顆粒應力模型阻止顆粒被過度堆積而超出允許的物理極限；
2	MP-PIC方法將顆粒相群體的特性映射到歐拉網格來計算連續顆粒應力場，而每個數值顆粒都具有獨立的速度和位置，并通過積分牛頓運動方程不斷更新其速度和位置，從而使其同時具有了TFM和CFD-DEM兩者的特點；
3	MP-PIC能計算大規模顆粒系統氣固流動，同時還能對固相在顆粒尺度進行描述。

圖1對不同的濃相氣固流動模型進行了對比，可以看到：DEM模型具有更高的解析度，但耗時最多；TFM計算耗時次之；MP-PIC雖能計算大規模顆粒系統氣固流動，同時能對固相在顆粒尺度進行描述，但由于使用的子模型最多，因而計算結果的準確性最差。在選擇模型，一定要根據自身應用的需求正確選擇。

圖1 濃相氣固多相流模型

2 氣相控制方程

氣相流動用基于體積平均的連續性和動量方程描述為：

• 方程（1）

• 方程（2）

方程（2）中氣相應力可以用方程（3）計算，其中氣相有效動力粘度包括了層流粘度和湍流粘度：

• 方程（3）

氣相壓力和密度間的關系可通過理想氣體狀態方程表示：

• 方程（4）

流場中任意點的氣相和固相的體積分數滿足：

• 方程（5）

3 固相控制方程

顆粒相在空間中分布用概率分布函數（PDF）描述f (x,u_g,ρ_p ,Vp,t )，其表示t時刻在位置x到具有速度為u_g，密度為ρ_p 和體積為Vp的顆粒的概率，f 隨時間t 的演變規律[1]表示為：

• 方程（5）

在實際的數值求解過程中，顆粒相的空間分布通過將概率分布分解為有限數量的數值粒子，并追蹤這些粒子的運動過程來獲得，而每個數值粒子表示一定數量的具有相同密度、體積、速度和位置的真實顆粒。由于使用數值粒子，從而極大地減少了計算需求，同時在大多數條件下不會降低計算結果的準確度。此外，通過對單個粒子賦予屬性從而可以使MP-PIC方法很自然地模擬具有粒徑分布的顆粒及多類型多粒徑的顆粒。

圖2 數值顆粒示意圖

通過顆粒的概率分布函數f，可以很方便地對其積分而得到顆粒的體積分數分布：

• 方程（7）：

顆粒相的其它整體屬性如平均顆粒速度和平均密度也可以通過相同的方式求得。

在每個時間步，每個計算粒子的位置和速度更新基于牛頓運動方程計算得到：

• 方程（8）：

• 方程（9）：

類似于雙流體模型，MP-PIC中的碰撞力也可以分為粘性應力和摩擦應力。粘性應力表示稀相時的二元碰撞而摩擦應力表示高濃度下的長時間接觸力。MP-PIC中只使用摩擦應力主對角元素的梯度，其作用是為了防止顆粒過度堆積而超出顆粒體積濃度所允許的物理極限值。

為了正確表征顆粒碰撞過程，可以在顆粒輸運方程中加入一個BGK項來表示碰撞作用對顆粒波動速度的阻尼作用。MP-PIC通過在顆粒分布函數的演變方程中加入了阻尼模型和各向同性模型來表示顆粒的稀相碰撞作用，這兩項通過對顆粒速度進行修正從而更加準確地描述顆粒間相互作用過程。顆粒連續方程右邊兩項分別表示碰撞趨于各向同性力[2]（Return-to-isotropy force）和碰撞阻尼力[3]（Damping force）對顆粒分布的影響。前者表示各向異性速度分布趨向與高斯各向同性分布的弛豫過程，后者則表示相對運動顆粒間的阻尼過程，即顆粒速度趨向于質量平均速度的弛豫過程。

• 方程（10）：

由于在濃相氣固流中，顆粒碰撞極為頻繁，因此采用一個考慮顆粒恢復系數和速度非平衡分布對阻尼時間作用的改進模型。最終的模型是一個平衡模型，且關鍵在于阻尼時間的求解。

時間尺度模型

• 方程（11）：

其中：

為平均碰撞頻率。

因而碰撞對顆粒速度的作用，

• 方程（12）：

其中網格平均顆粒速度在顆粒點位置的插值。

各向同性模型假設顆粒分布函數收斂到最終分布的時間正比于兩個弛豫時間的倒數。

• 方程（13）

• 方程（14）

不同于阻尼模型將弛豫時間直接用于速度修正，各向同性模型對速度修正分兩步進行。

01 速度取樣

對顆粒云中的每個粒子生成一個隨機數

如果

• 方程（15）：

否則，

• 方程（16）：

或

02 速度修正

由于只對有限數量粒子取樣從而會導致網格內平均速度和速度變量（顆粒溫度）不守恒。需要進行守恒修正：

• 方程（17）：

或

除了顆粒間的法向作用力，在一些非流化的區域，顆粒相互擠壓并存在相對運動。此時，顆粒間還存在著切向摩擦力，基于網格的平均靜和動摩擦可以用以下方程表示，

• 方程（18）：

式中c_s為摩擦系數，通常由實驗來確定。

本文為MP-PIC模型簡介，后續文章將詳細介紹 固相應力模型和MP-PIC的數值實現、 MP-PIC在OpenFOAM中的具體實現。敬請關注！

【參考文獻】

[1] Snider, D.M., An Incompressible Three-Dimensional Multiphase Particle-in-Cell Model for Dense Particle Flows. Journal of Computational Physics, 2001. 170(2): p. 523-549.

[2] O'Rourke, P.J. and D.M. Snider, Inclusion of collisional return-to-isotropy in the MP-PIC method. Chemical Engineering Science, 2012. 80: p. 39-54.

[3] O’Rourke, P.J. and D.M. Snider, An improved collision damping time for MP-PIC calculations of dense particle flows with applications to polydisperse sedimenting beds and colliding particle jets. Chemical Engineering Science, 2010. 65(22): p. 6014-6028.

來源：多相流在線