廣義守恒相場簡化多相流格子Boltzmann 方法

章詩婷，肖鴻威，周錦翔，牛小東，*

(1. 汕頭大學 工學院，汕頭 515063；2. 汕頭大學 智能制造技術教育部重點實驗室，汕頭 515063)

0 引 言

多相流動現(xiàn)象不僅與自然界、人們的日常生活息息相關，也普遍存在于工業(yè)生產領域[1-2]。對于計算流體力學（computational fluid dynamic，CFD）研究人員來說，由于多相流固有的非線性、物理性質（如密度和黏度）的巨大差異，以及因質量、動量和能量同時進行交換而產生移動界面的復雜拓撲變化[3]，因此開發(fā)一個簡單有效的數值模型來研究這些復雜現(xiàn)象在實踐中面臨著巨大的挑戰(zhàn)。到目前為止，已經有許多研究學者對多相流進行了數值仿真，并針對多相流問題發(fā)展了一些數值模擬方法，包括相場方法（phase-field method， PFM）[4-6]、 流 體 體 積 法（ volume of fluid，VOF）[7-8]、水平集方法（level set method，LSM）[9-10]、前沿跟蹤法（front-tracking method，F(xiàn)TM）[11-12]、光滑粒子流 體 動 力 學 方 法（ smooth particle hydrodynamics method，SPHM）[13]。其中，相場方法是研究這類流動最常用的方法之一。在相場方法中，尖銳界面被視為有限厚度的過渡區(qū)，物理參數在其上平滑變化。

相場方法作為擴散界面方法中研究復雜界面動力學的一種，已被許多學者應用在模擬多相流問題上。相場方法主要分為兩種：一種是用包含空間四階導數的Cahn-Hilliard（C-H）方程[14]來捕捉界面變化，另一種是具有二階精度的Allen-Cahn（A-C）界面控制方程[15]。當前，大多數描述多相流的相場方法都是基于C-H 理論，其方程源自梯度流法[16]。Lee 等[17-19]開發(fā)了一種基于總自由能假設的非整體相場方法；Kim[20]進一步改進了該方法，通過引入新的離散delta 函數提出了廣義連續(xù)表面張力公式。與非整體相場方法不同，Elliort 和Luckhaus[21]提出了一種具有修正自由能公式的整體相場方法，對其引入了對稱正定矩陣，其中的系數與不同流體之間的成對表面張力有關。Dong[2,22-24]使用整體相場方法進行了一系列應用。然而，我們觀察到，由于總自由能隨時間減少，C-H 方程在保持各相體積方面的性能較差[25]。此外，當構造非線性勢時，在兩兩表面張力不均勻的情況下很難滿足還原一致性條件。

格 子Boltzmann 方 法（lattice Boltzmann method，LBM）作為一種介觀方法，自20 世紀80 年代末以來，已發(fā)展成為一種流行的CFD 替代工具。因為它的簡單性、在并行計算機上的可擴展性以及易于處理復雜幾何圖形，被廣泛應用于模擬流體動力學[26-27]和非線性方程組[28-29]。LBM 使用基于粒子的密度分布函數來演化宏觀系統(tǒng)。最近，基于相場理論，許多研究者結合LBM 來仿真研究多相流問題。Liang 等[30]提出一個LB 模型，通過巧妙地構建流場LB 方程中的力項來進一步發(fā)展C-H 動力學。Reza 等[31]開發(fā)了另一個PF-LB 模型，使用C-H 方程跟蹤三種不同流體之間的界面。然而，上述兩種方法仍受到體積守恒和質量守恒的影響。在A-C 方程的動力學建模中，Reza 等[32]提出了一種描述三相流的守恒相場LB 模型，該模型可保持良好的質量守恒性，并且滿足還原一致性條件。Zheng 等[33]針對不混溶多組分流體開發(fā)了另一種滿足還原一致性條件的守恒PF-LB 方法。Hu 等[34]提出一種廣義守恒PF-LB 模型，其中每一相相變量滿足廣義eikonal 方程和體積分數約束。

在本文中，我們在單松弛LB 框架下提出一種簡化的多相流LB 方法，用于仿真具有大密度比、大黏度比的多相流問題。該模型基于Hu 等[34]提出的廣義守恒PF 理論，其優(yōu)點在于可以處理不可壓、不混溶多相流，并具有良好的還原一致性以及質量和體積守恒性。此外，該方法保持了我們早期開發(fā)的用于兩相流大密度比的SMLBM 的良好數值穩(wěn)定性[35]。與用于模擬多相流的常用單松弛和多松弛LB 方法不同，簡化多相流格子Boltzmann 方法（simplified multiphase lattice Boltzmann method，SMLBM）在單松弛LB 方法的框架內通過預測-校正策略來模擬流體系統(tǒng)和跟蹤界面演化，其計算過程中僅需要考慮平衡分布函數的演化。因此，SMLBM 大大減少了計算中對虛擬內存的需求，簡化了物理邊界條件的實施，而這兩點往往是限制傳統(tǒng)單松弛和多松弛LB 方法應用的主要因素。最近，Li 等[36]構造了一種統(tǒng)一的SMLBM 方法來模擬大密度比、大黏度比的磁流體兩相流問題，驗證了此方法模擬多相磁流體流動的效率性和準確性。Khan 等[37]運用SMLBM模擬研究了均勻磁場下磁流體液滴變形和氣泡在磁流體中合并的動力學過程。隨后，Chen 等[38]發(fā)展了一種三相流SMLBM 方法，并將其應用于中等雷諾數剪切流中的復合液滴數值研究，但研究中的各相流體之間的密度比和黏度比依舊較小。因此，在模擬多相流系統(tǒng)中，解決大密度比、大黏度比多相流問題還較為困難，尤其是三相流及以上的多相流問題。但在現(xiàn)實生活中，氣-液多相流中的各相流體之間的密度比普遍都達到了1 000 以上。所以為了提高模型對大密度比多相流問題數值仿真的魯棒性，并克服單松弛模型對黏度的敏感性，需要提高相場方程界面捕捉的穩(wěn)定性和精確性。因此，我們將廣義守恒PF 方程中擴散項構造的源分布函數吸收到相場平衡分布函數中，這使得該方法的實現(xiàn)更加簡單。

通過拉普拉斯定律和具有不同表面張力比的液滴透鏡擴展兩個基準問題的數值仿真，我們首先驗證了本方法的穩(wěn)定性和準確性。進一步，通過模擬不同黏度比的三相泊肅葉流問題以及不同密度比復合液滴在平板上的擴散算例，將數值結果與解析解進行比較，展示了本模型的精確性和處理大密度比、大黏度比復雜多相流問題的能力。

1 數值方法

1.1 廣義守恒相場模型

相場模型實質上就是一類特殊的擴散界面模型。為了準確地描述整個流體系統(tǒng)中多種非混溶流體共存的狀態(tài)以及他們在相互作用下的運動，在相場模型中將任意位置上多種非混溶流體所占的體積分數定義為該點的相變量φi，顯然存在如下約束：

1.2 廣義守恒相場簡化格子Boltzmann 模型

顯然，上述SMLBM 的整個實施過程只涉及平衡分布函數，而平衡分布函數又可以直接從宏觀變量計算得到。與傳統(tǒng)的單松弛和多松弛LB 方法相比，不再需要計算流場的碰撞演化過程，提高了模型的計算效率。此外，由廣義守恒相場方程（4）中擴散項構造的源分布函數被吸收到式（20）中的相場平衡分布函數中，這使得該方法更加簡單和穩(wěn)定。

2 模擬算例

2.1 拉普拉斯定律

本節(jié)首先進行拉普拉斯定律的基準試驗，以此來驗證模型的準確性。我們將兩個不同密度的液滴（流體1、流體2）同時浸入流體3 中，如圖1 所示。在表面張力的影響下，液滴內部的壓力大于液滴外部的壓力，在達到平衡時，液滴會以圓形保持靜止。同時液滴內外界面壓力差Δpij與表面張力σij之間滿足以下公式：

圖1 兩個液滴浸沒到另一種流體中示意圖Fig. 1 Diagram of two droplets immersed in another fluid

首先為了驗證網格的收斂率，分別在計算域為x×y= 100×50、200×100、400×200 和600×300 時進行仿真。我們設定參考長度Lc= 2Rij，并將其分別設置為30、60、120 和180。根據下述公式計算了三種相變量在不同網格大小下的數值誤差：

三種相變量的數值誤差如圖2 所示。可以發(fā)現(xiàn)三種相變量的誤差呈現(xiàn)相同的趨勢，即網格越小，相變量的誤差越小，說明模型的收斂精度與網格的大小有關。當計算域為100×50 時，網格最大，其收斂精度為1×10-3；當計算域為600×300 時，網格最小，收斂精度達到1×10-4，均在二階精度以內。從圖3 中可以看出，400×200 和600×300 兩種網格的誤差非常接近，可以判斷在這個網格數量區(qū)間內，模型的精度不會隨著網格數量的增加而有較大變化，進一步證明了網格獨立性。為了保證計算精度，同時盡可能地減少計算時間，我們選擇計算域為400×200 的網格大小。

圖2 當前模型的收斂誤差Fig. 2 Convergence errors of the present model

圖3 在不同網格大小下 E (φ3)隨時間變化Fig. 3 The time histories of E (φ3) under different grid sizes

圖4 中列出了界面壓力差在兩個液滴不同的表面張力（σ13、σ23）下和不同半徑（Rij= 20、30、40、50、60）下模擬值與理論值的對比。圖中可以清楚地看出模擬值和理論值吻合得非常好。我們對比了兩個液滴在不同半徑下相變量φ在水平中央線的模擬值和初始值，如圖5 所示。從圖中可以看出，在液滴半徑為Ri j=20、60時，相變量的模擬值和理論值依舊重合得很好。顯然，模型能夠準確地追蹤界面的變化且能夠保持較高的穩(wěn)定性。

圖4 液滴1 和液滴2 界面壓差的模擬值和理論值對比Fig. 4 Comparison of the pressure difference across the interfaces of droplet 1 and droplet 2 between the numerical results and theoretical solutions

圖5 相變量在水平中央界面分布對比圖（R ij=20、60）Fig. 5 Comparison of the simulated profile of the order parameter with the theoretical solutions (R ij=20, 60)

2.2 液滴透鏡

本節(jié)研究了一個透鏡在兩層不相溶流體間的擴散運動這一經典算例，液滴透鏡擴散問題在文獻[10,32,41,42]中被廣泛用于驗證模擬多相流的數值方法。液滴透鏡如果一開始存在于兩相不混溶的流體之間，在表面張力的作用下透鏡會發(fā)生形變，其最終的平衡狀態(tài)是上下兩層流體間相互作用的結果。透鏡的平衡狀態(tài)由三個表面張力系數決定，根據諾伊曼定律可以得出以下關系[43]：

圖6 液滴透鏡參數示意圖Fig. 6 Diagram of the spreading of a liquid lens

表1 液滴透鏡算例中表面張力系數的設定Table 1 Surface tension coefficients used in the simulation of spreading of a liquid lens

4 種情況下的透鏡平衡狀態(tài)如圖7 所示，顯然，不同的表面張力系數比會導致透鏡最后的平衡形狀不同。當σ23的值增加時流體2 的作用力會增大；同理，當σ13的值減小時，流體1 的作用力也會相應減小。正因為上下作用力的改變形成壓力差導致了透鏡形變的不同。我們將平衡時透鏡兩個三相點之間的距離d，以及h1和h2的模擬值和理論值進行了對比，如表2 所示。從表中可知，模擬值和理論值吻合得較好，最大的相對誤差值為5.6%。這種微小的誤差是由于在界面厚度為零的極限情況下計算了理論值。

圖7 不同表面張力系數下液滴透鏡的平衡狀態(tài)Fig. 7 The equilibrium shape of the liquid lens for different surface tension coefficients

圖8 展示了四組不同表面張力系數下液滴透鏡質量隨時間的變化。在圖中，以液滴透鏡的初始質量（m0）為評判標準，若質量守恒，液滴透鏡的初始質量和實時計算質量（m1）的比值為1。從圖中可以看出，四種方案的液滴透鏡在演變過程中的質量比曲線與準確質量比1 的直線基本重合，質量耗散誤差都在1%以內，表明本方法能夠保持良好的質量守恒性。

圖8 四種方案液滴透鏡的初始質量和計算質量比值隨時間的變化Fig. 8 Variations of liquid lens mass ratio with time for four cases

2.3 三相泊肅葉流

泊肅葉流動問題主要用于驗證黏度比對模型計算精度的影響。這里我們模擬了三相泊肅葉流動問題。假設各相流體之間的界面是光滑的，對于相場的相變量初始化使用了如下函數：

其中，h為矩形槽道的高度，槽道長度為2h，水平槽道內充滿了三層互不混溶的流體，如圖9 所示。

圖9 三相泊肅葉流示意圖Fig. 9 Diagram of Poiseuille flow with three phases

在三相泊肅葉流動中，流體由一個體力G=ρg驅動，g是x方向的恒定加速度。因為各相的流體密度為常量，整個系統(tǒng)的N-S 方程表述如下：

表3 三相泊肅葉流算例中密度和運動黏度的設定Table 3 Densities and kinematic viscosities used in three-phase Poiseuille flow example

我們繪制了速度分量ux在y軸上的分布曲線，并將當前結果與Hu 等的結果進行對比，如圖10 所示，發(fā)現(xiàn)三種情況下兩者的曲線吻合較好。圖10（a）中，由于三種流體之間的黏度比為1，沿y軸的速度剖面較為光滑。而在方案2 和方案3 中，黏度比增大，圖10（b、c）非常好地捕捉到了因黏度變化而引起的相界面處明顯的速度分量ux的跳躍。

圖10 當前結果與Hu 等結果[34]之間速度剖面的比較Fig. 10 Comparison of velocity profiles between present results and Hu et al.[34]

2.4 復合液滴鋪展

為了驗證模型是否能夠克服大密度比多相流交互界面處的較大壓力梯度問題，這里采用一個二維復合液滴在光滑平板上的鋪展問題來進行數值模擬。在這個算例中充分考慮了濕潤邊界條件，圖11 中流體1、流體2 和流體3 為非混溶的三種流體，與固體壁面從左至右形成的接觸角分別為θ13、θ12和θ23，由幾何關系可知，θ21= π -θ12。為了滿足三相表面張力系數平衡條件，流體i和流體j與光滑平面形成的接觸角θij可根據Young’s 法則定義如下：

圖11 三種流體與固體壁面接觸角示意圖Fig. 11 Diagram of a compound droplet on a solid substrate

為了解決三相流系統(tǒng)與固體平面接觸面上接觸角的不連續(xù)性問題，在濕潤邊界處理上設置加權接觸角θ1和θ2來判定濕潤度，保證界面均勻光滑地過渡。加權接觸角的定義如下：

濕潤邊界的實施還需要人為在固體壁面下側添加一層虛擬網格，如圖12 所示。使用中心有限差分格式離散單元格（i，1）后，通過相鄰網格的數值近似獲得虛擬網格（i，0）的值[45]：

圖12 濕潤邊界條件實施示意圖Fig. 12 Diagram of the implementation of wetting boundary conditions

復合液滴的初始形態(tài)如圖13 所示，流體1 和流體2 分別為一個四分之一圓，它們共同組成一個半圓復合液滴，周圍區(qū)域為流體3。此算例中，三相流體的運動黏度 υ相同。為了驗證模型在不同密度比下的準確性和穩(wěn)定性，設置了四組不同密度比的復合液滴鋪展算例進行對比，如表4 所示。復合液滴初始直徑2R選作特征長度，計算域設定為4R×2R，網格劃分成400×200 的 均 勻 網 格，界 面 厚 度 ξ = 4，參 考 速 度Uc= 0.01，流體區(qū)域四周邊界均采用無滑移邊界條件。

表4 復合液滴鋪展中三種流體的密度Table 4 Densities used in the simulation of spreading of the compound droplet

圖13 復合液滴在平面上鋪展初始形態(tài)Fig. 13 Initial shape of the compound droplet

在該算例中將濕潤邊界條件接觸角設置為：θ13=90°、θ12= 120°、θ23= 60°。假設任意兩相界面間的表面張力系數相同，其值為0.01，則有三相角φ1=φ2=φ3= 120°。整個復合液滴不考慮重力因素，主要在表面張力和黏性力作用下推動液滴形態(tài)變化。這里用無量綱參數Re和We來表征液滴的動力學特性，設置Re=30，Weij=2。

圖14 表示了在方案1 情況下，復合液滴在光滑平板上隨時間發(fā)展的狀態(tài)圖，不同顏色的線分別表示不同時間節(jié)點復合液滴所處狀態(tài)，這里時間節(jié)點分別取t= 0、0.10、0.25、0.50、0.75、2.00。從圖中可看到復合液滴隨著時間的演變，慢慢地向外鋪展，當復合液滴到達平衡狀態(tài)時其形態(tài)會保持不變，平衡時復合液滴的形狀主要由初始設置的接觸角和三相角給定。這里，我們令流體1 和流體2 與固體平面接觸的潤濕長度分別為L1和L2，其總潤濕長度為L，如圖15 所示。其中潤濕長度的理論值是根據Zhang 等[46]計算的結果，并與模擬值進行了對比，如表5 所示。從表中可以發(fā)現(xiàn)誤差值均小于1%，理論值與數值模擬結果非常吻合，驗證了該模型的精確性。

圖14 方案1 的復合液滴隨時間演化過程示意圖Fig. 14 Time evolution of the compound droplet of Case 1

圖15 方案1 的復合液滴平衡狀態(tài)Fig. 15 Equilibrium state of a compound droplet of Case 1

表5 方案1 的潤濕長度理論值與模擬值對比Table 5 Comparison of the numerical results and analytical solutions of wetting length in Case 1

方案2～方案4（如圖16 所示）的復合液滴的平衡形狀與方案1 相似。但是，隨著密度比的增加，復合液滴達到平衡狀態(tài)所需的時間逐漸增加。這是因為密度比的增加導致壓力梯度和界面上黏性力的增加，從而進一步延遲復合液滴擴散到平衡狀態(tài)。然而，即使對于密度比高達1 200 的方案4，本模型仍能給出穩(wěn)定解。四種情況下接觸角的數值結果與表6 中的理論值進行了比較，其中 Δθij表示各接觸角的相對誤差。結果表明，誤差均小于2%，進一步證明了模型的準確性。

表6 不同密度比復合液滴平衡狀態(tài)下接觸角的模擬值與理論值對比Table 6 Comparison of the numerical results and analytical solutions of the compound droplet for different cases under equilibrium state

圖16 不同情況下復合液滴的平衡形態(tài)Fig. 16 The equilibrium shape of the compound droplet for different cases

為了更清楚地觀察不同密度比下復合液滴的平衡形態(tài)，將四種情況下復合液滴平衡態(tài)圖繪制在一張圖中，如圖17 所示。我們發(fā)現(xiàn)了一個有趣的現(xiàn)象，隨著密度比的增加，流體1 和流體2 之間的交界線逐漸向左移動。此外，在較大的密度比下，潤濕邊界總長度L也增加，三相邊界點不僅向左移動，而且隨著密度比的增加向上移動。顯然，在方案4 中，復合液滴中流體1 和流體2 的形狀不像方案1 那樣規(guī)則，從而增加了總潤濕長度。有趣的是，四種情況下設置的接觸角是相同的，并且由于密度比的增加，復合液滴的形狀并不規(guī)則，無法計算潤濕長度的理論值，但是最終平衡狀態(tài)下接觸角仍與預設的平衡接觸角一致。

圖17 不同密度比下復合液滴平衡狀態(tài)圖Fig. 17 Equilibrium shape of compound droplets at different density ratios

為了更加直觀地觀察三相交界點的變化規(guī)律，繪制了在總潤濕長度L的變化過程中，不同密度比下三相交界點高度H的變化，如圖18 所示。首先在圖中可以明顯地看出，密度比增大的同時，復合液滴三相交界點的高度隨之增高，并在復合液滴演化的過程中一直保持這種趨勢。其次，在復合液滴到達平衡狀態(tài)時，密度比越大，總濕潤長度也會越長，其中方案4 的L是增長最顯著的。這一現(xiàn)象可歸因于兩種流體之間的作用力和壓力差，它們隨著密度比的增加而增加，因此復合液滴擴散得更多。

圖18 不同密度比下復合液滴三相交界點高度隨總潤濕長度的變化曲線Fig. 18 Curves of H varying with L for different cases

3 結 論

本文提出了一種廣義守恒相場簡化多相流格子Boltzmann 方法。該方法繼承了Hu 等推導出的廣義守恒相場方程的還原一致性和質量體積守恒的優(yōu)點，并結合早前提出的SMLBM 用于處理大密度比、大黏度比問題，其中廣義守恒相場方程的源分布函數被吸收入相場平衡分布函數中，以保持方法的良好穩(wěn)定性。此外，本方法大大降低了內存的需求，這通常是限制傳統(tǒng)單松弛或多松弛LB 方法應用的一個主要因素。為了驗證本方法在處理多相流時的網格獨立性和精確性，我們首先模擬了拉普拉斯定律和液體透鏡兩個基準算例，并獲得了令人滿意的結果。其次，對不同黏度比下的三相泊肅葉流進行了仿真，驗證了此方法能夠模擬大黏度比（最高可達500）問題。最后，對不同密度比下的復合液滴在潤濕基底上的擴散進行了模擬，數值結果表明，該方法能夠模擬大密度比多相流問題，其中最大密度比可達1 200。對此還發(fā)現(xiàn)在較大的密度比下，復合液滴輪廓變得不規(guī)則，導致總潤濕長度增加，并且還找到了復合液滴在不同密度比下三相交界點的變化規(guī)律。以此我們可以得出結論，本文所提出的方法在模擬多相流的準靜態(tài)問題上有較好的表現(xiàn)。