基于歐拉單元的水下氣泡高效數(shù)值仿真算法研究

蘇怡然，唐文勇，鄭紹文，李德聰

(1.上海交通大學(xué) 海洋工程國家重點實驗室，上海 200240;2.中國艦船研究設(shè)計中心，武漢 430064)

在實際工程中，氣泡脈動引發(fā)沖擊作用的實例有很多:如水下爆炸氣態(tài)產(chǎn)物對船體結(jié)構(gòu)造成沖擊、螺旋槳工作時產(chǎn)生的毫米量級的空化氣泡對螺旋槳的剝蝕，以及超聲輻照產(chǎn)生的微米量級氣泡對微生物的沖擊作用等。當(dāng)氣泡距離結(jié)構(gòu)足夠近時，其脈動過程往往伴隨著高速水射流［1］。水射流產(chǎn)生后會穿透氣泡并對附近結(jié)構(gòu)造成沖擊，此過程為氣泡破壞效應(yīng)的主要機理之一。

關(guān)于水下氣泡的運動特性，早期研究［1－2］主要基于一些理想假設(shè)，如氣泡呈球狀變化等。近年來，隨著計算機技術(shù)的發(fā)展，氣泡非球狀演化過程的數(shù)值模擬已成為可能。邊界積分方法假設(shè)流場無粘、無旋且不可壓縮，使用低維模型模擬高維空間，有效減小了計算量，在模擬水下爆炸氣泡射流現(xiàn)象中已有許多成功的應(yīng)用［3－7］。但上述假設(shè)限制了邊界積分法的應(yīng)用范圍:如螺旋槳空泡或超聲激發(fā)的氣泡，氣泡尺寸往往較小，此時粘性對氣泡脈動的影響顯著;再如水下爆炸氣泡回彈階段，激波引起的能量耗散將對氣泡演化規(guī)律造成影響，在該過程的數(shù)值模擬中流體可壓縮性將不可忽略;此外，在射流穿透氣泡后，氣泡表面速度勢變?yōu)槎嘀岛瘮?shù)，計算將被終止，需引入環(huán)狀氣泡模型才能對該氣泡進行模擬。因此，使用更具一般性的Navier-Stokes方程(以下簡稱NS方程)模擬氣泡脈動具有實際意義。

針對水下氣泡的非球狀脈動現(xiàn)象，Popinet等［8］使用MAC法成功模擬了毫米量級氣泡的射流過程，Hao等［9］將該法拓展到三維直角坐標(biāo)系中，Johnsen、McKee等［10－11］實現(xiàn)了可壓縮流體中射流現(xiàn)象的數(shù)值模擬，張凌新等［14］使用VOF法模擬了水下爆炸氣泡的潰滅現(xiàn)象。這些方法均需要將流體邊界取在遠離氣泡的位置，導(dǎo)致計算效率低下，尤其對于一些高強度氣泡，巨大的計算耗時將使數(shù)值模擬缺乏實用性。

本文在MAC法［13］基礎(chǔ)上提出一種適用于模擬氣泡脈動問題的近場壓力邊界條件，通過合理減小計算區(qū)域，提高了MAC法的計算效率，解決了高強度氣泡數(shù)值模擬中計算耗時較大的問題。通過引進Popinet［8］和Hirt［12］等對自由液面的處理方法，保證了自由液面邊界條件的計算精度和數(shù)值穩(wěn)定。

1 數(shù)學(xué)模型

1.1 控制方程

對于水下氣泡首次周期內(nèi)的脈動問題，使用軸對稱不可壓縮NS方程進行描述:

其中:Φ=p/ρ代表壓力，ρ為液體密度，g為重力加速度，γ為運動粘性系數(shù)，剪應(yīng)力分量計算方法如下:

軸對稱模型如圖1所示，方程需滿足的邊界條件為:

對歐拉域外圍邊界Sp使用壓力邊界條件，邊界上的壓力值由遠端流場確定。在求解遠端壓力分布時，本文將脈動氣泡假想為空間源匯，使用勢流理論計算遠場壓力，詳見本文第2章。

圖1 水下氣泡數(shù)值模型Fig.1 Numerical model for the simulation of underwater bubbles

1.2 離散方法

使用交錯網(wǎng)格對NS方程進行離散，壓力節(jié)點和速度節(jié)點分別位于網(wǎng)格中心和網(wǎng)格邊長中點。NS方程在空間上采用中心差分離散格式，在時間上使用一階隱式歐拉法進行離散。

使用預(yù)測－校正法解決不可壓縮流動的壓力－速度耦合問題［15－16］:在每個時間步內(nèi)，首先忽略壓力梯度項，通過動量方程求出預(yù)測速度，進而得到壓力泊松方程，使用逐次超松弛迭代法對該方程進行求解，最后用壓力梯度對預(yù)測速度進行校正。

在處理自由液面問題時，本文采用改進MAC法的思想，通過標(biāo)記點追蹤自由液面運動，實現(xiàn)了對自由液面的精確捕捉。首先在自由液面上布置一串具有先后順序的無質(zhì)量標(biāo)記點，并在每個時步內(nèi)使標(biāo)記點跟隨流場運動，從而實現(xiàn)自由液面形狀和位置的追蹤。標(biāo)記點的運動速度由四個相鄰節(jié)點速度進行雙線性插值得出。在計算過程中，需要不斷進行標(biāo)記點位置的重分配，以保證標(biāo)記點沿液面均勻分布，從而降低計算誤差。自由液面將歐拉網(wǎng)格劃分為流體網(wǎng)格、液面網(wǎng)格和空網(wǎng)格三類，其中空網(wǎng)格不參與計算。

液面法向邊界條件見式(11)，式中σ為表面張力系數(shù)，κ為液面曲率，ΦInt為自由液面附近流場壓力，ΦB為氣泡內(nèi)壓力。傳統(tǒng)處理方法簡單地認(rèn)為邊界單元壓力等于液面壓力，這樣處理在邊界法向壓力梯度較小時可滿足精度要求，但對于本文研究的高壓力梯度問題，會造成計算誤差與數(shù)值波動。為使計算誤差得到控制，采用Hirt等［14］提出的邊界單元壓力計算方法:認(rèn)為壓力在邊界附近呈線性分布，邊界單元壓力由相鄰流體單元壓力與液面壓力插值得到。

液面切向邊界條件見式(12)。根據(jù)Popinet等［8］提出的方法實現(xiàn)該條件，即認(rèn)為自由液面運動速度是由周邊流體網(wǎng)格處流速和空網(wǎng)格虛擬速度插值得到，通過人為指定虛擬速度的大小，可使流場滿足切向邊界條件。于是，該問題轉(zhuǎn)化為條件極值問題:

2 近場壓力邊界

對無窮流域問題進行數(shù)值模擬時，普遍需要通過較大的模型范圍來逼近無窮遠邊界，計算效率較低。尤其對于高強度氣泡，脈動過程中最大半徑可達初始半徑的30余倍，巨大的計算量將使數(shù)值模擬失去實用性。本文根據(jù)Best等［17］得到的規(guī)律，提出一種由近場壓力邊界描述無窮流域的數(shù)值算法，有效縮小了模型尺寸，提高了計算效率。

2.1 邊界壓力算法

Best等［17］提出，距氣泡一定距離之外的流場可近似看作發(fā)散狀。因此，將氣泡假想為空間點源(匯)，可認(rèn)為遠端流場由該點源(匯)激發(fā)。設(shè)V為氣泡體積，Q為點源(匯)強度，r為空間點到氣泡中心的距離，遠端流場速度勢為:

當(dāng)氣泡靠近壁面時，流場會受到鏡面效應(yīng)影響。此時，需在氣泡鏡像點處施加一個等強度的點源(匯)。設(shè)r1和r2分別表示空間點至氣泡中心和鏡像位置的距離，遠端流場速度勢為:

在忽略遠端粘性影響后，流場壓力可由伯努利方程求出。設(shè)Φ∞為無窮遠處壓力(靜水壓)，近場邊界壓力ΦBC可由式(16)求出:

將ΦBC作為近場邊界的壓力值，即可減小數(shù)值模擬所需的模型尺寸。

2.2 穩(wěn)定性證明

首先，對源強度施加小量誤差δQ，該誤差引起的近場邊界壓力誤差為:

由于時間步長Δt足夠小，式中第二項可以忽略。設(shè)R與·R·分別為氣泡等效半徑R對時間的一階和二階導(dǎo)數(shù)，則有 δ(·R·)? δ (R·)? δ(R)。根據(jù)Rayleigh-Plesset方程的積分形式:

忽略無窮小量后可得:

忽略無窮小量后可得:

可見，當(dāng)近場邊界距氣泡中心距離大于兩倍氣泡最大半徑時，近場壓力邊界條件已具有穩(wěn)定性。

2.3 邊界合理性判斷與數(shù)值震蕩抑制

當(dāng)標(biāo)記點移動至新網(wǎng)格時，其速度將由新的一組節(jié)點速度插值得到，這會引起標(biāo)記點速度的微幅波動。由于時間步長很短，該波動對于氣泡尺寸的影響可以忽略，但卻會使產(chǎn)生較大誤差，從而直接影響近場邊界壓力的準(zhǔn)確性。若不加控制，還易引發(fā)數(shù)值震蕩，導(dǎo)致計算失敗。

3 可靠性驗證

氣泡的球狀演化規(guī)律由 Rayleigh-Plesset方程給出:

其中:R為氣泡等效半徑，γ為運動粘性系數(shù)，S為表面張力系數(shù)，p∞為靜水壓力，pB為氣泡壓力。氣體狀態(tài)方程為:

式中:k為氣體絕熱指數(shù)，取1.4。

從本質(zhì)上講，Rayleigh-Plesset方程是NS方程在理想條件下的簡化形式，因此，算法可靠性的最基本要求即為可逼近Rayleigh-Plesset經(jīng)典解。對忽略重力作用與壁面作用的氣泡脈動過程進行數(shù)值模擬，考慮不同強度的氣泡(設(shè)氣泡強度系數(shù)為其壓力峰值與環(huán)境壓力的比值，這里取 50，100，150，200)，模型參數(shù)見表 1，計算結(jié)果如圖2所示。在強度系數(shù)為50時，氣泡半徑誤差為0.2%;隨著強度增高，半徑誤差略微變大，這是因為在模型尺寸不變的情況下，隨著氣泡體積變大，近場邊界距氣泡距離不斷減小。計算結(jié)果表明，各項誤差均在可接受范圍內(nèi)，其中半徑誤差小于1%，脈動周期誤差約為0.5% ～0.6%。

表1 數(shù)值模擬主要參數(shù)Tab.1 Parameters of numerical model

圖2 近場邊界方法與R－P方程計算結(jié)果對比Fig.2 Result comparisons between the near-field boundary method and the R-P equation

無窮流域問題數(shù)值模擬的傳統(tǒng)方法需使用較大的模型范圍，且模型邊界條件設(shè)為恒定靜水壓力。為驗證本文所述邊界條件在提高計算效率、降低計算誤差中的效果，表2將兩種邊界條件的誤差進行了對比。對于不同模型尺寸，恒壓邊界均會產(chǎn)生更大的誤差，尤其是脈動周期誤差。隨著模型尺寸的增大，兩種算法的誤差均有所減小，其中近場壓力邊界的誤差維持在0.5%左右，而恒壓邊界的誤差則從36%快速下降，但即使模型范圍達到10倍氣泡半徑，恒壓邊界誤差依然大于近場壓力邊界誤差。這說明本文所述邊界壓力算法在提高計算效率的同時也降低了計算誤差。

表2 近場邊界與恒壓邊界誤差比較Tab.1 Error comparison between near-field boundary and constant pressure boundary

然而，上述模型中自由液面均未發(fā)生變形。為驗證算法處理液面變形問題的可靠性，本文對重力作用下的氣泡脈動過程進行數(shù)值模擬，并與實驗數(shù)據(jù)進行了對比。Shima等［18］采用脈沖放電的方法產(chǎn)生水下氣泡，并提出由氣泡平衡半徑Re來確定氣泡壓力:

對于該實驗涉及的氣泡尺寸，流體粘性效應(yīng)影響可以忽略，但由于粘性作用顯著的微米量級氣泡的脈動實驗難以進行，所以這里僅將上述實驗數(shù)據(jù)同數(shù)值仿真結(jié)果進行對比，作為證明算法可靠性的依據(jù)，計算結(jié)果如圖3所示。在首次脈動周期內(nèi)，數(shù)值模擬結(jié)果與實驗數(shù)據(jù)吻合良好，氣泡半徑平均誤差為氣泡最大半徑的3.91%，脈動周期誤差為2.68%;當(dāng)氣泡進入回彈階段時，數(shù)值模擬得到的氣泡等效半徑明顯大于實驗結(jié)果。這是因為在回彈階段，多種物理過程產(chǎn)生的能量耗散效應(yīng)逐漸顯著，如氣泡界面上熱交換、可壓縮流體激波等。

Klaseboer等［4］在實驗水池中進行了不同裝藥量的爆炸實驗，并進行高速攝影，圖4給出了爆炸氣泡脈動過程與數(shù)值模擬結(jié)果的對比，兩者在趨勢上吻合較好。從數(shù)值模擬結(jié)果可以看出，在潰滅階段，受重力的作用，射流產(chǎn)生于氣泡底端并逐漸發(fā)展至穿透氣泡，射流方向與重力方向相反。實驗中監(jiān)測了爆炸氣泡附近壓力場隨時間的變化，由近到遠三個監(jiān)測點分別記為P1、P3、P5。圖5對比了三個監(jiān)測點壓力的實驗數(shù)據(jù)和數(shù)值計算結(jié)果?？梢钥闯?，兩者趨勢一致，誤差保持在15%以內(nèi)。產(chǎn)生該誤差的主要原因是試驗水池邊界和數(shù)值模型邊界之間的差異。數(shù)值模擬中沒有考慮池壁、水面等邊界影響。此外，由于算法沒有考慮流體壓縮性，所以不同監(jiān)測點同時達到壓力峰值，而非試驗中先后達到峰值。

圖3 重力作用下氣泡等效半徑隨時間變化Fig.3 Time history of bubble equivalent radius under the influence of gravity

圖4 水下爆炸氣泡脈動過程(虛線為數(shù)值模擬結(jié)果)Fig.4 Pulsation of underwater explosive bubbles(The dash line represents numerical results)

圖5 爆炸氣泡附近壓力隨時間變化Fig.5 Time evolution of the pressure near explosive bubble

4 近壁氣泡脈動過程分析

近壁氣泡脈動過程往往伴隨著射流現(xiàn)象，由于國內(nèi)外公開的相關(guān)實驗結(jié)果多為高速攝影，定量分析較為匱乏，故本文僅對典型近壁氣泡脈動過程進行定性分析。Tomita等［19］使用脈沖放電的方法制造出近壁氣泡，并進行高速攝影，如圖6(a)所示，圖6(b)為該過程的數(shù)值模擬計算結(jié)果。數(shù)值模型中，氣泡位于7.5 m水深，強度系數(shù)取50，初始半徑為0.15 m，氣泡距壁0.7 m。

在膨脹階段，氣泡始終保持近似球形，并被固壁略微排斥，膨脹階段停止于約0.5倍脈動周期;隨后進入潰滅階段，該階段內(nèi)氣泡受固壁的強烈吸引而產(chǎn)生大幅位移和變形，在潰滅階段尾期，射流形成于氣泡遠離固壁端;在124 ms時刻射流擊穿氣泡，使氣泡變?yōu)榄h(huán)狀。

射流形成時刻流線圖與壓力云圖的計算結(jié)果如圖7所示。此時，氣泡壓力大于靜水壓力，氣泡處于減速收縮階段，且射流后方存在局部高壓區(qū)域;流場邊界處流動近似發(fā)散狀，這與Best［17］提出的流速規(guī)律是相吻合的。

圖6 近壁氣泡脈動過程Fig.6 Bubble pulsation in the vicinity of a rigid boundary

圖7 射流形成時刻流線圖與壓力云圖Fig.7 Streamlines and pressure profiles at jet formation

5 結(jié)論

本文基于NS方程，提出一種近場壓力邊界算法，有效提高了計算效率。該方法適用于水下爆炸氣泡脈動、螺旋槳空泡潰滅等多種過程的數(shù)值模擬和相應(yīng)射流沖擊載荷的預(yù)報。本文主要結(jié)論如下:

(1)本文所述近場壓力邊界允許通過較小區(qū)域的流體模型模擬無窮流域問題，大幅提高了計算效率;

(2)MAC方法適于追蹤具有較大曲率的自由液面，允許人為的指定自由液面拓?fù)湫螤畹母淖儯谒職馀莘乔驙蠲}動過程的模擬中表現(xiàn)出色;

(3)本文旨在闡明近場壓力邊界，所以文中主要模擬了氣泡首次脈動周期內(nèi)的近似不可壓過程，關(guān)于氣泡回彈時流體壓縮效應(yīng)和液面熱交換作用的影響有待進一步研究。

［1］Blake J R，Gibson D C.Cavitation bubbles near boundaries［J］.Annual Review of Fluid Mechanics，1987，19(1):99－123.

［2］Plesset M S，Prosperetti A.Bubble dynamics and cavitation［J］.Annual Review of Fluid Mechanics，1977，9(1):145－185.

［3］Zhang Y L，Yeo K S，Khoo B C，et al.3D jet impact and toroidal bubbles［J］.Journal of Computational Physics，2001，166(2):336－360.

［4］Klaseboer E，Hung K C，Wang C，et al.Experimental and numerical investigation of the dynamics of an underwater explosion bubble near a resilient/rigid structure［J］.Journal of Fluid Mechanics，2005，537:387－413.

［5］ Zhang A M，Yao X L，Yu X B.The dynamics of threedimensional underwater explosion bubble［J］.Journal of Sound and Vibration，2008，311(4):1196 －1212.

［6］張阿漫，姚熊亮.水下爆炸氣泡與海底相互作用研究［J］.振動與沖擊，2008，27(3):92 －98.ZHANG A-man， YAO Xiong-liang. Interaction between underwater explosion bubble and seabed［J］.Journal of Vibration and Shock，2008，27(3):92－98.

［7］宗 智，何 亮，張恩國.水中結(jié)構(gòu)物附近三維爆炸氣泡的數(shù)值模擬［J］.水動力學(xué)研究與進展 A輯，2007，22(5):592－602.ZONG Zhi， HE Liang， ZHANG En-guo. Numerical simulation of a three-dimensional underwater explosion bubble near a structure［J］.Journal of Hydrodynamics Ser A，2007，22(5):592－602.

［8］ Popinet S，Zaleski S.Bubble collapse near a solid boundary:a numerical study of the influence of viscosity［J］.Journal of Fluid Mechanics，2002，464:137 －163.

［9］ Hao Y，Prosperetti A.A numerical method for threedimensional gas liquid flow computations［J］.Journal of Computational Physics，2004，196(1):126 －144.

［10］ Johnsen E， Colonius T. Numerical simulations of nonspherical bubble collapse［J］.Journal of Fluid Mechanics，2009，629:231 －262.

［11］ McKee S，Tome M F，F(xiàn)erreira V G，et al.The MACmethod［J］.Computers＆ Fluids，2008，37(8):907－930.

［12］ Nichols B D，Hirt C W.Improved free surface boundary conditions for numerical incompressible-flow calculations［J］.Journal of Computational Physics，1971，8(3):434 －448.

［13］Hawker N A，Ventikos Y.Interaction of a strong shockwave with a gas bubble in a liquid medium:A numerical study［J］.Journal of Fluid Mechanics，2012，701:59 －97.

［14］張凌新，尹 琴，邵雪明.水中氣泡潰滅的理論與數(shù)值研究［J］.水動力學(xué)研究與進展 A輯，2012，27(1):68－73.ZHANG Ling-xin，YIN Qin，SHAO Xue-ming.Theoretical and numerical studies on the bubble collapse in water［J］.Chinese Journal of Hydrodynamics Ser.A，2012，27(1):68－73.

［15］ Griffith B E.An accurate and efficient method for the incompressible Navier Stokes equations using the projection method as a preconditioner［J］.Journal of Computational Physics，2009，228:7565－7595.

［16］ Brown D L，Cortez R，Minion M L.Accurate projection methods for the incompressible Navier Stokes equations［J］.Journal of Computational Physics，2001，168(2):464－499.

［17］ Best J P.The dynamics of underwater explosions［D］.PhD dissertation. The University of Wollongong， New South Wales，1991.

［18］ Shima A，Tomita Y，Sendai.The behavior of a spherical bubble near a solid wall in a compressible liquid［J］.Ingenieur-Archiv，1981，51:243－255.

［19］Tomita Y，Shima A.High-speed photographic observations of laser-induced cavitation bubbles in water［J］.Acustica，1990，71(3):161－171.