變密度流體層中氣泡浮升運動的數值模擬

李 乾， 樓映中， 賀治國

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

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

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

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

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

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

1 數值模型

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

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

1.1 多水平集函數的構建及多物質流動系統的平滑化

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

(1)

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

(2)

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

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

H(φm)=

(3)

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

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

(4)

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

(5)

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

據此，先對氣液相界面進行平滑化，而后按液體序號依此完成各液體表面的平滑化，最終實現多物質流動系統的平滑化.

(6)

(7)

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

1.2 控制方程及水平集演化

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

(8)

和動量方程

(9)

(10)

(11)

(12)

(13)

(14)

(15)

(16)

(17)

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

(18)

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

(19)

(20)

(21)

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

1.3 數值計算步驟

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

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

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

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

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

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

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

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

(22)

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

2 模擬驗證

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

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

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

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

針對該算例，選用計算網格數為110×110×220，采用無滑移邊界條件，取Δt=0.005Δx，Re=474，Fr=0.64，We=1.0，ρg=0.001ρw，μg=0.01μw.

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

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

2.2 水中氣泡上浮算例

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

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

2.3 質量守恒性分析

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

(23)

(24)

(25)

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

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

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

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

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

MF1和MF2的模擬結果分別如圖6和7所示，MF3和MF4的模擬結果分別如圖8和9所示.由圖6～9可知，無論是均勻流體層中的氣泡運動，還是變密度流體層中的氣泡運動，整個流場的運動過程均可定性地分為3個階段：水下變形階段、融合階段以及慣性階段.

圖6 MF1的模擬結果Fig.6 Simulation results of MF1

圖7 MF2的模擬結果Fig.7 Simulation results of MF2

圖8 MF3的模擬結果Fig.8 Simulation results of MF3

圖9 MF4的模擬結果Fig.9 Simulation results of MF4

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

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

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

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

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

4 結論

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

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

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