近自由面通氣空泡誘導的飛濺水層閉合行為數值模擬

王廣航，王靜竹，杜 巖，王志英，王一偉，3

(1.中國科學院 力學研究所，北京 100190；2.中國科學院大學 未來技術學院，北京 100190；3.中國科學院大學 工程科學學院，北京 100190；4.廣東空天科技研究院，廣州 511458)

0 引 言

空泡與自由液面之間的相互作用是空泡動力學的經典課題之一，最早應用于研究水下爆炸，主要關注的是空泡非球形潰滅行為[1-8]。無量綱距離參數γ定義為空泡中心距離自由面之間的距離h和其在無限水域振蕩的最大半徑Rmax的比值，是影響空泡與自由液面相互作用的最重要參數。隨著 γ逐漸減小，空泡潰滅誘導的自由面變形主要表現為丘型、沖天型、皇冠型、濺射型、破碎型等[9-10]。前四種形式中，自由面與空泡通常不會發生界面的直接相互作用，表現為以界面間液態水為間接媒介的弱耦合作用。當空泡和自由面距離近的情況下的破碎型行為演化更為復雜，空泡膨脹誘導自由面向上飛濺形成透明水層，二者界面交迭形成氣流熱量和質量交換的通路。與此同時，通路內高壓氣體排出，壓力逐漸下降，外部氣流沿通路侵入，導致通路內部的壓力劇烈變化，驅動水層產生閉合等復雜行為。因此，空泡與自由面強耦合作用特點主要體現在空泡通氣、自由面飛濺及閉合等行為。

近年來，強耦合作用下自由面大變形與飛濺水層演化的研究也得到更多關注。Wang[11]發現，空泡通氣過程是由于空泡壁面發生了Rayleigh-Taylor不穩定，并根據飛濺水層是否閉合給出了兩種流態分類，與實驗結果吻合較好。Li[12]采用實驗和仿真方法研究了空泡與自由面相互作用，并給出了三種空泡破碎形態。Tian[13]基于歐拉有限元方法和流體體積法（volume of fluid, VOF）探究了近自由面空泡的通氣和飛濺水層閉合過程，發現飛濺水層閉合后產生的向上的射流比空泡不破裂時產生的自由面射流演化更高。Rosselló[14]采用伯努利方程建立模型模擬了飛濺水層閉合過程，將空泡膨脹過程視為圓柱形活塞入水過程，通過空泡膨脹速度間接得到空泡內與外界壓強差，發現飛濺水層閉合主要是由于空腔內外壓差導致的，表面張力對飛濺水層演化過程的影響是可以忽略的，他們還給出了空泡在 γ = 0.075～0.25之間會產生“子彈”狀射流。Marston[15]實驗研究了不同環境壓力的飛濺水層演化，指出空泡的快速膨脹所產生的吸力是導致飛濺水層閉合的原因?；谏鲜鲅芯?，空泡破碎與通氣、自由面飛濺與閉合等過程具有強非線性、非定常等特點。

空泡誘導自由面液面飛濺與閉合現象，與結構物入水[16-18]產生的表面閉合現象十分相似。Aristoff[19]發現，低邦德數（Bond number）下小球入水產生的空腔的潰滅過程主要是由于表面張力驅動的，并提出了高邦德數下耦合空腔動力學的飛濺水層演化的理論模型。Eshraghi[20]基于實驗與理論方法，提出了一套描述飛濺水層軌跡的物理模型，探究了水層內外壓差、表面張力與重力等對水層閉合的影響。Wang[21]通過理論與實驗方法對入水過程的表面閉合機制進行了研究，討論了重力、阻力、表面張力與空氣動壓對表面閉合水層的影響，發現空氣動壓主導了飛濺水層的表面閉合過程。上述入水誘導的飛濺水層建模方法、閉合機制討論及主控參數影響規律等研究過程，對本文空泡和自由面之間強耦合作用下自由面飛濺閉合過程的水層建模、受力分析、主控參數確定提供了重要參考。

以上主要通過定性實驗觀測解釋了飛濺水層閉合過程，對于其主要控制參數缺少定量測量，本文基于數值模擬獲得的飛濺水層圍成的通路內的速度場和壓力場信息，直接提取外部氣體侵入空泡內部導致的壓強差變化，并對水層受力進行分析，直接計算了各項作用力隨時間的演化，定量揭示了空泡與自由面強耦合作用誘導的水層閉合機制，并獲得了水層閉合時間與環境壓力的冪次律關系。

1 數值模擬方法

本文的數值模擬參考Zeng[22]和 Wang[23-24]對液滴內空泡動力學行為的研究，空泡初始化為一個振蕩的空化泡。軸對稱模型如圖1（a）所示，模型中采用柱坐標系O-rθz表達，原點O位于初始空泡中心，r軸沿著徑向方向，向外側為正，z軸代表軸向方向，向上為正。計算域長為60 mm、高為120 mm。邊界條件如圖1（b）所示，最右側為軸對稱邊界，其他為出口邊界。網格劃分如圖1（c）所示，采用正交笛卡爾網格進行劃分。對空泡和自由面附近進行局部加密，初始空泡在徑向方向包含60個網格，局部加密區域網格尺寸為5 μm，總網格量為384萬。

圖1 計算模型Fig.1 Computational modelling

本文基于OpenFOAM平臺，采用有限體積法直接求解N-S方程，并采用VOF來捕捉氣液界面[25-27]?；趬毫?隱式分裂算子（pressure-implicit with splitting of operators， PISO） 的多相可壓縮求解器（compressible-InterFoam） 模擬空泡和自由面之間的強耦合作用。氣體與液體均視為牛頓流體，連續性方程、動量方程、能量方程與輸運方程是：

式中： ?為梯度算子；t是時間；U是速度矢量；Ur是氣液界面相對速度；p是 壓力； ρm是混合物密度，ρm=ρlα+ρg(1-α)； ρl和 ρg分別是液體和氣體的密度；α是體積分數；g是重力加速度； σ是表面張力系數； μm是混合動力學黏度， μm=μlα+μg(1-α)； μl和 μg分別為液體和氣體的動力學黏度；κ是氣-液界面曲率，κ=?·n；T為溫度；Cv=Cvlα+Cvg(1-α)，Cvl和Cvg分別為液體和氣體的比熱容；λ為熱傳導系數；k為單位質量的動能。

考慮到空泡和外界空氣發生連通，空泡內氣體和外界大氣之間會發生熱交換等，仿真中氣體和液體分別采用perfectGas 和perfectFluid狀態方程[23,28]進行建模，如式（5、6）所示：

式中 ：R為狀態參數， ρl0為T=0時液體的密度。本文采用隱式歐拉格式和二階Gauss vanLeer格式來離散時間項和對流項。

圖2（a）展示的是不同網格分辨率下自由面飛濺與空泡隨時間的演化過程，其中粗、中及細網格分別對應網格總量為320萬、384萬、440萬。從對比結果可以看出，3種網格下的空泡輪廓基本吻合。對于飛濺水層模擬，中網格與細網格的結果基本吻合，粗網格結果與二者存在少許差異。圖2（b）展示的是3種網格下空泡正下方10Rmax處的壓強演化對比。從對比結果可知，中網格與細網格結果吻合較好，粗網格的結果在第二次反射的壓力值稍高于前兩者?？紤]到計算資源與計算精度的平衡，本文采用中網格來探究空泡和自由面之間的強耦合作用。

圖2 網格無關性驗證Fig.2 Grid independence verification

2 結果與討論

2.1 數值模擬方法驗證

本文數值模擬結果與文獻[11]的實驗結果進行了對比驗證。實驗中，脈沖激光束通過離軸反射鏡匯聚到水中，基于光學擊穿原理，產生一個振蕩的空化泡，通過改變自由面的位置調節無量綱距離 γ的大小。在數值模擬中，通過對比無限水域工況下實驗結果與數值模擬結果，獲得了數值模擬中空泡的初始壓力和半徑分別為30 MPa和0.3 mm，如圖3所示。

圖3 空泡在無限水域振蕩半徑的數值與實驗結果對比Fig.3 Comparison of bubble equivalent radius between experiment and numerical simulation in an unbounded liquid

以初始壓力30 MPa、半徑0.3 mm為初始條件，數值模擬了空泡在靠近自由面時發生的通氣、水層飛濺等過程。圖4分別是 γ = 0.15、0.20時實驗觀測和數值模擬的結果對比。圖中0時刻為高速相機開始采樣時刻。 γ = 0.15工況下，空泡在t= 10 μs時刻發生通氣，飛濺水層初始呈“V”型向外擴張，隨后飛濺水層逐漸向中心收縮，并在t= 124 μs附近完成閉合。γ = 0.20工況下，在空泡膨脹作用下，自由面在t= 10 μs向上飛濺并在t= 97 μs附近閉合。整體而言，針對空泡與飛濺水層的形態結果，數值模擬與實驗結果吻合較好，證明當前數值方法可用來探究空泡通氣、自由面液面飛濺及閉合強耦合作用。

圖4 實驗與數值模擬結果對比（左側為實驗結果，右側為仿真結果，深藍色為液體，白色為氣體，時間單位為μs ）Fig.4 Comparisons between experimental (left) and numerical simulation (right) results, γ = 0.15, 0.20(The deep blue and white present the phases of liquid and gas,respectively.The time unit is μs)

圖5展示的是 γ = 0.15和 γ = 0.20時空泡發生通氣、自由液面飛濺以及再閉合過程的壓力場與垂向速度場的演化過程。 γ = 0.15時，空泡在t= 10 μs發生通氣，空泡內氣體沿通道快速向外噴出，氣體噴出速度可達到1300 m/s，并在流場中產生了明顯的沖擊波。在氣體噴出約20 μs，空泡內壓力降低至20 kPa，隨著空泡繼續膨脹，外界空氣通過水層圍成的通道快速向內侵入，最高速度可達670 m/s（t= 40 μs左右），同時，通道受到內外壓強差作用快速增大。另一方面，在飛濺水層外側產生了明顯的旋渦，并逐漸向上方移動，與圖4（a）中實驗中觀察到的現象吻合。 γ = 0.20時的空泡和飛濺水層演化規律與 γ = 0.15基本相似，在t= 30 μs外界氣體開始侵入空泡內部，最大速度可達到420 m/s，通道內外壓強差達到70 kPa左右。

圖5 自由液面飛濺及閉合過程（左側為壓力場，右側為垂向速度場結果，時間單位為μs）Fig.5 Pressure fields (left) and vertical velocity fields (right) during free surface splash and closure processes when γ = 0.15, 0.20 (The time unit is μs)

2.2 閉合飛濺水層受力分析

由圖5可知，空泡內氣體首先向外部排出，由于空泡膨脹，空泡內部壓力低于外界壓力，外界氣體逐漸進入空泡內部。為了分析水層閉合機制，參考Eshraghi等[20]采用貼體坐標系對于小球入水誘導的飛濺水層進行的受力分析，本文采用固定坐標系對水層橫截面進行受力分析。如圖6（a）所示，對飛濺水層建立軸對稱坐標系，將坐標原點設在自由面上，垂直向上為^z軸正半軸，由原點向空泡向外膨脹方向為r正向。z?=表示無量綱化后的垂向高度，hc為飛濺水層閉合位置到自由面的高度，所以z?的范圍是0～1。

圖6 飛濺水層模型示意圖Fig.6 Schematic of the splash induced from the free surface

假設飛濺水層是軸對稱分布，不考慮其在 θ方向的分布，故在無量綱垂向高度z?，可將飛濺水層視為厚 度 為2a、 質量為關于軸呈軸對稱的環狀結構，如圖6（b）所示。取環狀結構上的微元進行受力分析，如圖6（c）所示，此微元主要受到空氣阻力（Fd） 、 表面張力 （Fσ1、Fσ2）、水層兩側的壓差力（F?p）以及重力（Fg）的作用。

環狀水層受到空氣阻力為：

其中：CD為阻力系數， ρg為空氣密度，a為水環厚度的一半，r(t) 為水環距離軸的距離，v為速度，阻力方向與速度v方向相反。水環受到兩個方向的表面張力[29-30]，作用于水環周向的表面張力，方向為徑向并指向空泡內，可表示為：

式中 σ=0.072 N/m為表面張力系數。水層對該水環的拉力，方向與速度方向垂直，可表示為：

水環兩側壓差帶來的壓差力，方向與速度v方向垂直，促使水層擴張或閉合，壓差力可表示為：

其中 ?p(t)為水環兩側內外壓差。為分析飛濺水層在閉合過程中受到的主要驅動力，以 γ = 0.20工況為例，根據本文數值模擬獲取的速度場與壓力場結果計算上述各力。其中，飛濺水層閉合的控制條件為：

式中vr(t)為水環徑向速度。

考慮到飛濺水層閉合是由于各力在r方向的共同作用，故對各力在徑向方向的分量進行分析，如圖7所示。Fdr、Fσ1r、Fσ2r、F?pr分別表示水環受到的空氣阻力、表面張力以及壓差力在徑向方向的分量，正負號代表各力的方向，其中沿徑向向外為正方向。由圖5可知，在t= 40 μs左右，外界氣體沿通道均勻向內流動，外界氣壓大于通道內部壓力，飛濺水層開始向中心收縮，從此開始統計各力在徑向方向分量。

圖7 γ = 0.20時不同無量綱垂向高度 z?水環所受各力在徑向方向的分量Fig.7 The radial components of forces acting on the water ring at different dimensionless vertical heights normalized by hc at γ = 0.20

由圖7可知，飛濺水層閉合的過程中，阻力、表面張力在徑向方向的分量基本保持不變，穩定在0 N（o(10-3)）附近，壓差力在水平方向的分量遠大于阻力與表面張力。當無量綱垂向高度z?在0.4～0.6附近，水環受到的壓差力隨著時間演化呈現逐漸增加的趨勢。結合圖5通道內速度場與壓力場可知，外界空氣不斷侵入通道內，且氣流逐漸穩定，使得通道內壓力逐漸增加。根據公式（10）可知，r(t)對于壓差力的影響更大，空泡的繼續膨脹導致了圖7（a、b、c）中水環內外壓差力逐漸增加。當z?在0.7左右，水環受到的壓差力出現較大的波動，結合速度場可知，z?=0.7處通道內氣體流動不均勻，導致了水環兩側壓差力的波動?？傊?，不同無量綱垂向高度z?上的水環受力分布規律基本相似，在飛濺水層閉合過程中，壓差力始終占據主導地位，表面張力與空氣阻力是可以忽略的。

2.3 飛濺水層閉合行為

考慮到飛濺水層在閉合過程中主要受到內外壓差的影響，故本文探究了不同 γ和環境壓力下的飛濺水層閉合特性。圖8所示是不同 γ空泡內平均壓強隨時間的演化關系，其中虛線表示空泡通氣過程?？张輧葰怏w的平均壓強是通過將空泡內單元網格的體積與壓強乘積積分后與空泡總體積的比值獲取的。這里定義空泡在初始與外界空氣發生連通的時刻作為空泡通氣開始的時刻tstart，飛濺水層閉合的時刻作為通氣結束的時刻tend，兩者之差作為空泡發生通氣的時間長度 ?t=tend-tstart，飛濺水層閉合的位置與初始時刻自由面之間的高度差作為飛濺水層的閉合高度hc, 如圖9（b）所示。

圖8 不同 γ情況下的空泡內部平均壓強隨時間的演化Fig.8 Evolution of averaged pressure in the bubble with time for differentγ

圖9 不同 γ時空泡通氣時間、飛濺水層閉合高度、水層閉合形態對比Fig.9 Bubble ventilation time, splash closure height, and the morphology of splash closure for differentγ

由圖8中可知，空泡在初始發生通氣時，空泡內部平均壓力大于外界環境壓力。結合圖5速度場可知，針對 γ = 0.15的工況，空泡氣體通過通道向外快速噴出，持續時間約20 μs，由于通道內外壓差，外界空氣逐漸向空泡內侵入，直至水層閉合。

圖9展示的是不同 γ情況下空泡通氣時間、飛濺水層閉合高度以及水層閉合時刻的形態。

隨著 γ的增加，空泡發生通氣的時間逐漸減少，空泡在 γ = 0.15附近通氣時間急速下降，在 γ = 0.20之后基本保持穩定，維持在80 μs左右。當 γ逐漸增大時，空泡開始發生通氣的時間逐漸延后，空泡的勢能更多地傳遞給上方水層，使得水層動能增加，獲得較大的垂向速度，使得水層閉合的高度隨 γ的增加呈現不斷增加的趨勢。

同時，在 γ = 0.15下，對比分析不同環境壓力對飛濺水層閉合行為的影響，如圖10所示。隨著環境壓力的增加，空泡形態基本相似，飛濺水層閉合高度、閉合的時刻逐漸減小。圖11定量給出了 γ =0.15時，空泡通氣時間、飛濺水層閉合高度隨不同環境壓力的演化?？张萃鈺r間和飛濺水層閉合高度隨著環境壓力的增加逐漸減少。當環境壓力小于1×105Pa時，通氣時間迅速增加；當環境壓力大于1×105Pa時，通氣時間逐漸減少并基本穩定在20 μs左右。這是因為當環境壓力增加，空泡初始時刻內外壓差減少，導致了空泡振蕩強度減弱，削弱了空泡誘導自由面飛濺水層的發展和空泡通氣過程，使得空泡較晚發生通氣，飛濺水層較早發生閉合，縮短了空泡通氣時間。

圖10 γ = 0.15時不同環境壓力下的飛濺水層閉合形態Fig.10 Morphology of splash closure for different ambient pressure conditions at γ = 0.15

圖11 γ = 0.15時空泡通氣時間和飛濺水層閉合高度隨環境壓力的演變Fig.11 The evolutions of bubble ventilation time and splash closure height under different ambient pressure at γ = 0.15

圖12展示了 γ = 0.15工況下飛濺水層閉合時間隨環境壓力的演變。圖中實心圓點來自本文數值模擬，虛線為Lee等[31]研究中獲得的冪次律關系。Lee發現在小球入水后，水層閉合時間與環境壓力之間近似滿足t～p-1關系。本文發現，在環境壓力為5×104～ 1×106Pa時，水層閉合時間與環境壓力之間近似滿足t～p-0.92關系。這是因為通氣空泡誘導與入水問題形成的飛濺水層閉合過程的受力類型和主要驅動力相同，均主要由水層兩側的壓差力主導。因此，二者冪次律關系相似，差異可能來源于實驗觀測和數值模擬中水層閉合時間的測量誤差。

圖12 γ = 0.15時飛濺水層閉合時間隨環境壓力的演變Fig.12 The evolutions of splash closure time under different ambient pressure at γ = 0.15

3 結 論

本文對空泡與自由面強耦合作用誘導的飛濺水層閉合行為進行了數值模擬研究?；陂_源的OpenFOAM平臺，采用有限體積法直接求解N-S方程，并使用VOF方法捕捉氣液界面。得到如下結論：

1）飛濺水層在閉合過程中的主要驅動力是水層兩側壓差力，重力、表面張力與空氣阻力可以忽略。

2）空泡通氣時間隨著無量綱距離 γ和環境壓力的增加而減少，飛濺水層閉合時的高度隨著無量綱距離γ的增加而增加，隨著環境壓力的增加而不斷減小，水層閉合時間與環境壓力之間近似滿足t～p-0.92關系，該冪次律與入水問題的t～p-1相似。

本文數值模擬了空泡通氣、水層飛濺與閉合等現象，旨在為二者的強耦合作用研究提供參考。在未來的研究中，需要引入粒子圖像測速等技術，定量測量通氣過程的速度場演化，為數值模擬研究提供重要數據對比。同時，在理論建模方面，考慮水層演化對空泡行為的影響，建立二者雙向強耦合物理模型。