本期闡述MI-PIC模型的數值實現，重點介紹顆粒相運動求解過程，對固相應力模型進行調整的固相速度修正模型以及氣固間插值和平均方法。

歐拉-拉格朗日方法中，顆粒相和連續相在計算尺度和描述方式的不同導致不能用與單一流體相同的方式進行求解，而是在各自的相空間分別求解。

MP-PIC方法氣固場求解過程的四個階段（圖1）

a	顆粒相信息集中 - 平均過程；
b	流體場和顆粒相連續場求解；
c	流體信息和顆粒連續信息插值到離散顆粒位置；
d	顆粒速度和位置更新。

圖1 顆粒相求解的四個階段

從圖1可以看到：由于流體相的信息存儲在計算網格體心，而顆粒相的計算粒子隨時間變化在空間不斷更換位置，為了在拉格朗日相空間和歐拉網格間傳遞顆粒特性來實現相間耦合，必須使用插值方法來完成此項操作。

通常插值權重函數為三個方向上插值算子的乘積，且這三個算子需要滿足相互正交和有界。如x方向的算子[1]：

• 方程1

式中，x_ξ表示網格節點位置，x_p表示計算顆粒位置。

一般來說，任何滿足此條件的函數都可以用來計算權重，當使用笛卡爾網格-切割體網格，為了減少計算復雜度，通常采用三線性插值函數。例如，網格上的固相體積分數可以通過如下插值操作得到，

• 方程2

其中，V_p,i 是顆粒體積，n_p,i 是數值粒子中包含的顆粒數，N_p是當前網格ξ內的數值粒子數量，V_ξ為網格體積。此外，在更新顆粒速度時需要計算顆粒應力的梯度，對于梯度的計算如下式所示：

• 方程3

式中，N_ξ 表示數值粒子周圍的網格，即粒子位置的場梯度為周圍網格場梯度的插值的和。插值函數的乘積和梯度的詳細定義和運算可以參考文獻[2]。

在MP-PIC中，將顆粒相信息分配到每個顆粒，讓它們單獨攜帶并在一個時間步內以彼此獨立的特定屬性運動。當進入下一個時間步后，先將網格內屬于此網格或影響此網格的所有粒子信息收集，在歐拉網格上重新構建出連續的顆粒相場以便計算顆粒間作用力。因而可以將顆粒相既當做連續相又可以當做離散相。

圖2為MP-PIC求解流程圖，其中（a）為氣固相求解，（b）為固相求解。

圖2 MP-PIC求解流程圖

新時間層的顆粒位置可以通過對方程積分得到更新：

• 方程4

顆粒速度則通過對加速度方程積分得到，對所有項都采用新時間層的量后有隱式形式

• 方程5

其中：

為顆粒位置的插值隱式速度；

為顆粒位置的插值隱式氣相壓力梯度；

為顆粒位置的插值隱式固相壓力梯度。

、和為從歐拉網格中心點上計算得到的平均場量在數值粒子所處位置的插值。

圖3為考慮了顆粒阻尼耗散和各項同性松弛對顆粒運動影響的示意圖。在這其中顆粒應力模型通常是一個分段函數，在顆粒體積分數接近最大堆積濃度時，會得到間斷解，容易導致數值不穩定。

圖3 數值顆粒運動過程

由于固相應力的計算相對于其他幾項對顆粒濃度的變化極為敏感，為了維持求解過程的穩定性，速度方程中的顆粒應力項的計算需要特殊處理。此外，當顆粒應力占主導時，由其計算的修正速度非常大而與實際不符，因而不能直接在方程5中使用。通常方程5中顆粒速度的計算被分割成兩個部分，即速度為固相應力和其他所有作用力時間積分的總和，這樣能為計算帶來更大的穩定性。

對于某一特定粒子，

• 方程6

因此，排除了固相應力后計算的中間速度為：

• 方程7

從連續的顆粒固相應力梯度計算的顆粒速度為：

• 方程8

為保障計算的穩定性，需要忽略分母中的隱式項，因而重新寫為：

• 方程9

1

固相速度修正

通過對固相應力模型的分析（詳見：MP-PIC顆粒間作用力：固相應力）可以知道，顆粒應力模型計算出的速度將遠大于顆粒運動速度，如果直接用于顆粒運動，將會違背守恒原理。而實際過程中，顆粒接近大量堆積顆粒時是逐漸減速后反射的。因此，應力模型的作用應該是防止顆粒進入已經處于堆積狀態的顆粒群而出現過度堆積。同時，應力模型給出了顆粒所受應力的方向，即加速度方向。通過加速度方向以及顆粒的當前運動速度以及恢復系數等參數就能夠對顆粒速度進行修正。

根據對濃度計算方式的不同，速度修正方法分為顯式和隱式方法。

1.1 顯式形式

在更新顆粒速度的計算中，為了保持計算的穩定性，減小顆粒應力項對顆粒體積分數的敏感性，將顆粒速度的計算分割成了兩部分。通過對前一部分的計算，可以對速度進行更新，并通過此中間速度得到顆粒云的中間位置

• 方程10

根據得到的中間顆粒位置，可以計算出中間顆粒濃度，并用此濃度近似下一時刻濃度從而計算顆粒應力和移動顆粒，并最終得到新的顆粒濃度分布

• 方程11

由前面分析可知，由于顆粒應力僅在顆粒接近堆積濃度時才變得重要，因而對于絕大部分區域，中間顆粒濃度基本上等于下一時間步的顆粒濃度。并在此實現中將顯式地用于顆粒應力的計算進而計算修正速度。為了更為直觀的介紹速度修正過程，假設顆粒僅在垂直方向運動，如圖4所示。

圖4 顯式速度修正

圖4展示了單個粒子相對于堆積區的三種典型的運動狀態：

A為顆粒處于堆積區內部，運動受到周圍粒子或壁面的限制，運動速度基本與堆積區顆粒的整體速度一致；

B為顆粒朝堆積區運動，將受到顆粒應力作用而逐漸減速，甚至會反彈；

C為顆粒遠離堆積區，不再受顆粒應力作用。

方程12表示當應力模型計算的速度朝上時（僅在垂直方向），用于限制顆粒速度的速度修正方程。圖4右側的四種情況則更為詳細地描述了如何確定修正速度，具體的實現過程可以參考Garg[3]。

• 方程12

1.2 隱式形式

隱式形式中不再用中間濃度來計算顆粒應力，而是直接計算出下一時間步的顆粒濃度。類似于雙流體模型，對固相建立連續方程或者用顆粒控制方程得到：

• 方程13

將第二項用中間時間層的量近似，并分裂成兩項：

• 方程14

之后認為第二項滿足中間時間層的連續方程，并將第三項中的速度用速度應力方程代替：

• 方程15

最后對最后一項用鏈式法則，得到關于固相濃度的方程：

• 方程16
  

求得新時間層的固相濃度后就可以計算出單元顆粒速度。

2

插值和平均方法

相比于通常的歐拉-拉格朗日方法，MP-PIC方法顆粒相與流體網格結合得更為緊密，需要頻繁地在兩者之間交換信息。精確而高效的插值和平均算法是保證MP-PIC計算魯棒性和一致性的前提，特別是在非結構網格中，由于丟失了與坐標軸線相關的拓撲信息，插值和平均算法需要更具有通用性。存儲流體信息的歐拉網格是按一定規則組織的離散的空間點集，插值計算即意味著基于這些空間點的坐標和場值，以及假定的場變化來構建連續的場分布。對于結構化且網格線與坐標軸線共線的網格，有許多高階、有效的插值格式，這些高精度插值格式通過在三個維度上進行重復操作即可實現將數據從歐拉網格傳遞到空間中任意點。在實際的氣固流動模擬中，絕大多數的設備裝置都具有復雜的幾何結構。雖然也可以通過臺階網格（Step Grid）來近似邊界，但對計算域仍會引入較大的誤差。

相對來說，非結構網格對流場幾何邊界的保真程度相當高，能很好地用于流體相的數值模擬。然而，當涉及到氣-固兩相間插值操作時則會對計算過程以及程序開發帶來了很大的困難，很多在結構網格上充分發展，相當成熟的高精度方法不能在非結構網格上使用。此外，在非結構網格上選擇插值方法時，有幾點要求需要確保，如全面覆蓋流場、插值量在空間中連續變化、具有能恢復構造離散點值的能力、僅與幾何相關以及強連通性。由于涉及到空間插值，就會伴隨有空間分解，將復雜的多面體網格單元分解為基礎的四面體能很方便的進行插值計算。雖然也可以在六面體等復雜的單元內插值，但無疑會極大地增加計算難度。本文涉及的非結構網格主要針對多面體網格（Polyhedronmesh）。在接下來的內容中將介紹單元分解方法以及相關的四面體內插值算法。

2.1 插值方法

• Cell-Vertex Method
  

Cell-Vertex方法主要基于將每個多面體網格單元分解成多個四面體單元，并在四面體單元內進行插值計算。

分解策略如下：

1

對于給定的多面體網格單元，初步分解成多個小多面體，每個小多面體由單元中心（Cell Centre）和每個面的所有頂點組成；　

2

對于每個小多面體，基于原多面體的面再次進行分解，根據使多面體質量最優的原則，在面的多個頂點中選擇一個點作為分解基點，將面上與基點不相鄰點間連線，從而分割原來的面，最后用組成分割成的三角面的三個頂點與單元中心點組成四面體。　

如圖5所示，分解后得到的一個單元體由面基礎點Base，構成一條邊的A、B點以及網格中心點C構成。

圖5 計算網格分解為四面體單元

將計算網格分解成由四面體組成的子網格系統后，此網格系統由單元中心點和各個單元頂點組成，其中單元中心點存儲了所有的流體信息（同位網格）。由于插值過程在每個四面體內進行，因此在插值前需要知道單元頂點的流場值。這一步由單元中心到單元頂點的平均過程完成。單元頂點處的值通過對環繞此頂點的所有單元值進行平均，由于假設流場量在空間中線性變化，而各單元中心距此頂點距離不同，因而其影響也不盡相同，需要通過權重來表示。在這里使用中心距頂點的距離的倒數作為權重。計算表達式為：

• 方程17
  

• 方程18

其中Φ_v 和Φ_i 分別表示單元頂點和第個單元中心的場量，w_i 為歸一化的權重值，r_i 為從單元中心到頂點的距離。采用距離的倒數作為權重計算是眾多計算方法中一種應用較廣且比較簡單的方式，此外，還可以通過拉格朗日乘子法（Lagrangemultipliers）來計算權重。

• 權重計算

在將整個計算域用上述任意一種方法分解后，會得到眾多的四面體單元，這些單元相互不重疊，且完全充滿整個空間。計算域中任意一點必然會位于其中某一四面體內或者單元中心、單元頂點甚至單元面上。這樣，由于四面體結構的固定性，任意一個特定點上的流場值都可以通過四面體單元四個頂點上的值插值得到。這其中涉及到四面體內權重的計算，在這里使用基于重心的權重算法。

假設流場在空間連續線性地變化，因此就可以根據空間點的位置計算權重，這些權重值在顆粒位置的下一次更新前可以一直使用。考慮一個由四個頂點r1，r2，r3，r4組成的四面體，其內部任意點r都可以用四個頂點的組合表示，即

• 方程19

通過將r用分量形式表示，加上λ的和為1的條件，可以建立關于λ的代數方程組：

• 方程20

其中T 為3×3的矩陣，其具體的表達式為[4]：

• 方程21

計算得到的作為當前點r的權重值。

2.2 平均方法

通過積分牛頓運動方程后，顆粒移動到新位置，因而使顆粒云的空間分布得到更新。從連續顆粒應力模型計算速度的方程可知，計算顆粒中間速度需要使用顆粒群整體信息的平均值，即拉格朗日場的平均表述形式。其位于氣相流場的網格節點上，包括顆粒濃度、顆粒應力梯度和顆粒整體的平均速度以及這些場的梯度場。而粒子是離散的，因此需要將這些離散的信息收集起來建立連續的場。平均方法主要有單元平均，雙重網格平均以及矩平均[5]，這里主要介紹雙重網格平均。

雙重網格平均在單元分解的基礎上判定顆粒歸屬的四面體，對于單元中心點，顆粒將單元內所有顆粒信息都映射到單元中心上，而對于單元頂點則需要建立雙重網格，在整個映射過程中，仍然使用基于四面體重心計算的權重系數。

• 方程22

• 方程23

其中Φ為需要重構的平均顆粒場，N為與中心或頂點相關的顆粒數，V_tet為顆粒所在四面體的體積，除以4是因為雙重網格的原因。

由于假設構建的場在四面體內線性變化，因此，在計算過程中可以認為梯度在四面體內為常量，梯度計算采用有限元中的一階梯度估算，即線性梯度重構。假設一個偏離四面體頂點x₁的任意點x₁+h，其值可以通過泰勒展開式近似為[6]：

• 方程24

因而，對其他三個頂點都用此式近似，得到如下3×3系統：

• 方程25

從而可以計算出當前四面體的顆粒平均場梯度▽Φ。

來源：多相流在線