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

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

0 引 言

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

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

1 基本理論

1.1 控制方程

在任意歐拉拉格朗日算法（ALE）中，除拉格朗日坐標和歐拉坐標之外，引入?yún)⒖甲鴺讼颠M行控制方程的描述。若物質(zhì)的速度為， 網(wǎng)格速度為，則相對速度。所以在ALE方法中，可以得到質(zhì)量守恒、動量守恒、能量守恒的表達式分別為：

1.2 材料模型與狀態(tài)方程

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

1.3 耦合算法

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

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

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

2 數(shù)值模型

2.1 網(wǎng)格域

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

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

2.2 邊界條件

圖2 網(wǎng)格密度示意圖Fig.2 Schematic diagram of mesh density

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

3 數(shù)值模擬中的主要影響因素

3.1 網(wǎng)格密度

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

3.2 接觸剛度

接觸剛度（PFAC）是用來描述耦合系統(tǒng)之間剛度的參數(shù)，它用來計算流體和結(jié)構(gòu)之間的耦合力。接觸剛度的大小會影響穿透現(xiàn)象的發(fā)生，如果接觸剛度太大，會使物質(zhì)產(chǎn)生鋼化現(xiàn)象，物質(zhì)表面過分剛硬會影響程序計算的穩(wěn)定。

3.3 二次砰擊

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

4 算例分析

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

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

4.1 數(shù)值方法驗證

針對落體砰擊試驗中的外飄砰擊力和砰擊壓力載荷進行預報，由試驗模型尺度確定仿真模型區(qū)域1的x，y，z方向尺度分別為0.22 m，0.25 m和0.005 m，其他區(qū)域按照圖4給定的比例進行確定。其中，細網(wǎng)格區(qū)域1、區(qū)域4的網(wǎng)格在x、y方向的尺度分別為0.001 1 m和0.001 25 m。接觸剛度相關(guān)的參數(shù)dmax=0.000 3 m和pmax=17 000 Pa。接觸阻尼取為0.9倍臨界阻尼。通過與公開發(fā)表的文獻[6]進行對比，如圖5所示?？梢钥闯隼碚撝蹬c試驗值吻合較好，驗證本文采用的數(shù)值模型和計算參數(shù)的正確性。

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

4.2 不同網(wǎng)格密度仿真模型計算結(jié)果

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

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

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

表1 細網(wǎng)格區(qū)域的網(wǎng)格尺度Tab.1 Mesh size in fine grid region

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

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

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

4.3 不同接觸剛度仿真模型計算結(jié)果

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

表2 不同網(wǎng)格密度下砰擊載荷峰值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 接觸剛度參數(shù)Tab.3 Parameters of the contact stiffness

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

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

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

4.4 不同剖面形狀計算結(jié)果

在接觸剛度為k4、網(wǎng)格密度為mesh3時，對2種不同剖面砰擊載荷的對比如圖12和圖13所示?？梢钥闯觯?種剖面形式在時間從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 砰擊壓力系數(shù)與規(guī)范值的比較

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

表5 砰擊壓力系數(shù)的數(shù)值預報結(jié)果Tab.5 Numerical results of the slamming pressure coefficient

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

6 結(jié) 語

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

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

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

3）在相同網(wǎng)格密度和接觸剛度參數(shù)下，型線1在二次砰擊發(fā)生時（1.0～1.2 s），砰擊力和砰擊壓力P1、P2明顯要大于型線2的值，1.2 s之后2種型線砰擊載荷變化幾乎一致?？梢?，二次砰擊現(xiàn)象會導致局部位置砰擊載荷的增加，應該引起結(jié)構(gòu)設計的關(guān)注。

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

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

[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]陳震, 肖熙.平底結(jié)構(gòu)砰擊壓力峰值分析[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]張曉波.船底結(jié)構(gòu)砰擊時的流固耦合數(shù)值模擬[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]王輝.船體結(jié)構(gòu)局部強度設計中的砰擊載荷確定方法[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.