彈體出水引起的非穩態流場數值計算方法

姚熊亮，張成，楊衡，張忠宇

(1.哈爾濱工程大學船舶工程學院，黑龍江哈爾濱150001;2.中國艦船研究設計中心，湖北武漢430064)

流場處于非穩定狀態往往是由于砰擊作用引起的，自然界中存在各種各樣的砰擊現象.人們往往對某些顯而易見的砰擊現象感興趣，如波浪載荷的砰擊作用，但事實上很多常見的物理現象也可能產生瞬間的砰擊作用，如彈體出水時，彈頭與自由面撞擊，就會產生一種砰擊現象.砰擊現象的特點是流場瞬間變化劇烈，瞬間的砰擊作用引起的水動力非常大.這種砰擊作用使結構周圍流場在初始時刻附近處于不穩定狀態，而這種非穩態流場的動壓力值比穩態流場的動壓力值大幾倍甚至幾十倍.另外水下爆炸沖擊波載荷作用下的船體結構流固耦合問題也是以流場瞬間劇烈變化為特點，對船體結構造成非常強烈的沖擊，進而導致非常嚴重的毀傷.但無論是強烈的沖擊問題還是弱砰擊問題，均是以流場瞬間劇烈變化為特點，解決這類問題的關鍵在于必須考慮流場可壓縮性.

如今，解決如水下爆炸這樣的強沖擊問題有很多方法，著名的有延遲勢法［1-2］、雙漸近法(DAA)［3-6］等.這些方法能較為成功的應用于水下爆炸領域，因此一直被沿用至今.由于彈體出水砰擊問題和水下爆炸問題均須考慮流場可壓縮性，因此采用擅長解決水下爆炸問題的雙漸近法，將其應用于彈體出水砰擊問題.但該法存在局限性，即是基于線性化假設提出的一種理論［7］，從理論上來說是不能用于分析航行態或做大幅運動結構的瞬態運動問題.彈體出水恰是典型的具有初始航速的瞬態運動問題，這樣的流固耦合問題不但以“砰擊”為主要特點，還有一重要特點即“非線性”.所以對彈體出水問題的分析，采用傳統的雙漸近法和不考慮可壓縮性的勢流理論方法均不適合.

針對雙漸近法只能用于解決準靜態條件下的流固作用的局限性，在DAA2法的基礎上，將運動非線性引入雙漸近法，提出改進方法，即非線性雙漸近法(NDAA).使雙漸近法可應用于航行態結構或做大幅運動結構的流固耦合問題.并以具有初始航速的彈體為例，采用兩種方法(DAA2以及NDAA)捕捉彈體出水時可能會發生的砰擊現象，模擬彈體砰擊引起的非穩態流場.將NDAA與DAA2法的數值計算結果進行對比，分析考慮非線性效應的流場和彈體出水時可能會發生的物理現象.

1 非線性雙漸近法(NDAA)

浸在水中的結構物與流場的瞬態耦合問題有這樣的特點:早期瞬態(高頻)響應主要是聲輻射問題，即所謂平面波或曲面波近似解問題;后期瞬態(低頻)響應主要是"虛質量"問題.用延遲勢原理可分別得出其早期(高頻)近似式(E.T.A)和后期(低頻)近似式(L.T.A).Geers等(1975～1980)，根據這些特點提出了雙漸近法(DAA)，就是用漸近展開匹配法導出能用于中間頻率的微分方程式.當趨于高頻時其解漸近于早期近似解(E.T.A);在低頻時就漸近于后期近似解(L.T.A)［3］.

雙漸近法(DAA)之所以能廣泛應用于水下爆炸領域，其中一個重要原因是它考慮了流體可壓縮性，考慮可壓縮性的勢函數微分方程式:

DAA法就是從上述控制方程的基礎上經推導得到.DAA法的基本方程如下［4］:

一階雙漸近法方程(DAA1):

二階雙漸近法方程(DAA2):

結構動力方程:

式中:pS為流體中的散射壓力，pi為入射波壓力，Mf為流體質量矩陣，Ωf為流體頻率矩陣，Af為流體單元的面積矩陣，uⅠ為入射波速度，x為結構位移，G為坐標轉換矩陣.

Geer［5］在他的研究中指出，一階雙漸近法DAA1在高頻時相當于平面波近似式(PWA)，在低頻時相當于虛質量近似式(VMA).二階雙漸近法DAA2在高頻時相當于曲面波近似式(VWA)，在低頻時相當于修正的虛質量近似式(VMA)［5］.

但無論對DAA1還是DAA2，在其理論推導的過程中，均引入了一個線性假設，即

式中:P=pS+pi，即總的流場動壓力值，φ為流場速度勢.

式(3)中壓力P和速度勢φ呈線性關系.對高頻動載作用下結構做線性小擾動時的周圍流場進行動壓力計算，這種假設是合理的.因為速度勢的空間導數本身就很小，其二次項就更小了，壓力的線性項即速度勢的時間導數比其空間導數大得多.因此，忽略速度勢隨空間的變化，認為速度勢是時間的函數，這種情況下，對壓力作線性假設是合理的.但水中結構物如果處于航行態或做大幅運動，其非線性效應對水動力值影響比較顯著，速度勢的空間導數就不能忽略了.為此，對DAA2法進行改進，將速度勢的空間導數引入到DAA2法的數值計算，從而將非線性的影響考慮到水中結構物的瞬態流固耦合計算當中.

DAA2法是從延遲勢法基礎上得到，延遲勢的速度勢方程如下:

將式(3)代入式(4)，可得延遲勢的壓力方程:

DAA1方程和DAA2方程均是在延遲勢壓力方程式(5)的基礎上得到，也就是說，DAA1方程和DAA2方程是在經過式(3)線性變換的情況下得到的.

求解DAA方程式(1)，相當于求解式(5)，在通過DAA2方程式(1)的求解可以得到流場散射壓力ps，即延遲勢方程式(5)中的ps.

由式(3)可得

式(6)中的流場散射壓力pS通過求解DAA方程式(1)得到，入射波壓力pi作為已知條件，那流場總速度勢φ即可得解.

在得到流場速度勢的基礎上，可以通過伯努利方程得到流場動壓力計為Pd:

式中:V為結構運動速度.

解式(8)得到的流場動壓力值與解DAA2方程得到的流場動壓力值的本質區別在于式(8)中的第2項和第3項，式(8)中的第2項仍是線性項，但考慮了結構的航速效應，第3項是非線性項，考慮了結構的運動非線性.

求解式(8)得到考慮非線性效應的流場動壓力，將該考慮非線性效應的流場動壓力ps+pi代入結構運動方程中，于是式(2)可寫成如下形式:

將式(6)、(8)、(9)與DAA2方程(1)聯立求解的方法稱為NDAA法.

采用三角形單元對結構進行離散，如圖1所示.在對式(8)的求解中，需要對速度勢方向導數▽φ進行求解，求解思路如下［8］

式中:符號j表示第j個三角形單元，A、B、C表示一個三角形單元的3個節點，x、y、z表示直角坐標系下的坐標，n表示單位法向矢量.

圖1 節點i及周圍單元，j表示其中的一個單元Fig.1 Node i and its surrounding elements，j is element number

通過對式(10)的求解，就可以得到單元j的速度勢梯度▽φ了，其他單元類同.注意這里用(uj，vj，wj)代替▽φ.

2 彈體出水引起的非穩態流場數值模擬

考慮導彈從自由面以下一定深度處以一軸向初速度做出水運動，在彈體橫向存在一定橫向初速度，除受流場動壓力以外，不受其他任何外力作用.研究彈體在出水過程中其周圍流場各節點動壓力值的變化，并將流場動壓力P、時間T、結構運動速度V無量綱化.

無量綱結點運動速度為

無量綱流場動壓力為

無量綱時間為

式中:V0為彈體初始運動速度，l0為初始時刻彈體頂點距自由液面的距離.

圖2 物理模型Fig.2 Physical model

2.1 非穩態流場時間分布

本節從時間的角度分析彈體各點流場動壓力，分析彈體出水過程中流場動壓力的變化、非線性效應不能被忽略的流場狀態以及彈體出水時可能會發生的物理現象.

圖3給出了2幅相互對應的圖，第1幅表示彈體頂點流場動壓力，第2幅表示彈體頂點垂向速度.2幅圖在各個時間段都一一對應.彈體頂點出水前，其流場動壓力以及垂向速度均是從高頻非穩態區線性過渡到低頻穩態區;彈體頂點出水后，彈體頂點作為結構的一部分，仍然存在垂向速度，但此時該點已位于空氣中，流場動壓力為0.

另外在如圖3所示的非穩態區，流場動壓力極大.非線性雙漸近法NDAA和二階雙漸近法DAA2的數值計算結果很接近，即對于高頻非穩態區，非線性效應不明顯.在穩態區，流場動壓力相對較小，非線性雙漸近法NDAA比二階雙漸近法DAA2的數值計算結果偏大，即對于低頻穩態區的流場，非線性效應顯著.另外從圖3中還可以看出，t=0.4是流場非穩態區和穩態區的分界線.

圖3 彈體頂點壓力及其垂向速度時歷曲線Fig.3 The pressure-time curve and the velocity-time curve of the top node of the missile

穩態區所對應的結構垂向速度也趨于穩定狀態，垂向速度V平均值大概0.94左右.如果按不考慮可壓縮性的勢流理論的觀點，彈體以垂向速度在水中做勻速運動(忽略橫向速度對彈體頂點壓力的影響)，彈體頂點的理論壓力值大概在0.896附近.值得注意的是，采用NDAA法，穩態區的壓力平均值大概在0.95左右，采用DAA2法，穩定后的壓力平均值在0附近.可見，DAA2法雖然理論上在結構運動后期趨于勢流理論中的虛質量假設，但是由于它沒有考慮運動非線性，在后期處理有初始航速的結構時，其數值計算結果是與勢流理論解不符.而NDAA法不但在前期與DAA2法結果符合良好，能較好的反映非穩態流場的本質特征，而且在后期趨于勢流理論中的虛質量假設，能較好的反映穩態流場的物理特性.

該節點位于彈體頭部與中部交界處，彈體頂點于t=1.1時開始出水.而該節點此時還未露出水面，但也明顯受到彈體頂點出水的影響，在其壓力時歷曲線上可以觀察到2個較為明顯的壓力峰值，即雙峰值.該節點在t=1.5以后出水，出水后流場壓力降為0.如圖4所示.

圖4 彈體頭部與中部交界處某點壓力時歷曲線Fig.4 The pressure-time curve on the shoulder of the missile

該節點反映了彈體頭部的壓力變化特征，經歷非穩態、穩態、再到非穩態這樣一個過程.第1個非穩態是由瞬間初速度產生的砰擊而引起的，第2個非穩態是由彈頭與自由面砰擊而引起的.可以從該點的壓力曲線上觀察到第1個非穩態對應的穩態過程，無法觀察到第2個非穩態對應的穩態過程.說明彈頭出水時會使流場一直處于非穩定狀態，直至整個彈頭均露出水面.

另外，和彈體頂點類似，非線性效應對非穩態區流場動壓力值影響微弱，對穩態流場動壓力值影響顯著.穩態區流場動壓力平均值的NDAA法數值計算結果在0.85左右，DAA2法的數值計算結果在0左右.

圖5所示節點位于彈體中部距彈體肩部1/4彈體長度的地方.該節點雖然沒有布置在彈體頭部，但其壓力值明顯受到了彈頭出水的影響.彈頭出水時，和圖4類似，在其壓力時歷曲線上也可以觀察到雙峰值現象.由于沒有對雙峰值的產生機理做深入探討，初步判斷，雙峰值可能是由彈頭形狀引起的.因為彈頭出水的時間恰好是雙峰值產生的時間，并且彈頭形狀近似三角形，三角形結構屬斜壁出水而非直壁出水，斜壁浮體的大幅運動問題現在也正是水動力學的重要難點.

圖5 彈體中部壓力時歷曲線Fig.5 The pressure-time curve on the middle of the missile

該節點流場動壓力主要經歷4個階段，即非穩態、穩態、非穩態、穩態.2個非穩態區的產生機理和圖4類似，不再贅述.第2個穩態區對應的是彈體中間的圓柱部分出水，而彈體中間圓柱部分出水屬直壁出水，這意味著結構直壁出水時不會產生砰擊現象.

由圖6可以看出，彈體尾部節點壓力值比其他節點都稍大，穩定值在0.98左右.尾部節點壓力曲線存在有趣而且特殊的物理現象，即其曲線存在顯著的拍振特征.大概是3個周期1個拍，拍振發生的前提條件存在2種相近的頻率特征.對于彈體來說，這2種頻率可能來源于彈體尾部的局部振動以及流場動壓力頻率，而這2種頻率又接近，導致拍振的發生.這就意味著彈體某些局部結構可能會發拍振現象.拍振現象的發生可能導致彈體局部結構振幅很大，振動強烈，進而導致結構破壞.

圖6 彈體尾部壓力時歷曲線Fig.6 the pressure-time curve on the bottom of the missile

總結圖3～6，彈體出水需經歷4個階段，即非穩態、穩態、非穩態、穩態這4個過程，如圖7所示.

圖7 彈體出水時流場物理特性變化過程及其成因Fig.7 Changing process of the fluid characters andthe causes at missiles flying out of the water surface

2.2 非穩態流場空間分布

針對不同階段的流場動態特性，取分別表征4個流場狀態(非穩態、穩態、非穩態、穩態)的4個特征時間點(初始時刻附近，t=0.45;出水前，t=0.9;出水瞬間，t=1.27;出水后，t=3).采用非線性雙漸近法(NDAA)對彈體周圍流場動壓力的空間分布做詳細分析.

將彈體中間圓柱部分沿軸線從中間切開，展開成矩形，圖8、9給出了彈體中間圓柱部分周圍流場動壓力沿彈體軸向和周向的分布情況.用(X、Y、Z)分別表示彈體各點軸向、橫向、法向坐標，坐標原點位于彈體頂點，并且軸向正方向指向彈體尾部，L表示彈體長度，R表示彈體中間圓柱部分的最大半徑，θ表示彈體中間圓柱部分的周向角度，θ=0對應彈體橫向迎流面的中線.

圖8 非穩態流場(t=0.45)與穩態流場(t=0.9)對比Fig.8 Comparison between unsteady flow field(t= 0.03 s)and the steady flow field(t=0.06 s)

圖9 非穩態流場(t=1.27)與穩態流場(t=3)對比Fig.9 Comparison between unsteady flow field(t= 1.27)and the steady flow field(t=3)

無論是彈體初始速度引起的非穩態流場，還是彈頭出水時與自由面砰擊而產生的非穩態流場，水動壓力值均比對應的穩態流場壓力值大得多.非穩態流場動壓力值沿空間的分布也相對復雜，尤其是沿彈體軸向.彈體屬細長體，軸向受力不均容易使彈體產生較大的彎矩，進而導致彈體折斷，說明非穩態流場數值模擬十分必要.

另外在初始時刻附近，可以觀察到較為明顯的低壓區(見圖10)，分布在彈體肩部附近.低壓區的壓力極小，很可能在彈體肩部形成空化現象［9］，即產生空泡，空泡區域的大小跟彈體周圍壓力的分布以及空化壓力有關［10-11］.在出水瞬間，低壓區開始下移即空泡區下移，也就是說，彈體出水可能導致空泡形成的區域以及位置發生改變.如果空化壓力確定以后，可以通過流動壓力的空間分布來預測空泡區域的大小及其分布位置.這里不對彈體的空化現象做深入討論.

圖10 流場動壓力沿彈體軸向分布(θ(2π)=0)Fig.10 Dynamic pressure along the axial direction (θ(2π)=0)

圖11 流場動壓力沿彈體周向分布(X/L=0.5)Fig.11 Dynamic pressure along the circumferential direction(X/L=0.5)

低壓區即空泡區下移也可以從流場動壓力沿軸向分布圖中看出來(見圖11).初始時刻的流場動壓力比出水瞬間流場動壓力稍大，也就是說，彈體瞬間初速度引起的砰擊比彈頭出水引起的砰擊要強烈.但二者均比各自對應的穩態流場壓力大得多，出水前的穩態流場壓力比出水后稍大.彈體壓力沿周向分布較規律，接近正余弦分布，尤其是初始時刻附近，流場動壓力沿軸向基本呈均值為0的余弦分布.

由圖12可以看出，初始時刻附近彈體頭部流場動壓力關于軸線呈對稱分布.橫向運動幾乎不對頭部壓力分布造成影響，對頭部壓力主要造成影響的是垂向運動.因為在彈體垂向，彈體突然以某一初速度出水.而彈體頭部恰好是迎流面，從而產生劇烈的砰擊，流場處于不穩定狀態，導致彈體頭部流場動壓力值極大.t=0.45的彈頭壓力系數最大值可達到5.5.如果按不考慮可壓縮性的勢流理論的觀點，彈體以某一初速度在水中做勻速運動，彈體頂點的壓力系數大概在1左右.可見，如果采用不考慮可壓縮性的勢流理論方法，是無法對這種高頻高幅非穩態流場進行數值模擬的.

圖12 初始時刻附近(t=0.45)彈體頭部流場動壓力Fig.12 The hydrodynamic on the head of the missile at the beginning(t=0.45)

圖13 出水前(t=0.9)彈體頭部流場動壓力Fig.13 The hydrodynamic on the front part of the missile before flying out of the water surface(t=0.9)

由圖13可以看出，出水前彈體頭部穩態流場關于軸線呈非對稱分布，關于橫向對稱.這是因為流場處于穩定狀態時，彈體頭部的砰擊影響越來越微弱，可壓縮性影響越來越小，流場動壓力分布越來越趨于勢流理論解.從圖中還可以看出，彈體頭部壓力場分布是橫向速度和垂向速度共同作用的結果，這與勢流理論解是一致的.這也說明非線性雙漸近法在后期趨于勢流理論的虛質量假設.

圖14 出水瞬間(t=1.27)和出水后(t=3)彈體頭部流場動壓力Fig.14 The hydrodynamic on the front part of the missile when flying out of the water surface(t=3)

圖14為部分彈頭出水，出水部分壓力為0，未出水部分壓力很大.而部分彈頭出水，正是砰擊較為強烈的時刻，整個彈頭出水后，砰擊現象也就消失了，說明只有彈頭部分與自由面產生砰擊作用.彈體中間的圓柱部分出水時，不會產生砰擊現象.

3 結論

在DAA2法的基礎上，將運動非線性引入雙漸近法，即非線性雙漸進法(NDAA)，比較了兩種方法的數值計算結果，分析非線性對彈體出水流固耦合動態特性的影響.簡單探討了彈體出水時雙峰值的產生機理以及彈體出水時其局部結構可能會發生拍振的現象.探討了非穩態流場的成因，給出了不同階段的流場動態特性，以及各個特征時間點流場動壓力的空間分布情況.得到主要結論如下:

1)在非穩態區，流場動壓力極大，非線性效應不明顯;在穩態區，非線性效應顯著.對于具有一定初始航速的結構，如彈體出水的流固耦合問題，應該考慮非線性效應.

2)彈體出水過程中可觀察到砰擊現象，第1次砰擊發生在彈體以某一初速度瞬間啟動時，第2次砰擊發生在彈頭出水時.砰擊時，彈體周圍流場處于不穩定狀態.彈體整個出水過程可以概括成4個物理狀態，即非穩態、穩態、非穩態、穩態.非穩態流場動壓力值比穩態流場大得多.

3)NDAA法在前期與雙漸近法結果符合良好，能較好的反映非穩態流場的本質特征，在后期趨于勢流理論中的虛質量假設，能較好的反映穩態流場的物理特性.

4)在研究彈體出水過程中還觀察到了2個重要的物理現象:一是雙峰值現象，二是拍振現象.

5)彈體未出水時，彈體周圍低壓區位于彈體肩部，這可能是肩空泡產生的原因.隨著彈體頭部開始出水，低壓區開始下移，這可能導致空泡區下移.

［1］HUANG H.Transient interaction of plane acoustic waves with a spherical elastic shell［J］.J Acoust Soc Amer， 1969，45:661-670.

［2］DYKA C T，INGEL R P.Transient fluid-structure interaction in naval applications using the retarded potential method［J］.Engineering Analysis with Boundary Elements，1998，21:245-251.

［3］GEERS T L.Residual potential and approximation methods for three dimensional fluid-structure interaction problems［J］.J Acoust Soc Amer，1971，49:1505-1510.

［4］GEERS T L.Doubly asympotic approximation for transient motions of submerged structures［J］.J Acoust Soc Ame，1978，64:1500-1508.

［5］GEERS T L，FELIPPA C A.Doubly asymptotic approximations for vibration analysis of submerged structures［J］.J Acoust Soc Amer，1980，73:1152-1159.

［6］GEERS T L，ZHANG Peizhane.Advanced DAA methods for shock response analysis［R］.Denver:University of Colorado，1992:252-696.

［7］SPRAGUE M A，GEERS T L.Response of empty and fluidfilled submerged spherical shells to plane and spherical，step-exponential acoustic waves［J］.Shock and Vibration，1999(6):147-157.

［8］ZHANG Aman，YAO Xiongliang.The interaction between multiple bubbles and free surface［J］.Chinese Physics B，2008，17(3):927-938.

［9］權曉波，魏海鵬.潛射導彈大攻角空化流動特性計算研究［J］.宇航學報，2008，29(6):1701-1705.

QUAN Xiaobo，WEI Haipeng.Numerical simulation on cavitation of submarine launched missile＇s surface at large angles of attack［J］.Journal of Astronautics，2008，29 (6):1701-1705.

［10］羅金玲，毛鴻羽.導彈出水過程中氣/水動力學的研究［J］.戰術導彈技術，2004，7(4):23-25.

LUO Jinling，MAO Hongyu.Research on aerodynamics and hydrodynamics in exiting water process of missile［J］.Tactical Missile Technology，2004，7(4):23-25.

［11］REN Bing，LI Xuelin，WANG Yongxue.Experimental investigation of instantaneous properties of wave slamming on the plate［J］.China Ocean Engineering，2007，21(3): 533-540.