于龍超，閆發(fā)鎖，趙九龍，王 淇

(1.哈爾濱工程大學(xué)，黑龍江 哈爾濱 150001; 2.滬東中華造船(集團)有限公司，上海 200129)

不同曲率球面體入水砰擊載荷數(shù)值計算

于龍超1，閆發(fā)鎖1，趙九龍2，王 淇1

(1.哈爾濱工程大學(xué)，黑龍江 哈爾濱 150001; 2.滬東中華造船(集團)有限公司，上海 200129)

從三維邊界元方法出發(fā)，基于Wagner的自由液面抬升理論，采用數(shù)值仿真的方法研究了圓球入水的問題。通過數(shù)值模擬與軸對稱體入水試驗結(jié)果的正確性及適用性進行驗證的基礎(chǔ)上，并分析了不同半徑的球體在不同速度下的入水總體受力和表面壓力的變化規(guī)律。計算結(jié)果表明隨著球體半徑和入水速度的增加，峰值壓力迅速增加，受力峰值與球體半徑的平方和入水速度平方成正比關(guān)系；球體表面壓力系數(shù)的分布計算結(jié)果表明，在測試點與水開始接觸時壓力系數(shù)最大，然后迅速減?。粔毫ο禂?shù)和入水速度無關(guān)，但和入水深度相關(guān)。

球體砰擊；邊界元法；受力峰值；壓力系數(shù)

Abstract:Based on Wagner’s theory,the numerical simulation by 3D BEM method was used to study the slamming problem of spheres with different curvatures.The numerical results were validated by the test data from an axisymmetric body slamming test[11].The force acting on the body and the pressure at 5 specified locations are analyzed for the spheres with different radiuses and entry speeds.The result shows that the peak pressure increases rapidly with the increase of the radius and velocity,and the peak pressure has a linear relationship with the square of the radius of sphere and velocity.To calculate the distribution of the pressure coefficient along the surface of the sphere,we selected some test points in the sphere at different angles.The result shows the pressure coefficient has the maximum value when the test point touches water,and then decreases rapidly，and the pressure coefficient is sensative to the entry depth.

Keywords:sphere water entry； BEM； peak pressure； pressure coefficient

結(jié)構(gòu)物入水是一個瞬態(tài)過程，在砰擊入水過程中會受到巨大的沖擊載荷，甚至?xí)?dǎo)致結(jié)構(gòu)物的破壞，因此砰擊問題的研究對于入水結(jié)構(gòu)物的設(shè)計具有重要意義。如船舶在惡劣海況中航行的時候，由于船舶的劇烈起伏運動，首尾會頻繁的入水，入水過程中所產(chǎn)生的砰擊壓力會造成船體結(jié)構(gòu)的局部變形，嚴重的會造成船舶結(jié)構(gòu)的破壞。隨著海洋軍事和民用航海領(lǐng)域的不斷發(fā)展，砰擊問題廣泛存在于救生艇的落放、宇宙飛船的海上回收、海上飛機的降落，以及空投魚雷入水等研究領(lǐng)域之中。入水砰擊研究具有重要的工程應(yīng)用背景。

砰擊是一個十分復(fù)雜的物理現(xiàn)象，要涉及到固、液、氣三種介質(zhì)的耦合，理論分析比較困難，目前仍沒有一種可以適合各種情況的理論模型。為了問題的研究，通常使用結(jié)構(gòu)物入水模型作為力學(xué)模型來分析。自20世紀(jì)30年代以來，結(jié)構(gòu)入水問題的研究取得了很大的進展。Von Karman[1]最早采用附加質(zhì)量代替流體的作用分析入水沖擊的問題，提出了附加質(zhì)量方法來計算入水的沖擊載荷，將水上飛機降落在水面過程理想為二維楔形體入水的過程，并使用動量守恒定理推導(dǎo)出了入水沖擊載荷計算公式。Wagner[2]將Von Karman的方法理論化,同時考慮到物體入水時有液面升高的現(xiàn)象，提出了小斜升角模型的近似平板理論，為入水問題的研究提供了基礎(chǔ)。在Wagner的研究中，引入了水波影響因子，利用伯努利方程算出了沖擊壓力在結(jié)構(gòu)物濕表面上的分布情況，使理論更加適合實際問題。

近些年來，人們對于入水壓力的問題還是以二維問題的方面進行研究。Wu[3]等提出了二維楔形體的迭代求解方法，并且利用了邊界元方法，使用柯西積分對楔形體入水近似解進行了求解。但是在工程實際中使用的結(jié)構(gòu)形式多樣性，二維方法有一定的局限性，因此研究三維物體的數(shù)值算法對于工程應(yīng)用具有很強的實用價值。Faltisen等[4]以非軸對稱的船體模型入水為研究對象，基于邊界元方法、發(fā)展了Wagner理論，研究了三維入水砰擊的數(shù)值方法。Korobkin等[5-6]從逆推的Wagner問題及線性Wagner問題兩方面出發(fā)研究了三維入水砰擊理論，解決了流場流動、壓力分布及流域形狀三方面的問題。

針對三維物體入水問題，國內(nèi)學(xué)者針對實際應(yīng)用，陳震等[7]采用基于有限體積法的數(shù)值仿真方法對三維球鼻艏的入水砰擊問題進行了研究。王珂等[8]采用數(shù)值仿真的方法研究了三維剛性回轉(zhuǎn)體的入水問題，分析了不同回轉(zhuǎn)角度的回轉(zhuǎn)體的砰擊壓力峰值的變化規(guī)律，并研究了以一定初速度入水時，結(jié)構(gòu)質(zhì)量對砰擊壓力的影響。

基于三維邊界元方法進行了球面體三維入水方面的研究，在數(shù)值上方面從三維模型出發(fā)，繼承了Wagner的自由液面抬升理論，引入了浸深因子CW[9]以確定等效自由液面的高度，并考慮網(wǎng)格運動，將物體的入水過程分成若干時間步，一般根據(jù)入水速度及入水時長進行劃分，而且還要和網(wǎng)格的垂向尺寸相協(xié)調(diào)，計算每一個時間步對應(yīng)的濕網(wǎng)格中心點速度勢，直到時間步長的全部結(jié)束。然后再根據(jù)每一時間步記錄的數(shù)據(jù)進行后續(xù)的流場計算,實現(xiàn)對入水問題的求解。文中所使用的算法不同于以往的二維算法：首先，建立了實際的物面模型，并進行了網(wǎng)格劃分，使網(wǎng)格數(shù)足夠多，形成和實際結(jié)構(gòu)物面極為相似的三維邊界元模型；同時考慮了三向流場，在實際意義上更接近三維算法。

1 數(shù)值求解方法

1.1Wagner方法和自由液面的處理

問題的假設(shè)是二維物體或者三維軸對稱體豎直落入半無界水域，入水模型被視為剛性結(jié)構(gòu)，并且認為物體在落水的短時間內(nèi)速度保持定常。在入水過程中，作用于模型上的流體慣性力占主導(dǎo)地位[10]，忽略流體黏性、表面張力、可壓縮性、重力以及入水模型表面與流體間的氣墊效應(yīng)，并且認為流體是無旋的。Von Karman最早采用附加質(zhì)量的方法計算物體入水的受力問題。物體表面的邊界條件施加于流體自由表面的初始位置，通過引入線性化的伯努利方程，對自由面邊界條件進行簡化。這里在原始Wagner方法的基礎(chǔ)上進行了改進，引入了邊界元方法，考慮具體浸水物面，采用非線性伯努利方程對壓力進行求解。

假設(shè)水是理想不可壓縮的無旋流體。這樣，流場速度勢φ滿足Laplace方程

在物面和無窮遠處滿足

式中：Vn是物面外法線n方向上(指向流場)的速度。

將濕物面劃分為N個面元，可以得到離散化的方程組：

式中：aij是j面元對i面元的影響系數(shù)，σj為j面元上的分布源密度，其中aij為

通過式(3)和(4)便可以求解出各面元中心點處的分布源密度，進一步可得到個面元中心點的速度勢

進一步可以求出各面元中心點處的誘導(dǎo)速度

忽略大氣壓力和重力影響，由Bernoulli方程可求得物面沖擊壓力

物體在流場中所受到的力為

圖1 發(fā)展的Wagner方法Fig.1 Generalized Wagner method

在VonKarman的方法中，忽略了自由表面的升高，物面與水面的交線位于靜水面，很容易通過物體入水的深度和物面形狀得到。如圖1所示，在Wagner方法中，自由表面條件應(yīng)用在升高的自由表面上，也就是考慮了自由液面升高的影響。在笛卡爾坐標(biāo)系o-xyz中，物體形狀可以表示為Z=f(x,y)。如果t表示物體以速度V(t)入水的時間，則由于物體的運動導(dǎo)致物面形狀方程為

t時刻自由液面的高度使用z=ζ(x,y,t)表示，在歐拉系統(tǒng)中滿足dζ/dt=φz，所以點(x,y)處自由表面的升高為

應(yīng)注意到垂向速度φz是從升高的自由表面，而不是原靜水面z=0。如果定義物面和自由水面的交界線為

由式(9)和(10)得

與Von-Karman的理論相比，可以看到濕表面和升高的自由水面是預(yù)先未知的。所以，它們的交界線，即物體的水線也是未知的。所有這些都需要求解后才知道。以上的過程也稱為發(fā)展的Wagner方法。

文中所使用的算法在對自由液面進行線性化處理的時候，根據(jù)Wagner理論，考慮到了液面抬升。通常情況下，物體在入水過程中，自由液面與物體相交的位置均要高于原始的自由液面，流體有沿物面向上爬升的趨勢，嚴格來說在計算的過程中必須要追蹤自由液面的位置，這樣才能使求解更加準(zhǔn)確，但是對三維入水模型來講，采用式(12)進行自由液面的追蹤數(shù)值處理上還存在困難。例如文獻[4]在三維船體入水砰擊的自由液面追蹤中，采用的方法是沿著船體與自由液面交界的一圈均勻標(biāo)記若干點，之后分別以這些標(biāo)記點為起點，平行于原始的自由液面畫出射線，最終由作出的許多射線近似的擬合出一個抬升后的自由液面。本文算法也沿用了這種“射線延伸”的思想。采用線性化的自由液面，采用經(jīng)驗公式換算出物面與自由液面每一時刻的觸點高度，以確定該觸點的位置，并以之為起點平行于水平面的平面作為新的自由液面，即是“自由液面”。為了確定等效自由液面的位置，同時引入了浸深因子Cw=η/b。其中b=vt為物體的入水深度，η為自由液面升高加上b。

對于圓球體等軸對稱鈍體結(jié)構(gòu)，Cw的取值如下：

其中，D是圓球的直徑。

圖2 自由液面處的網(wǎng)格處理方法Fig.2 Processing method for the grid at free surface

1.2網(wǎng)格入水過程中的處理

文中的計算程序是把三維邊界元模型入水的時間歷時分為若干個時間步進行了分析，在每一時間步時刻都要統(tǒng)計等效自由液面以下的面元以及節(jié)點的信息。但是在實際的計算模擬中，只有在水面下離等效自由面較遠的網(wǎng)格能夠保持完整性，在等效自由液面附近，沿著模型表面與等效自由面的交界線有一圈網(wǎng)格被自由液面所截斷，生成新的面元網(wǎng)格；由于在程序中所使用的網(wǎng)格為長寬比值適中的四邊形面元，并且為了提高程序的準(zhǔn)確性，不能單純的將這部分網(wǎng)格忽略，因此在程序中對于這部分網(wǎng)格進行了截斷重生處理，如圖2所示。

2 砰擊載荷的數(shù)值計算結(jié)果與試驗比較

為了研究高速物體結(jié)構(gòu)入水過程中所受到的水動力壓力大小，Adrian等[11]完成了一系列的入水受力測試試驗，主要測試了若干鋁合金制軸對稱體以定常速豎直入水過程中的垂向受力數(shù)據(jù)，以供研究探討。試驗中所使用的模型如圖3所示。

本文使用了三維邊界元方法對文獻[11]試驗進行了數(shù)值模擬，并與模型實驗的結(jié)果進行了對比。圖4中T1為Wagner方法對應(yīng)的半球體完全浸水時刻，通過球體模型的數(shù)值模擬可以發(fā)現(xiàn)，簡化的半球體加密封蓋模型在完全入水之前的壓力歷時和球體模型試驗的結(jié)果吻合較好；在入水之后，由于半球缺少球體的一部分的表面積，所以在入水之后半球的總體受力迅速減小，與試驗相比有誤差。實際上，我們重點關(guān)注的是球體入水的最大沖擊力，出現(xiàn)在球體入水的初期。所以，為了研究球體入水的峰值壓力問題，采用了簡化的半球體加密封蓋模型，能夠準(zhǔn)確計算出半球體入水的最大受力，也是表面最大壓力的大小和發(fā)生的時間。雖然T1之后的時段結(jié)構(gòu)有偏差，但是這種情況下大大減少了計算量，提高了計算效率。

圖3 試驗中使用的模型Fig.3 Model utilized in the test

圖4 半球以20 m/s速度入水垂向受力計算值與試驗值比較Fig.4 Comparison of the pressures from calculation and model test(V =20 m/s)

3 數(shù)值模擬及結(jié)果分析

在驗證了算法的正確性之后，使用Abaqus/CAE程序建立了相應(yīng)的邊界元模型分別建立了半徑為0.5、0.6、0.8、1.0 m半球型邊界元模型，半球縱深H=D/4，并進行了網(wǎng)格的劃分和節(jié)點的提取，網(wǎng)格和節(jié)點數(shù)見表1所示。

表1 模型的網(wǎng)格劃分Tab.1 Division of the model grid

圖5 邊界元模型Fig.5 Mesh of boundary elements

所建立的邊界元模型的外表面為物體的濕表面，同時也是面元法線的正方向一側(cè)。圖5(a)是實際入水試件表面，在計算入水壓力時，只考慮其表面上的壓力值。在半球完全入水之前，參與計算的僅僅是其中的一部分面元，這樣可以與抬升后的自由液面對稱可得到一封閉的面元體結(jié)構(gòu)，這樣就可以根據(jù)面元法求解方程，進行面元源強的求解；隨著時間的推移，當(dāng)所有面元全部位于抬升后的自由液面以下時，這樣與等效自由液面得到的面元體就不是閉合結(jié)構(gòu)，不符合勢流理論的求解原理，方程不成立。因此根據(jù)實際情況邊界元模型加上了一個上蓋，如圖5(b)所示的結(jié)構(gòu)。這樣在圖5(a)中面元完全入水之后，圖5(a)和圖5(b)共同組成了一個閉合的邊界元模型，經(jīng)與抬升后的自由液面對稱便可以形成兩個閉合的邊界元模型。在計算的過程中，認為在半球入水之后液體立即將半球吞沒。

為了研究球面上不同位置的表面壓力分布，在計算的過程中針對性的選取了5個位置(圖6)，計算了這些位置處的壓力，并統(tǒng)計了壓力峰值。在砰擊試驗中，物面各點的壓力常用p=0.5CpρV2表示。其中Cp為無量綱沖擊壓力系數(shù)，反映物面各點的壓力水平；ρ為流體密度；V指物體入水速度。采用無量綱的表示形式，還計算了半球入水過程中壓力系數(shù)沿球體表面的分布情況，選取的測點如圖6所示。這些點在下文中用其所在位置的徑向射線與球體對稱中心線之間的夾角表示。這些點的標(biāo)號在圖6中由左至右分別為：0°單元、10°單元、20°單元、30°單元和40°單元。

圖6 半球體表面壓力計算點分布Fig.6 Location of the points chosen for pressure calculation

圖7給出了不同半徑在不同速度下入水壓力隨時間的變化曲線，通過計算可以發(fā)現(xiàn)在球體在入水之后，壓力迅速達到峰值，然后迅速回落。針對同一半徑的球體，受力峰值對入水速度非常敏感。例如圖7(c)，半徑為0.8 m的球體在10 m/s時的受力峰值約為100 kN，而速度為20 m/s時受力的峰值約為400 kN。不同半徑的球體以不同速度入水時的受力峰值見圖8所示，其中橫坐標(biāo)用速度的平方表示?？梢娗蝮w的峰值受力與入水速度的平方呈線性關(guān)系。針對不同的入水速度，受力峰值與半徑的關(guān)系統(tǒng)計見圖9所示，同樣存在線性的比例關(guān)系。

圖7 球體以不同速度入水時垂向受力的歷時曲線Fig.7 Time history of vertical forces acting on the sphere at different entry speeds

圖8 受力峰值隨入水速度平方變化Fig.8 Peak force with t square of entry velocity

圖9 受力峰值隨半球半徑平方變化Fig.9 Peak force with square of sphere radius

圖10為半徑1 m的球體在5個位置點的壓力系數(shù)隨時間的變化。因為球體入水過程中，0°單元先著水，依次是10°、20°、30°和40°單元。從時歷變化可見，各點的壓力在觸水時刻達到最大值，然后呈級數(shù)形式下降。先觸水的幾點壓力在初期的數(shù)毫秒內(nèi)下降迅速，然后降速變緩。點的位置與球體中心線越接近，在入水時間內(nèi)的壓力水平越高。圖11為不同半徑的球體以20 m/s的速度入水時，0°單元上壓力的時間變化。半徑越大，其壓力水平越高。

圖10 半徑1 m半球時壓力系數(shù)(V=10 m/s)Fig.10 Time history of Cp at 10 m/s entry speed

圖11 不同半徑球體0°面元壓力系數(shù)(V =20 m/s)Fig.11 Time history of Cp on 0° element

半徑為0.5 m的球體上的0°單元以不同速度入水時壓力與入水深度之間的關(guān)系見圖12所示?？梢钥闯觯煌娜胨俣惹闆r下，無量綱的壓力系數(shù)幾乎重合。說明對同一球體而言，只球體上各點的Cp與入水速度無關(guān)。而對于不同球體而言，綜合考慮壓力與浸入深度b和球體尺度r的關(guān)系，可以得到一個統(tǒng)一的壓力系數(shù)，如圖13所示。

圖13 不同半徑球體0°面元壓力系數(shù)(V =20 m/s)Fig.13 Cp on 0° element of different spheres

4 結(jié) 語

以三維邊界元方法編寫了半球入水的砰擊載荷計算程序，與軸對稱體入水試驗的結(jié)果進行對比，驗證了程序的正確性，在此基礎(chǔ)上可以對方法進行改進，對任意三維物體入水砰擊載荷進行計算。

1)通過數(shù)值計算可以發(fā)現(xiàn)半球在入水的過程中，最大砰擊力出現(xiàn)在入水初期，并在峰值出現(xiàn)后迅速減?。蛔杂梢好嬖谌胨跗趯η蝮w的影響較小，在球體入水的后期影響較大。

2)通過不同半徑和不同速度的半球模型對比可以發(fā)現(xiàn)，入水的最大砰擊力與半徑的平方成正比例關(guān)系，與入水速度的平方成正比關(guān)系。

3)通過在球體表面選取的測試點壓力系數(shù)的計算，結(jié)果表明壓力系數(shù)在測試點和水接觸時達到最大值；球體壓力系數(shù)隨著球體在水中的深度而變化，局部壓力系數(shù)只取決于水中形狀，和入水的速度無關(guān)。

[1] VON KARMAN T.The impact of seaplane floats during landing[R].Washington D C,USA:National Advisory Committee for Aeronautics,NACA Technical Notes 321,1929.

[2] WAGNER V H.Phenomena associated with impact s an d sliding on liquid surfaces[J].Z Angew Math Mech,1932,12(4):193-215.

[3] 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 Structures,2004,19:277-289.

[4] FALTINSEN O M,CHEZHIAN M .A generalized wagner method for three-dimensional slamming[J].Journal of Ship Research,2005,49(4):279-287.

[5] SCOLAN Y M,KOROBKIN A A.Three-dimensional theory of water impact:Part 1.inverse wagner problem[J].Journal of Fluid Mechanics,2001,440:293-326.

[6] KOROBKIN A A,SCOLAN Y M.Three-dimensional theory of water impact:Part 2.linearized wagner problem[J].Journal of Fluid Mechanics,2006,549:343-373.

[7] 陳震,肖熙.三維球鼻艏入水砰擊研究[J].船舶工程,2007,29(4):61-64.(CHEN Zhen,XIAO Xi.Impacting study on the water entry of 3D bulbous bows[J].Ship Engineering,2007,29(4):61-64.(in Chinese))

[8] 王珂，陳剛，袁洪濤.三維回轉(zhuǎn)體入水砰擊載荷預(yù)報[J].船舶工程,2012,34(1):12-15.(WANG Ke,CHEN Gang,YUAN Hongtao.Prediction of the slamming pressure on a 3-D axisymmetric structure[J].Ship Engineering,2012,34(1):12-15.(in Chinese))

[9] JR A B W,MORRISON A M,BALDWIN J L.Prediction of impact pressures forces and moments during vertical and oblique water entry[M].1977:SWC/WOL/TR 77-16.

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

[11] CONSTANTINESCU A,ALAOUI A E M,NEME A,et al.Numerical and experimental studies of simple geometries in slamming[J].International Journal of Offshore and Polar Engineering,2011,21(3):216-224.

Numerical calculation of slamming load for different curvature of the sphere

YU Longchao1,YAN Fasuo1,ZHAO Jiulong2,WANG Qi1

(1.Harbin Engineering University,Harbin 150001,China; 2.Hudong-Zhonghua Shipbuilding (Group) Co.,Ltd.,Shanghai 200129,China)

O353.4

A

10.16483/j.issn.1005-9865.2016.01.005

1005-9865(2016)01-0033-07

2015-01-08

于龍超(1990-)，男，河南鹿邑人，碩士生，主要從事海洋工程結(jié)構(gòu)物設(shè)計分析。E-mail:yanfasuo@163.com