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

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

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

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

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

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

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

1 數學模型

1.1 控制方程

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

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

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

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

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

1.2 離散方法

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

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

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

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

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

2 近場壓力邊界

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

2.1 邊界壓力算法

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

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

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

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

2.2 穩定性證明

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

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

忽略無窮小量后可得:

忽略無窮小量后可得:

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

2.3 邊界合理性判斷與數值震蕩抑制

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

3 可靠性驗證

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

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

式中:k為氣體絕熱指數，取1.4。

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

表1 數值模擬主要參數Tab.1 Parameters of numerical model

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

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

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

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

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

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

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

圖4 水下爆炸氣泡脈動過程(虛線為數值模擬結果)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 近壁氣泡脈動過程分析

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

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

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

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

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

5 結論

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

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

(2)MAC方法適于追蹤具有較大曲率的自由液面，允許人為的指定自由液面拓撲形狀的改變，在水下氣泡非球狀脈動過程的模擬中表現出色;

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

［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］宗 智，何 亮，張恩國.水中結構物附近三維爆炸氣泡的數值模擬［J］.水動力學研究與進展 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，Ferreira 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］張凌新，尹 琴，邵雪明.水中氣泡潰滅的理論與數值研究［J］.水動力學研究與進展 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.