船首外飄砰擊載荷的數值預報

于鵬垚，盧 雷，甄春博，王天霖

(1.大連海事大學 交通運輸裝備與海洋工程學院, 遼寧 大連 116026；2.中船重工船舶設計研究中心有限公司, 遼寧 大連 116001)

0 引 言

船舶在高海況下航行時，船體首部的外飄區域往往承受劇烈的砰擊載荷。準確地預報外飄砰擊載荷是合理地設計船首抗砰擊結構的前提。Von Karman[1]采用動量守恒方法計算了楔形體進入靜水面時所受的砰擊載荷。Wagner[2]在Von Karman研究的基礎上，考慮了楔形體入水時引起的波面升高，并給出了楔形體表面的壓力分布公式。Dobrovol'skaya[3]基于解析公式給出楔形體恒速入水時的自相似解。Zhao，Sun等[4–6]采用非線性邊界元法計算了二維剖面的入水砰擊載荷，并進一步分析了發生流動分離后剖面砰擊載荷的變化。Sun[6]采用非線性邊界元法和模態疊加法對圓殼入水過程中的砰擊載荷與結構響應進行預報，結果與試驗值吻合較好。陳震等[7–8]采用MSC.Dytran軟件仿真了楔形體和平底結構的入水砰擊載荷。 Aquelet等[9]討論了在采用任意歐拉拉格朗日算法預報砰擊載荷時，接觸阻尼對數值噪聲的影響。Stenius等[10]通過改變參數對楔形體入水問題進行了大量計算，分析了網格密度、接觸剛度、接觸阻尼、時間步長對局部砰擊壓力載荷的影響。張曉波[11]，駱寒冰[12]，Wang[13]等針對小尺度的楔形體模型或外飄模型在自由落體時所受的砰擊載荷進行預報，并將預報值與試驗值進行比較。

船體外飄剖面隨船體運動時，剖面底部入水后流體與剖面可能發生分離，分離后的液面再次砰擊到外飄區域，這一過程中會捕捉一定的氣體從而導致二次砰擊問題更為復雜。過往的研究普遍針對小尺度的楔形體模型和外飄模型，而且模型的運動狀態為自由落體運動或勻速運動。由于模型型線和模型運動速度的影響，很難反映出實船外飄剖面的二次砰擊現象。本文利用耐波性理論確定的船波相對運動規律，對實船尺度的外飄剖面砰擊載荷進行預報，分析二次砰擊現象、網格密度和接觸剛度參數對砰擊載荷的影響。

1 基本理論

1.1 控制方程

在任意歐拉拉格朗日算法（ALE）中，除拉格朗日坐標和歐拉坐標之外，引入參考坐標系進行控制方程的描述。若物質的速度為， 網格速度為，則相對速度。所以在ALE方法中，可以得到質量守恒、動量守恒、能量守恒的表達式分別為：

1.2 材料模型與狀態方程

船體結構砰擊現象中涉及空氣、水和船體3種材料。其中，本文對于空氣和水采用NULL材料模型進行模擬，采用剛體模型對船體的外飄型線進行模擬。對于流體材料需要結合相應的流體狀態方程來實現控制方程的求解。狀態方程描述了體積變形與壓力的關系。本文采用線性多項式狀態方程來表達空氣，采用Gruneisen方程來表達水，具體的參數取值參照文獻[11]。

1.3 耦合算法

采用罰函數耦合方法來模擬砰擊現象中結構與流體的耦合現象。罰函數耦合方法類似于一個彈簧系統，彈簧的一端與結構網格（主物質）的節點相連，彈簧的另一端與流體網格（從物質）的節點相連。當流體節點浸入結構網格表面時，在界面處產生耦合力會分布到歐拉流體的節點上來保證結構表面的物面條件。耦合力與結構網格進入流體網格中的深度成正比，可以表達為：

式中：P為耦合力；k為接觸剛度；d為浸入深度。

與文獻[11–13]不同，本文給出耦合力與浸入深度的變化關系。當浸入深度d=0時，P=0；當d=dmax時，P=pmax。dmax和pmax分別為所允許的最大浸入深度和砰擊時達到的最大砰擊壓力。本文后續均采用上述2點給出耦合力與浸入深度的線性關系，此時接觸剛度相當于pmax與dmax的比值。

2 數值模型

2.1 網格域

建立如圖1所示的砰擊載荷數值預報模型，圖中流體域的尺度代表各自之間的長度比例，可根據結構剖面的大小確定流體域的大小。區域1到區域3為水域，區域4–區域7為空氣域，曲線為外飄剖面結構。區域1、區域4為結構與流體的主要作用區域，采用細網格進行劃分。區域2–區域3，5–區域7采用尺度漸變的網格進行過渡。XOY面內結構網格的長度是區域1網格在x方向上網格長度的1.5倍。模型在z方向（垂直于紙面向外）的尺度為1個網格的長度。網格密度的示意如圖2所示。

圖1 數值模型Fig.1 Numerical model

2.2 邊界條件

圖2 網格密度示意圖Fig.2 Schematic diagram of mesh density

在add’a’面施加對稱條件，采用一半模型進行船體外飄剖面砰擊入水的模擬。在abb’a’，bcc’b’，cdd’c’平面添加無反射邊界條件，來模擬無限尺度的流域。約束abcd，a’b’c’d’平面內節點z方向的線位移和繞x，y方向的角位移，模擬二維砰擊現象。通過添加負y方向的加速度，實現重力加速度。給定剖面初速度或強制位移，實現不同運動規律下砰擊載荷的預報。

3 數值模擬中的主要影響因素

3.1 網格密度

網格密度是影響數值仿真計算的重要參數，網格尺寸過大會導致仿真結果不真實，不能很好地表達結構與流體之間的耦合作用；網格尺寸過小會影響計算效率，導致工程應用性差。因此網格密度的選取至關重要。

3.2 接觸剛度

接觸剛度（PFAC）是用來描述耦合系統之間剛度的參數，它用來計算流體和結構之間的耦合力。接觸剛度的大小會影響穿透現象的發生，如果接觸剛度太大，會使物質產生鋼化現象，物質表面過分剛硬會影響程序計算的穩定。

3.3 二次砰擊

由于船體首部各剖面形狀均不同，在曲率變化大的剖面可能會有二次砰擊現象發生，而且發生二次砰擊的剖面計算所得的砰擊壓力峰值要遠大于首次砰擊產生的砰擊壓力峰值，因此在船首一些特殊剖面必須考慮二次砰擊的現象并且對二次砰擊的砰擊壓力峰值進行計算，才能夠完整的模擬整個砰擊過程。

4 算例分析

對2種外飄型線，進行船舶與波浪相對運動規律下的砰擊載荷的預報，討論了網格密度、接觸剛度和剖面形式對砰擊載荷的影響，并將數值預報得到的砰擊壓力系數與規范值進行對比。型線1與型線2僅在P1以下形狀不同，外飄區域的形狀完全相同。依照剖面的尺度確定仿真模型區域1的x，y，z方向尺度分別為22 m，25 m，0.5 m，其他區域按照圖3給定的比例進行確定。接觸阻尼取為0.9倍的臨界阻尼。

圖3 實船尺度剖面Fig.3 Section of an actual ship

4.1 數值方法驗證

針對落體砰擊試驗中的外飄砰擊力和砰擊壓力載荷進行預報，由試驗模型尺度確定仿真模型區域1的x，y，z方向尺度分別為0.22 m，0.25 m和0.005 m，其他區域按照圖4給定的比例進行確定。其中，細網格區域1、區域4的網格在x、y方向的尺度分別為0.001 1 m和0.001 25 m。接觸剛度相關的參數dmax=0.000 3 m和pmax=17 000 Pa。接觸阻尼取為0.9倍臨界阻尼。通過與公開發表的文獻[6]進行對比，如圖5所示。可以看出理論值與試驗值吻合較好，驗證本文采用的數值模型和計算參數的正確性。

圖4 落體模型剖面Fig.4 Section of the free-drop model

4.2 不同網格密度仿真模型計算結果

圖5 理論值與試驗值對比Fig.5 Comparison between numerical results and experimental results

圖6 船舶與波浪相對運動Fig.6 Relative motion between the ship and the wave

在與接觸剛度相關的參數dmax=0.0 4 m和pmax=2.40×106 Pa時，改變網格的密度，針對2種型線進行船舶相對運動關系下的砰擊載荷預報。不同網格密度下對應細網格區域在x，y方向的尺度dx，dy如表1所示。砰擊力與砰擊壓力的對比如圖7和圖8所示，為了節省空間僅給出P1位置的砰擊壓力隨時間變化曲線。砰擊力與砰擊壓力的峰值如表2所示。

表1 細網格區域的網格尺度Tab.1 Mesh size in fine grid region

可以看出，對于型線1和型線2，不同網格下的砰擊力和砰擊壓力變化趨勢一致；當網格尺度從mesh1到mesh3逐漸減小時，型線1的砰擊力和砰擊壓力的峰值有增大的趨勢，而型線2砰擊力和砰擊壓力的峰值變化不大。由圖9可以看出型線1出現了液面與剖面分離后二次砰擊現象，分析正是由于此種現象，導致不同網格下，型線1的砰擊力與砰擊壓力的峰值更難達到收斂。

圖7 不同網格密度下砰擊力Fig.7 Slamming forces with different mesh densities

圖8 不同網格密度下P1處砰擊壓力Fig.8 Slamming pressure of P1 with different mesh densities

4.3 不同接觸剛度仿真模型計算結果

在網格密度為mesh3時，改變與接觸剛度相關的參數dmax和pmax，針對2種型線進行船舶相對運動關系下的砰擊載荷的預報，參數dmax和pmax如表3所示。砰擊力與砰擊壓力的對比如圖10和圖11所示。砰擊力與砰擊壓力的峰值如表4所示。

表2 不同網格密度下砰擊載荷峰值Tab.2 Peak values of slamming loads with different mesh densities

圖9 t=0.8 s時流場液面Fig.9 Water elevation of the fluid at t=0.8 s

表3 接觸剛度參數Tab.3 Parameters of the contact stiffness

可以看出，對于型線1和型線2，不同接觸剛度下的砰擊力和砰擊壓力變化趨勢一致；對于型線2，當接觸剛度從k1到k3逐漸增大時，砰擊力和砰擊壓力的峰值變化不大，已經收斂，但在載荷下降時間段內（1.2 s以后），數值噪聲有增大的趨勢；由型線2的砰擊峰值和接觸剛度收斂范圍，推斷當接觸剛度從k3到k5時，型線1砰擊力和砰擊壓力峰值也收斂。但計算結果表明，當接觸剛度從k3到k5時，砰擊力和P1，P2砰擊壓力峰值有增大的趨勢, 分析原因為二次砰擊后，砰擊力與砰擊壓力的峰值受接觸剛度的影響更為敏感。文獻[14]在平底結構入水問題中說明持續時間較短的砰擊壓力峰值的差別并不會導致應力響應的差別，在理論預報的砰擊壓力峰值與試驗值存在差別的情況下，結構應變與試驗依然吻合較好，因而后續可從結構應變響應的角度對型線1的接觸剛度的收斂性進行研究。

圖10 不同接觸剛度下砰擊力Fig.10 Slamming pressures with different contact stiffness

圖11 不同接觸剛度下的砰擊壓力Fig.11 Slamming forces with different contact stiffness

4.4 不同剖面形狀計算結果

在接觸剛度為k4、網格密度為mesh3時，對2種不同剖面砰擊載荷的對比如圖12和圖13所示。可以看出，2種剖面形式在時間從1.0 s到1.2 s之間，由于二次砰擊的影響，導致型線1的砰擊力和P1，P2處砰擊壓力比型線2多一個很大的峰值。在1.2 s之后，2種型線的砰擊載荷相差不大。

表4 不同接觸剛度下砰擊載荷峰值Tab.4 Peak values of slamming loads with different contact stiffness

圖12 不同型線下的砰擊力Fig.12 Slamming forces of different sections

圖13 不同型線下的砰擊壓力Fig.13 Slamming pressures of different sections

5 砰擊壓力系數與規范值的比較

在接觸剛度為k2、網格密度為mesh3時，利用2種不同剖面確定砰擊壓力和砰擊速度換算得到砰擊壓力系數，如表5所示。P1恰好位于曲率變化比較大的位置，考慮到剖面入水后的液面升高，單純依據P1局部位置的底升角確定砰擊壓力是不合理的，這里選取P1位置上部曲率變化平穩位置確定P1的底升角。由圖14可以看出，型線2各點對應的砰擊壓力系數恰好位于2種規范曲線[15]之間，而型線1的P1，P2對應的砰擊壓力系數則明顯大于規范值，可見目前規范曲線尚未考慮到二次砰擊引起砰擊壓力峰值的增大。

表5 砰擊壓力系數的數值預報結果Tab.5 Numerical results of the slamming pressure coefficient

圖14 砰擊壓力系數理論值與規范值Fig.14 Numerical and rule results of the slamming pressure coefficient

6 結 語

本文通過對船首外飄砰擊載荷的預報得到如下結論：

1）通過對比數值預報結果與落體砰擊模型試驗結果，驗證本文在網格劃分和接觸剛度等參數選取的合理性。

2）通過研究不同網格密度、接觸剛度下2種型線的砰擊力與砰擊壓力載荷的數值預報，可以看出對于型線2，在網格密度相對稀疏，dmax約為模型寬度的0.002，pmax約為剖面可達到的最大砰擊壓力時，即可得到穩定的砰擊載荷數值結果；由于二次砰擊現象，導致型線1砰擊力和砰擊壓力載荷的峰值受網格密度和接觸剛度的參數影響更加敏感，但不同參數下，砰擊力與砰擊壓力載荷的變化趨勢一致，后續可從結構響應的角度對二次砰擊現象中的網格密度和接觸剛度的收斂性進行研究。

3）在相同網格密度和接觸剛度參數下，型線1在二次砰擊發生時（1.0～1.2 s），砰擊力和砰擊壓力P1、P2明顯要大于型線2的值，1.2 s之后2種型線砰擊載荷變化幾乎一致。可見，二次砰擊現象會導致局部位置砰擊載荷的增加，應該引起結構設計的關注。

4）由砰擊壓力系數與規范值的比較看出，對于未發生二次砰擊（型線2的P1，P2，P3）和不受二次砰擊影響的位置（型線1的P3），砰擊壓力系數與規范值相差不大。但對于受二次砰擊壓力較大的位置（型線1的P1，P2），砰擊壓力系數明顯大于規范值，可見當前規范仍需更合理地考慮二次砰擊現象對砰擊載荷的影響。

畢淑敏在《提醒幸福》里說過：我們從小就習慣了在提醒中過日子，提醒注意跌倒、提醒注意路滑、提醒受騙上當、提醒榮辱不驚……先哲們提醒了我們一萬零一次，卻從不提醒我們幸福。現實生活中，我們的物質生活越來越富裕，社會福利越來越好，可是幸福卻常常被人遺忘。回顧這一年來自己所經歷的一切，真實地感受到幸福也需要被提醒。

[1]VON KARMAN T.The impact on seaplane floats during landing[R].NACA

參考文獻：

TN321, oct.1929.

[2]WAGNER H.Uber stass-und gleitvorgange und der oberflache von flussigkeiten[J].ZAMM, 1932, 12(4): 193–215.

[3]DOBROVOL’SKAYA Z N.On some problems of similarity flow of fluids with a free surface[J].Journal of Fluid Mechanics, 1969, 36: 805–829.

[4]ZHAO R, FALTINSEN O M.Water entry of two-dimensional bodies[J].Journal of Fluid Mechanics, 1993, 246: 593–612.

[5]ZHAO R, FALTINSEN O M, AARSNES J V.Water entry of arbitrary two-dimensional sections with and without flow separation [C]// Proc.21st Symposium on Naval Hydrodynamics, 1996: 408–423.

[6]SUN Hui.A boundary element method applied to strongly nonlinear wave-body interaction problems [D].Norway:Norwegian University of Science and Technology, 2007:15–100.

[7]陳震, 肖熙.二維楔形體入水砰擊仿真研究[J].上海交通大學學報, 2007, 41(9): 1425–1428.

[8]CHEN Zhen, XIAO Xi.The simulation study on water entry of 2D wedge bodies[J].Journal of Shanghai Jiaotong University,2007, 41(9): 1425–1428.

[9]陳震, 肖熙.平底結構砰擊壓力峰值分析[J].上海交通大學學報, 2006, 40(6): 983–987.

[10]CHEN Zhen, XIAO Xi.Analysis about the slamming pressure peak value on a flat bottom structure[J].Journal of Shanghai Jiaotong University, 2006, 40(6): 983–987.

[11]AQUELET N, SOULI M, OLOVSSON L.Euler-lagrange coupling with damping effects: application to slamming problems[J].Computer Methods in Applied Mechanics and Engineering, 2006, 195: 110–132.

[12]STENIUS I, ROSN A, KUTTENKEULER J.Explicit FE-modeling of fluid-structure interaction in hull-water impacts[J].Int Shipbuild Prog, 2006, 53, 1031–1121.

[13]張曉波.船底結構砰擊時的流固耦合數值模擬[D].大連: 大連理工大學, 2007.

[14]駱寒冰, 吳景健, 王珊, 等.基于顯式有限元方法的二維楔形剛體入水砰擊載荷并行計算預報[J].船舶力學, 2012, 16(8):907–913.

[15]LUO Hang-bing, WU Jing-jian, WANG Shan, et al.Parallel computing simulation of water entry of a 2D rigid wedge using an explicit finite element method [J].Journal of Ship Mechanics, 2012, 16(8): 907–913.

[16]WANG S, GUEDES S C.Slam induced loads on bow-flared sections with various roll angles[J].Ocean Eng, 2013, 67:45–57.

[17]FALTINSEN O M.Hydrodynamics of high-speed marine vehicles[M].New York: Cambridge University Press, 2005.

[18]王輝.船體結構局部強度設計中的砰擊載荷確定方法[J].中國造船, 2010, 51(2): 68–77.

[19]WANG Hui.Slamming load determination in local structure design of ships[J].Shipbuilding of China, 2010, 51(2): 68–77.