空氣壓強對垂直入水空泡影響的數值研究

王 聰，何春濤，權曉波，魏英杰

(1.哈爾濱工業大學航天學院，150001哈爾濱;2.中國運載火箭技術研究院總體設計部，100076北京)

物體高速入水時，在水面打開的1個開口空腔，稱為入水空泡.入水空泡的生成、發展以及潰滅過程，涉及流固耦合、介質突變、自然空化以及強非線性等，是1個非常復雜的過程，尤其會對空投魚雷、導彈、超空泡射彈等入水瞬間的流體動力特性和水下彈道的穩定性產生較大的影響.對入水問題的研究最早可以追溯到十九世紀末A.M.Worthington和R.S.Cole［1］所開展的關于液滴和球體等物體入水時引起的自由液面噴濺和入水空泡的研究等.1929年Von K［2］在不考慮重力、水的粘性力以及空氣慣性力的基礎之上完成了飛機水面降落沖擊載荷的研究.上世紀50年代Albert May等［3-7］在美國海軍的支持下開展了大量關于入水問題的基礎研究，主要研究內容包括垂直和傾斜入水空泡、流體動力、入水撞擊和入水彈道等.1968年J.G.Waugh［8］等人在美國水下海軍作戰中心的支持下基于弗勞德數和空氣密度等參數的相似性，完成了不同頭型的Mk25空投魚雷縮比模型傾斜入水的實驗研究.Albert May［7］和David Gilbarg等［9］提出將空氣壓強與標準大氣壓強的比值，或密度與標準大氣壓強下空氣密度的比值作為入水空泡閉合前的空泡動力相似參數，并對不同壓強和密度條件下的球體低速入水空泡形態進行了試驗研究.研究表明空氣壓強對入水空泡的運動特性和空泡閉合方式等具有不可忽略的影響，因此空氣壓強常被作為最主要的入水空泡流動特性的相似參數之一.根據以上情況，本文將采用數值研究的方法對上述特性開展進一步的研究.隨著測試設備和手段的發展，越來越多的研究人員對入水問題的機理開展進一步的研究.Jeffrey M.Aristoff等［10］利用PIV等對低密度的球體垂直入水過程中球體和空泡壁面速度變化規律進行了試驗研究;Cyril Duez等［11］不同親水性的球體垂直入水過程所產生的噴濺進行了對照試驗研究;Hongmei Yan等［12］對低Re數(o(10))條件下球體垂直入水空泡動力學進行了試驗研究.

綜上所述，過去對入水問題的研究大多都是通過試驗的方法，近年來人們(在大量寶貴的實驗數據基礎之上)開始發展入水空泡的數值計算方法［13-14］，由于數值計算方法大多是半經驗型的，因此，對于一些復雜影響因素難以考慮周全.本文在RANS方程基礎之上引入VOF多相流模型，并通過引入水和水蒸汽之間的輸運關系，在考慮自由液面的同時，還考慮了由入水所引起的自然空化現象，對帶150°錐角圓柱體勻速垂直入水空泡形態與空氣域壓強之間的變化規律進行了數值研究.在計算過程中采用動網格技術實現圓柱體入水過程的絕對運動與流場運動之間的耦合求解過程.本文假設空氣域為理想氣體，空氣密度隨空氣域壓強變化而發生相應的變化.分別對1/8～4倍標準大氣壓強(P0=101 325 Pa)條件下，150°錐角頭型圓柱體以10～50 m/s速度勻速垂直入水的空泡進行了數值研究.

1 控制方程

1.1 連續性方程和動量方程

混合相連續性方程為

空氣相連續性方程為

水蒸汽連續性方程為

混合相動量方程為

式中:ui表示混合物速度分量(m/s);ρm為混合物密度(kg/m3)，表達式如下:

μm是混合物動力粘性系數(m/s)，

下標m、l、g和v分別表示混合相、水、空氣和水蒸汽;ag、av和1-ag-av分別為空氣、水蒸汽和水的體積分數;Re、Rc表示水蒸汽產生和凝結過程的源項.

1.2 空化模型

高速入水物體在入水過程中產生空泡，僅有物體頭部與水接觸時，在物體頭部邊緣分離處會出現較大的壓強梯度變化，當壓強降低到飽和蒸汽壓強之下時，可能會出現自然空化現象.自然空化是1個非常復雜的相變過程.本文通過式(3)中的源項來表示該過程(即Zwart-Gerber-Belamri空化模型)，表達式如下:

式中:RB=10-6m為氣泡半徑;anuc=5×10-4為不可凝結氣體體積分數;Fvap、Fcond為經驗常數，分別取50和0.001.

1.3 湍流模型

本文計算中采用標準k-ε湍流模型，湍動能k和湍能耗散率ε的輸運方程分別為

式中:μt=ρmCμk2/ε為湍流粘性系數;σk=1.0和σε=1.3分別為k和ε的Prandtl數;C3ε=tan｜v/u｜，v為平行于重力場矢量的流速分量，u為垂直于重力場矢量的速度分量;Cμ、C1ε和C2ε為經驗常數，分別取0.09、1.44，和1.92;Gk=為由平均速度梯度所引起的湍動能生成項.

2 數值求解

2.1 計算域和邊界條件

本文計算模型采用150°錐角頭型圓柱體，直徑D=10 mm，長度L=40 mm.計算域設置和模型表面網格分布如圖1所示，空氣域長度為60D，水域為400D，徑向為80D.文獻［1］中的研究表明，軸對稱物體垂直入水所產生的空泡具有軸對稱性質.因此，本文數值計算均采用二維軸對稱計算模型.使用四邊形網格，網格數為324 750.同時還采用了107 000、168 800和515 000三種不同數量的網格進行對照計算，計算結果表明網格數較大的3種計算模型的空泡形態和空泡表面閉合時間等參數之間差值小于5%.綜合考慮計算結果精度和計算量等因素，最后選取324 750數量的網格作為本文計算網格.

圖1 計算域和圓柱表面網格分布

2.2 數值方法

本文采用有限體積法對控制方程進行空間離散，并應用PISO方法對動量方程和連續性方程聯立求解建立速度場和壓力場的耦合關系.計算過程中對流項的離散采用二階迎風格式;水蒸汽相的離散采用QUICK格式;湍流輸運方程的離散采用二階迎風格式;對時間的離散采用一階顯示格式.計算時間步長的選取根據網格尺度、運動速度等因素綜合確定.并通過動網格技術實現圓柱體的運動過程與計算域內的網格更新.由于本文計算采用結構化網格，網格動態更新過程中主要是被拉伸或被壓縮變形，因此，采用動態層法進行網格更新.動態層法的核心是當網格被拉伸到大于設定的理想最大尺度((1+as)hideal)時，該網格分離成兩層新的網格，當網格被壓縮到小于設定的理想最小尺度(achideal)時網格與其相鄰的網格合并成新的網格，即

引入動態網格后，相當于網格按一定速度作對流運動，在實際計算中應減掉網格的對流效應.因而引入動網格后控制體V對變量φ的統一形式的守恒方程可寫為

式中:ρm為混合物密度;為速度矢量;g為動網格的運動速度;Γ為擴散系數;Sφ為標量φ的源項;?V表示控制體積V的邊界.

3 數值計算方法對比與驗證

本文通過對Albert May等［7］提出的勻速垂直入水理想空泡模型進行對比計算，驗證本文所采用的數值計算方法的可信性.為了便于分析，引入如圖2所示的球體垂直入水空泡柱坐標剖面示意圖.

圖2 垂直入水空泡坐標

坐標原點建立在未擾動水平自由水面與球體對稱軸交點處，徑向半徑與未擾動水平自由水面重合，豎直向下為軸向正方向.在該二維軸對稱坐標系下，垂直入水早期空泡形態的理想空泡模型為［7］

式中CD為入水物體的阻力系數，d為入水物體的特征尺寸.該入水空泡模型主要是用于描述入水過程中，空泡閉合之前發展較為充分時期的外部輪廓.

根據以上理想空泡模型邊界條件要求，對球體以30 m/s從空氣域勻速垂直入水過程進行數值計算，空氣域環境壓強為101 325 Pa，空氣、水和水蒸汽的密度分別為1.225、998.2和0.554 kg/m3.將以上計算結果選取入水空泡充分發展的某一時刻的空泡形態，與相同時刻下由式(12)理想空泡模型的擬合結果進行對比.阻力系數取CD=0.20［15］，并用球體直徑作為特征尺寸分別對空泡直徑Da和長度La進行無量綱化.數值計算與理論公式計算結果對比如圖3所示.

圖3 球體入水空泡輪廓計算結果與經驗公式對比

由圖3可以看出，數值計算結果與理想空泡模型擬合結果具有較好的一致性.在空泡輪廓兩端具有一定的誤差，產生該誤差的主要原因是由于理想空泡模型本身忽略了自由表面處和分離點處的空泡輪廓的不規則性.通過以上對球體垂直入水空泡形態的計算方法的驗證，說明了本文所采用的數值計算方法是合理可信的.并在此基礎上展開空氣域壓強(或密度)變化對帶150°錐角頭型圓柱體勻速垂直入水空泡形態的影響規律的數值計算研究.

4 空氣壓強對空泡影響分析

以下對帶150°錐角頭型圓柱體分別以10、20、30、40和50 m/s的速度勻速垂直入水的空泡生成過程、自然空化以及不同空氣壓強對空泡閉合影響規律等進行數值計算，并對結果進行分析.

4.1 自然空化對空泡形態的影響分析

圖4給出了圓柱體以30 m/s勻速垂直入水過程中不同時刻空泡的形態輪廓圖.零時刻圓柱頭部與自由水面接觸.

由圖4可以看出，隨著入水時間的推移，圓柱體逐漸打開1個與空氣相同的空腔，即入水空泡.空泡的生成過程實質上是1個圓柱體與水之間的能量轉換過程.圓柱體撞擊自由水面將動能傳遞給水體，并通過阻力作用持續地將動能傳遞給水體，動能將驅使水體運動，同時水體又受到由重力所產生的靜壓的阻礙作用.空泡壁面附近的水介質在動能和勢能的共同作用下運動.由于圓柱體撞擊自由液面時，在自由液面附近的水介質會產生1個繞圓柱體肩部的軸向速度分量，使得部分水介質以該速度相上運動，形成噴濺.并由于空泡壁面的運動，使得自由液面附近的水面上升，并超過了未擾動自由液面高度.在空泡入口處，由于噴濺的作用，使得入口逐漸變小以至于空泡表面閉合，阻止了外界空氣進入到空泡內部，而由于慣性的作用空泡體積繼續增加，使得空泡內部壓強迅速降低，分布變得更加復雜.圖5為標準大氣壓下，圓柱體以50 m/s勻速垂直入水空泡閉合瞬間肩部壓力梯度云圖.

圖4 圓柱體以30 m/s勻速垂直入水時的空泡云圖

圖5 圓柱體以50 m/s入水圓柱肩部壓力云圖

由圖5可以看出，由于頭部繞流和肩部銳邊使得流體在流過銳邊時產生較大的順壓梯度.根據伯努利方程可知，當dP/ds=0時，壓強降低到最低值.根據柯乃普［16］對空化的定義，如果壓強值低于水的飽和蒸汽壓強(本文取3 510 Pa)，并在時間和空間上持續一定長度，則會產生自然空化.由于此時空泡處于剛剛閉合或快要閉合狀態，空泡內部壓強接近空氣壓強，遠遠高于水的飽和蒸汽壓強，其將阻礙空化的發展.上述壓強梯度變化所對應的肩部自然空化的水蒸汽體積分數如圖6所示.

圖6 以50 m/s速度入水的圓柱肩部水蒸汽體積分數云圖

由圖6可以看出，只是在肩部后小范圍內出現局部的自然空化，且水蒸汽的含量較低.通過對10、20、30和40 m/s速度垂直入水的計算結果的自然空化狀態對比分析表明，隨著入水速度的增加，空化區域呈增大趨勢.若以10 m/s和20 m/s速度垂直入水，在圓柱體肩部幾乎不發生空化.在10～50 m/s速度范圍內，帶150°錐角頭型圓柱體勻速垂直入水，自然空化的范圍較小，只是出現在肩部，對入水空泡形態的影響較小.

4.2 空氣壓強對空泡表面閉合影響分析

圖7 給出了圓柱體在不同速度條件下勻速垂直入水時空泡無量綱閉合時間Ts/(D/V)與空氣壓強之間的關系曲線.

圖7 入水空泡閉合時間與空氣壓強之間的關系曲線

其中，Ts為空泡表面閉合時間，D為圓柱體直徑，V為圓柱體入水速度;P0為標準大氣壓強(本文取101325 Pa)，P為實際空氣壓強.計算壓強分別為4P0、2P0、P0、(1/2)P0、(1/4)P0和(1/8)P0.由圖7可以看出，不同速度條件下入水空泡無量綱閉合時間隨壓強的變化趨勢基本一致，當空氣壓強高于1個標準大氣壓時，隨著壓強的升高，空泡在自由液面上閉合即表面閉合時間縮短，基本呈線性趨勢變化;當壓強低于1個標準大氣壓強時，隨著壓強的降低，表面閉合時間迅速增加，呈非線性趨勢變化.同時，可以看出在相同壓強條件下不同速度入水空泡表面閉合時間趨于常數，當空泡閉合方式由表面閉合轉變為在自由液面之下某一點閉合即深閉合后，空泡閉合時間變化規律變得較為復雜.圖中P/P0=1/4時，10 m/s速度入水空泡閉合方式由表面閉合轉變成為深閉合.隨著空氣壓強的降低，低速度入水空泡最先出現空泡閉合方式的轉變，即由表面閉合轉變成深閉合.圖8給出了圓柱體以10 m/s速度分別在2P0、P0、(1/2)P0和(1/4)P0時，空泡第一次閉合瞬間空泡輪廓圖.

圖8 不同壓強條件下入水空泡閉合瞬間輪廓圖

由圖8可以看出，在壓強高于(1/2)P0時空泡在表面發生閉合，當壓強為(1/4)P0時空泡在表面以下某點發生閉合.并且隨著空氣壓強的降低空泡閉合時間延后，當降低到一定程度后空泡閉合方式發生轉變，即由表面閉合轉變為自由液面以下某一點的深閉合.根據本文對不同速度條件下圓柱體入水計算結果分析表明，隨著速度的增加，空泡發生閉合方式轉變的臨界壓強降低.同時可以看出，在空泡閉合時刻，入水點附近自由液面的高度隨著空氣壓強的降低而增加，空泡直徑隨著壓強的降低也呈增加趨勢.圖9為在不同空氣壓強條件下，不同速度垂直入水空泡閉合瞬間空泡最大直徑變化規律.

由圖9可以看出，在同一壓強條件下，以不同速度入水，其入水空泡最大無量綱直徑幾乎趨于某一定值;同時，隨著空氣壓強的降低入水空泡閉合瞬間最大無量綱直徑呈遞增趨勢，在空泡發生閉合方式轉變之前，空泡無量綱直徑隨壓強變化的趨勢較為穩定，當壓強降低到一定程度后，空泡無量綱直徑出現較大幅度的增加，并出現不穩定現象.

圖9 空泡最大直徑變化規律

5 結論

本文對帶150°錐角頭型圓柱體勻速垂直入水空泡形態進行了數值計算研究.計算工況為:入水大氣壓強從(1/8)P0～4P0，入水速度從10 m/s到50 m/s.通過對以上計算結果分析得出以下結論.

1)通過對不同速度和不同壓強條件下入水自然空化的計算結果分析表明，在本文所計算入水速度范圍內，在圓柱體肩部出現小范圍內的局部自然空化，自然空化對整個空泡形態影響較小.

2)隨著空氣域壓強的降低，空泡表面閉合無量綱時間呈增加趨勢變化;當空氣壓強高于1個標準大氣壓強時，空泡無量綱閉合時間隨壓強的變化趨勢較為穩定;當壓強低于1個標準大氣壓強時，空泡無量綱閉合時間變化較為復雜，呈現出較強的非線性.

3)當空氣域壓強降低到一定程度時，空泡閉合方式將由表面閉合向深閉合轉變.

4)空泡閉合瞬間最大無量綱直徑在同一空氣域壓強條件下趨于某一穩定值;并隨著壓強的降低，出現一定的不穩定性.

［1］WORTHINGTON A M，COLE R S.A study of splashes［M］.New York:Longmans Green and Company，1908:13-119.

［2］VON K.The impact on seaplane floats during landing［R］.Washington:National Advisory Committee for Aeronautics，1929:1-8.

［3］MAY A.Review of water-entry theory and date［J］.Journal of hydrodynamics，1970，4(4):140-142.

［4］MAY A，WOODHULL J.C.Drag coefficients of steel spheres entering water vertically［J］.Journal of Applied Physics，1948，19(12):1109-1121.

［5］MAY A.Effect of surface condition of a sphere on its water-entry cavity［J］.Journal of Applied Physics，1951，22:1219-1222.

［6］MAY A.Vertical entry of missile into water［J］.Journal of Applied Physics，1952，23(12):1362-1372.

［7］MAY A.Water entry and the cavity-running behavior of missiles［R］.Virginiai Nowal Surface Weapons Center，1975:1-37.

［8］WAUGH J G.Water-entry pitch modeling［J］.Journal of hydrodynamics，1968，2(2):87-92.

［9］GILBARG D，ANDERSON R A.Influence of atmospheric pressure on the phenomena accompanying the entry of spheres into water［J］.Journal of applied physics，1948，19(2):127-139.

［10］ARISTOFF J M，TRUSCOTT T T，TECHET A H，et al.The water entry of decelerating spheres［J］.Physics of fluids，2010，22(3)，032102:1-7.

［11］DUEZ C，YBERT C，CLANET C，et al.Making a splash with water repellency［J］.Nature physics，2007，3:180-183.

［12］YAN H，LIU Y，JAKUB K，et al.Cavity dynamics in water entry at low Froude numbers［J］.Journal Fluid Mech，2009，641:441-461.

［13］陳學農，何友聲.平頭物體三維帶空泡入水的數值模擬［J］.力學學報，1990，22(2):129-138.

［14］ARISTOFF J M，BUSH J W M.Water entry of small hydrophobic spheres［J］.Journal Fluid Mech，2009，619:45-78.

［15］李玉柱，苑明順.流體力學［M］，北京:高等教育出版社，2003:173-179.

［16］柯乃普，戴利，哈密脫.空化與空蝕［M］.水利水電科學研究院，譯.北京:水利出版社，1981:1-12.