楊 衡，孫龍泉，龔小超，姚熊亮

（哈爾濱工程大學 船舶工程學院，哈爾濱 150001）

結構入水砰擊是典型的結構與流體瞬態(tài)流固耦合作用過程，在船舶工程、航空領域等中十分常見，如船艏砰擊、海上救生艇拋落、空投魚雷入水、水上飛機降落著水問題等。入水初期砰擊載荷以“瞬態(tài)”、“大幅”為特點，其砰擊壓力的脈寬往往在毫秒量級，瞬態(tài)沖擊載荷作用下往往使結構產生高頻的彈性振動，甚至對局部結構產生塑性破壞，因此，入水砰擊問題越來越引起人們的重視，同時其在軍事和民用上有著很現(xiàn)實的工程應用背景。結構自空中進入水面的過程中所承受的外部載荷大小外部載荷主要來自流體動壓力即流體沖擊壓力。當沖擊發(fā)生時，流體動壓力沿物體濕表面的分布特征及其大小取決于入水速度、入水角、結構形式、流體粘性、可壓縮性、結構彈性以及其它一些因素［1］，張阿漫等［2］通過實驗方法得到不同環(huán)境下條件氣泡脈動特性，得到剛性平面、45°斜面及舷側曲面條件下氣泡射流砰擊規(guī)律，對不同入水角和不同底部形式入水砰擊問題提供研究思路。由于結構物入水沖擊過程涉及結構運動非線性、自由表面運動非線性、氣穴效應、流固耦合效應、重力效應等，特別是伴隨自由面的抬升、飛濺，入水空泡現(xiàn)象等，在數(shù)學處理上有很大困難，張阿漫等［3］基于邊界積分開展近自由面水下爆炸氣泡作用下的運動規(guī)律進行研究，對自由面運動進行較好模擬。自Von Karman［4］提出用附加質量法計算物體入水砰擊載荷以來，結構的入水砰擊問題得到了廣泛而深入的研究，研究對象為簡單二維楔形剛性體，基于勢流理論進行求解。在彈性體入水砰擊研究中，流體與結構的耦合效應成為研究重點，盧熾華等［1］通過非線性伯努利方程，將流體域邊界元方程與結構有限元方程耦合求解，并考慮了自由表面的非線性邊界條件。陳鐵云等［5］采用混合歐拉－拉格朗日方法對高速船舶砰擊現(xiàn)象的沖擊載荷進行數(shù)值計算。顧慈祥等［6］計算了平頭旋轉殼體垂直入水時的水彈性問題。在國外對水彈性在入水砰擊問題上進行分析計算的有Sharov等［7－9］。在入水問題的非線性方面，Vorus［10］提出了一種改進的wagner模型，在其計算中，考慮了自由液面條件和伯努利條件的非線性，而物面條件是利用相當平板來模擬的。Korobkin［11］發(fā)展了Logvinovich提出的模型，原始Logvinovich方法（OLM）利用相當平板理論模擬物面條件，對自由面條件進行了線性化處理，而計算壓力的伯努利方程則是非線性的。Wu等［12］采用邊界元法討論了全非線性自由面條件下楔形體自由下落的模型。此外，SPH法，CIP法，VOF方法等數(shù)值方法都已經應用在入水沖擊問題的求解上。本文基于邊界元方法，結構邊界為流固耦合交界面，計算效率較SPH、VOF方法明顯提高，計算中考慮了結構運動非線性、流固耦合效應等，對結構入水問題進行了數(shù)值模擬；但基于邊界元方法計算結構入水砰擊問題時，無法精確模擬流體飛濺、自由液面大變形等問題。

本文的基于雙漸進法［13－14］（DAA），推導了考慮流體可壓縮性［15］的結構入水沖擊過程非線性雙漸進法（NDAA）流固耦合運動方程，研究分析彈性體結構入水沖擊過程的流固耦合效應。流體域速度勢采用邊界元方法計算，自由面速度勢滿足自由表面的非線性條件。結構為彈性體，運用有限元方法獲得結構的彈性響應，通過非線性伯努利方程將流體域速度勢和結構彈性響應耦合求解，通過與模型試驗值進行對比，驗證本文方法可靠性基礎上，討論彈性體與剛性體入水砰擊過程載荷特性差異，著重探討了入水速度對彈性結構入水砰擊載荷與彈性動力響應特性的影響。

1 結構入水NDAA數(shù)值計算方法

結構入水沖擊過程是一個流體和結構瞬態(tài)、強非線性的相互作用過程，為了將問題簡化，假設流場為各向同性、無粘、無旋但可壓縮的理想流體，計算中考慮結構的重力效應，計算模型如圖1所示。本文將結構視為彈塑性結構，在重力作用下以一定初速度入水，結構在外界水動力的作用下做流固耦合運動，同時結構的動力響應又反作用于流場，實現(xiàn)結構入水過程中流固耦合運動的數(shù)值模擬。

考慮可壓縮性的流體速度勢滿足如下控制方程：

假設可壓縮性對自由液面速度勢影響可忽略，計算模型需要滿足如下邊界條件［16］：

自由面上SF運動學邊界條件：

圖1 數(shù)值模型Fig．1 Numerical model

由面上SF動力學邊界條件：

在物面SB上：

在無窮遠邊界S∞邊界條件：

初始邊界條件：

式中，g為重力加速度；n是法向矢量坐標；V為結構運動速度，ω為彈性結構的振動速度。

DAA法就是從控制方程（1）的基礎上經推導得到。二階DAA基本方程［17］如下：

結構動力方程

式中：ps為流體中的散射壓力，Mf為流體質量矩陣，Ωf為流體頻率矩陣，Af為流體單元的面積矩陣，G為坐標轉換矩陣。

DAA方程在其理論推導的過程中，均引入了一個線性假設，即

二階DAA法是從延遲勢法基礎上得到，延遲勢的速度勢方程如下：

延遲勢的壓力方程

方程式（11）的求解可以得到流場散射壓力ps，流場總速度勢即可由得解。

對于結構入水砰擊問題，以“瞬態(tài)”，“大幅”為特點的非線性運動，其非線性效應對水動力值影響比較顯著，速度勢的空間導數(shù)就不能忽略了。本文在二階DAA法基礎上，將速度勢的空間導數(shù)引入到二階DAA法的數(shù)值計算，從而將非線性的影響考慮到結構入水砰擊的瞬態(tài)流固耦合計算當中。

在得到流場速度勢的基礎上，可以通過非線性伯努利方程［16］得到流場動壓力Pd。

式中：V為結構運動速度。

式（12）得到的流場動壓力值與二階DAA方程得到的流場動壓力值的本質區(qū)別在于式（12）中考慮了結構的速度效應項及結構的運動非線性項

求解式（12）得到考慮非線性效應的流場動壓力，將該考慮非線性效應的流場動壓力Pd代替結構運動方程中ps，結構運動方程變?yōu)椋?/p>

聯(lián)合求解式二階DAA方程、式（12）、（13）即為非線性雙漸進法（NDAA法）。

2 圓柱殼入水砰擊載荷特性分析

2.1 數(shù)值驗證

應用本文方法計算編制計算程序，計算結果與文獻［16］某縮比模型入水過程試驗數(shù)據(jù)進行對比分析。計算模型與試驗模型具有相同質量分布、慣性矩、材料特性，計算初始參數(shù)如下：入水角78°，入水速度Vx＝1．341 3 m／s、Vy＝0．835 8 m／s、Vz＝－10．068 9 m／s。下圖給出結構入水過程中三方向速度及砰擊壓力計算值與試驗值對比結果。

圖2 速度曲線對比Fig．2 The compare of the velocities between the numerical result and the experiment

圖3 砰擊載荷對比Fig．3 The compare of the slamming loads

圖2是縮比模型入水過程中三個方向的速度變化曲線，V1是橫向速度、V2是法向速度、V3是垂向速度，從圖中可以看出，圓柱入水過程中受到的垂向載荷變化劇烈，圓柱的垂向速度變化劇烈，相比于垂向速度和法向速度，橫向速度的變化曲線最為平緩。

圖3給出距離模型尾部4．6 m處測點入水砰擊載荷對比，從圖中可以看出，應用本文計算方法得到的砰擊載荷的計算值與試驗值峰值、脈寬相近，說明本文方法的可靠性。5．但本文計算模型由于未考慮結構自身阻尼特性，流場的飛濺、自由液面變形引起的水動力變化，僅考慮結構自身在水動力作用下的彈性振動，計算得到的壓力曲線與試驗壓力曲線在0．5 MPa、0．25 MPa附近的震蕩形式存在一定差別。

2.2 彈性圓柱體傾斜入水模擬

入水角度為30°，垂向速度為30 m／s，水平速度為0 m／s。取圓柱底部與側部迎水面特征節(jié)點進行分析，特征節(jié)點分布如圖4所示，其中節(jié)點70為入水點。

圖4 節(jié)點分布Fig．4 The distribution of notes

圖5 底部節(jié)點砰擊壓力對比Fig．5 The compare of the slamming loads on button

采用本文計算方法計算圓柱體傾斜入水過程各特征節(jié)點砰擊壓力分布及總體受力。圖5、6給出各節(jié)點砰擊壓力對比。

由圖可知，各節(jié)點砰擊壓力呈現(xiàn)相同的規(guī)律，即特征點接觸到流體時，會產生較大的砰擊峰值，砰擊壓力迅速下降，形成二次砰擊峰值，二次砰擊峰值約為第一次砰擊峰值的1／3。由于結構的彈性效應，入水初期，圓柱特征節(jié)點砰擊壓力呈現(xiàn)波動，隨入水深度的增加，結構表面的壓力主要為靜水壓力，此時，節(jié)點的壓力波動情況減弱，隨侵水深度的增加，壓力略有增加。

比較七個特征節(jié)點的壓力變化時歷曲線，由圖5可以看出，底部三個節(jié)點中，節(jié)點70壓力峰值最大，越靠近底部中心壓力峰值越小；而比較圖6中各節(jié)點壓力曲線，遠離初始入水點的特征節(jié)點的壓力峰值較大。

圖6 側部節(jié)點砰擊壓力對比Fig．6 The compare of the slamming loads on side

圖7 水平方向合力變化曲線Fig．7 The history curve of the horizontal force

圖8 垂直方向合力變化曲線Fig．8 The history curve of the vertical force

圖9 俯仰力矩變化曲線Fig．9 The history curve of the moment

圖10 俯仰角變化曲線Fig．10 The history curveof the pitch angle

圖11 俯仰角速度變化曲線Fig．11 The history curve of the pitch angle velocity

圖7、8給出圓柱入水過程水平方向和垂直方向合力時歷曲線。如圖6所示，水平及垂向合力存在雙峰值現(xiàn)象，分別對應圓柱體尾部入水和頭部入水兩個時間點，尾入水時刻水平合力較頭入水時刻略小，垂向合力較頭部入水時刻減小一半。

結構入水過程中，由于流體載荷的作用不均勻，會產生繞質心轉動的俯仰力矩，俯仰力矩的變化是結構入水姿態(tài)變化的主要因素。下圖給出俯仰力矩變化曲線。

由圖9可以看出，在圓柱頭尾入水時刻，俯仰力矩也會出現(xiàn)兩個峰值，與圓柱體所受到的合外力的變化相關。下圖給出在俯仰力矩作用下，圓柱入水過程俯仰角變化及圓柱質心位置的運動軌跡。

由俯仰角曲線可以看出，在整個入水過程中俯仰角逐漸減小，在圓柱質心入水前，俯仰角速度先反向增大，質心入水后，俯仰角速度減小，俯仰角速度的變化趨勢與俯仰力矩的變化曲線相同。

圖12 圓柱質心運動軌跡Fig．12 The motion track of the center of mass

圖12給出圓柱體入水過程質心位置變化，由圖可以看出，在入水初期，圓柱質心先向水平方向正向移動，隨入水深度的增加，質心向水平方向負向移動，質心運動曲線變化與圓柱水平方向合力時歷曲線變化密切相關，質心運動變化率隨水平方向合力變化而變化。

3 艦船入水砰擊載荷特性分析

在船舶與海洋工程領域，砰擊現(xiàn)象廣泛存在，如救生艇拋落、船舶從船臺下水、高速滑行艇艇底砰擊、甲板上浪等，砰擊的存在使船體局部產生應力集中區(qū)，船體產生局部變形甚至破壞。船舶在海上航行時，砰擊可能發(fā)生在船艏、船底、船艉等位置，尤其以底部砰擊最為重要，船底入水過程中受到巨大的砰擊載荷，會改變船體的垂向加速度，船底局部產生高頻振動，第一次世界大戰(zhàn)期間，發(fā)生了多起因船底砰擊而損壞的事故，船底砰擊越來越受到人們重視。

本節(jié)在前文基礎上，將雙漸近法應用于大型水面艦船的砰擊載荷計算中，相對于前文的回轉體，在計算大型水面艦船的入水砰擊流固耦合特性時，計算域變大，結構外形更加復雜，在砰擊載荷作用下彈性變形更顯著。同時，根據(jù)圖13船體入水砰擊試驗，船體入水過程中無空泡產生，且本文在計算中忽略了自由液面飛濺及變形的影響，但考慮了自由液面效應、流場可壓縮性、船體運動的非線性及伯努利方程的非線性。

圖13 船體入水砰擊Fig．13 Slamming Character on Water Entry of Hull

本節(jié)采用雙漸近法對某型艦船入水砰擊流固耦合運動進行數(shù)值模擬，研究船體垂直與傾斜入水砰擊過程中載特性及船體變形。某型艦船計算模型如圖13所示，船長為38．6 m，型深為 3．78 m。

圖14 船體型線Fig．14 The shiplines

3.1 艦船垂直入水砰擊

本節(jié)討論該水面艦船以初速度（30 m／s）入水，除受到因船底砰擊而產生的水動力作用外，不受其他任何外力的作用，研究該艦船的入水砰擊過程中受力狀態(tài)（流場動壓力）、Mises應力變化及船底變形。

在該工況計算中，靜止的流場在船底入水瞬間受到劇烈擾動，產生較大的砰擊載荷作用于船體表面，由于船底型線不均勻性，首先接觸水面的局部產生較大的砰擊載荷，隨入水深度的增加，應力逐步向船尾船底。下圖給出船體垂直入水過程中側部和底部Mises應力變化云圖。

圖15船體入水初期Mises應力變化云圖，距船艏部0．177 m處首先觸水，Mises應力逐漸增大，在船底，Mises應力向船艉逐漸傳遞，在船體側部，Mises應力向船甲板逐漸傳遞。由于船體對稱性，Mises應力分布呈對稱分布。在整個入水過程中，船體表面產生微幅振動。

取船底三個節(jié)點進行受力分析，如圖16所示，節(jié)點26為首先觸水點。下圖給出三個節(jié)點入水砰擊過程中垂向加速度和砰擊壓力對比。

圖15 船體Mises應力變化云圖Fig．15 The nephogram of Mises stress changes

圖16 特征節(jié)點分布Fig．16 The distribution of notes

由圖17、18可以看出，在船底砰擊過程中，越靠近船艉處，節(jié)點的垂向加速度及砰擊載荷越大。當船體開始入水后，由于砰擊載荷的作用范圍較小，船體整體仍處于加速運動狀態(tài)，砰擊載荷增加。

圖19，20給出船底三個節(jié)點彈性振動位移曲線及船體變形圖。

由圖19、20可以看出，船底有向內凹的趨勢。由于船底左右底板相交的區(qū)域為尖狀設計，如節(jié)點26處，并未產生明顯的變形，而船底型線之間的部分及船側部，如節(jié)點15，變性較大，局部振動明顯，而尾部節(jié)點6，由于船艉板的影響，其彈性變形要小于船中部。因此，在船底局部設計過程中應對應力集中區(qū)及船底中部進行局部加強設計。

圖17 垂向加速度對比Fig．17 The compare of the vertical accalarations

圖18 砰擊壓力對比Fig．18 The compare of the slamming loads

圖19 彈性振動位移對比Fig．19 The compare of the elastic vibration displacement

圖20 船底變形前后對比（1∶2）Fig．20 The compare of the hull deformation（1∶2）

3.2 艦船傾斜入水砰擊

本節(jié)討論該水面艦船橫向傾斜10°入水砰擊，以初速度（垂向速度20 m／s、法向速度5 m／s）入水，下圖給出船體傾斜入水過程中側部和底部Mises應力變化云圖。

由圖21船體入水初期Mises應力變化云圖可以看出，船體尾部首先觸水，Mises應力逐漸增大，在船底，Mises應力呈扇形逐漸傳遞，在船體側部，Mises應力向船甲板逐漸傳遞。

圖22 特征節(jié)點布置Fig．22 The distribution of notes

取船底三個節(jié)點進行受力分析，如圖22所示。圖23，24給出該船傾斜入水過程中三個節(jié)點入水砰擊過程中垂向加速度和砰擊壓力對比。

圖23 垂向加速度對比Fig．23 The compare of the vertical accalarations

圖24 砰擊壓力對比Fig．24 The compare of the slamming loads

圖25 船體變形前后對比（1∶2）Fig．25 The compare of the hull deformation（1∶2）

由圖23，24可以看出，在船底砰擊過程中，先觸水的節(jié)點的垂向加速度及砰擊載荷較后觸水的小。

從船底變形圖中可以看出，船底傾斜入水時，船底及側部變形較垂直入水砰擊時產生較大不同，船體先觸水一側的船底與側部產生較大變形。

4 結 論

本文基于波動理論，采用雙漸近方法，建立了結構入水砰擊流固耦合數(shù)值模型，計算中考慮了自由液面非線性條件、運動非線性及重力效應，對彈性圓柱體以及某型艦船入水砰擊過程（以某一初速度運動、瞬態(tài)砰擊載荷作用下的受迫振動以及后期自由振動）進行流固耦合數(shù)值模擬，實現(xiàn)對本文所開發(fā)程序的工程化應用，結果表明：

（1）結構入水過程中受到明顯的砰擊載荷，對于彈性體入水，砰擊載荷作用后，由于結構自身的彈性，結構底部產生明顯的彈性振動，砰擊壓力逐漸衰減。

（2）圓柱體傾斜入水時，剛彈耦合效應明顯，結構入水過程中水平、垂向合力在頭部入水和尾部入水時刻出現(xiàn)兩個峰值。

（3）船體垂直入水過程中，由于船底結構形式的影響，Misses應力分布沿船底中線對稱分布，船底在砰擊載荷的作用下出現(xiàn)較大變形。船底傾斜入水時，舷側也出現(xiàn)較大變形，可見，控制入水砰擊的角度，對船舶入水砰擊安全性至關重要。

［1］盧熾華，何友聲．二維彈性結構入水沖擊過程中的流固耦合效應［J］．力學學報，2010（3）：129－140．LU Chi-hua，HE You-sheng．Coupled analysis of nonlinear interaction between fluid and structure during impact［J］．Acta Mechanica Sinica，2010（3）：129－140．

［2］張阿漫，王詩平，白兆宏，等．不同環(huán)境下氣泡脈動特性實驗研究［J］．力學學報，2011，43（1）：71－83．ZHANG A-man，WANG Shi-ping，BAI Zhao-hong，et al．Experimental study on bubble pulse features under different circumstances［J］．Chinese Journal of Theoretical and Applied Mechanics，2011，43（1）：71－83．

［3］張阿漫，姚熊亮．近自由面水下爆炸氣泡的運動規(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．

［4］Von Karman T．The impact of seaplane floats during landing［R］．Washingt on DC，USA：National Advisory Committee for Aeronautics，NACA Technical Notes，1929，321．

［5］陳鐵云，王剛．高速船舶結構設計中流體沖擊載荷的數(shù)值計算［J］．上海交通大學學報，1998，32（11）：30－33．CHEN Tie-yun，WANG Gang． Numerical calculation of hydrodynamic impulsive loading in structural design of high speed ship［J］．Journal of Shanghai Jiaotong University，1998，32（11）：30－33．

［6］顧懋祥，程貫一，張效慈．平頭旋轉殼撞水水彈性效應的研究［J］．水動力學研究與進展，1991，6（1）：41－51．GU Mao-xiang， CHENG Guan-yi， ZHANG Xiao-ci． The study of hydroelastic effect on the water-entry of a revolutionary shell with a flat nose． ［J］． Journal of Hydrodynamics，1991，6（1）：41－51．

［7］Sharov V F．Ship bottom impact upon wave ［J］．Sudostroyeniye，1958，4：5－9．

［8］Meyerhoff E K．Die berechnung hzdroelastis cherstobe［J］．Schiffstechnik，1965，12：19－30．

［9］Zhao R，F(xiàn)altinsen O．Water entry of two-dimensional bodies［J］．Fluid Mech，1993，（246）：593－612．

［10］Vorus W．A flat cylinder theory for vessel impact and steady planin resistance［J］．Journal of Ship Research，1996，40（2）：89－106．

［11］Korobkin A A，Scolan Y M．Three-dimensional theory of water impact：Part2．Linearized Wagner problem ［J］．Journal of Fluid Mechanics，2006，549：343－373．

［12］Wu G X，Sun H，He Y S．Numerical simulation and experimental study of water entry of a wedge in free fall motion［J］．Journal of Fluids and Strucures，2004，19：277－289．

［13］Geers T L，Residual potential and approximation methods for three dimensional fluid-structure interaction problems［J］．Acoust．Soc．Amer，1971，49：1505－1510．

［14］ Geers T L，Doubly asympotic approximation for transient motions of submerged structures［J］．Acoust．Soc．Ame，1978，64：1500－1508．

［15］王詩平，孫士麗，張阿漫，等．可壓縮流場中氣泡脈動數(shù)值模擬［J］．力學學報，2012，44（3）：513－519．WANG Shi-ping，SUN Shi-li， ZHANG A-man，et al．Numerical simulation of bubble dynamics in compressible fluid［J］．Acta Mechanica Sinica，2012，44（3）：513－519．

［16］孫士麗．具有自由液面效應的高頻流固耦合問題［D］．哈爾濱：哈爾濱工程大學，2010．