一種基于化學勢LBM多相流模型的改進方法

,,,

(廣西師范大學廣西多源信息挖掘與安全重點實驗室，廣西桂林541004)

晶格Boltzmann方法(lattice Boltzmann method，LBM)由于物理背景清晰，算法簡單，計算效率高和適合大規模并行計算等諸多優點，已在多相流領域得到眾多學者們的關注和應用[1-5]。目前，比較流行的LBM多相流模型主要有：①1991年，Gunstensen等[6]提出的顏色梯度模型；②1993年，Shan等[7]提出的直接刻畫粒子間相互作用的偽勢模型；③1995年，Swift等[8-9]基于自由能理論提出的自由能模型。其中偽勢模型和自由能模型應用最為廣泛。

偽勢模型，又叫Shan-Chen模型，是通過一種假設的相互作用勢來刻畫流體微觀粒子之間的相互作用。該模型雖然算法簡單且容易捕獲和跟蹤相界面，但也有一些局限性：一是只有當有效密度滿足φ(x)=φ0exp(-ρ/ρ0)(其中ρ表示密度，ρ0和φ0表示參考密度和質量)時，該模型才符合熱力學一致性。但在實際應用中，往往有效密度并不取這一特殊形式；二是只適用于密度比較小的多相流系統；三是當模擬達到平衡時，相界面處存在較大的虛速度。多年來，眾多學者對其進行了研究和改進：2003年，Zhang等[10]采用勢函數的形式計算粒子間的相互作用力，使得偽勢模型的穩定性得以提高；2006年，Yuan等[11]通過將各類狀態方程引入模型，擴展了原始偽勢模型的應用范圍，實現了更高的密度比，但卻依然存在熱力學不一致的缺點；2009年，Kupershtokh等[12-13]采用折合變量的思想，提出了一種以加權的方式將有效密度形式和勢函數形式結合起來計算粒子之間的相互作用的模型，當取合適的權值時，模型的穩定性得到大大提高；2013年，Hu等[14]通過在狀態方程前添加一個可調的參數，能夠很好地提高模擬的穩定性；2018年，Qin等[15]引入高階差分來計算非理想流體的相互作用勢，有效地提高了模型的精確度，使得計算結果更加可靠。

基于熱力學自由能理論，1995年，Swift等[8-9]將壓力張量耦合到平衡分布函數，提出了自由能模型。雖然在理論上具有熱力學一致性，但是原始的自由能模型并不滿足伽利略不變性這一物理基本規律[9，16]，構建的新平衡態分布函數也較為復雜。針對這些問題，2000年，Inamuro等[16]通過調節平衡態分布函數中各項系數的方式來減少模型不具有伽利略不變性所帶來的影響；2015年，Wen等[17]通過分析偽勢模型和自由能模型各自的缺陷和特點，將基于熱力學計算得到的非理想力通過偽勢模型的力項技術合并到晶格Boltzmann方程中，提出了同時滿足熱力學一致性和伽利略不變性的多相流模型；2017年，Wen等[18-19]又從化學勢(chemical potential，CP)的角度考慮，構建了一種基于化學勢的晶格Boltzmann方法多相流模型，該模型相比于基于壓力張量的模型[17]避免了計算具有各向異性的壓力張量及其散度，只需要計算化學勢的梯度，從而大大提高了計算效率和模型的穩定性。

值得注意的是，在Kupershtokh等[12-13]的工作中，通過引入折合變量pr=p/pc,ρr=ρ/ρc,Tr=T/Tc，式中p表示壓強，ρ表示密度，T表示溫度，pc、ρc、Tc分別表示臨界壓強、臨界密度和臨界溫度，pr、ρr、Tr分別表示折合壓強、折合密度和折合溫度，對van der Waals狀態方程(VDW EOS)進行了重新定義，發現基于折合變量的狀態方程相比原方程沒有任何與具體物質相關的參數，因此可以定性地描述真實流體的相關特性，使得方程更具有普適性。

受到以上工作的啟發，本文擬將折合變量引入到化學勢多相流模型中，構建一種基于折合的化學勢LBM多相流模型，并對該模型在性能和實用性上進行數值模擬驗證。

1 改進化學勢模型

1.1 化學勢LBM多相流模型

在非理想流體系統中，含有梯度平方近似的自由能泛函可以表示為[9,17,20]：

(1)

其中積分第一項是自由能密度函數，第二項是密度梯度對自由能的貢獻，κ是與界面厚度和表面張力有關的系數?；瘜W勢由自由能密度函數決定：

(2)

其中,

(3)

由方程(2)和(3)可得：

μ=ψ′(ρ)-κ2ρ。

(4)

對于各類流體系統，其狀態方程的通用形式為：

p0=ρψ′(ρ)-ψ(ρ)。

(5)

系統的非理想力可以由化學勢求得[18]，即

F(x)=-ρμ+

(6)

另外，本文使用標準單馳豫晶格Boltzmann演化方程：

(7)

(8)

其中ωi表示權重系數。在不同的離散速度集模型DnQm(D表示維度，Q表示離散速度，n表示空間維度數，m表示離散速度數)中，ωi、cs與ei的取值是不一樣的，下面給出本文使用的D2Q9離散速度模型的相關參數：

(9)

(10)

在D2Q9模型下，流體的宏觀密度和宏觀速度可由分布函數求得：

(11)

(12)

平衡狀態下，方程(11)和(12)則是系統質量守恒和動量守恒的基本條件。

本文采用Shan-Doolen力項模型，通過修正平衡分布函數中的速度，將非理想力耦合到LBM演化方程中[21]。非理想力對碰撞過程的影響使得平衡態粒子的動量增加為

ρueq=ρu+τF。

(13)

(14)

1.2 折合表示的化學勢

根據方程(5)可以求得自由能密度與狀態方程之間的關系：

(15)

其中，C表示積分的常數項。因此，由方程(4)和(13)可知，只要給定狀態方程的形式，即可求得各類流體對應的化學勢。

以van der Waals流體為例，其狀態方程和對應化學勢的形式為：

(16)

(17)

引入折合變量(以下折合后的變量下標均標注為r)：

(18)

則得到折合后的狀態方程：

(18)

可以看出，相比原方程，折合狀態方程的形式更加簡潔，并且與參數a、b無關，避免了與具體物質有關的參數影響。因為折合后的方程是無量綱的，所以此形式下的狀態方程更具有普適性[12]。

在格子單位下，以格子時間Δt和格子步長Δx為特征量，將化學勢進行折合，得到折合化學勢：

由圖1可知，鹵汁1在復鹵過程中總酸含量總體上比鹵汁2要高，鹵汁1復鹵過程中總酸最低含量為0.91 g/kg，最高含量為3.88 g/kg；鹵汁2復鹵過程中總酸最低含量為1.13 g/kg，最高含量為3.79 g/kg。鹵汁中總酸變化的原因可能是，在豆干鹵制過程中，工廠每天要進行補料和補水，所以總酸會在一定的范圍內波動。鹵汁1在鹵制過程中，總酸含量在2.71 g/kg上下波動；鹵汁2在鹵制過程中，總酸含量在2.31 g/kg上下波動。而變化不大的原因是工廠每天會對鹵汁總酸進行監控。所以總酸只在一定的范圍內波動[9]。

(19)

則van der Waals流體對應的折合表示的化學勢為：

(20)

引入折合表示的流體速度和粒子分布函數：

(21)

相應地，折合表示的晶格Boltzmann演化方程變為：

(22)

平衡態分布函數變為：

(23)

(24)

(25)

這里給出折合表示的計算非理想力的分量形式：

(26)

力項模型依然采用Shan-Doolen的平衡態速度修正模型[21]，折合形式下平衡態速度和流體宏觀速度被重新定義為：

(27)

2 模型驗證

2.1 兩相共存密度

在氣液相變的數值模擬中，Maxwell等面積重構方法解出的兩相共存密度被廣泛作為驗證模型熱力學一致性的基準數據。將計算域劃分為100×200的長方形均勻網格流場，流場中部為液體，其余部分為氣體，氣液兩相界面設為水平的，四周采用周期性邊界條件，初始流體速度為零，為方便與原化學勢多相流模型作對比，根據文獻[18-19]，設弛豫時間τ=1.3，參數κ=0.01。

圖1 兩相共存密度曲線與Maxwell方法理論解的比較Fig.1 Comparison of two-phase coexistence density curves and theoretical solution of Maxwell’s method

為了能夠體現改進模型的擴展性，我們同時對VDW EOS、Carnahan-Starling狀態方程(CS EOS)、Redlich-Kwong狀態方程(RK EOS)、Peng-Robinson狀態方程(PR EOS)進行了模擬，并與Maxwell方法理論結果進行比較。如圖1所示，在不同狀態方程下，當前模型的數值計算結果與Maxwell方法的理論結果都能很好地吻合，這說明當前模型是符合熱力學一致性的。我們也對不同狀態方程模擬的溫度范圍進行了考察，如表1所示，與原始化學勢模型相比，使用當前模型模擬的溫度更低，溫度范圍更廣，這說明當前模型具有很好的穩定性。此外，使用當前模型顯著提高了模擬的兩相密度比，例如使用PR EOS和RK EOS能夠模擬的最低折合溫度都為0.34，最大兩相密度比分別為3.07×108和5.10×108，均超過了108。

2.2 Laplace定律驗證

Laplace定律通常用來評價多相流系統中表面張力的影響，也是驗證多相流模型可行性的重要依據。根據Laplace定律，當液滴與周圍氣體達到穩定狀態時，液滴內外壓差ΔP與液滴半徑R滿足：

ΔP=Pin-Pout=σ/R，

(28)

其中，σ表示液滴的表面張力,Pin和Pout分別表示液滴的內部壓強和外部壓強。

表1 采用原始化學勢模型與當前模型能夠模擬的溫度范圍和最大兩相密度比

由方程(28)可知，液滴內外壓差ΔP與半徑R存在線性關系，也就是說要想判斷模型是否符合Laplace定律，只需考察ΔP與R的倒數是否成一次函數關系即可。選取VDW EOS進行模擬，在網格大小為256×256的流場中央設置一個靜止液滴，四周采用周期性邊界條件，弛豫時間τ=1.3，參數κ=0.01?？疾煺酆蠝囟萒r分別為0.6、0.7、0.8時，半徑分別為30、35、40、45、50、55、60(格子單位)的液滴內外壓差與半徑之間的關系，液滴周圍為當前溫度下飽和的氣體。如圖2所示，在不同溫度下靜止液滴內外壓差ΔP隨半徑R的倒數的增加而線性增加，并且兩者成正比關系，此結果與Laplace定律相一致。這里每條線都有一個穩定的斜率，代表在當前溫度下表面張力的大小。由此驗證了當前模型滿足Laplace定律。

圖2 VDW流體在不同溫度下ΔP與1/R的關系Fig.2 Relationship between ΔP and 1/R of VDW fluid at different temperatures

2.3 降低計算中的虛速度

在使用LBM多相流模型進行模擬時，虛速度是衡量模型性能的一個重要指標。如果虛速度過大，則模擬中很難將其與真實流體速度區分開，使得計算結果失去意義。有效降低虛速度的影響一直是多相流領域研究的熱門話題。為方便比較，我們仍然以VDW EOS為例，計算區域設置為256×256大小的均勻網格，弛豫時間τ=1.3，κ=0.01，網格中央設置一個靜止液滴，液滴周圍為當前溫度下飽和的氣體，四周采用周期性邊界條件?？疾飚斦酆蠝囟萒r分別為0.82和0.92時，不同半徑液滴兩相界面處虛速度的大小。

如圖3所示，當Tr分別等于0.82和0.92時，當前模型的虛速度都明顯低于原始化學勢(CP)模型的，由此可以說明使用當前模型可以有效降低虛速度對流體行為的影響。

圖3 原始化學勢(CP)模型與VDW EOS的虛速度比較Fig.3 Comparing the spurious current of the original chemical-potential(CP) model and the present model

3 模型應用：液滴撞擊液膜

作為模型的測試案例和初步應用，我們模擬了單液滴在無重力影響的情況下，以一定初始速度撞擊液膜的過程。同樣選取VDW EOS，在網格區域大小為2 000×600的二維平面內設置一個液滴，計算域的底部鋪展一層液膜，高度為整個計算域高度的1/12，液滴半徑R=100，折合溫度Tr=0.6，弛豫時間τ=1.25，參數κ=0.01，流場上下邊界采用半程反彈邊界條件，左右邊界采用周期性邊界條件。在液滴撞擊液膜的過程中，雷諾數Re具有至關重要的作用[22-23]，通常Re=UD/νl，其中U為液滴撞擊液膜的初始速度，D為液滴的直徑，νl為液體的粘性。

圖5 鋪展半徑r與無量綱時間t*的函數關系Fig.5 Spread radius r as a function of the nondimensional time t*

4 結論

本文基于化學勢晶格Boltzmann方法多相流模型理論，通過引入折合變量，構建了一種新的化學勢多相流模型。對該模型從性能和實用性進行驗證，得到以下結論：

①使用新模型計算幾種常用狀態方程的兩相共存密度曲線，在更低的相對溫度下與Maxwell方法的理論結果符合得非常好，說明新模型仍具有原始化學勢模型的熱力學一致性。

②從穩定性、Laplace定律和虛速度等方面對新模型進行檢驗發現：與原化學勢模型相比，新模型在滿足Laplace定律的同時，顯著地擴展了模擬的溫度范圍，模型穩定性被大大提高，虛速度的影響也被有效降低，幾種常用的狀態方程能夠模擬的兩相密度比最高超過了108。

③利用新模型再現了液滴撞擊液膜的過程，并且撞擊后液滴飛濺水花的鋪展半徑能夠很好地服從文獻中所提出的冪律函數關系，說明新模型可以可靠地應用于實際模擬中，具有較好的應用前景。