張晨星，谷立祥，權(quán)曉波，魏海鵬，程少華

（北京宇航系統(tǒng)工程研究所，北京，100076）

0 引 言

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

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

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

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

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

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

1.1 基本控制方程

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

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

動量守恒方程為

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

在連續(xù)性方程和動量守恒方程的基礎(chǔ)上，VOF模型增加了第n相體積分?jǐn)?shù)的輸運方程：

通過求解該方程得到控制單元內(nèi)第n相的體積分?jǐn)?shù)αn來表征該流體單元內(nèi)的流體相分布。

1.2 湍流模型

本文采用標(biāo)準(zhǔn)k-ε湍流模型，其湍流動能方程為

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

1.3 網(wǎng)格劃分及邊界條件設(shè)置

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

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

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

1.4 數(shù)值計算方法

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

1.5 仿真結(jié)果驗證

如圖2所示，本文通過對仿真得到的入水空泡形態(tài)和文獻(xiàn)［2］的試驗結(jié)果進(jìn)行對比來驗證數(shù)值仿真的正確性?？梢钥吹剑抡孑^好地模擬出了球體入水空泡的擴(kuò)張、收縮以及閉合等演化過程。仿真得到的空泡形態(tài)在球體入水初期和試驗結(jié)果的一致性較好。仿真中空泡的閉合時間相比試驗提前了約6.7 ms，存在10.14%的時間誤差。出現(xiàn)誤差的原因可能包括網(wǎng)格離散引起計算誤差、試驗中對高速攝像進(jìn)行時間標(biāo)定時存在測量誤差以及仿真模型中球體的物理特性與試驗中不完全相同等。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖7為沖擊及流動形成階段流場的速度云圖和氣水界面圖?？梢园l(fā)現(xiàn)，在球體臨近入水時刻，球體下方的空氣被高速排出。隨著球體的入水，固液接觸線處的水體在沖擊高壓產(chǎn)生的壓差力的作用下沿著球體表面迅速向外延伸，形成高速射流。同時，球體的動能不斷轉(zhuǎn)化為水的動能，當(dāng)球體表面的水體獲得足夠大的徑向速度時，水體脫離球體表面開始形成入水空泡，流動也逐漸穩(wěn)定。

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

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

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

由于空泡某一截面的徑向位置是徑向速度對時間的積分，因此空泡壁面的徑向速度是影響入水空泡形態(tài)的重要參數(shù)?？张荼诿嫔狭黧w的徑向速度隨軸向位置和時間的變化曲線如圖9所示。

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

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

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

在空泡發(fā)展階段，空泡壁面的徑向速度均為正值，空泡壁面各點均保持?jǐn)U張。而一旦空泡壁面的徑向速度出現(xiàn)負(fù)值，則表明空泡出現(xiàn)了部分收縮，空泡演化進(jìn)入深度閉合階段。

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

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

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

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

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

圖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可知，空泡閉合之后，閉合點上下方的空泡壁面繼續(xù)收縮，閉合點擴(kuò)展為閉合線。徑向運動的水體在閉合位置繼續(xù)相撞形成高壓，高壓區(qū)沿著閉合線向上、向下擴(kuò)展，壓力等值線呈現(xiàn)“8”字形。提取t=59.75 ms 時刻空泡壁面的徑向速度，發(fā)現(xiàn)閉合線上方的徑向速度的量值要比下方的大，因此敞口空泡底部的沖擊更為劇烈，因而敞口空泡底部的沖擊高壓區(qū)要比尾空泡上方的高壓區(qū)大。

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

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

3 結(jié) 論

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

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

b）球體入水空泡形成之后，空泡壁面各處的徑向速度經(jīng)歷了先快速增大而后緩慢減小的階段?？拷嫖恢玫目张荼诿妫捎谑艿降乃膲毫^小，因此其徑向速度衰減得最慢。

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

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