基于數值仿真的球體垂直入水空泡演化過程研究

張晨星，谷立祥，權曉波，魏海鵬，程少華

（北京宇航系統工程研究所，北京，100076）

0 引 言

物體由空氣介質穿越氣水界面運動至水中的過程稱為入水過程。最早對物體入水開展研究的是英國學者Worthington A M和Cole R S［1］，他們通過高速閃光相機對球體落入液體中所產生的噴濺進行了拍攝，得到了大量的球體入水圖像。此后，隨著魚雷空投入水、水上飛機水上著陸、火箭助推器以及航天器海面回收等工程問題的出現，大量的學者對入水問題開展了相應的研究。

作為最早得到研究的入水問題，球體入水一直受到學者們的關注。球體入水的過程往往伴隨著入水空泡的產生和演化。入水空泡的演化過程涉及到固液沖擊和多相流動，是球體入水問題的一項重要研究內容，國內外學者均對球體入水空泡開展了相應的研究工作。

早期對球體入水空泡的研究主要以試驗和理論為主，Gilbarg D 和Anderson R A［2］通過試驗研究了液面壓力和入水速度對球體入水空泡閉合時間的影響，發現液面壓力越高，入水速度越快，球體入水空泡表面閉合的速度就越快；May A［3］研究了球體表面狀態對入水空泡的影響，發現球體表面污染越嚴重、表面沾有的液體黏性越強，球體入水產生空泡所需的速度越小；此外，May A［4］在理論和試驗的基礎上，提出了適用于不同外形物體入水空泡早期形態的計算公式；Duez 等［5］提出了用于判斷物體入水是否產生空泡的依據，指出物體入水空泡的產生受到入水速度和物體表面親/疏水性兩個因素的影響；在Duez 研究的基礎上，Truscott T T［6］開展了旋轉球體的垂直入水試驗，發現除了入水速度和球體表面親/疏水性外，旋轉也會對空泡的形成產生影響；馬慶鵬等［7］對球體垂直入水開展了試驗研究，得到了球體入水空泡的演化圖像，分析了入水速度和表面沾濕狀態對入水空泡流動的影響。

在20 世紀80 年代以后，隨著計算機技術的發展和流體力學數值計算方法的成熟，數值仿真成為研究球體入水問題新的有效方法。Lee M等［8］為驗證其提出的高速入水空泡動力學模型，對球體的高速入水過程進行了數值模擬，得到了球體高速入水時的空泡形態；張偉偉等［9］采用多物質ALE 算法對球體的入水空泡進行了計算，發現該算法能夠較好地對球體入水過程進行仿真模擬；余彧等［10］研究了不同的計算參數對親/疏水性球體入水過程仿真結果的影響；孫釗等［11］對不同親/疏水性球體的入水空泡進行了仿真計算，得到了4種不同的空泡形態。

綜上所述，在對球體入水空泡的研究方面，國內外學者在試驗、理論和數值仿真方面均取得了一定的成果，但目前對球體入水空泡的數值仿真研究多集中于對算法的研究以及對不同因素對入水空泡的影響規律的研究，對入水空泡演化過程的詳細分析研究較少。而試驗研究會受到試驗條件和測量設備的限制，獲取的球體入水流場信息以圖像和壓力數據為主，難以對入水空泡的演化進行細致的分析。因此，本文采用數值仿真的方法對球體垂直入水過程進行了仿真計算，基于仿真結果對球體入水空泡的演化過程進行全面系統地分析。

1 數學模型

1.1 基本控制方程

球體入水涉及到多相流動，因此本文采用均質多相流的N-S 方程作為流動的基本控制方程。多相流模型選取VOF 模型，即將多相流體視為單一的流體介質混合物，各相具有相同的壓力場和速度場。αl和αg分別表示液相和氣相的體積分數，兩者滿足關系式：

式中ρm為混合介質的密度，且ρm=αlρl+αgρg，ρl和ρg分別為液相和氣相的密度；ui為速度分量。

動量守恒方程為

式中μm為混合介質的黏性系數，且μm=αlμl+αgμg，μl和μg分別為液相和氣相的黏性系數；μt為湍流黏性系數。

在連續性方程和動量守恒方程的基礎上，VOF模型增加了第n相體積分數的輸運方程：

通過求解該方程得到控制單元內第n相的體積分數αn來表征該流體單元內的流體相分布。

1.2 湍流模型

本文采用標準k-ε湍流模型，其湍流動能方程為

式中k為湍流動能；ε為湍流能量耗散率；Gk和Gb為湍流生成項；YM為湍流耗散項；Sk和Sε為源項，C1ε=1.44，C2ε=1.92，σk= 1.0，σε=1.3。

1.3 網格劃分及邊界條件設置

Aristoff J M 等［12］對直徑為25.4 mm，不同密度比的小球以2.17 m/s 左右速度的入水過程開展試驗研究，本文采用二維軸對稱模型對其試驗中密度比為2.3的球體的入水過程展開仿真計算。如圖1所示，空氣域的高度L2=160 mm，水域高度L3=400 mm，徑向直徑L1=400 mm。球體壁面第一層網格高度取為0.002 mm，網格均采用結構化網格且在球體附近加密，網格數量為158 400。入口條件為壓力入口，入口壓力設置為標準大氣壓強p0=101 325 Pa，出口條件為壓力出口的壓力為p=p0+ρlgh，其中g為重力加速度，h為水深。本文中壓力出口處水深h=L3=400 mm，因此壓力出口處的壓力設置為105 242 Pa。球體表面和計算域側面定義為壁面條件。計算初始時刻，球體距水面1 mm，設置球體的初始速度為2.17 m/s。

圖1 計算域、邊界條件及網格劃分示意Fig.1 Schematic of computational area, boundary conditions and mesh

為方便對球體的入水過程進行描述，如圖1 所示，記小球下落的方向為軸向x，記垂直于下落的方向為徑向y，對稱軸和自由水平液面的交點為坐標原點。定義角度θ為球心和球面一點連線與軸向的夾角，逆時針為正。此外，入水深度為球心距水面的距離，取球心到達水面的時刻為零時刻。

1.4 數值計算方法

本文應用流體計算軟件對球體入水過程進行計算。采用有限體積法對控制方程進行離散，采用SⅠMPLE算法對壓力-速度場的耦合進行求解。壓力場的離散選用PRESTO格式，體積分數的離散采用Geo-Reconstruct格式，動量方程和能量方程的離散采用二階迎風格式。由于入水過程中球體在計算域中發生運動，因此采用動網格技術來實現網格的移動和更新。同時，編寫用戶自定義函數來計算球體的運動速度和位移。

1.5 仿真結果驗證

如圖2所示，本文通過對仿真得到的入水空泡形態和文獻［2］的試驗結果進行對比來驗證數值仿真的正確性。可以看到，仿真較好地模擬出了球體入水空泡的擴張、收縮以及閉合等演化過程。仿真得到的空泡形態在球體入水初期和試驗結果的一致性較好。仿真中空泡的閉合時間相比試驗提前了約6.7 ms，存在10.14%的時間誤差。出現誤差的原因可能包括網格離散引起計算誤差、試驗中對高速攝像進行時間標定時存在測量誤差以及仿真模型中球體的物理特性與試驗中不完全相同等。

圖2 球體入水空泡形態試驗和仿真結果對比Fig.2 Comparision of experimental and numerical results of cavity's shape which was formed during sphere's waterentry

球體的表面特性是球體的一項重要物理特性，球體表面出現沾濕、沾灰等污染現象都會影響球體的表面特性，即會改變球體表面的潤濕性［11］。表面接觸角是表征球體表面潤濕性的一個常用物理量，在Aristoff J M 試驗中， 小球的表面接觸角為（120±5）°，考慮球體表面受到污染，取球體表面接觸角為150°，對球體的入水過程進行仿真計算，得到球體入水空泡的閉合時間為62.7 ms，空泡閉合時間的誤差降至5.14%。可見球體表面的污染特性會對球體入水空泡的閉合時間產生一定的影響。

針對Aristoff J M 的球體入水試驗，張偉偉等［9］采用多物質ALE仿真算法進行了仿真模擬，重點分析了球體密度和入水速度對入水空泡閉合時間的影響。得到的空泡閉合時間與試驗值也存在約10%的偏差，可見實現與試驗結果高度吻合的仿真較為困難。

2 球體入水空泡演化過程分析

2.1 入水空泡演化過程階段劃分

圖3 給出了球體入水過程中空氣體積分數的云圖。可以看到，在球體撞擊水面之后，球體底部的水體被排開，球體底部被沾濕。部分水體沿著球體表面的切向斜向上涌起，形成高于液面的射流。隨著球體的下落，球體的沾濕面積迅速增大，直至水體從球體表面分離。水體從球體表面分離之后，球體后方的空氣源源不斷地進入水面以下，形成一個開口的空泡，通常稱之為入水空泡。空泡形成之后，首先隨著球體的下落被不斷地拉伸，同時向徑向擴張，空泡呈現為倒圓臺的形態。隨著球體入水深度的不斷增加，水面和球體中間的部分空泡截面開始由擴張變為收縮，收縮最快的空泡截面率先收縮為一點，將入水空泡分割為兩個部分：上方的敞口空泡和下方附著在球體上的尾空泡。這種現象通常稱為空泡的深度閉合，空泡閉合時，上方的敞口空泡呈現為倒錐形，下方的尾空泡呈現為錐形，且尾空泡的體積明顯小于敞口空泡。空泡的深度閉合阻斷了空氣進入下方的尾空泡，同時產生上下兩股高速射流。向上射流的上涌使得敞口空泡逐漸消失，向下的射流則擊打在球體表面，使得水體從球體的上表面沾濕球體。隨著球體進一步下落，向上的射流涌出水面，而尾空泡則不斷被壓縮，泡內氣體不斷泄漏，空泡潰滅。

圖3 球體入水空泡演化過程示意Fig.3 Schematic of the evolution of cavity formed during sphere's water-entry

根據球體沾濕面積的變化規律對球體入水初期的空泡演化過程進行階段劃分。球體的沾濕角度隨時間的變化曲線如圖4所示，這里的沾濕角度θw是指固液接觸線處的角度值。記球體入水至沾濕角度快速增大的階段為沖擊及流動形成階段，可以看到在此階段球體沾濕角度在6.6 ms的時間內由0°增加至約78.5°。流動形成之后，球體沾濕面積增大的速度明顯放緩，空泡壁面各處沿徑向均在不斷擴張，這個階段即為空泡的發展階段。而在空泡深度閉合階段，空泡壁面則部分沿徑向擴張，部分向內收縮，直至壁面某處在空間收縮為一點。這兩個階段是空泡發展演化最重要的階段，即空泡發展及深度閉合階段。空泡閉合后，敞口空泡的長度不斷減小并最終消失。尾空泡內則不斷被壓縮，氣體不斷泄漏，直至球體表面完全浸濕，這個階段即為空泡的潰滅階段。

圖4 沾濕角度θw隨時間變化曲線Fig.4 Curve of wet angle θw versus time

2.2 沖擊及流動形成階段空泡演化分析

入水沖擊是指球體與水接觸，球體壓縮水體產生一個高于水中聲速的沖擊波的過程，這個過程時間很短，通常為微秒級［14］，球體會受到極大的沖擊載荷。沖擊及流動形成階段流場的壓力云圖和氣水界面圖如圖5所示。可以發現，球體經歷入水沖擊后，球體沾濕部分下方的流場形成高壓區。隨著時間的推移，球體的沾濕面積不斷增大，流場的高壓區在不斷地擴展，同時高壓區的壓力在逐漸降低。

圖5 沖擊及流動形成階段流場壓力云圖及氣水界面圖Fig.5 Contours of pressure and air-water interface of the flow field during impact and flow formation stage

進一步對作用在球體表面上的壓力進行分析，不同時刻球體表面的壓力隨θw的變化曲線如圖6 所示，圖中虛線依次表示不同時刻球體的沾濕角度。可以看到，在入水初期的某一固定時刻，球體表面的壓力分布規律如下：球體正下方出現高壓，且高壓向球體外側不斷拓展、增大，在固液接觸線附近達到了最大值。壓力出現上述分布的原因在于，入水初期，沾濕部分的水體在沖擊作用下開始流動，而在固液接觸線處的水體流動則尚未形成，因此與球體的沖擊作用最為劇烈，從而壓力在此處最高。隨著時間推移，球體的沾濕面積不斷增大，球體附近水體的流動逐漸穩定，球體受到的沖擊力逐漸轉變為穩定的流體阻力，因此球體下方的壓力不斷減小，同時最大壓力的位置逐漸回歸到球體的正下方。

圖6 不同時刻球體表面無量綱壓力隨角度變化曲線Fig.6 Curve of dimensionless pressure versus θw at different times

圖7為沖擊及流動形成階段流場的速度云圖和氣水界面圖。可以發現，在球體臨近入水時刻，球體下方的空氣被高速排出。隨著球體的入水，固液接觸線處的水體在沖擊高壓產生的壓差力的作用下沿著球體表面迅速向外延伸，形成高速射流。同時，球體的動能不斷轉化為水的動能，當球體表面的水體獲得足夠大的徑向速度時，水體脫離球體表面開始形成入水空泡，流動也逐漸穩定。

2.3 空泡發展及深度閉合階段空泡演化分析

空泡發展階段t=17.15 ms 時流場的壓力云圖、速度云圖、流線圖和氣水界面圖如圖8所示。由壓力云圖可知，由于開口空泡和水面大氣相連通，因此泡內的壓力接近大氣壓力。而空泡外部的壓力則沿徑向逐漸增加至水體靜壓，因此在空泡發展階段，空泡壁面會受到周圍水的壓力，且水深越深，壓力就越大。同時可以看到在固液接觸線的下方出現了低壓區。而由速度云圖及流線圖可以發現，隨著球體的下落，球體后方的空氣會高速進入空泡內部，泡內空氣的流動速度要明顯大于空泡周圍水體的流動速度。

圖8 空泡發展階段流場示意Fig.8 Schematic of the flow field during the cavity development stage

由于空泡某一截面的徑向位置是徑向速度對時間的積分，因此空泡壁面的徑向速度是影響入水空泡形態的重要參數。空泡壁面上流體的徑向速度隨軸向位置和時間的變化曲線如圖9所示。

圖9 空泡發展階段空泡壁面徑向速度變化Fig.9 Curve of radial velocity of the cavity's wall versus time and radial position during the cavity development stage

由圖9a可知，對于某一固定空泡截面，空泡壁面的徑向速度在球體經過時迅速增大，而后空泡壁面的徑向速度在壓差力的作用下經歷了一個緩慢的衰減過程。軸向位置越大，即水深越深的截面處，由于球體經過該位置時的速度更低，且壓差力也更大，因此該處空泡壁面能達到的最大徑向速度越小，且速度衰減得越快。

由圖9b可知，對于某一固定時刻，在空泡發展階段的早期，空泡壁面的徑向速度沿著軸向降低，球體表面處空泡的徑向速度最小。但在空泡發展階段的后期，空泡壁面最小徑向速度的位置出現在了水面和球體之間的一點。這是由于該點上方的空泡壁面受到的壓差力相對較小，因此徑向速度衰減較慢。而該點下方的空泡壁面雖然受到較大的壓差力，但由于其仍受到球體運動的影響，因此徑向速度衰減也相對較慢。

在空泡發展階段，空泡壁面的徑向速度均為正值，空泡壁面各點均保持擴張。而一旦空泡壁面的徑向速度出現負值，則表明空泡出現了部分收縮，空泡演化進入深度閉合階段。

空泡深度閉合階段流場的壓力云圖及氣水界面圖如圖10 所示。可以看出，在深度閉合階段的初期，流場的壓力場特征和空泡發展階段相類似，空泡內的壓力和大氣壓接近。隨著空泡的進一步收縮，收縮位置下方的空泡受到壓縮，泡內的壓力開始上升。空泡閉合的瞬間，徑向運動的水體在閉合點處撞擊，形成高壓區。

圖10 空泡深度閉合階段流場的壓力云圖及氣水界面圖Fig.10 Contour of pressure and air-water interface of the flow field during the cavity deep closure stage

圖11為空泡在深度閉合階段的速度云圖、流線圖及氣水界面圖。由圖11b和11c可知，空泡開始收縮后，收縮位置的空氣在水體的擠壓下，向上、向下排出。由于空氣向下排出時空間有限，會造成下方空泡的壓力上升，流動受阻，因此下方空泡內的空氣流速較低。而空氣向上排出時，由于無空間限制，因此流動較為順暢，流速也更高。同時，在空泡收縮的位置，存在速度滯止的區域。空泡閉合時如圖11d所示，沖擊產生的高壓驅使水體向上、向下運動，形成高速射流。

圖11 空泡深度閉合階段流場速度云圖、流線圖及氣水界面圖Fig.11 Contour of velocity, streamlines and air-water interface of the flow field during the cavity deep closure stage

對空泡形態的演化進行進一步分析，提取深度閉合階段空泡壁面的徑向速度隨軸向位置的變化曲線如圖12 所示。可以看出，在深度閉合階段初期，空泡中間偏向球體位置處的壁面徑向速度在周圍水體壓力的作用下開始減小為負值，空泡在此處收縮。但收縮位置上下方的空泡壁面仍保持擴張，這種現象通常稱為空泡的“頸縮”。隨著空泡閉合的進行，空泡壁面越來越多的部分開始處于收縮狀態，甚至剛剛脫離球體表面的空泡壁面也在收縮。這表明相比于球體向水體轉移的動能，壓力已在空泡的發展中起主要的影響作用。此外還可以發現，水面處的空泡壁面一直在保持擴張，這是由于水面的壓力過小。

圖12 空泡深度閉合階段空泡壁面徑向速度隨軸向位置變化曲線Fig.12 Curve of radial velocity of the cavity's wall versus radial position during the cavity deep closure stage

2.4 空泡潰滅階段空泡演化分析

空泡潰滅階段流場的壓力云圖和氣水界面圖如圖13所示。由圖13a可知，空泡閉合之后，閉合點上下方的空泡壁面繼續收縮，閉合點擴展為閉合線。徑向運動的水體在閉合位置繼續相撞形成高壓，高壓區沿著閉合線向上、向下擴展，壓力等值線呈現“8”字形。提取t=59.75 ms 時刻空泡壁面的徑向速度，發現閉合線上方的徑向速度的量值要比下方的大，因此敞口空泡底部的沖擊更為劇烈，因而敞口空泡底部的沖擊高壓區要比尾空泡上方的高壓區大。

圖13 空泡潰滅階段流場壓力云圖及氣水界面圖Fig.13 Contour of pressure and air-water interface of the flow field during the cavity collapse stage

同時，射流在壓差力驅使下繼續運動。向上的射流使得上方的敞口空泡的體積快速減小直至空泡消失，球體上方流場的壓力逐漸穩定。向下的射流則擊打在球體的上表面對球體產生沖擊。尾空泡內的壓力隨著球體的下落而增加，因此空泡被壓縮，導致泡內含氣量逐漸減小直至空泡完全潰滅。

3 結 論

本文對球體低速垂直入水進行了仿真計算，通過和試驗結果的對比驗證了仿真的正確性，重點對球體入水空泡的演化過程分階段進行了分析，研究得到以下結論：

a）球體在入水初期，下方流場會形成高壓區。壓差力引起水體沿球體表面延伸形成高速射流。同時，球體表面上固液接觸線處的壓力最高。流動穩定之后，壓力最高點回歸到球體的正下方。

b）球體入水空泡形成之后，空泡壁面各處的徑向速度經歷了先快速增大而后緩慢減小的階段。靠近水面位置的空泡壁面，由于受到的水的壓力較小，因此其徑向速度衰減得最慢。

c）空泡發生深度閉合，原因在于水和空泡之間壓差力使得空泡壁面的徑向速度衰減為負，徑向速度衰減快的空泡壁面率先收縮為一點，空泡閉合。

d）空泡閉合時會在閉合位置形成高壓，進而形成兩股高速射流。由于水和閉合點上方敞口空泡的壓差力更大，因此向上射流的速度比向下射流的高。