耦合界面力的兩相流相場格子Boltzmann模型?

李洋 蘇婷 梁宏 徐江榮

(杭州電子科技大學理學院,杭州 310018)

(2018年6月25日收到;2018年9月14日收到修改稿)

1 引 言

格子Boltzmann(LB)方法[1,2]是近三十年來發展起來的一種流體系統建模和模擬方法.它不再基于宏觀連續介質模型的Navier-Stokes方程,而是直接從微觀模型出發通過描述流體粒子分布函數的演化再現復雜流動的宏觀行為.因此,LB方法相比傳統數值方法有一些獨特優勢,而受到了眾多學者的關注,并已廣泛應用于多個學科領域,如多孔介質中傳熱傳質、湍流的直接數值模擬、多相流體流動與傳熱、顆粒流[1,2]等.特別是LB方法的介觀屬性,可以使其直接描述多相流體系統中流體間與流固間的相互作用,因而可以方便地模擬多相流體流動,這也被認為LB方法區別傳統數值方法的最大優勢之一.

目前,許多學者已經根據流體間微觀相互作用力的不同物理背景,提出了多種不同類別的兩相流LB模型[1,2],包括顏色模型、偽勢模型、自由能模型、基于相場理論的模型.在這些提出的LB模型中,基于相場理論的兩相流LB模型[3?9]因具有堅實界面追蹤的物理機制,近年來受到了廣泛關注,并且已成功地應用于復雜多相流動問題的研究[10?13].He等[3]基于相場理論提出了第一個不可壓兩相流LB模型,其基本原理是通過引入序參數分布函數來追蹤流體界面,而用另一個壓力分布函數來求解流場.然而,Zheng等[4]指出He模型存在所恢復界面追蹤方程與相場理論中Cahn-Hilliard方程不一致的問題,并通過在界面演化方程中引入序參數分布函數的差分項,提出了與相場方程相一致的LB模型.另外,針對He模型所模擬大密度比的限制,Lee和Liu[5]提出了二階混合差分格式來離散流體間的相互作用力,從而提高了模擬大密度多相流動問題的數值穩定性.然而,一些學者指出混合差分格式的使用會導致Lee模型不滿足全局質量守恒[14].受Zheng等模型的啟發,Zu和He[6]在界面演化方程中引入一種平衡態分布函數的空間差分項,從而也能夠準確恢復到Cahn-Hilliard方程.另外,為了恢復正確的不可壓Navier-Stokes方程,Zu和He在流場的LB算法中引進一些復雜的源項,這導致模型中速度和壓力滿足一組隱式方程,需要通過預估與校正方法來求解.最近,Liang等[7]也從相場理論出發構造了一類界面求解精度更高且數值穩定性更好的兩相流LB模型,并且流場中宏觀速度和壓力的計算更為簡便.

綜上所述,兩相流的相場LB模型得到了較大的發展.然而,在相場LB模型中,存在著一類特殊的界面力[15].通過Chapman-Enskog理論分析發現,相場LB模型中均存在著界面力的尺度與理論分析結果不一致的問題.另外,為了描述流場的信息,相場LB模型中均含有一些復雜的外力項.針對這些問題,本文考慮了界面力的尺度效應,并提出了更加簡單的外力項分布函數.通過多尺度分析,本模型也能恢復到正確的宏觀控制方程.

2 控制方程和數學模型

2.1 控制方程

在相場理論中,用于描述多相系統中相界面運動的宏觀控制方程是著名的Cahn-Hilliard(CH)對流擴散方程[6,7],

這里,?是標識不同流體的序參數,u是流體速度,M是遷移率,μ是自由能的化學勢,在相場理論中可以表示為

其中?l和?g分別表示液相和氣相的序參數,而相界面可以表征為?=(?l+?g)/2的等值線,參數β和κ由界面厚度D和表面張力σ決定,

另外,為了描述多相流體流動,界面追蹤CH方程需要耦合流體動力學方程.根據文獻[7],含有外力項的不可壓縮黏性流體的Navier-Stokes(NS)方程可表述為:

其中,ρ是流體密度,p是壓力,ν是流體黏性系數,F為總外力,包含外力項G和表面張力Fs,在本文的數值模擬中,我們取常用的勢形式Fs=???μ或Fs=μ??[6,7,15]來減少流體界面處的虛假速度.

2.2 兩相流格子Boltzmann模型

為了求解CH方程和NS方程,本文引入了兩個獨立的分布函數fi和gi,其對應單松弛LB模型的演化方程可分別表示為:

其中τf和τg是分別與遷移率M和運動學黏性系數ν有關的松弛因子,fi(x,t)是粒子在t時刻x位置的序參數分布函數,gi(x,t)是密度分布函數,相應的平衡態分布函數和分別被定義為[7]:

這里

其中η是調節遷移率大小的參數,ωi是權系數,cs為格子聲速.格子速度的離散模型決定著η和ωi的取值,對于D2Q9模型,相應的權系數√ω0=4/9,ω1?4=1/9, ω5?8=1/36,cs=離散速度ci為

其中,c= δx/δt,δx是單位格子大小,δt為時間步長.本文中,我們取c=δx=δt=1.

為了準確恢復到CH方程,Liang等[7]在演化方程引入了可以局部計算的時間導數項,其對應的源項分布函數可定義為

在不同的兩相流相場LB模型中,外力項分布函數具有不同的形式[6,7,15].在原始的相場LB模型中,外力項分布函數可定義為[7]

基于Chapman-Enskog多尺度分析,我們發現速度的計算表達式中外力應為o(ε)尺度,然而Fa為o(ε2)尺度,這就導致速度計算式出現了尺度不一致的問題.界面力尺度與理論分析不一致問題同樣存在于其他的相場LB模型[6,15]中.另外,我們注意到為了恢復正確的不可壓Navier-Stokes方程,相場LB模型在流場的演化方程中均引入了較為復雜的外力項分布函數,這影響了模型的實現難度和計算效率.針對上述問題,我們考慮了尺度效應,并提出了一種新的更加簡單的外力項分布函數

本模型中序參數?的計算如下:

宏觀量密度ρ可以看作序參數?的線性函數,

其中ρl與ρg分布代表液相和氣相的密度.修正后的速度表達式及壓力計算式可表示為:

通過Chapman-Enskog多尺度分析[7],可以證明本模型中界面追蹤的LB演化方程(6a)可正確地恢復到Cahn-Hilliard控制方程,并且遷移率M與松弛因子關系可表示為

接下來,我們簡要地論述本文提出新的外力項分布函數的LB演化方程(6b)如何恢復到不可壓流體的Navier-Stokes方程.根據外力項分布函數的表達式(14),其對應的矩條件為

為了恢復正確的宏觀控制方程,我們引入如下的多尺度展開式:

其中?是非常小的展開參數.將LB演化方程(6b)Taylor展開到δ2t量級,并應用上述的多尺度展開式可得

其中,D1i= ?t1+ciα?1α. 由分布函數和Gi的矩條件,gi的零階矩可定義為[10]

進一步,應用多尺度展開式到(17)和(24)式,可得分布函數gi滿足如下性質

利用上述結果,對(23b)式分別求零階和一階速度矩,可以得到t1尺度上的宏觀方程

同理,t2尺度上的宏觀方程可從(23c)式的零階和一階速度矩得到,

其中,階數為O(δtMa2)的一些項已被略去.代入上述結果可將(30)式改寫為

對t1和t2尺度上的方程進行尺度黏合,最終發現可以正確地恢復到不可壓縮黏性流體的Navier-Stokes方程組(5a)和(5b),并且運動學黏性系數ν(23b)式和各階矩條件,經系列計算最終可得到為

在古詩文的詞匯中，實詞的意義相對固定，隨著時代的更迭，會產生相應的拓展和引申，但總的來說學生易于掌握，難的莫過于虛詞。劉彥和將常見的虛詞作了如下歸類：“至于夫惟蓋故者，發端之首唱；之而于以者，乃扎句之舊體；乎哉矣也，亦送末之?？啤！辈⒅该魈撛~在作者回環婉轉的巧妙運用下，文辭就會變得嚴密。因此，作者在創作過程中會努力避免虛詞的運用出現謬誤，那么我們在閱讀過程中同樣要將虛詞理解準確。否則，即使理解了文意也是不嚴密的。當然，在教學實踐中不能完全按劉彥和的方法來辨析，畢竟《文心雕龍》成書于南朝，此后的詞匯在詞性和詞義上都產生了變化與演進。

在數值模擬中,需要對時間導數項和空間導數項進行離散計算.本文中選用顯式的歐拉差分格式來計算時間導數項,

針對空間導數項,為了保證多相系統的全局質量守恒,我們采用如下的各向同性的二階中心差分格式[14],

其中χ表示任意變量.

3 數值結果和討論

在本節中,首先通過模擬幾個典型的兩相流問題,包括靜態液滴、液滴合并、亞穩態分解和單模瑞利-泰勒不穩定性問題,來驗證所提出的兩相流相場LB模型的準確性和穩定性.接著,利用本模型來研究隨機擾動分布的多模瑞利-泰勒不穩定性問題,著重考察雷諾數對相界面的影響規律.

3.1 靜態液滴

首先模擬靜態液滴這一基本的兩相流問題用于驗證本模型,并由此測量本模型產生的虛假速度.該問題的初始條件設置如下:半徑R=25的圓形液滴坐落在100×100的計算區域內,邊界四周采用周期性邊界條件,其他的相關物理參數設置為ρl=100,ρg=1,?l= ??g=1,D=5,σ =0.001,τf=0.8,M=0.1,ν=0.1.序參數的初始分布給定為

在這里(xc,yc)=(50,50)是液滴中心的坐標.根據方程(16),可以得到相應的初始密度分布.圖1(a)給出平衡狀態下整個計算區域的速度分布,以及黑色的實線表示液滴的初始形狀,紅色的虛線表示平衡狀態下的液滴形狀.從圖中可以發現它們吻合得很好,這驗證了本模型能夠準確的求解相界面.此外,圖1(b)給出了沿著線x=Ny/2的密度分布,可以發現密度場的解析解與通過相場LB方法模擬獲取的數值結果相一致.相界面處存在虛假速度一直是數值方法普遍存在的非物理現象.我們進一步統計了本模型產生的虛假速度,結果表明其虛假速度的最大值約為3.57×10?9.在前人的相場LB模型中,Zu和He[6]的模型產生的虛假速度為10?8量級,而Liang等[7]提出的模型在單松弛情況下產生的虛假速度為7.8×10?6.對比前人相場LB模型,可以發現本文修正的LB模型可以獲得相對較小的虛假速度.

圖1 靜態液滴測試 (a)平衡狀態下的整個區域的速度分布圖,實線和虛線分別代表液滴的初始形狀和在平衡狀態下的形狀;(b)沿著x=Ny/2的密度分布的LB模擬結果和解析解Fig.1.Static droplet test:(a)The velocity distribution of the whole domain at the equilibrium state and the solid and dashed lines respectively represent the initial shape of the droplet and its equilibrium shape;(b)density profile across the interface obtained from LB simulation and corresponding analytical solution along x=Ny/2.

3.2 液滴合并

圖2 兩個液滴的合并過程(d=5,D=2.4) (a)t=0;(b)t=20000;(c)t=40000;(d)t=60000Fig.2.The combination of two static droplets(d=5,D=2.4):(a)t=0;(b)t=20000;(c)t=40000;(d)t=60000.

圖3 兩個液滴的合并過程(d=5,D=2.6) (a)t=0;(b)t=20000;(c)t=23000;(d)t=50000Fig.3.The combination of two static droplets(d=5,D=2.6):(a)t=0;(b)t=20000;(c)t=23000;(d)t=50000.

3.3 亞穩態分解

進一步,我們模擬了更為復雜的相分離問題[6,17]來驗證本文提出的兩相流相場LB模型.亞穩態分解,也稱相或組分分離現象,是混合流體的一種基本特征.對處于亞穩態狀態的均勻混合流體施加一個很小的擾動,此時流體系統是不穩定的,相分離現象將會發生.在我們的模擬過程中,計算區域設置為NX×NY=100×100,四周采用周期性邊界條件.初始的序參數分布設定為

根據密度與序參數之間的線性關系,可以得到相應的初始密度分布.模擬中其他的物理參數設置為:液相密度ρl=10,氣相密度ρg=1,?l=1,?g=0,τf= τg=1,界面厚度D=4,遷移率M=0.5.圖4展示了表面張力σ=0.001時兩相分離現象的演化過程,從圖中可以觀察到,兩種流體初始時刻處于混合狀態,隨著演化時間的增加,流體在流體間的相互作用力下漸漸地聚集,從而系統中形成了小液滴且尺寸隨時間不斷地增大,最終在t=160000時亞穩態分解形成兩相流體,上述相分離過程與文獻結果[17]定性地相一致.另外,我們進一步研究了表面張力對相分離過程的影響,將表面張力調整到σ=0.01.圖5展示了表面張力σ=0.01時兩相分離現象的演化過程.從圖中可以發現,增大表面張力不會影響相界面的演化圖案,但可以加速相分離過程,相比σ=0.001情形可以更快的達到相分離的穩定狀態.

3.4 瑞利-泰勒不穩定性

圖4 σ=0.001時相分離的時間演化過程 (a)t=0;(b)t=28000;(c)t=38000;(d)t=49000;(e)t=113000;(f)t=160000Fig.4.Time evolution during the phase separating process at σ=0.001:(a)t=0;(b)t=28000;(c)t=38000;(d)t=49000;(e)t=113000;(f)t=160000.

圖5 σ=0.01時相分離的時間演化過程 (a)t=0;(b)t=3000;(c)t=4000;(d)t=7000;(e)t=11000;(f)t=16000Fig.5.Time evolution during the phase separating process at σ=0.01:(a)t=0;(b)t=3000;(c)t=4000;(d)t=7000;(e)t=11000;(f)t=16000.

在本小節,我們模擬了兩相流體間的經典界面不穩定性現象,即瑞利-泰勒(Rayleigh-Taylor,RT)不穩定性問題[18,19].當密度較大的流體放置于密度較小的流體上層時,給流體界面處施加一個微小擾動,在重力的作用下,上層流體流入下層,下層流體從兩側向上涌入,最終達到混沌混合狀態,這種現象稱之為RT不穩定現象.最早對RT不穩定性問題的研究可追溯到著名學者Rayleigh[18]對云層形成的研究,后來學者Taylor[19]又在原子核聚變中發現了RT不穩定性現象.自此以后,許多學者對RT不穩定性問題開展了理論、實驗或者數值研究,有興趣的讀者可以參考最近關于RT不穩定性的文獻綜述[20].

我們考慮一個矩形微管道,其高度L和寬度W的比例為4,對兩相界面處施加一個微小擾動,其擾動函數給定為

初始的序參數分布則設定為

雷諾數(Re)和Atoods數(At)是影響RT不穩定性的兩個重要無量綱參數,分別定義為

其中,g是重力加速度,為了耦合重力效應,在模擬中對計算區域的所有流體點施加一個重力,其定義為

在數值模擬中,液相密度ρl=3,氣相密度ρg=1,相應的At數為0.5,其他物理參數設置為:?l= ??g=1,W=128,=0.04,D=4,τf=0.8,表面張力σ=0.0001.佩克萊數(Pe)的定義如同于參考文獻[6,7],并取為50.另外,我們討論雷諾數對相界面的影響,并通過改變雷諾數的大小可以確定流體黏性系數,進而確定松弛因子τg的值.上下壁面均采用無滑移的邊界條件,左右邊界應用周期邊界條件.圖6給出了Re=256時非混相RT不穩定性中相界面隨時間的演化圖.從圖中可以看到,在初始時刻,由于兩相界面處受到微小擾動的作用,液相重流體在中間的凹界面處往下運動,而氣相輕流體沿著兩側的壁面向上升起,從而分別形成了尖釘和氣泡.隨著演化的不斷進行,尖釘和氣泡繼續相互滲透,伴隨著尖釘兩側向上卷起形成了一對旋轉方向相反的漩渦.此后,在流體間的剪切力下,漩渦不斷地發展,其尺寸不斷地增大.可以發現,各個時刻相界面的演化過程與He等[3]的結果在定性上相一致.我們進一步地考察了Re數對RT不穩定性演化特性的影響,并且模擬了Re=1024較大的情形.圖7給出了Re=1024時非混相RT不穩定性中相界面演化圖.從圖中可以發現,在演化前期,不穩定性的界面演化圖案與Re=256情形相類似:液相流體與氣相流體相互滲透形成尖釘和氣泡,隨著時間演化,尖釘向上卷起形成兩個關于中間軸對稱的漩渦.接下來,相界面呈現出與低Re數情形顯著不同的形態.一對漩渦繼續發展,其尺寸不斷地增大,伴隨著在漩渦尾端形成二級漩渦,這一現象并未在低Re數的情形中觀察到.隨后,流體間界面不穩定性繼續發展,在高流體間剪切力作用下,形成了多對漩渦的復雜界面結構.當Re較大時,系統中相界面的非線性程度和流場較為劇烈,此時界面力效應比較顯著.因此,我們通過對比本模型和Li等[15]的模型所模擬Re=1024條件下RT不穩定性的結果,來體現界面力的效應.為了消除界面追蹤LB模型所帶來的影響,本文采用的界面追蹤LB演化方程(6a)應用于Li等的模型中.圖8給出了Li等的模型計算獲取的非混相RT不穩定性中相界面隨時間的演化.對比圖7的結果,可以發現兩種LB模型獲得的相界面圖案在不穩定性前期差別不大,這是由于在演化前期,不穩定性處于線性增長階段,流場較為溫和,界面力的效應不顯著.而在演化后期,不穩定性進入非線性增長階段,界面力的效應開始顯著.從圖8可以看出,相界面相比圖7變得光滑,在尖釘的尾端未觀察到顯著的二級漩渦.進一步,我們也定量地統計了不同Re數下隨時間演化的尖釘和氣泡在y方向上的位置,并將結果展示在圖9,其中時間t已經被特征時間無量綱化.從圖中可以看到,氣泡振幅和尖釘振幅都隨時間不斷地增大,并且尖釘振幅的增長速度大于氣泡的增長速度.此外,還可以發現,隨著Re數的增大,氣泡振幅和尖釘振幅在演化前期幾乎不受影響,而在演化后期,氣泡振幅和尖釘振幅隨著Re數的增大而增大,這是由于在演化前期,不穩定性較弱,流體與流體剪切力較小,而在演化后期,不穩定性較劇烈,流體與流體間的剪切力較大.作為對比,我們還將He等[3]的結果和Li等[15]的模型的計算結果列在圖9.可以看到,本文數值模擬獲得的尖釘和氣泡位置在演化前期與文獻模型計算的結果相吻合,而在演化后期,本模型計算結果與文獻模型結果有一定的差異.這是由于He等未考慮界面力的影響以及Li等未準確地處理界面力的尺度效應所造成的.

圖6 Re=256時,本模型計算獲取的RT不穩定性演化圖 (a)t=0;(b)t=1;(c)t=2;(d)t=3;(e)t=4;(f)t=5Fig.6.At Re=256,time evolution of RT instability obtained by the present model:(a)t=0;(b)t=1;(c)t=2;(d)t=3;(e)t=4;(f)t=5.

圖7 Re=1024時,本模型計算獲取的RT不穩定性演化圖 (a)t=0;(b)t=1;(c)t=2;(d)t=3;(e)t=4;(f)t=5Fig.7.At Re=1024,time evolution of RT instability obtained by the present model:(a)t=0;(b)t=1;(c)t=2;(d)t=3;(e)t=4;(f)t=5.

圖8 Re=1024時,Li等[15]模型計算獲取的RT不穩定性演化圖 (a)t=0;(b)t=1;(c)t=2;(d)t=3;(e)t=4;(f)t=5Fig.8.At Re=1024,time evolution of RT instability obtained by model of Li et al.[15]:(a)t=0;(b)t=1;(c)t=2;(d)t=3;(e)t=4;(f)t=5.

圖9 不同Re數下隨時間演化的(a)氣泡和(b)尖釘在y方向上的位置Fig.9.Time evolution of y coordinate of(a)bubble and(b)spike poistions at different Re.

通過文獻調研表明,對RT不穩定性問題的研究絕大多數為單模情形,針對隨機擾動的多模RT不穩定性的研究則相對較少[20].本文的最后,我們將利用修正的兩相流相場LB模型對隨機擾動的多模RT不穩定性進行模擬,并研究Re數對相界面演化的影響規律.物理問題設定為一個L×L的計算區域,對兩相界面處施加一個初始的微小擾動,其中擾動函數給定為

其中,an和bn是滿足高斯正態分布的隨機數,序參數的初始分布則給定為

在數值模擬中,一些物理參數設置如下:ρl=3,ρg=1,?l= ??g=1,L=512,=0.08,D=4,σ=0.0001.

首先模擬Re=4096的情形,其對應的流體黏性為0.01.圖10給出了Re=4096時多模RT不穩定性中密度界面隨時間演化的圖案,其中時間t是被特征時間無量綱化的時間.從圖中可以看到:初始階段,擾動隨時間逐漸發展:輕流體和重流體相互滲透,并在t=1時刻形成一排“蘑菇”形狀;緊接著,一排“蘑菇”繼續增長,變得纖細而修長(如圖t=2);當t>2時,流體與流體間的剪切力變得非常大,非線性程度在加劇,造成了流體界面的混合程度在不斷地加強,最終形成了非常復雜的混沌界面結構.作為對比,我們進一步模擬了Re=1024的情形.圖11給出了Re=1024時多模RT不穩定性中密度界面隨時間的演化圖.從圖中可以看出:在初始時刻,不穩定性發展的比較緩慢,界面處于抑制狀態,未觀察到顯著的“蘑菇”形狀;接著不穩定性進入穩定期,流體混合程度在降低;而在演化后期,在兩相系統中也觀察到較為混亂的相界面圖案.最后,還模擬了Re=700的多模RT不穩定性問題.圖12給出了Re=700時多模RT不穩定性中密度界面隨時間的演化圖.從圖中可以發現初始的不穩定性發展極其緩慢,最終形成的的流體界面相對光滑,未觀察到明顯的混沌現象.

圖10 Re=4096時,多模RT不穩定性的演化 (a)t=0;(b)t=1;(c)t=2;(d)t=3;(e)t=4;(f)t=5Fig.10.Time evolution of a multiple-mode RT instability at Re=4096:(a)t=0;(b)t=1;(c)t=2;(d)t=3;(e)t=4;(f)t=5.

圖11 Re=1024時,多模RT不穩定性的演化 (a)t=0;(b)t=1;(c)t=2;(d)t=3;(e)t=4;(f)t=5Fig.11.Time evolution of a multiple-mode RT instability at Re=1024:(a)t=0;(b)t=1;(c)t=2;(d)t=3;(e)t=4;(f)t=5.

圖12 Re=700時,多模RT不穩定性的演化 (a)t=0;(b)t=1;(c)t=2;(d)t=3;(e)t=4;(f)t=5Fig.12.Time evolution of a multiple-mode RT instability at Re=700:(a)t=0;(b)t=1;(c)t=2;(d)t=3;(e)t=4;(f)t=5.

4 結 論

本文基于Cahn-Hilliard相場理論提出了一種修正的兩相流LB模型,通過引入新的簡化的外力項分布函數,解決了界面力的尺度與多尺度分析結果不一致的問題,并且通過Chapman-Enskog理論分析,可以證明本模型能夠準確地恢復到CH方程和NS方程.通過模擬一系列經典兩相流的數值算例,包括靜態液、液滴合并、亞穩態分解和RT不穩定性問題,用于測試本模型的準確性和可行性.對于靜態液滴算例,發現當前的相場LB模型能準確地求解相界面,并且可以獲得量級為10?9極小的虛假速度.還模擬了兩個液滴在表面張力作用下的合并問題,獲得了與前人文獻相一致的數值結果:當液滴間距小于兩倍的界面厚度時,液滴在表面張力的作用下自動發生合并;而當液滴間距大于兩倍界面厚度,液滴不會發生合并現象.針對亞穩態分解問題,我們利用本文相場LB模型可以成功地捕獲相界面的分離過程,并發現增加表面張力σ不影響界面的演化圖案,但可以加速相分離進展.最后,利用本文修正的兩相流相場LB模型模擬了單模的RT不穩定性問題和隨機擾動的多模RT不穩定性問題,并分析了Re數對相界面的演化規律的影響,結果表明本模型可以準確地追蹤單模RT不穩定性中尖釘與氣泡的演化規律.另外,針對多模情形,發現高Re數時,流體界面在演化前期形成一排“蘑菇”形狀,而在演化后期,流體界面則形成非常獨特的混沌拓撲結構;而當Re數較低時,流體界面則相對比較光滑,演化后期未觀察到非?；煦绲幕旌犀F象.最后,本文對界面力的分析思路可以推廣到基于Allen-Cahn相場理論的LB方法中.相比界面追蹤的Cahn-Hilliard方程,基于Allen-Cahn方程的LB方法[9,21]在求解相界面方面具有較小的數值耗散,從而可以應用于大密度比兩相流問題.