倪寶玉，孫士麗，孫龍泉，張 成

(1.哈爾濱工程大學 船舶工程學院，哈爾濱 150001;2.上海交通大學 船舶海洋與建筑工程學院，上海 200240;3.中國艦船研究設(shè)計中心，武漢 430064)

水下航行器出、入水過程涉及到流固耦合問題，自70年代以來，在數(shù)值研究流體結(jié)構(gòu)相互作用的方法上取得了很大的進展，并深入到許多復(fù)雜的流體結(jié)構(gòu)系統(tǒng)。其研究方法可分為弱耦合法和強耦合法［1］。強耦合方法則是對流體和結(jié)構(gòu)分別建立運動方程，采用數(shù)值方法將其耦合求解。弱耦合法即對流體作簡化，將流體對結(jié)構(gòu)的作用表達為流體力的形式，如附加質(zhì)量力等。該方法形式簡單，計算量相對較少，在工程中得到廣泛的應(yīng)用［2］。

用附加質(zhì)量度量流體慣性大小，它是影響水下航行器的水動力特性和動載荷特性的重要參數(shù)之一［3－4］，附加質(zhì)量和阻尼力的計算是直接影響到水下結(jié)構(gòu)運動的關(guān)鍵性問題［5－7］。球體等簡單形狀物體的附加質(zhì)量可通過解析解求得，復(fù)雜形狀物體的附加質(zhì)量通常通過估算法或Hess－Smith［8］方法求得。水下航行器入水過程中，由于流場的作用引起的附加質(zhì)量會改變結(jié)構(gòu)的動力學和運動學特性［9］。賈亮等［10］把結(jié)構(gòu)簡化為梁單元，主要對附加質(zhì)量對結(jié)構(gòu)彎矩的影響進行了研究。李明等［9］將航行體簡化為梁模型，對出水過程中附加質(zhì)量的變化航行體橫向振動的影響進行了研究。王科等［11］應(yīng)用邊界積分方程的方法研究了FPSO型采油平臺的附加質(zhì)量和阻尼力問題。袁振偉等［12］分析了各種因素對動力附加質(zhì)量和動力附加阻尼的影響。

飛行器入水過程中附加質(zhì)量是隨入水深度發(fā)生變化的，在這個過程中附加質(zhì)量的改變有多大，附加質(zhì)量隨時間改變的快慢如何，其對結(jié)構(gòu)的運動參數(shù)有多大影響，目前相關(guān)文獻研究較少。為此，本文首先在分析國外文獻的基礎(chǔ)上，給出含有附加質(zhì)量變化率的入水砰擊力簡約表達式。其次，給出考慮附加質(zhì)量變化的結(jié)構(gòu)運動方程，考慮由于附加質(zhì)量的變化產(chǎn)生的阻尼系數(shù)的改變項。最后以彈體入水為例，分析彈體入水過程中剛體運動、彈性振動引起的附加質(zhì)量的變化情況，以及相應(yīng)砰擊力的變化。

1 落水砰擊力的解析解

為了分析方便，本文首先研究物體以某一定常速度垂直落水情況下，其對應(yīng)的水動力解析解。如圖1所示，定義一個迪卡爾坐標系O-XYZ，原點O位于未擾動的自由液面上，Z軸豎直向上。物體以一豎直速度W(W為負值時表示物體向下運動)運動的情況。假定流體滿足不可壓縮、無旋、無粘的假設(shè)，可應(yīng)用勢流理論進行求解。假設(shè)在自由液面SF和無窮遠邊界S∞速度勢φ=0。

圖1 物體垂直入水運動示意圖Fig.1 The sketch map of vertical entry in the water of a body

為了能夠?qū)⑽矬w所受的水動力表達成含有附加質(zhì)量變化率的形式，本文從能量觀點［13－14］獲得砰擊力。流場中的動能E可寫成如下形式:

式中:ρ為流體密度，φ為速度勢，▽為哈密爾頓算子，R為S包圍的整個流域，S為包含物體濕表面Sb(t)、自由面SF和無窮遠邊界面S∞在內(nèi)的整個流域的邊界面，n為流域邊界的法向(指向流域外為正)?？紤]到自由面邊界和無窮遠處邊界條件，可得:

式中:ψ=φ/W為單位速度勢，Ma是附加質(zhì)量，根據(jù)定義為:

將輸運公式應(yīng)用于式(1)，經(jīng)過變化［14］:

考慮到物體所受的流體動壓力滿足:

式中未考慮靜水壓的影響。將式(5)、式(2)代入式(4)中，即可獲得含有附加質(zhì)量變化率的流體動壓力計算表達式為:

由上式可見，在考慮附加質(zhì)量變化率后，物體所受的流體動壓力不僅和加速度d W/d t有關(guān)，還和附加質(zhì)量變化率d Ma/d t有關(guān)。有趣的是，方程(6)貌似和著名的動量定律導(dǎo)出的公式)相矛盾。而實際上，對于動量定理中，F(xiàn)e是剛性物體受到的所有合外力，M是剛體本身的質(zhì)量。對于式(6)中F僅僅是物體受到的水動力，Ma為物體擾動的流體的質(zhì)量。換句話說，根據(jù)牛頓第三定律來看，流體本身除了受到來自物體的擾動力F，在其他邊界(如無窮遠邊界)還將受到其他的力。關(guān)于這點，文獻［14］已經(jīng)有了深入的討論，并通過數(shù)學演繹推導(dǎo)出的式(6)實質(zhì)上和勢流理論中著名的近場公式是一致的。

關(guān)于瞬態(tài)濕表面積物體附加質(zhì)量的求解，本文采用時域格林函數(shù)法，其對應(yīng)的具體相關(guān)理論見文獻［8］，主要的思想是引入時域格林函數(shù)求解如下定解條件:

式中:vn為物體運動的法向速度。自由液面的初始條件為:

速度勢的積分方程可表示為:

式中:p為固定點，q為動點，?rpq為點p，q的距離為點p，q關(guān)于自由面的鏡像點的距離，ωl(τ)表示τ時刻彈體瞬時表面與xoy平面的交線，N是平面曲線ωl(τ)的單位法向量，指向彈體內(nèi)部，VN是周線ωl(τ)沿法向N的法向運動速度為時域格林函數(shù)的興波項。通過邊界積分法求得速度勢φ后，代入式(3)即對應(yīng)時刻和對應(yīng)運動方向的附加質(zhì)量。邊界積分方程的求解方法見文獻［8］或文獻［15］。

2 考慮附加質(zhì)量變化率的結(jié)構(gòu)運動方程

本節(jié)在式(6)的基礎(chǔ)上，將附加質(zhì)量變化率考慮到結(jié)構(gòu)運動方程中，推導(dǎo)針對彈體模型的砰擊力計算公式。首先將式(6)中d Ma/d t項時空分離，改成如下形式:

式中:l表示彈體軸向浸濕長度，對于上文考慮的彈體垂直入水，有d l/d t=－W，d Ma/d l表示單位入水長度引起的附加質(zhì)量的變化。將d Ma/d t時空分離之后的好處在于d l/d t將直接與彈體運動速度有關(guān)，而附加質(zhì)量的空間變化率d Ma/d l則和彈體運動狀態(tài)(橫向剛體運動、軸向剛體運動、彈性振動)有關(guān)。

設(shè)結(jié)構(gòu)的運動方程可離散成多自由度系統(tǒng)的運動方程式，如下所示:

式中:M為結(jié)構(gòu)的質(zhì)量矩陣;C為結(jié)構(gòu)的內(nèi)阻尼矩陣;K為結(jié)構(gòu)的剛度矩陣;F為結(jié)構(gòu)運動引起的水動力合力;f0為除此之外的外界激振力。將式(13)代入式(14)得到:

設(shè)x=φq，φ為干模態(tài)振型，q為彈體真空中的主坐標矢量:

式中:M'e為主質(zhì)量矩陣，Me(t)為主附加質(zhì)量矩陣，Ce為主阻尼矩陣，Ke為主剛度矩陣。由上式可見，附加質(zhì)量不僅會影響質(zhì)量陣，還會影響結(jié)構(gòu)阻尼陣，對應(yīng)主阻尼陣增量為所以，彈體在出入水過程中，由于附加質(zhì)量的影響，其對應(yīng)的振動頻率、阻尼系數(shù)等會劇烈變化，對應(yīng)的結(jié)構(gòu)振動特性較之水中有很大變化。

3 算例分析

3.1 方法驗證

本文基于勢流理論，采用邊界元方法，開發(fā)三維計算程序，計算彈體各階濕模態(tài)的附加質(zhì)量隨入水長度變化規(guī)律。在探討規(guī)律之前，先通過簡單的彈體運動驗證本文程序的有效性。取圓柱體水中繞流模型作為剛體附加質(zhì)量的計算模型，參數(shù)見表1。

無限長直圓柱水中繞流模型單位長度附加質(zhì)量的理論解是πa2，a是圓柱體半徑。

從表1可以看出，用數(shù)值方法計算附加質(zhì)量存在一定的誤差，并且長度直徑比越大，計算所得的附加質(zhì)量越接近于理論值.這是因為選取的圓柱體模型為有限長，而圓柱體水中繞流模型附加質(zhì)量理論解是基于無限長直圓柱假設(shè)而得到的結(jié)論。因此，本文數(shù)值解有較高的精度，驗證了用數(shù)值方法計算剛體附加質(zhì)量的正確性和有效性。

表1 剛體運動附加質(zhì)量Tab.1 Added mass of the rigid body motion

3.2 彈體斜入水算例

選擇某一工況進行物體入水附加質(zhì)量變化研究。模型入水示意圖如圖2所示，其中O－XYZ為大地坐標系，o－xyz為彈體隨體坐標系。計算工況為:彈體全長L=1.0 m，初始俯仰角 θ=65°，浸濕長度 l選取為物體水下部分在垂向的投影長度，則最大的浸濕長度為0.9 m。初始彈體運動水平速度方向向右，大小為0.6 m/s，垂向速度方向向下，大小為16 m/s。

圖2 彈體斜入水運動示意圖Fig.2 The sketch map of oblique entry in the water of a body

圖3 彈體落水過程數(shù)值模擬Fig.3 The simulation of a missile entering into water

彈體入水運動的姿態(tài)變化示意圖如圖3所示。

從圖3中可見，在t=0.01 s時，彈體頭部基本完全入水，流場在物體擠壓下開始上涌;0.03 s時彈體左側(cè)流體下陷，右側(cè)流體上涌，這是因為彈體在向下運動的同時，還存在一向右的速度，使彈體擠壓其右側(cè)流場，偏離其左側(cè)流場;t=0.06 s時，彈體尾部開始入水，自由液面變形較大;t=0.08 s時，流體完全淹沒彈體，自由面產(chǎn)生大變形，此后彈體將在流場中自由運動，液面在較長時間后將趨于平靜。本文的數(shù)值模擬中沒有考慮彈體入水攜帶的肩空泡問題，故將在一定程度上與真實物理過程有所差異。

本文最主要關(guān)心彈體入水過程的中瞬態(tài)附加質(zhì)量Ma以及附加質(zhì)量變化率d Ma/d l問題。Ma及d Ma/d l本身和結(jié)構(gòu)的落水運動沒有關(guān)系，僅和濕表面形狀以及模型所在水下位置有關(guān)。圖4～圖8中給出圖3入水模型相關(guān)方向的Ma及d Ma/d l的變化曲線，各圖中橫坐標Z0為彈體頭部頂點在大地坐標系下的垂向坐標。

圖4 X方向附加質(zhì)量Mxx及其變化率d Mxx/d lFig.4 The added mass in X direction Mxx and its variation d Mxx/d l

從圖4可見，X方向附加質(zhì)量Mxx隨彈體入水深度增大而增大，當Z0超過0.9 m以后，附加質(zhì)量隨浸水長度變化微小。因為Z0小于0.9 m時，濕表面隨浸水長度的增加而增加，附加質(zhì)量也隨之增加。Z0超過0.9 m后，彈體完全入水濕表面不再改變，只不過模型所在位置的水深發(fā)生變化，附加質(zhì)量Mxx變化很小。

圖5 Y方向附加質(zhì)量Myy及其變化率d Myy/d lFig.5 The added mass in Y direction Myy and its variation d Myy/d l

Z0小于0.9 m時，附加質(zhì)量變化率d Mxx/d l隨浸水長度l的增加而增加。Z0超過0.9 m，附加質(zhì)量變化率d Mxx/d l隨入水長度的增加而迅速減小，并趨于0。這是因為Z0小于0.9 m，主要影響附加質(zhì)量的是濕表面積，Z0大于0.9 m，主要影響附加質(zhì)量的是模型所在位置的水深。濕表面對附加質(zhì)量的影響要比水深影響大得多。

從圖5可見，Y方向附加質(zhì)量 Myy及其變化率d Myy/d l的變化趨勢和X方向基本一致，只不過Y方向附加質(zhì)量比X方向附加質(zhì)量稍小，X方向附加質(zhì)量的最大值達到32.6kg，而X方向附加質(zhì)量的最大值只有28kg左右。這是因為彈體OYZ平面上存在傾角，因此在Y方向的投影面積比在X軸的投影面積小，致使Y方向附加質(zhì)量比X方向附加質(zhì)量略小。

圖6 Z方向附加質(zhì)量Mzz及其變化率d Mzz/d lFig.6 The added mass in Z direction Mzz and its variation d Mzz/d l

圖7 一階彈性振動附加質(zhì)量M11及其變化率d M11/d lFig.7 The added mass of the first vibration M11 and its variation d M11/d l

圖8 Y方向和一階彈性振動耦合附加質(zhì)量My1及其變化率d My1/d lFig.8 The coupling added mass My1 and its variation d My1/d l

從圖6可見，軸向附加質(zhì)量變化率d M/d l存在兩個明顯的峰值，峰值出現(xiàn)時刻分別對應(yīng)物體頭部和尾部入水?？梢姖癖砻娣eZ向投影的劇烈變化導(dǎo)致了附加質(zhì)量變化率的峰值存在。同時這兩種峰值在垂向沖擊力的變化中也可以觀察到。另外，在第一個和第二個峰值之間，附加質(zhì)量變化率dMzz/dl隨浸水長度的變化很小，這與Mxx，Myy的變化趨勢是不同的，這主要是因為結(jié)構(gòu)在Z方向投影面積的變化很小，Z0小于0.9 m時，Z方向附加質(zhì)量Mzz近于線性增加。Z0大于0.9 m后，Z方向附加質(zhì)量Mzz基本維持某一定值附近。

圖7為一階彈性振動的附加質(zhì)量及其變化率隨入水長度的變化曲線。從中可見，在物體尾部浸入水中時變化特別大。這是由結(jié)構(gòu)本身的一階振型決定的。本文研究的結(jié)構(gòu)對應(yīng)的一階振動在尾部振動最為劇烈，其它部分振動很小。致使在彈體尾部浸入水中時，附加質(zhì)量變化特別大。Z0大于0.9 m后，一階彈性振動附加質(zhì)量M11隨入水長度的變化同樣不明顯。

圖8為彈體一階彈性振動與Y方向剛體運動的耦合附加質(zhì)量變化趨勢。對于本文研究的彈體結(jié)構(gòu)，其一階彈性振動方向主要是沿Y方向，在X方向振動量很小，故其與X方向剛體運動耦合附加質(zhì)量非常小，此處不再給出。從圖8中可見，一階彈性振動與Y方向耦合附加質(zhì)量卻非常大。因此，彈體在做剛體運動的同時伴隨彈性振動，耦合附加質(zhì)量是不可忽略的。同時，耦合附加質(zhì)量呈現(xiàn)波動的趨勢，這和振型密切相關(guān)的。

3.3 彈體直入水算例

在研究了剛體運動、彈性振動誘發(fā)的附加質(zhì)量隨浸濕長度變化的趨勢后，本文接著研究彈體落水過程中垂向砰擊力的變化趨勢。選取的計算工況示意圖如圖1所示，其中彈體長1.0 m，邊界元模型與3.2節(jié)一致。垂直入水速度為10 m/s，與3.2節(jié)算例中自由落水有所不同的是，本節(jié)中強制彈體落水速度保持10 m/s勻速不變，對應(yīng)的彈體姿態(tài)也保持垂直不變。

這里采用兩種數(shù)值方法計算彈體遭受的砰擊力，一是通過公式(5)求解物面上每點壓力，然后直接積分求得合力(稱為“數(shù)值解”);二是應(yīng)用公式(6)，即先求得附加質(zhì)量以及變化率，然后直接獲得合力。由于附加質(zhì)量需要數(shù)值計算獲得，而一旦得到附加質(zhì)量變化率，則可應(yīng)用式(6)計算合力，故本文稱此方法為“半數(shù)值半解析解”。

圖9 “數(shù)值解”與“半數(shù)值半解析解”對比Fig.9 The compare between the numerical results and the semi-analytical ones

圖9中F3為計算得到的垂向力。從圖中可見無論是直接積分法還是通過附加質(zhì)量變化求得的砰擊力均存在兩個峰值，第一個峰值是彈體頭部觸水時，第二個峰值是整個彈體完全浸沒水中時。對比圖9和圖6可發(fā)現(xiàn)，垂向砰擊力的兩個峰值與軸向附加質(zhì)量變化率峰值出現(xiàn)的時刻相同，所以彈體在頭部、尾部入水瞬間對應(yīng)的砰擊力主要來源于附加質(zhì)量變化率一項，即源于濕表面積的劇烈變化。此外，在誤差允許范圍內(nèi)，公式(7)計算得到的力與公式(5)計算得到的力十分接近，再次表明應(yīng)用附加質(zhì)量及其變化率計算砰擊力的可行性。

4 結(jié)論

基于勢流理論，本文對彈體落水問題進行研究，基于能量法求得入水砰擊力與附加質(zhì)量及變化率的簡約關(guān)系。并在此關(guān)系基礎(chǔ)上研究了附加質(zhì)量變化對彈體結(jié)構(gòu)振動的影響，以及入水過程中剛體運動、彈性振動引起附加質(zhì)量的變化。結(jié)論如下:

(1)彈體入水過程中垂向砰擊力存在兩個峰值，分別對應(yīng)彈體頭、尾部接觸水面時刻，該砰擊力主要源于這兩個時刻濕表面積的劇烈變化;

(2)彈體入水過程中附加質(zhì)量的變化主要取決于彈體形狀和濕表面積的變化，水深對其影響相對較小;

(3)附加質(zhì)量的變化不僅影響結(jié)構(gòu)系統(tǒng)振動的質(zhì)量陣，還對阻尼系數(shù)陣產(chǎn)生的影響。即附加質(zhì)量不僅改變彈體振動頻率，也改變阻尼系數(shù)。

［1］蔣 莉，沈孟育.求解流體與結(jié)構(gòu)相互作用問題的ALE有限體積方法［J］.水動力學研究與進展，2000，15(2):148－155.JIANG Li，SHEN Meng-yu.ALE finite volume computations of fluid-structure interaction problems［J］.Journal of Hydrodynamics，2000，15(2):148 －155.

［2］王基盛，楊慶山.流體環(huán)境中結(jié)構(gòu)附加質(zhì)量的計算［J］.北方交通大學學報，2003，27(1):40 －43.WANG Ji-sheng，YANG Qing-shan.Calculation on added mass of structures in fluid environments［J］.Journal of Northern Jiaotong University，2003，27(1):40 －43.

［3］弓三偉，陸宏志，鄒正平，等.彈性體和帶空泡體的附加質(zhì)量動態(tài)數(shù)值計算［J］.計算機仿真，2010，27(4):348－354.GONG San-wei， LU Hong-zhi， ZOU Zheng-ping， et al.Hydrodynamic numerical computation of the added mass of elastic body ＆ cavitation bubble［J］.Computer Simulation，2010，27(4):348 －354.

［4］Newman J N，Sclavounos P D.The computation of wave loads on large offshore structures［C］. Trondheim:Tapir Publishers，1998，2:605 －622.

［5］Buchner B，Vogt A.The effect of bow flare angle on FPSO green water loading［C］.Proceedings of 19th OMAE.New Orleans，2000:1 －8.

［6］Kashiwagi M，Endo K.Yamaguchi H.Wave drift forces and moments on two ships arranged side by side in waves［J］.Ocean Engineering，2005，32:529－555.

［7］Kenneth K C.Underwater shock analysis of a missile launch tube［R］.Research and Technology Department，Naval Surface Weapons Center，1982 －05 －20.

［8］戴遺山，段文洋.船舶在波浪中運動的勢流理論［M］.北京:國防工業(yè)出版社，2008.

［9］李 明，尹云玉.水下航行體動態(tài)響應(yīng)計算的附加質(zhì)量探討［J］.導(dǎo)彈與航天運載技術(shù)，2008，4:16－18.LI Ming，YIN Yun-yu.Effect of the added mass on transverse response calculation of submarine vehicle in exiting-water process［J］.Missiles and Space Vehicles，2008，4:16 －18.

［10］賈 亮，榮克林.附加質(zhì)量對結(jié)構(gòu)出水載荷的影響［C］.全國結(jié)構(gòu)動力學學術(shù)研討會論文集，2007:37－144.

［11］王 科，張 犀，張志強，等.FPSO型采油平臺附加質(zhì)量與阻尼力的分析研究［J］.船舶力學，2011，15(2):40－47.WANG Ke，ZHANG Xi，ZHANG Zhi-qing，et al.Numerical analysis of added mass and damping of FPSO type platform，Journal of Ship Mechanics，2011，15(2):40 －47.

［12］袁振偉，褚福磊，王三保.橫向振動圓柱體在流體中的動力學特性［J］.工程力學，2006，23(8):49－52.YUAN Zhen-wei，CHU Fu-lei，WANG San-bao.Dynamic characteristics of a lateral vibrating cylinder immersed in unbounded fluids［J］.Engineering Mechanics，2006，23(8):49－52.

［13］Molin B，Cointe R，F(xiàn)ontaine E.On energy arguments applied to the slamming force［C］.11st International Workshop on Water Waves and Floating Bodies，Hamburg，1996.

［14］Wu G X.Hydrodynamic force on a rigid body during impact with liquid［J］.Journal of Fluids and Structures.1998，12(5):549－559.

［15］張阿漫，姚熊亮.近自由面水下爆炸氣泡的運動規(guī)律研究［J］.物理學報，2008，57(1):339－353.ZHANG A-man， YAO Xiong-liang. The law of the underwater explosion bubble motion near free surface［J］.Acta Physica Sinica，2008，57(1):339－353.