李 乾， 樓映中， 賀治國

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

氣體在液體中的運動廣泛存在于自然界和工程領域中，如發(fā)動機水下排氣、水下爆炸引起的氣泡運動、石油開采.其中的動力學過程在一個多世紀以來一直受到人們的關注[1-2].

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

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

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

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

1 數(shù)值模型

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

圖1 模型設置示意圖Fig.1 Schematic of model setting

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

首先，以氣液相界面和各液體表面構建多個水平集函數(shù)φm(X)(下文簡寫為φm).其中：m為流體序號；X為點的空間坐標，X=(x,y,z).當m=1時，流體的種類為氣體，則有

(1)

當m≥2時，流體的種類為液體，則有

(2)

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

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

H(φm)=

(3)

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

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

(4)

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

(5)

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

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

(6)

(7)

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

1.2 控制方程及水平集演化

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

(8)

和動量方程

(9)

(10)

(11)

(12)

(13)

(14)

(15)

(16)

(17)

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

(18)

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

(19)

(20)

(21)

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

1.3 數(shù)值計算步驟

MLS法的主要計算流程如下：

(1) 構建多個水平集函數(shù)；

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

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

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

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

計算中的各流體間互不混溶，并可近似為不可壓縮流體[16].在此條件下，氣泡內(nèi)的初始壓力為1個標準大氣壓(1.013×105Pa)，對應的氣體初始密度為1 kg/m3，水、氣間的密度比為 1 000∶1.

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

(22)

式中：u、v、w為空間直角坐標系的無量綱速度分量.

2 模擬驗證

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

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

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

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

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

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

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

2.2 水中氣泡上浮算例

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

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

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

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

(23)

(24)

(25)

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

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

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

氣泡在自由表面爆裂時，液氣交界面的變化十分劇烈，流場中的水氣混合過程相當復雜.為對該過程進行研究，設置了4組算例(MF1～MF4)，各算例的主要參數(shù)如表1所示.4組算例中，氣泡的半徑為1，球心為(x,y,z)=(0,0,-3)，水深區(qū)間為[-6,-2.0].4組算例的網(wǎng)格數(shù)，無量綱數(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可知，無論是均勻流體層中的氣泡運動，還是變密度流體層中的氣泡運動，整個流場的運動過程均可定性地分為3個階段：水下變形階段、融合階段以及慣性階段.

圖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

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

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

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

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

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

4 結(jié)論

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

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

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