跨介質(zhì)航行器高速斜入水跳彈現(xiàn)象研究

周可，黃振貴，陳志華，劉想炎，王浩

（南京理工大學(xué) 瞬態(tài)物理國家重點實驗室，南京 210094）

隨著跨介質(zhì)航行器的發(fā)展，近來越來越多的國家開始研究跳彈現(xiàn)象，能在飛行過程中反復(fù)出入水面，對作戰(zhàn)計劃的隱蔽性有很大的提升。跳彈現(xiàn)象也已經(jīng)在軍事上得到應(yīng)用，也就是“跳彈桅頂轟炸”。早在二戰(zhàn)時期，W.Barnes就利用跳彈轟炸德國魯爾水壩，隨后在太平洋戰(zhàn)爭中，美國也效仿英國的跳彈，對日軍艦隊進(jìn)行轟炸。相比于高空投射，低空的跳彈有著更好的隱蔽性和精準(zhǔn)性，因此航行器的跳彈現(xiàn)象也得到了更大的關(guān)注，如果能較好地掌握這一現(xiàn)象的規(guī)律，就能開發(fā)出一種突防能力更好武器。

入水問題的研究最早開始于19世紀(jì)末，Worthington等利用閃光照相機對剛性球體垂直入水過程中發(fā)生的一系列現(xiàn)象進(jìn)行了試驗研究。對于物體入水穩(wěn)定性方面的研究，侯宇等通過試驗探究了射彈小角度高速入水時產(chǎn)生的超空泡對射彈的彈道軌跡的影響，以及高負(fù)載對彈體的損傷。胡青青探討了4種頭型、3種長徑比射彈高速入水時的穩(wěn)定性和空泡形態(tài)，并分析了其對入水過程的影響。黃鴻鑫等通過數(shù)值模擬方法，從空泡形態(tài)、速度衰減以及俯仰角變化探討了射彈頭部形狀和質(zhì)心位置對高速入水穩(wěn)定性的影響。鄭磊基于STAR-CCM+和ABAQUS聯(lián)合仿真對彈體進(jìn)行了垂直和傾斜入水雙向耦合計算，分析了在入水過程中彈體的動力學(xué)特性、彈體表面變量以及應(yīng)力應(yīng)變情況。

跳彈現(xiàn)象則是物體入水失穩(wěn)方面的一種表現(xiàn)，對于跳彈問題的系統(tǒng)性研究最開始以剛性球體為研究對象。Richardson研究了圓球傾斜入水跳彈的問題，討論了當(dāng)球體傾斜入水后彈跳時所涉及的機制和力。Miloh等通過數(shù)學(xué)方程的推導(dǎo)得出了橢圓球體入水發(fā)生跳彈的臨界角的經(jīng)驗表達(dá)式。Soliman等用鋼球和硬鋁球射入淺水深度和干沙中，總結(jié)了不同介質(zhì)下跳彈發(fā)生的臨界角，臨界角隨著速度的增加而減小，但存在一個截止角，在任何速度下都不會發(fā)生彈跳。Johnson探究了旋轉(zhuǎn)和非旋轉(zhuǎn)球形彈丸對跳彈行為的影響。隨著研究的進(jìn)一步深入，圓柱體以及射彈等回轉(zhuǎn)體的入水跳彈現(xiàn)象開始被關(guān)注。Nguyen等通過數(shù)值模擬計算，探究了圓柱體入水時的角度和速度對跳彈行為的影響，總結(jié)出了每個入水角都有一個臨界速度發(fā)生跳彈。Hutchings提出了一種彈跳理論，該理論允許考慮彈丸自旋的影響，計算了球體和旋轉(zhuǎn)圓柱體的彈跳臨界角，并將旋轉(zhuǎn)圓柱體的理論應(yīng)用于第二次世界大戰(zhàn)期間開發(fā)的“巴恩斯·沃利斯”跳彈炸彈。Moxnes等通過試驗分析和數(shù)值模擬計算總結(jié)了不同類型的射彈在水中發(fā)生跳彈的規(guī)律。Gold等利用UZI沖鋒槍發(fā)射9 mm子彈，探究了相同速度下不同角度子彈高速入水下的跳彈規(guī)律。

綜上所述，目前對于物體跳彈現(xiàn)象的研究著重于球體、圓柱體和射彈，而對航行器入水時發(fā)生跳彈現(xiàn)象的研究還較少，特別是高速入水。公開的文獻(xiàn)僅有Chen等通過數(shù)值模擬計算了不同頭型的航行器入水后的跳彈行為，得出錐鼻形頭部更容易改變軌跡。袁緒龍等提出了跨介質(zhì)航行器高速入水過程中產(chǎn)生的巨大法向過載和彎矩。李永利等利用數(shù)值模擬技術(shù)研究了航行器在不同頂角下小角度低速入水時的跳彈現(xiàn)象。陳國明等通過試驗和數(shù)值仿真研究了航行器發(fā)生跳彈現(xiàn)象的條件和規(guī)律。本文通過數(shù)值模擬方法，模擬了航行器在高速斜入水后產(chǎn)生跳彈現(xiàn)象的運動過程，討論并分析了不同入水角、不同速度對高速傾斜入水時航行器跳彈行為的影響，揭示產(chǎn)生跳彈的機理，為航行器高速入水跳彈的設(shè)計提供參考。

1 數(shù)值計算理論

1.1 流體基本控制方程

VOF多相流模型是建立在歐拉網(wǎng)格下的界面追蹤方法。該方法中，互不相容的流體組分共用一套動量方程，并通過引入相體積分?jǐn)?shù)這一變量來實現(xiàn)對計算域內(nèi)相間界面的追蹤，更加適用于捕捉流體自由面的變化和求解氣液相的比例。其中表示相的體積占所在網(wǎng)格體積的比值，0表示為空氣相，1表示為水相，0<<1表示為交界面。VOF方法就是通過求出整個計算域內(nèi)各網(wǎng)格單元的相分?jǐn)?shù)，從而構(gòu)建出界面，如圖1所示。

圖1 VOF交界面示意圖Fig.1 Interface of the VOF: a) phase fractional spatial distribution; b) interface reconstruction

對于每個控制單元，其連續(xù)性方程為：

式中：為整體網(wǎng)格單元內(nèi)的平均數(shù)學(xué)密度。VOF多相流模型中，密度的定義為=+，和分別為空氣和水的密度，和為水和空氣在單元格的密度。

動量守恒方程為：

式中：為網(wǎng)格單元的動力黏性系數(shù)；為網(wǎng)格單元的湍流黏性系數(shù)，=+。

兩相之間采用的質(zhì)量傳輸方程為Schnerr-Sauer空化模型。該空化模型的計算效率高，穩(wěn)定性強，模擬效果好。相關(guān)計算表達(dá)式為：

式中：分別代表蒸發(fā)率和冷凝率；、分別為蒸汽相的體積分?jǐn)?shù)和密度；為氣泡半徑；為單位體積空泡數(shù)。

1.2 湍流模型

采用Realizable-模型，它是基于標(biāo)準(zhǔn)-模型發(fā)展而來的，基本思想和-模型一致。其用某種數(shù)學(xué)方式建立湍動黏度和湍動能以及湍流耗散率的數(shù)學(xué)聯(lián)系，并分別以湍動能和湍動耗散率建立運輸方程，最終實現(xiàn)方程的封閉。

湍動黏度的表達(dá)式為：

湍動能和湍動耗散率的輸運方程為：

1.3 六自由度彈道模型

本文利用DFBI模塊耦合重疊網(wǎng)格模擬剛體復(fù)雜的六自由度運動。在全局慣性坐標(biāo)系下，剛體繞質(zhì)心的運動方程為：

在局部坐標(biāo)系下，DFBI旋轉(zhuǎn)運動方程為：

1.4 重疊網(wǎng)格技術(shù)

重疊網(wǎng)格技術(shù)是對計算區(qū)域和部件區(qū)域單獨劃分設(shè)計網(wǎng)格，主要由背景網(wǎng)格和組件網(wǎng)格相互重疊而成，如圖2所示。兩者之間會有一部分重疊區(qū)域，但不存在連通關(guān)系，在初始化的時候會進(jìn)行布爾運算，計算域外的網(wǎng)格會被刪除并排除在外，使該部分網(wǎng)格不進(jìn)行數(shù)據(jù)交換。根據(jù)布爾運算剩余的重疊網(wǎng)格區(qū)域，建立插值關(guān)系，并且建立相應(yīng)的數(shù)學(xué)關(guān)系進(jìn)行數(shù)據(jù)交換。相比動網(wǎng)格技術(shù)中持續(xù)的網(wǎng)格更新，重疊組件中組件網(wǎng)格隨物體一起運動，不需要網(wǎng)格變形，也不需要重新生成網(wǎng)格，只需要定義組件運動規(guī)律和網(wǎng)格重疊邊界條件，即可實現(xiàn)對物體在流域中運動的流體力學(xué)問題的模擬，避免因網(wǎng)格變形、網(wǎng)格單元生成或消失帶來的計算誤差，具有較大的優(yōu)越性。

圖2 重疊網(wǎng)格插值示意圖Fig.2 Overset interpolation diagram

2 數(shù)值模擬方法

2.1 計算模型

計算模型的幾何尺寸如圖3所示，全長=1.28 m，彈徑=0.128 m，質(zhì)量為16.991 kg，分布均勻，頭部為30°的拱形，尾部為線性截斷體。

圖3 航行器物理模型Fig.3 The physical model of the vechicle

采用的計算軟件為STAR-CCM+，基于N-S方程，選用Realizable-湍流，歐拉多相流模型，用VOF流體域體積模型捕捉自由液面的變化和Schnerr and Sauer空化模型模擬空化過程，引用重疊網(wǎng)格技術(shù)耦合六自由度方程來捕捉航行器的入水運動。為了減少計算量，本文預(yù)估了航行器的大致軌跡，并在軌跡附近的網(wǎng)格進(jìn)行加密作為計算域。網(wǎng)格截面如圖4所示，軸由紙面向外。計算域如圖5所示，長10，寬3.9，水域深2.7，空氣域為1.2。計算域四周以及底部為壓力出口，壓力根據(jù)水深依據(jù)公式=設(shè)置，頂部設(shè)置為速度進(jìn)口，速度幅值為0。湍流運輸方程采用二階迎風(fēng)格式，-湍流求解選用STAR-CCM+內(nèi)嵌的AMG線性求解器，松馳方案選用高斯–賽德爾。時間上采用隱式非穩(wěn)態(tài)，時間步長為2×10s。分別對航行器在初始入水角為15°和入水速度為100、150、200 m/s，以及初始速度為150 m/s，入水角為8°、15°、25°一共5個工況進(jìn)行數(shù)值模擬計算。

圖4 Z=0截面網(wǎng)格Fig.4 Schematic diagram of Z=0 section mesh

圖5 計算域及邊界條件設(shè)置Fig.5 Computational domain and boundary condition settings

2.2 數(shù)值計算方法驗證

為了驗證本次數(shù)值模擬方法的有效性，本文選取參考文獻(xiàn)[22]中的試驗結(jié)果作為對照。本次驗證模型采用文獻(xiàn)中直徑為6 mm、長為24 mm的回轉(zhuǎn)圓柱體，質(zhì)量為4.88 g，入水速度為106.6 m/s，水箱尺寸為310 mm×310 mm×660 mm，水深為600 mm。在數(shù)值模擬中，模型質(zhì)量、幾何尺寸、入水速度與試驗保持一致。通過對比圓柱體的速度衰減和位移變化以及不同時刻的空泡演化來驗證模擬的有效性。

速度曲線的數(shù)值模擬和試驗結(jié)果對比如圖6a所示，最大誤差約在5%，并且在入水時刻速度的劇烈變化，數(shù)值模擬能很好地捕捉。圓柱體的位移變化曲線如圖6b所示，最大誤差也約在5%，在誤差允許范圍之內(nèi)。

圖6 試驗數(shù)據(jù)與模擬數(shù)據(jù)對比Fig.6 Comparison of experimental data and simulation data:a) velocity curve; b) displacement curve

圓柱體入水時在=0.8 ms和=2 ms空泡輪廓的對比如圖7所示。圖7清楚地展示了空泡的產(chǎn)生以及閉合2種狀態(tài)。在=0.8 ms時，空泡完全包裹了圓柱體。隨著圓柱體繼續(xù)的下落，空泡逐漸的被拉長。在=2 ms時，空泡出現(xiàn)了閉合。通過與試驗圖像的對比可以發(fā)現(xiàn)，模擬結(jié)果與試驗結(jié)果較為一致，進(jìn)一步驗證了模擬方法的有效性。

圖7 空泡演變圖Fig.7 Cavitation evolution: a) experiment; b) simulation

為了避免網(wǎng)格數(shù)量對計算結(jié)果產(chǎn)生的影響，在數(shù)值計算之前建立了3種不同的網(wǎng)格密度來對網(wǎng)格無關(guān)性進(jìn)行檢驗?；诓煌W(wǎng)格密度計算下航行器入水深度的變化曲線如圖8所示?？梢钥闯?，在下潛階段，網(wǎng)格密度的影響不大，在后期上升過程中，不同網(wǎng)格密度之間有一定的差距，隨著網(wǎng)格的加密，差距也越來越小。網(wǎng)格數(shù)量330萬以上的時候，曲線幾乎重合，相差較小。從計算成本和計算結(jié)果的準(zhǔn)確性綜合考慮，本文采用330萬網(wǎng)格數(shù)量進(jìn)行數(shù)值計算。

圖8 不同網(wǎng)格密度下航行器深度變化Fig.8 Vehicle depth variation at different mesh densities

4 結(jié)果及討論

4.1 航行器跳彈機理分析

通過數(shù)值模擬計算航行器在150 m/s，入水角為15°工況下的跳彈現(xiàn)象來闡述跳彈發(fā)生的機理，整個跳彈現(xiàn)象過程如圖9所示。首先，航行器的頭部接觸水面時，頭部會受到很大的沖擊載荷，同時產(chǎn)生噴濺和濺水。因為航行器頭部為卵圓形，使得頭部兩側(cè)受力不均勻，其中下半部遠(yuǎn)大于上半部，以此產(chǎn)生對航行器較大的俯仰力矩，如圖10所示。此力矩會使得航行器產(chǎn)生一個逆時針方向的俯仰角速度，它隨時間的變化曲線如圖11所示?？梢钥闯?，俯仰角速度在入水后呈現(xiàn)增加的趨勢，且方向一直為正方向，航行器的俯仰角增大，彈道發(fā)生彎曲。航行器頭部沾濕面積逐漸增加，使得航行器受到的流體動力變大，俯仰力矩持續(xù)增加。隨著航行器繼續(xù)前行，在=0.128 s時，尾部也開始出現(xiàn)沾濕。因為頭部和尾部的沾濕面出現(xiàn)在質(zhì)心兩端，所以兩端對質(zhì)心的力矩是反向的，使得航行器在該節(jié)點出現(xiàn)力矩減小，致使角速度也減小。最后航行器兩端都受力的情況下達(dá)到力矩平衡，維持在386 N·m左右。此后俯仰角速度大小呈平穩(wěn)增長趨勢，航行器以較為平穩(wěn)的速度抬升，然后沖出水面。當(dāng)航行器頭部離開水面瞬間，頭部受到的流體動力消失，俯仰力矩開始迅速減小。

圖9 航行器入水出水過程Fig.9 Schematic diagram of the process of entering and exiting water of the vehicle

圖10 俯仰力矩變化曲線（v=150 m/s，α=15°）Fig.10 Change of pitch moment curve (v=150 m/s, α=15°)

圖11 俯仰角速度變化曲線（v=150 m/s，α=15°）Fig.11 Change of pitch angular velocity curve(v=150 m/s, α=15°)

通過上述分析可將航行器的跳彈運動過程分為3個階段：第一階段為入水，航行器頭部撞擊水面，頭部下緣受力，給航行器一個逆時針方向的俯仰力矩；第二階段為浸水，此時航行器兩端的沾濕面積逐漸增加，使得航行器俯仰力矩出現(xiàn)微波動直至穩(wěn)定；第三階段為出水，此階段航行器俯仰力矩維持穩(wěn)定，航行器以較為平穩(wěn)的速度沖出水面。

從圖9可以看出，航行器的空泡經(jīng)歷了入水撞擊、空泡形成、開空泡和空泡閉合幾個階段。航行器頭部撞擊液面，形成了不對稱空泡。這是因為撞擊水面時，下表面的沾濕面積大于上表面，所以下表面的空泡形成快于上表面，同時尾部只有下表面沾濕，使得航行器兩側(cè)受力不同。隨著航行器的姿態(tài)和運動軌跡的變化，空泡變得越來越不對稱。同時，可以觀察出空泡呈現(xiàn)先擴張、后收縮的趨勢，空泡后半段在擴張到最大直徑后，水壓開始擠壓空泡壁，空泡收縮并開始閉合。

4.2 航行器初速度對跳彈行為的影響

航行器在不同入水速度下水相的變化如圖12所示。從圖12中可以直觀地看出3種工況下航行器在水中發(fā)生跳彈的全過程對比，不同速度下航行器的出水時間有著明顯區(qū)別，隨著速度的增加，出水時間節(jié)點就越早。其深度和水平位移如圖13所示（圖13b中箭頭標(biāo)注為航行器出水時間節(jié)點），也能看出3種入水速度下出水時間的差異，并且航行器最大下潛深度隨著速度的改變變化不顯著，但出水時的水平位移隨著速度的增加而減小。速度的變化趨勢如圖14所示，不同入水速度下速度變化的斜率有所不同，初始速度為200 m/s的速度變化最快。由表1可知，雖然三者速度改變的快慢不同，但是在出水時它們的出入水速度比幾乎一致，這是由于它們出水時間不一致所致。

表1 不同入水速度下入水出水速度、時間變化Tab.1 Change of water inflow and outflow speed and time under different water entry speeds

圖12 不同速度工況下的水相圖Fig.12 The water phase diagram of different speed conditions

圖13 不同速度工況的位移變化曲線Fig.13 Change of displacement curve under different speed conditions a) Vertical displacement; b) Horizontal displacement

圖14 不同速度工況的速度變化曲線Fig.14 Change of velocity curve under different speed conditions

俯仰力矩的變化如圖15所示?？梢钥闯?，入水階段的俯仰力矩峰值隨著速度的提升而明顯增加，在從浸水階段進(jìn)入出水階段的過渡時間明顯縮短。初始速度越快，在該階段越能快速穩(wěn)定俯仰力矩，這有利于縮短出水時間。最后在出水階段時，俯仰力矩隨著速度的增加而提升。俯仰角速度變化的如圖16所示。不同速度工況下，除了在浸水階段俯仰角速度有一小段下降外，在出水前的其他時間段呈增加趨勢。這是因為尾部受力增大，使得俯仰力矩為負(fù)，從而使得俯仰角速度減小。最后俯仰角速度發(fā)生驟降，標(biāo)志著航行器頭部已經(jīng)離開了水面。由此可知，入水速度對跳彈現(xiàn)象有著明顯的影響，并且可推測當(dāng)初始速度低于某一臨界點時，航行器在該角度就無法跳出水面。

圖15 不同速度工況下俯仰力矩的變化曲線Fig.15 Change of pitch moment curve under different speed conditions

圖16 不同速度工況下俯仰角速度的變化曲線Fig.16 Change of pitch angular velocity curve under different speed conditions

4.3 不同入水角對跳彈行為的影響

在航行器初始速度為150 m/s下研究入水角對跳彈行為的影響，其仿真結(jié)果水相圖如圖17所示，俯仰角角度大小變化如圖18所示?？梢钥闯?，隨著航行器入水角的增大，最后出水時的角度也有著明顯增大。由表2可知，入水角越大時，角度改變量也會明顯增加。

表2 不同入水角入水出水時角度的變化量Tab.2 Variation of angle when entering and exiting water at different entry angles

圖17 不同入水角度工況下的水相圖Fig.17 Water phase diagrams at different water entry angle angles

圖18 不同入水角度工況下俯仰角的變化曲線Fig.18 Change of pitch angle curve under different water entry angle conditions

不同入水角下航行器質(zhì)心垂直位移和水平位移變化趨勢如圖19所示。從圖19a可以看出，3種入水角中，航行器質(zhì)心垂直位移有著較大差異，其中25°入水角的最大垂直位移幾乎為8°入水角的8倍。圖19b中，3種工況在出水前的水平位移較為一致，但是因為出水時間的不同，使得它們的出水點位移有著較大差別，入水角為8°、15°、25°對應(yīng)的出水點位移分別為3.31、6.43、8.78 m。由此可知，入水角度對航行器下潛深度的影響比初始速度大。

圖19 不同入水角度工況下位移的變化曲線Fig.19 Change of displacement curve under different water entry angle conditions

進(jìn)一步分析不同入水角帶來的差異性。入水瞬間，25°和15°入水角工況的俯仰力矩峰值相差不大，但是都遠(yuǎn)大于8°的工況。在浸水階段，不同角度的俯仰力矩變化趨勢幾乎一致，在該階段不同角度的俯仰力矩幾乎都維持在386 N·m，而8°入水角的浸水階段時間較短，俯仰力矩還未穩(wěn)定就進(jìn)入了出水階段，離開水面。從圖21也能看出，這一階段俯仰角速度的增加較為勻速，航行器偏轉(zhuǎn)較為平穩(wěn)。航行器在剛?cè)胨畷r，角速度迅速增加，當(dāng)尾部觸碰空泡壁沾濕后，3個工況角速度均出現(xiàn)下降過程。隨著航行器姿態(tài)的改變，入水角為15°和25°的工況，航行器頭尾兩端作用力相互作用，使得力矩穩(wěn)定，角速度平穩(wěn)增加，直至最后離開水面；而入水角為8°的工況因為沒有經(jīng)過力矩平衡這一過程，在角速度下降后，又迅速增加，然后離開水面，整個過程中角速度大小隨入水角的減小而增大。

圖20 不同入水角度工況下俯仰力矩的變化曲線Fig.20 Change of pitch moment curve under different water entry angle conditions

圖21 不同入水角度工況下俯仰角速度的變化曲線Fig.21 Change of pitch angular velocity curve under different water entry angle conditions

5 結(jié)論

本文采用數(shù)值模擬的方法對跨介質(zhì)航行器的小角度入水的跳彈過程進(jìn)行了研究，分別研究了入水速度和入水角的變化對跳彈行為帶來的影響，分析對比了不同工況下航行器的位移、俯仰力矩和俯仰角速度變化以及出水時間的差異，得到以下結(jié)論：

1）航行器高速入水時，發(fā)生跳彈的主要原因是航行器上下以及頭尾兩側(cè)受力不均，形成了向上的偏轉(zhuǎn)力矩，改變了航行器的俯仰角，會形成不對稱空泡。整個跳彈過程可以分為入水階段、浸水階段和出水階段3個階段.

2）相同入水角，不同入水速度下，速度越大，越容易發(fā)生跳彈現(xiàn)象，出水時間也越早，航行器出水時的速度約為入水速度的70%。航行器在水中下潛的深度受速度的影響不大，但是在出水點位移隨著速度的增加而減小。入水速度越大，航行器所受到到的俯仰力矩越大，俯仰角速度越大。

3）相同入水速度，不同入水角工況下，入水角越小，越能促進(jìn)跳彈行為的發(fā)生，航行器出水時角度改變量越小。入水角的改變對航行器位移的變化更加明顯，隨著入水角的增加，航行器的下潛深度和出水點位置也隨之增大。