耦合界面力的兩相流相場(chǎng)格子Boltzmann模型?

李洋 蘇婷 梁宏 徐江榮

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

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

1 引 言

格子Boltzmann(LB)方法[1,2]是近三十年來(lái)發(fā)展起來(lái)的一種流體系統(tǒng)建模和模擬方法.它不再基于宏觀連續(xù)介質(zhì)模型的Navier-Stokes方程,而是直接從微觀模型出發(fā)通過描述流體粒子分布函數(shù)的演化再現(xiàn)復(fù)雜流動(dòng)的宏觀行為.因此,LB方法相比傳統(tǒng)數(shù)值方法有一些獨(dú)特優(yōu)勢(shì),而受到了眾多學(xué)者的關(guān)注,并已廣泛應(yīng)用于多個(gè)學(xué)科領(lǐng)域,如多孔介質(zhì)中傳熱傳質(zhì)、湍流的直接數(shù)值模擬、多相流體流動(dòng)與傳熱、顆粒流[1,2]等.特別是LB方法的介觀屬性,可以使其直接描述多相流體系統(tǒng)中流體間與流固間的相互作用,因而可以方便地模擬多相流體流動(dòng),這也被認(rèn)為L(zhǎng)B方法區(qū)別傳統(tǒng)數(shù)值方法的最大優(yōu)勢(shì)之一.

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

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

2 控制方程和數(shù)學(xué)模型

2.1 控制方程

在相場(chǎng)理論中,用于描述多相系統(tǒng)中相界面運(yùn)動(dòng)的宏觀控制方程是著名的Cahn-Hilliard(CH)對(duì)流擴(kuò)散方程[6,7],

這里,?是標(biāo)識(shí)不同流體的序參數(shù),u是流體速度,M是遷移率,μ是自由能的化學(xué)勢(shì),在相場(chǎng)理論中可以表示為

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

另外,為了描述多相流體流動(dòng),界面追蹤C(jī)H方程需要耦合流體動(dòng)力學(xué)方程.根據(jù)文獻(xiàn)[7],含有外力項(xiàng)的不可壓縮黏性流體的Navier-Stokes(NS)方程可表述為:

其中,ρ是流體密度,p是壓力,ν是流體黏性系數(shù),F為總外力,包含外力項(xiàng)G和表面張力Fs,在本文的數(shù)值模擬中,我們?nèi)〕Ｓ玫膭?shì)形式Fs=???μ或Fs=μ??[6,7,15]來(lái)減少流體界面處的虛假速度.

2.2 兩相流格子Boltzmann模型

為了求解CH方程和NS方程,本文引入了兩個(gè)獨(dú)立的分布函數(shù)fi和gi,其對(duì)應(yīng)單松弛LB模型的演化方程可分別表示為:

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

這里

其中η是調(diào)節(jié)遷移率大小的參數(shù),ωi是權(quán)系數(shù),cs為格子聲速.格子速度的離散模型決定著η和ωi的取值,對(duì)于D2Q9模型,相應(yīng)的權(quán)系數(shù)√ω0=4/9,ω1?4=1/9, ω5?8=1/36,cs=離散速度ci為

其中,c= δx/δt,δx是單位格子大小,δt為時(shí)間步長(zhǎng).本文中,我們?nèi)=δx=δt=1.

為了準(zhǔn)確恢復(fù)到CH方程,Liang等[7]在演化方程引入了可以局部計(jì)算的時(shí)間導(dǎo)數(shù)項(xiàng),其對(duì)應(yīng)的源項(xiàng)分布函數(shù)可定義為

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

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

本模型中序參數(shù)?的計(jì)算如下:

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

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

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

接下來(lái),我們簡(jiǎn)要地論述本文提出新的外力項(xiàng)分布函數(shù)的LB演化方程(6b)如何恢復(fù)到不可壓流體的Navier-Stokes方程.根據(jù)外力項(xiàng)分布函數(shù)的表達(dá)式(14),其對(duì)應(yīng)的矩條件為

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

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

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

進(jìn)一步,應(yīng)用多尺度展開式到(17)和(24)式,可得分布函數(shù)gi滿足如下性質(zhì)

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

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

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

對(duì)t1和t2尺度上的方程進(jìn)行尺度黏合,最終發(fā)現(xiàn)可以正確地恢復(fù)到不可壓縮黏性流體的Navier-Stokes方程組(5a)和(5b),并且運(yùn)動(dòng)學(xué)黏性系數(shù)ν(23b)式和各階矩條件,經(jīng)系列計(jì)算最終可得到為

在古詩(shī)文的詞匯中，實(shí)詞的意義相對(duì)固定，隨著時(shí)代的更迭，會(huì)產(chǎn)生相應(yīng)的拓展和引申，但總的來(lái)說學(xué)生易于掌握，難的莫過于虛詞。劉彥和將常見的虛詞作了如下歸類：“至于夫惟蓋故者，發(fā)端之首唱；之而于以者，乃扎句之舊體；乎哉矣也，亦送末之常科。”并指明虛詞在作者回環(huán)婉轉(zhuǎn)的巧妙運(yùn)用下，文辭就會(huì)變得嚴(yán)密。因此，作者在創(chuàng)作過程中會(huì)努力避免虛詞的運(yùn)用出現(xiàn)謬誤，那么我們?cè)陂喿x過程中同樣要將虛詞理解準(zhǔn)確。否則，即使理解了文意也是不嚴(yán)密的。當(dāng)然，在教學(xué)實(shí)踐中不能完全按劉彥和的方法來(lái)辨析，畢竟《文心雕龍》成書于南朝，此后的詞匯在詞性和詞義上都產(chǎn)生了變化與演進(jìn)。

在數(shù)值模擬中,需要對(duì)時(shí)間導(dǎo)數(shù)項(xiàng)和空間導(dǎo)數(shù)項(xiàng)進(jìn)行離散計(jì)算.本文中選用顯式的歐拉差分格式來(lái)計(jì)算時(shí)間導(dǎo)數(shù)項(xiàng),

針對(duì)空間導(dǎo)數(shù)項(xiàng),為了保證多相系統(tǒng)的全局質(zhì)量守恒,我們采用如下的各向同性的二階中心差分格式[14],

其中χ表示任意變量.

3 數(shù)值結(jié)果和討論

在本節(jié)中,首先通過模擬幾個(gè)典型的兩相流問題,包括靜態(tài)液滴、液滴合并、亞穩(wěn)態(tài)分解和單模瑞利-泰勒不穩(wěn)定性問題,來(lái)驗(yàn)證所提出的兩相流相場(chǎng)LB模型的準(zhǔn)確性和穩(wěn)定性.接著,利用本模型來(lái)研究隨機(jī)擾動(dòng)分布的多模瑞利-泰勒不穩(wěn)定性問題,著重考察雷諾數(shù)對(duì)相界面的影響規(guī)律.

3.1 靜態(tài)液滴

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

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

圖1 靜態(tài)液滴測(cè)試 (a)平衡狀態(tài)下的整個(gè)區(qū)域的速度分布圖,實(shí)線和虛線分別代表液滴的初始形狀和在平衡狀態(tài)下的形狀;(b)沿著x=Ny/2的密度分布的LB模擬結(jié)果和解析解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 兩個(gè)液滴的合并過程(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 兩個(gè)液滴的合并過程(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 亞穩(wěn)態(tài)分解

進(jìn)一步,我們模擬了更為復(fù)雜的相分離問題[6,17]來(lái)驗(yàn)證本文提出的兩相流相場(chǎng)LB模型.亞穩(wěn)態(tài)分解,也稱相或組分分離現(xiàn)象,是混合流體的一種基本特征.對(duì)處于亞穩(wěn)態(tài)狀態(tài)的均勻混合流體施加一個(gè)很小的擾動(dòng),此時(shí)流體系統(tǒng)是不穩(wěn)定的,相分離現(xiàn)象將會(huì)發(fā)生.在我們的模擬過程中,計(jì)算區(qū)域設(shè)置為NX×NY=100×100,四周采用周期性邊界條件.初始的序參數(shù)分布設(shè)定為

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

3.4 瑞利-泰勒不穩(wěn)定性

圖4 σ=0.001時(shí)相分離的時(shí)間演化過程 (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時(shí)相分離的時(shí)間演化過程 (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.

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

我們考慮一個(gè)矩形微管道,其高度L和寬度W的比例為4,對(duì)兩相界面處施加一個(gè)微小擾動(dòng),其擾動(dòng)函數(shù)給定為

初始的序參數(shù)分布則設(shè)定為

雷諾數(shù)(Re)和Atoods數(shù)(At)是影響RT不穩(wěn)定性的兩個(gè)重要無(wú)量綱參數(shù),分別定義為

其中,g是重力加速度,為了耦合重力效應(yīng),在模擬中對(duì)計(jì)算區(qū)域的所有流體點(diǎn)施加一個(gè)重力,其定義為

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

圖6 Re=256時(shí),本模型計(jì)算獲取的RT不穩(wěn)定性演化圖 (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時(shí),本模型計(jì)算獲取的RT不穩(wěn)定性演化圖 (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時(shí),Li等[15]模型計(jì)算獲取的RT不穩(wěn)定性演化圖 (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數(shù)下隨時(shí)間演化的(a)氣泡和(b)尖釘在y方向上的位置Fig.9.Time evolution of y coordinate of(a)bubble and(b)spike poistions at different Re.

通過文獻(xiàn)調(diào)研表明,對(duì)RT不穩(wěn)定性問題的研究絕大多數(shù)為單模情形,針對(duì)隨機(jī)擾動(dòng)的多模RT不穩(wěn)定性的研究則相對(duì)較少[20].本文的最后,我們將利用修正的兩相流相場(chǎng)LB模型對(duì)隨機(jī)擾動(dòng)的多模RT不穩(wěn)定性進(jìn)行模擬,并研究Re數(shù)對(duì)相界面演化的影響規(guī)律.物理問題設(shè)定為一個(gè)L×L的計(jì)算區(qū)域,對(duì)兩相界面處施加一個(gè)初始的微小擾動(dòng),其中擾動(dòng)函數(shù)給定為

其中,an和bn是滿足高斯正態(tài)分布的隨機(jī)數(shù),序參數(shù)的初始分布則給定為

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

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

圖10 Re=4096時(shí),多模RT不穩(wěn)定性的演化 (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時(shí),多模RT不穩(wěn)定性的演化 (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時(shí),多模RT不穩(wěn)定性的演化 (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 結(jié) 論

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