李 乾， 樓映中， 賀治國(guó)

(浙江大學(xué) 海洋學(xué)院；港口海岸與近海工程研究所，浙江 舟山 316021)

氣體在液體中的運(yùn)動(dòng)廣泛存在于自然界和工程領(lǐng)域中，如發(fā)動(dòng)機(jī)水下排氣、水下爆炸引起的氣泡運(yùn)動(dòng)、石油開(kāi)采.其中的動(dòng)力學(xué)過(guò)程在一個(gè)多世紀(jì)以來(lái)一直受到人們的關(guān)注[1-2].

數(shù)值模擬是研究氣體在液體中運(yùn)動(dòng)的一種有效方法[3-4].由于液體與氣體之間的密度比一般遠(yuǎn)大于重液體與輕液體之間的密度比，所以氣液相界面的捕捉一直是數(shù)值模擬的重點(diǎn)與難點(diǎn).目前，捕捉氣液相界面最有效且最具有代表性的方法有兩種，即流體體積(VOF)法和水平集法.VOF法由Hirt等[5]在1981年最先提出，其將供體-受體結(jié)合方程用于求解對(duì)流F函數(shù)，以確保數(shù)值解.但該方法若不進(jìn)行幾何重構(gòu)，可能會(huì)產(chǎn)生較大的數(shù)值耗散，進(jìn)而影響計(jì)算精度.

而在水平集法中，氣液兩相流中的相界面被設(shè)為零水平集.水平集函數(shù)是符號(hào)距離函數(shù)，相界面的一側(cè)符號(hào)為正，另一側(cè)符號(hào)為負(fù)，函數(shù)值的大小為空間位置與零水平集(相界面)的距離.據(jù)此，氣液界面的形狀和曲率可以較容易地通過(guò)求解水平集函數(shù)直接獲得.水平集法由Osher等[6]在1988年最先提出，并將之用于捕捉鋒面運(yùn)動(dòng).該方法在形狀拓?fù)渖暇哂休^大的優(yōu)勢(shì)，但可能伴有質(zhì)量損失.Sussman等[7]首次將水平集法運(yùn)用于計(jì)算不可壓縮二相流.在此之后，水平集法在計(jì)算流體力學(xué)中的應(yīng)用得到了很好的發(fā)展.Adalsteinsson等[8]提出水平集函數(shù)可以由一個(gè)擴(kuò)展速度推進(jìn)，而非流體速度.Jiang等[9]提出用高精度的加權(quán)本質(zhì)無(wú)振蕩(WENO)有限差分格式來(lái)求解流場(chǎng)可有效減少質(zhì)量損失.現(xiàn)階段的水平集法已經(jīng)能夠較好地捕捉變形劇烈的流體交界面[10-12].針對(duì)多相問(wèn)題，Starinshak等[12]提出了當(dāng)計(jì)算域中含多種材料時(shí)水平集函數(shù)的處理方法，但應(yīng)用該方法時(shí)，一旦各水平集函數(shù)內(nèi)部(正數(shù)部分)出現(xiàn)重疊，可能在重疊處出現(xiàn)嚴(yán)重的質(zhì)量損失問(wèn)題，甚至導(dǎo)致程序的崩潰.

當(dāng)氣體在層化的海洋環(huán)境等非均勻流體中運(yùn)動(dòng)時(shí)，由于周圍液體的密度不恒定，其運(yùn)動(dòng)過(guò)程應(yīng)當(dāng)被視為一種在變密度流體層中的運(yùn)動(dòng).此時(shí)，氣泡周圍的變密度流體層可能在此種復(fù)雜運(yùn)動(dòng)過(guò)程中扮演著十分重要的作用.

近年來(lái)，已有不少關(guān)于氣泡動(dòng)力學(xué)的研究論文，如梅登飛等[13]的黏性與非黏性顆粒在流化床中的氣泡行為模擬.然而，雖然有一些氣液耦合的研究，如李少白等[14]的非牛頓流體中線雙氣泡相互作用的數(shù)值模擬，但以上研究主要考慮兩種流體間的相互作用，并沒(méi)有考慮多種液體的存在，或液體內(nèi)部間的物質(zhì)交換對(duì)氣泡上升運(yùn)動(dòng)的影響.此外，目前已有一些研究分析了流體黏性對(duì)氣泡運(yùn)動(dòng)的影響[2]，但關(guān)于氣泡在變密度流體層中浮升運(yùn)動(dòng)的研究仍較少，尤其缺乏對(duì)氣體運(yùn)動(dòng)穿過(guò)多種液體并進(jìn)入大氣環(huán)境這種復(fù)雜過(guò)程下的氣泡浮升運(yùn)動(dòng)精細(xì)結(jié)構(gòu)的捕捉.因此，為實(shí)現(xiàn)氣泡在變密度流體層中浮升運(yùn)動(dòng)的有效模擬，并考慮到變密度流體層可以等效為若干雙層流體的疊加，本文開(kāi)發(fā)了一種耦合氣體和兩層不同密度液體的直接模擬模型，對(duì)氣泡在變密度流體層中的運(yùn)動(dòng)進(jìn)行計(jì)算分析.基于5階加權(quán)本質(zhì)無(wú)振蕩和3階Runge-Kutta高精度格式[4]提出一種多水平集函數(shù)耦合數(shù)值(MLS)法， 用以捕捉氣體表面以及各液體表面；使用投影法計(jì)算壓力，并通過(guò)求解多物質(zhì)流動(dòng)系統(tǒng)的平滑化方程計(jì)算交界面處的密度和黏度.通過(guò)與已有文獻(xiàn)結(jié)果的比較，驗(yàn)證了所建立的MLS法的有效性.采用該模型對(duì)水面薄油層的兩層流體中的氣泡運(yùn)動(dòng)進(jìn)行研究，實(shí)現(xiàn)了對(duì)復(fù)雜條件下氣泡浮升運(yùn)動(dòng)精細(xì)結(jié)構(gòu)的有效捕捉.

1 數(shù)值模型

耦合氣泡和不同密度流體的數(shù)值模型在計(jì)算中采用水和油作為典型變密度流體層，并將水作為底層流體，油作為上層流體，油的上方則為空氣.模型設(shè)置如圖1所示.其中：d1為氣泡上方水層的厚度；d2為油層的厚度.

圖1 模型設(shè)置示意圖Fig.1 Schematic of model setting

1.1 多水平集函數(shù)的構(gòu)建及多物質(zhì)流動(dòng)系統(tǒng)的平滑化

首先，以氣液相界面和各液體表面構(gòu)建多個(gè)水平集函數(shù)φm(X)(下文簡(jiǎn)寫(xiě)為φm).其中：m為流體序號(hào)；X為點(diǎn)的空間坐標(biāo)，X=(x,y,z).當(dāng)m=1時(shí)，流體的種類為氣體，則有

(1)

當(dāng)m≥2時(shí)，流體的種類為液體，則有

(2)

式中：Ωm為m號(hào)流體所在的區(qū)域；ψm為流體的表面微元集，即流體界面的水平零集；s為坐標(biāo)點(diǎn)到流體界面的距離.

由于對(duì)控制方程進(jìn)行離散化后，交界面處微元的密度和黏度不再能夠簡(jiǎn)單地考慮為某一流體的密度和黏度，所以需要進(jìn)行平滑化處理.因此，引入Heaviside函數(shù)H(φm)，則有

H(φm)=

(3)

式中：ε為界面厚度，即平滑層厚度的一半.引入H(φm)后，便可實(shí)現(xiàn)兩相流的平滑化，如水氣兩相流的密度及黏度可以表示為

ρ(X)=ρg+(ρw-ρg)H(φwg)

(4)

μ(X)=μg+(μw-μg)H(φwg)

(5)

式中：ρw為水的密度；ρg為氣體的密度；μw為水的黏度；μg為氣體的黏度；φwg為以水氣交界面為零水平集且水相內(nèi)符號(hào)為正的水平集函數(shù).

據(jù)此，先對(duì)氣液相界面進(jìn)行平滑化，而后按液體序號(hào)依此完成各液體表面的平滑化，最終實(shí)現(xiàn)多物質(zhì)流動(dòng)系統(tǒng)的平滑化.

(6)

(7)

式中：n為流體種類總數(shù)；ρm、μm分別為m號(hào)流體的密度和黏度.

1.2 控制方程及水平集演化

采用三維不可壓縮Navier-Stokes方程組作為控制方程.該方程組包括連續(xù)性方程

(8)

和動(dòng)量方程

(9)

(10)

(11)

(12)

(13)

(14)

(15)

(16)

(17)

式中：φ為以氣液相界面為零水平集的水平集函數(shù)，等效于φm=1；κ(φ)為曲率；δ(φ)為Diracδ函數(shù).使用投影法求解壓力項(xiàng)和下一時(shí)刻的速度.根據(jù)水平集演化方程[4]得到

(18)

以及加入質(zhì)量修正的重距離化方程[4]

(19)

(20)

(21)

求解下一時(shí)刻的水平集函數(shù)，并依此計(jì)算下一時(shí)刻流場(chǎng)的密度和黏度分布.式中：τ為人工時(shí)間；φm0為τ=0時(shí)即求解式(18)時(shí)的φm；S(φm0)為平滑符號(hào)函數(shù)；λ為質(zhì)量補(bǔ)償系數(shù)；Ω為計(jì)算域.

1.3 數(shù)值計(jì)算步驟

MLS法的主要計(jì)算流程如下：

(1) 構(gòu)建多個(gè)水平集函數(shù)；

(2) 耦合所有水平集函數(shù)，求解整個(gè)流場(chǎng)的密度及黏度分布；

(3) 根據(jù)整個(gè)流場(chǎng)的密度及黏度分布和現(xiàn)有速度場(chǎng)，計(jì)算壓力場(chǎng)和下一時(shí)刻的速度場(chǎng)；

(4) 根據(jù)新的速度場(chǎng)，計(jì)算新的各水平集函數(shù)；

(5) 重復(fù)步驟(2)～(4)直到模擬時(shí)間結(jié)束.

計(jì)算中的各流體間互不混溶，并可近似為不可壓縮流體[16].在此條件下，氣泡內(nèi)的初始?jí)毫?個(gè)標(biāo)準(zhǔn)大氣壓(1.013×105Pa)，對(duì)應(yīng)的氣體初始密度為1 kg/m3，水、氣間的密度比為 1 000∶1.

計(jì)算采用結(jié)構(gòu)網(wǎng)格，通過(guò)有限差分法對(duì)控制方程中的各項(xiàng)進(jìn)行離散.同時(shí)，對(duì)于水平集法，采用具有5階精度的WENO格式對(duì)空間項(xiàng)進(jìn)行離散，并采用具有3階精度的Runge-Kutta格式對(duì)時(shí)間項(xiàng)進(jìn)行離散.計(jì)算中需要滿足以下Courant-Friedrich-Levy(CFL)條件，即可保證模型計(jì)算時(shí)的收斂性[4]：

(22)

式中：u、v、w為空間直角坐標(biāo)系的無(wú)量綱速度分量.

2 模擬驗(yàn)證

通過(guò)與文獻(xiàn)的模擬氣泡在自由表面爆裂和水中氣泡上浮結(jié)果的對(duì)比，以驗(yàn)證MLS法的有效性.

2.1 氣泡在自由表面爆裂算例

以Gu等[4]的氣泡在自由表面爆裂算例作為驗(yàn)證算例1(C1)，該算例為在[-3,3]×[-3,3]×[-6,6]的計(jì)算域中，設(shè)置一個(gè)半徑為1，球心坐標(biāo) (x,y,z)=(0,0,-3.2)的氣泡，水深區(qū)間為[-6,-2.0]，水的上方為氣體.C1示意圖如圖2所示.

圖2 C1示意圖Fig.2 Schematic of C1

針對(duì)該算例，選用計(jì)算網(wǎng)格數(shù)為110×110×220，采用無(wú)滑移邊界條件，取Δt=0.005Δx，Re=474，F(xiàn)r=0.64，We=1.0，ρg=0.001ρw，μg=0.01μw.

本文模型(M1)的計(jì)算結(jié)果與C1的模擬結(jié)果對(duì)比如圖3所示.由圖3可知，M1與C1的結(jié)果一致表明整個(gè)流場(chǎng)的演化過(guò)程可定性地分為3個(gè)階段.第1階段為水下變形階段，由于氣泡各處的壓力不同，氣泡會(huì)在上升的同時(shí)發(fā)生形變；因底部受到的壓力較大，表面張力和黏性力難以維持其原有形態(tài)，氣泡底部將發(fā)生凹陷.第2階段為融合階段，氣泡在浮力的作用下，穿透液體表面，在液體中部形成開(kāi)口；此后隨氣泡繼續(xù)上升，開(kāi)口將進(jìn)一步擴(kuò)大.第3階段為慣性階段，此時(shí)氣泡完全與上方氣體融合，伴隨上升的液體形成液柱結(jié)構(gòu)；后續(xù)在重力和黏性力的主導(dǎo)作用下，由于液柱內(nèi)各處速度不同，液柱進(jìn)一步發(fā)生分離現(xiàn)象.另一方面，對(duì)比結(jié)果表明，M1精確地捕捉到了液柱二次分離過(guò)程中產(chǎn)生的細(xì)小流體結(jié)構(gòu)，具有更好的連續(xù)性和質(zhì)量守恒性.

圖3 水氣交界面的演變過(guò)程Fig.3 Evolution process of water-air interface

2.2 水中氣泡上浮算例

以Premlata等[2]的水中氣泡上浮算例作為驗(yàn)證算例2(C2)，該算例的參數(shù)為Re=100，F(xiàn)r=1.0，We=200，ρg=0.001ρw，μg=0.01μw.針對(duì)該算例選用的計(jì)算網(wǎng)格數(shù)為100×100×200，采用無(wú)滑移邊界條件，并取Δt=0.005Δx.本文模型(M2)的計(jì)算結(jié)果與C2的模擬結(jié)果的對(duì)比如圖4所示.由圖4可知，直到t=0.4時(shí)，氣泡都幾乎是球形的；當(dāng)t>0.4時(shí)，氣泡底部持續(xù)發(fā)生凹陷；當(dāng)t=1.6時(shí)，氣泡發(fā)生顯著的變形.此后，氣泡繼續(xù)向上移動(dòng)，形成環(huán)形的“甜甜圈狀結(jié)構(gòu)”.M2與和C2所獲得的氣泡形態(tài)演變過(guò)程的模擬結(jié)果具有良好的一致性.

圖4 氣泡形狀演變過(guò)程Fig.4 Time evolution of bubble shapes

2.3 質(zhì)量守恒性分析

為進(jìn)一步精確地檢驗(yàn)MLS法的有效性，計(jì)算每一時(shí)刻的相對(duì)質(zhì)量變化率ηmc，則有

(23)

(24)

(25)

式中：mc為計(jì)算質(zhì)量.M1和M2中的ηmc隨時(shí)間的變化如圖5所示.氣泡在自由表面爆裂時(shí)，水氣相界面的變化明顯比在水中上浮時(shí)水氣相界面的變化更劇烈.這一點(diǎn)與前者的質(zhì)量損失較大相符合.同時(shí)，兩算例中的最大相對(duì)質(zhì)量變化率均小于0.6%[4]，說(shuō)明了MLS法已具有良好的質(zhì)量守恒性，已能夠以較小的質(zhì)量損失較為精確地捕捉氣泡的運(yùn)動(dòng).

圖5 相對(duì)質(zhì)量變化率Fig.5 Relative change rate of mass

3 變密度層中氣泡浮升運(yùn)動(dòng)的數(shù)值模擬

氣泡在自由表面爆裂時(shí)，液氣交界面的變化十分劇烈，流場(chǎng)中的水氣混合過(guò)程相當(dāng)復(fù)雜.為對(duì)該過(guò)程進(jìn)行研究，設(shè)置了4組算例(MF1～MF4)，各算例的主要參數(shù)如表1所示.4組算例中，氣泡的半徑為1，球心為(x,y,z)=(0,0,-3)，水深區(qū)間為[-6,-2.0].4組算例的網(wǎng)格數(shù)，無(wú)量綱數(shù)以及離散方法等均與C1相同.在MF2和MF3中，取油的密度ρo=0.8ρw及黏度μo=2.0μw.

表1 各算例主要參數(shù)表Tab.1 Main parameters in simulation cases

MF1和MF2的模擬結(jié)果分別如圖6和7所示，MF3和MF4的模擬結(jié)果分別如圖8和9所示.由圖6～9可知，無(wú)論是均勻流體層中的氣泡運(yùn)動(dòng)，還是變密度流體層中的氣泡運(yùn)動(dòng)，整個(gè)流場(chǎng)的運(yùn)動(dòng)過(guò)程均可定性地分為3個(gè)階段：水下變形階段、融合階段以及慣性階段.

圖6 MF1的模擬結(jié)果Fig.6 Simulation results of MF1

圖7 MF2的模擬結(jié)果Fig.7 Simulation results of MF2

圖8 MF3的模擬結(jié)果Fig.8 Simulation results of MF3

圖9 MF4的模擬結(jié)果Fig.9 Simulation results of MF4

對(duì)比圖6(b)與圖8(b)、圖7(b)與圖9(b)可知，在水下變形階段時(shí)，增大氣泡上方液體層的厚度會(huì)增加氣泡處于該階段的時(shí)間.例如，當(dāng)氣泡上方的液體層厚度d(d=d1+d2)從0.2增大到0.4時(shí)，該階段的持續(xù)時(shí)間約從0.3增大到0.6.同時(shí)，厚度的增加會(huì)改變氣泡與上方氣體融合時(shí)的形狀和運(yùn)動(dòng)狀態(tài)，進(jìn)而影響到后續(xù)兩個(gè)階段的變化.另外，氣泡上方液體層的種類對(duì)水下變形階段的氣泡形態(tài)影響較小，這主要是由于油層的厚度較薄，變密度層中的氣泡和均勻流體中的氣泡在液體下所受到的壓力基本相似.

在融合階段，當(dāng)氣泡在上方較厚的液體層作用下，產(chǎn)生了較大的形變而變得足夠扁平時(shí)，氣泡上下的壓力差將會(huì)減少，而壓力作用的面積則會(huì)增大，使得壓差液柱的橫截面積增加.另一方面，雖然厚度相同的均勻水體與變密度流體層中的氣泡運(yùn)動(dòng)具有相似的發(fā)展過(guò)程，但是由于變密度層的非均勻性，各流體間的相互作用更為復(fù)雜.具體而言在該階段中，油和水體的上部分會(huì)因受氣泡給予的推力而向兩側(cè)邊界運(yùn)動(dòng)，而其下方的水體則會(huì)因受壓力而向上運(yùn)動(dòng)形成壓差液柱.之后，在四周邊界的約束下，部分液體回流.由于油與水互不混溶，油體將主要在水體的表面流動(dòng).

在慣性階段，壓差液柱的發(fā)展是流場(chǎng)變化的主要特點(diǎn).此時(shí)，流場(chǎng)主要受慣性、黏性力以及重力的控制.在該階段早期，雖然受到向下的有效重力，但在慣性的作用下，壓差液柱將繼續(xù)上升.然而，相較于均勻流體層中的運(yùn)動(dòng)，由于油層和水體是互不混溶的，當(dāng)回流的流量較大時(shí)，油層回流會(huì)沖擊壓差液柱的底部，對(duì)壓差液柱產(chǎn)生剪切和擠壓的作用，使得黏性力難以保持壓差液柱和下方水體的連續(xù)，最終產(chǎn)生分離現(xiàn)象.

4個(gè)算例(MF1～MF4)中液體能達(dá)到的最大高度隨著時(shí)間變化的結(jié)果如圖10所示，其中zh為液體最高點(diǎn)在z軸方向上的無(wú)量綱坐標(biāo)值.由圖10可知，當(dāng)d=0.2時(shí)，MF1和MF2中水下變形階段的時(shí)間均約為0.3，同時(shí)最大zh均大于3.0.在隨后的融合階段中，由于氣泡兩側(cè)的液體向氣泡底部流動(dòng)，液面的最大高度會(huì)有所下降.在慣性階段中，由于氣泡已經(jīng)與上方氣體融合，液體內(nèi)部間的相互作用成為影響液氣界面的主要因素.因此，MF1與MF2中液氣分布的區(qū)別將進(jìn)一步增大，兩者zh的最大差值可達(dá)0.5.另一方面，當(dāng)d=0.4時(shí)，由于氣泡處于水下變形階段的時(shí)間較長(zhǎng)，氣泡形態(tài)趨于扁平，氣泡上下兩側(cè)的壓力差較小，故最大zh將比d=0.2時(shí)減小約2.8.

圖10 zh演化圖Fig.10 Time evolutions of zh

4 結(jié)論

本文針對(duì)氣泡在變密度流體層中的運(yùn)動(dòng)開(kāi)發(fā)了一種耦合氣體和兩層不同密度液體的直接模擬模型，并基于水平集法提出了一種求解多水平集函數(shù)耦合的數(shù)值方法.通過(guò)與已有數(shù)值模擬結(jié)果的對(duì)比驗(yàn)證了MLS法的有效性.采用該模型對(duì)水、油兩層流體中氣泡的運(yùn)動(dòng)進(jìn)行了研究，有效地捕捉了復(fù)雜條件下氣泡浮升運(yùn)動(dòng)的精細(xì)結(jié)構(gòu).研究結(jié)果表明，當(dāng)液體和液體中的氣泡可視為不可壓縮流體，且不考慮各液體間的擴(kuò)散作用時(shí)，主要結(jié)論如下.

(1) 若邊界對(duì)流體的運(yùn)動(dòng)具有一定的約束，相較于均勻密度流體層中的運(yùn)動(dòng)，由于邊界的約束，油層將會(huì)回流，從而對(duì)水體產(chǎn)生剪切擠壓的作用，使得氣泡在變密度流體層中運(yùn)動(dòng)時(shí)，其下方的壓差液柱將更容易在底部發(fā)生斷裂分離.

(2) 增大氣泡上方初始液體厚度會(huì)使氣泡處于水下變形階段的時(shí)間更長(zhǎng)，同時(shí)氣泡底部的形態(tài)也將更為扁平.