基于改進Wagner方法的軸對稱體入水“后時段”砰擊載荷數(shù)值計算

趙九龍，閆發(fā)鎖，樊 磊

（1.哈爾濱工程大學，哈爾濱150001；2.交通運輸部 上海打撈局，上海 200090；3.德州農工大學，德克薩斯州 77840，美國）

基于改進Wagner方法的軸對稱體入水“后時段”砰擊載荷數(shù)值計算

趙九龍1，2，閆發(fā)鎖1，樊 磊3

（1.哈爾濱工程大學，哈爾濱150001；2.交通運輸部 上海打撈局，上海 200090；3.德州農工大學，德克薩斯州 77840，美國）

文章基于三維邊界元方法研究了三維軸對稱體入水砰擊載荷的數(shù)值算法。算法從三維力學模型出發(fā)，繼承Wagner自由液面抬升理論，引入浸深因子Cw以確定自由液面抬升高度，將自由液面線性化處理，同時考慮網格運動，在自由液面附近對網格進行截斷重構，以確保水下濕面積的精準。算法中使用考慮加權運動項的非線性伯努利方程計算得到入水結構的表面壓力，進一步積分可得入水結構的總體受力；另外，算法中引入虛擬面元的概念，將砰擊載荷計算時歷延長至液面高于墜物制高點之后，擴大了傳統(tǒng)入水載荷計算的時歷范圍，并且文中引用Alaoui半球入水試驗，對算法的正確性及適用性進行了驗證。

軸對稱體；砰擊載荷；液面升高；網格截斷；虛擬面元

0 引 言

結構物在短時間內迅速進入液體，隨之產生液體飛濺的現(xiàn)象稱之為砰擊。然而，砰擊涵蓋的面非常廣泛，廣義的砰擊概念是指結構物與流體發(fā)生的碰撞現(xiàn)象，一般都涉及到固、液、氣三種介質的相互耦合，是一種瞬態(tài)的強非線性、強耦合和極其復雜的物理現(xiàn)象。在工程實際中，砰擊問題非常易見，例如：水上飛機水面降落、船舶在惡劣海況中的底部砰擊、艦船上的救生艇落水、航天運載器的水面回收、魚雷空投入水等問題都屬于該范疇。

砰擊現(xiàn)象對應的的實際力學模型要涉及到固、液、氣三種介質的耦合，而且具體情況復雜多變，目前為止仍未能提出一種貼切真實情況、普適性強的理論模型。但是為了便于對問題展開研究，通常使用“結構物入水沖擊”模型作為具體的力學模型來分析。

許多工程結構物在高速入水時，所面臨的挑戰(zhàn)是入水過程瞬時所受到的巨大水動壓力，這往往會導致入水結構的局部大變形甚至破斷。近年來人們在計算入水壓力方面做了許多工作，但是大部分都是以二維視角展開研究，研究基于的入水模型也多見于二維楔形體，數(shù)值方面做得比較好的有Wu等[1-2]研究出了二維V形楔的迭代求解方法，并且運用二維邊界元法，使用柯西積分計算了楔形體入水的相似解，許國冬[3]討論了V形楔多種入水形式的相似解。但是由于實際工程應用中人們使用的結構形式多樣，二維入水算法存在著一定的應用局限性，因此對于三維入水算法的開發(fā)有著迫切的需求。

目前還沒有非常成熟的三維入水算法供于應用。近年來的Korobkin和Scolan[4-5]從逆推的Wagner問題及線性Wagner問題兩方面出發(fā)研究了三維入水砰擊理論，解決了流場流動、壓力分布及流域形狀三方面的問題。Faltisen和Chezhian[6]以非軸對稱的船體模型入水為研究對象，基于邊界元方法、發(fā)展的Wagner理論，研究了三維入水砰擊的數(shù)值方法。

基于此，本文也在三維入水方面進行了研究，主要研究了基于三維邊界元方法的砰擊理論算法及其在砰擊載荷計算中的應用。數(shù)值方面的研究從三維力學模型出發(fā)，基于三維邊界元方法，繼承Wagner自由液面抬升理論，引入了“浸深因子”Cw[7]以確定自由液面抬升的高度，在這一高度處建立新的自由液面并稱之為“等效自由液面”。算法中考慮到網格運動，對入水過程中每一時間步時刻對應的網格編號、節(jié)點坐標等進行追蹤計算，由此得出網格與等效自由液面的位置關系，進一步可精確追蹤到每一時間步時刻對應的濕網格數(shù)據(jù)，同時對等效自由液面附近的網格根據(jù)實際情況進行截斷重生，以確保水下濕面積的精確；通過時間步的循環(huán)計算，最終可得到入水結構的表面“壓力”和總體“受力”兩種載荷時歷。另外，提出了虛擬面元的概念，使算法可以應用于墜物浸水之后的砰擊時段（液面高于入水物體的最高點，即所謂的入水“后時段”，將傳統(tǒng)的砰擊載荷計算時歷延長。

1 數(shù)值算法

1.1 改進Wagner方法在砰擊載荷計算中的應用原理

Wagner在Von-Karman[7]線性砰擊理論的基礎上，提出了自由液面的抬升效應，采用平板擬合法，相對精確了浸濕半寬、附加質量、沖擊力及壓力分布的求解。在求解過程中所采用的相當平板是對具體物面的近似，而且采用線性伯努利方程對壓力分布進行求解，整個過程雖然計算量小，但是計算精度還有待進一步提高。本文在原始Wagner方法的基礎上進行了改進，引入了邊界元方法[8]，考慮具體浸水物面，采用非線性伯努利方程對壓力進行求解。

談及流體載荷的計算，近年來出現(xiàn)了多種計算方法，比較常用的有攝動法、有限元法、有限差分法、邊界元法。由于邊界元方法將空間的問題轉化為邊界表面的積分問題，將問題減小一維，大大減少了未知量的數(shù)目，另一方面它能夠處理無窮遠處的邊界條件，對外部問題有很好的適用性，正是由于這兩點才使得邊界元方法在流體載荷計算中得以迅速發(fā)展。本文所使用的邊界元方法適用于求解無界流場中的三維無升力繞流問題，分布奇點法的所有要點都包含在內。速度勢可以用分布源和分布偶極來描述，本文算法僅僅涉及到分布源的概念。

控制方程及邊界條件如公式（1）、（2）：

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

式中：aij是j面元對i面元的影響系數(shù)：

通過求解方程組（3）便可以求解出各面元中心點處的分布源密度，進一步可得到個面元中心點的速度勢：

緊接著可以求出各面元中心點處的誘導速度：

再進一步可根據(jù)各時刻各點的速度勢數(shù)值，對時間進行微分計算便可得到這樣就可以使用非線性的伯努利方程：

計算得到各時刻對應的繞流壓力。

從上面的推導可以計算出無界流中運動物體的表面繞流壓力值，但是人們通常研究的入水砰擊現(xiàn)象發(fā)生在半無界流場中，而且正是自由液面的附近。為了解決理論的適用問題，考慮線性的自由液面，由于自由液面上的速度勢φ=0，即可將自由液面下的濕表面沿著自由液面進行上下對稱處理，形成一個沿自由液面對稱的閉合曲面s+s′，這樣就將問題等同于考慮在無界流場中運動的s+s′，使問題適用于上述理論。

1.2 自由液面處理方法

本文算法在對自由液面進行線性化處理的時候，根據(jù)Wagner理論，考慮到了液面抬升，如圖1所示的入水模型。通常情況下，物體入水過程中，自由液面與物體相交點處的位置均要高于原始的自由液面，流體沿物面有向上爬升的趨勢，嚴格來講計算過程中必須要追蹤自由液面的運動形式，這樣一定會使求解更加精確，但是對三維入水模型來講，自由液面的追蹤非常的困難，而且鮮有適用性比較廣泛的方法，例如Faltinsen和Chezhian[6]，在三維船體入水砰擊的自由液面追蹤中也不能夠很精確地描述自由液面的具體形狀及變化，他們所采用的方法是沿著船體與自由液面交界的一圈均勻標記若干點，之后分別以這些標記點為起點，平行于原始的自由液面畫出射線，最終由做出的許多射線近似地擬合出一個抬升后的自由液面。該方法目前來講無疑是一種比較準確的追蹤方法，基于此，在自由液面的處理過程中，本文程序也沿用了這種“射線延伸”的近似，但是較之有所簡化。

從工程實用角度，為了簡化流固交界附近數(shù)值處理的難度和計算量。本文算法采用線性化的等效自由液面，使用經驗公式換算出物面與自由液面每一時刻的觸點高度，以確定該觸點位置，并以之為起點平行于水平面引出射線，認為該射線所在的水平面便是新的自由液面，在這里稱之為“等效自由液面”。為了確定等效自由液面的位置，同時引入了浸深因子Cw的概念，如圖2所示。

圖1 Faltinsen關于自由液面的追蹤Fig.1 The free surface tracking method by Faltinsen

圖2 本文采用的自由液面處理方法Fig.2 The free surface caculation by this paper

如果知道Cw值，便很容易可以確定等效自由液面的位置。本文參照了美國水面武器研究中心白橡樹實驗室的一份研究報告[9]，其中詳細介紹了軸對稱物體垂向入水與任意物體傾斜入水兩種類別的Cw值選取辦法。

表1 Cw選取方法Tab.1 Calculation of Cw

顯然使用該方法可以使砰擊載荷計算得到很大程度上簡化，提高了計算效率。

1.3 面元網格的運動及重構

本文將邊界元模型入水運動的過程劃分為若干時間步進行分析，對于每一時間步時刻都要統(tǒng)計出等效自由液面以下的面元及其節(jié)點信息，這是計算誘導速度系數(shù)的基本前提。很顯然，水面以下遠離等效自由液面的網格保持原有的完整性，因此不必對其幾何形狀進行改變；而在等效自由液面附近，沿著物體與自由液面的交界線有一圈網格被等效自由液面截斷，生成新的四邊形、五邊形和三角形面元；但是考慮到本文算法中要求提供的面元必須是長寬比比值適中的四邊形面元，而且為了提高程序的精確性，不能單純地將這部分網格進行忽略，因此在程序中對這部分網格進行了截斷重生處理，具體的處理方法如下圖3所示。

另外，考慮到網格運動，于是便不能直接使用常規(guī)的非線性伯努利方程來計算載荷，需要對其進行改進，添加運動項目。根據(jù)泰勒公式變換得到：

圖3 自由液面交界處的網格截斷及重構方法Fig.3 The gridding truncation method nearby the boundary of the free surface and the object

圖4 墜物砰擊后時段邊界元模型示意圖Fig.4 The BEM of the later entry time

1.4 砰擊后時段算法適用性研究

按照上述的方法僅僅可以計算墜物未完全浸水的砰擊時段，而對于墜物制高點低于等效液面之后的時段則無法進行計算。在物體完全浸水之后，短時間內物體的上表面事實上還沒有被液體覆蓋，因此在計算中物體的濕表面s與投影面s′將無法組成封閉面，本文所引用的無界流面元法則不能適用，在后續(xù)求解方程組的過程中會出現(xiàn)無解的情況。

如圖4所示，T1時刻之后s與s′分離，之后便是本文所指的“后砰擊”時段。

對此，本文提出了添加虛擬面元的方法，使入水“后時段”的砰擊載荷計算得以實現(xiàn)，下圖5是本文所設計的三種入水邊界元模型。圖中a是實際的入水邊界元模型，b是添加的虛擬邊界元模型，每個面都有兩種顏色，其中棕色代表物體的觸水濕表面。第一種模型根據(jù)實際情況對a圖邊界元模型加了一個上蓋，并且劃分邊界元網格，即b圖所示圓形結構。這樣在T1時刻之后，a、b共同組成了一個閉合的邊界元模型，經與抬升后的自由液面對稱便可以形成兩個閉合的邊界元模型，符合理論要求。如果使用這樣的水動力模型進行計算的話，實際上就是認為在半球入水之后液體立即將半球體吞沒，即半球上蓋之上被流體覆蓋。第二種模型，實際上就是模擬了在T1時刻之后a的上緣與抬升后的自由液面間的液面形狀，這一段液面上一定存在著分布源，它對流場性質有重要的影響，本文計算時將其視作濕表面的一部分，可以計算出其上的分布源強度，從而對a表面面元速度勢進行修正，最后在計算試件垂向壓力時并不是對a和b上的所有面元進行加權計算，僅僅考慮a上的面元壓力，b僅僅是作為一個虛擬的濕表面。第三種模型類似于第二種，只是假想T1時刻后a上緣的自由液面成圓筒狀，因此將b設計成了圓筒狀的一層面元。

圖5 3種虛擬邊界元模型的添加實例Fig.5 3 adding methods for the virtual surface elements

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

2.1 Alaoui剛性半球定常速垂向入水[15]

本文使用Abaqus/CAE程序建立了相應的邊界元模型。半球的縱深H=77.34 mm，根據(jù)Nisewanger提供的半球入水自由液面升高公式計算可得在T1時刻半球恰恰完全浸沒于水（半球底面恰與等效自由液面平齊）。對于18 m/s的情況來講，T1=2.61 ms；對于20 m/s的情況，T1=2.41 ms。試驗分別給出了6 ms和5 ms的時歷曲線，因此可以斷定T1之后的數(shù)值是完全浸沒之后的數(shù)值（試件最高點低于等效自由液面），但是由于入水情況的瞬時性，半球周圍的液體被拍擊而濺向四周，在浸沒之后的短時間內是不會將水下物體覆蓋，從以往的試驗中也可以觀察到這種現(xiàn)象。這就提出了新的問題：T1時刻之后的濕表面的構成問題，即要探尋一種合適的水動力模型（邊界元模型）來精確T1時刻之后的壓力計算。

根據(jù)圖5提出的三種添加邊界元模型的方法進行砰擊載荷求解，計算結果如圖6-8所示。

圖6 18 m/s半球入水受力時歷Fig.6 The force-time curve of the hemisphere for 18 m/s

圖7 20 m/s半球入水受力時歷Fig.7 The force-time curve of the hemisphere for 20 m/s

可見使用邊界元模型2的虛擬面元形式計算結果與實驗匹配良好，所以對于半球入水“后時段”的砰擊載荷計算可以采用該添加虛擬面元的方法進行計算。

3 結 語

（1）本文基于三維邊界元方法，編寫了可用于計算軸對稱結構入水砰擊載荷的計算程序，在自由液面處理方面繼承Wagner液面抬升理論，引入“浸深因子”Cw以確定液面抬升高度，簡化了自由液面處理方式，提高了程序計算效率。

（2）算法中考慮了網格運動，在自由液面附近對網格進行截斷重構處理，精準了每一時間步對應的水下濕面積。

（3）使用該程序算法計算了Alaoui圓球體入水載荷，載荷計算結果與試驗數(shù)據(jù)匹配良好，驗證了本文計算程序的可行性與準確性；并從另一方面說明了線性自由液面對于砰擊載荷計算有較好的適用性。

（4）提出了入水“后時段”的砰擊載荷計算方法，擴大了傳統(tǒng)的砰擊載荷時歷計算范圍。

圖8 半球阻力系數(shù)時歷Fig.8 The resistance-coefficients-time curve of the hemisphere

[1]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.

[2]Wu G X.Numerical simulation of water entry of twin wedge[J].Journal of Fluids and Structures,2006,22:99-108.

[3]Xu G D,Duan W Y,Wu G X.Numerical simulation of oblique water entry of an asymmetrical wedge[J].Ocean Engineering,2008,35(16):1597-1603.

[4]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.

[5]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.

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

[7]Von Karman T.The impact on seaplane floats during landing[R].NACA TN 321.Oct.1929.

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

[9]Prediction of impact pressure,forces,and moments during vertical and oblique water entry[R].NSWC.15,January,1977.

[10]Nisewanger C R.Experimental determination of pressure distribution on a sphere during water entry[J].NAVWEPS 7808, 1961.

[11]Faltinsen O M,Landrini M,Greco M.Slamming in marine applications[C].Journal of Engineering Mathematics,2004,48: 187-217.

[12]Krobkin A A,Pukhnachou V V.Initial stage of water impact[J].Annual Review Fluid Mechanics,1988,20:159-185.

[13]Takagi K,Dobashi J.Influence of trapped air on the slamming of a ship[J].Journal of Ship Research,2003,47(3):187-193.

[14]Sharov V F.Ship bottom impact upon wave[J].Sudostroyeniye,1958,4:5-9.

[15]Aboulghit EI Malki Alaoui,Alain Neme,Nicolas Jacques.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 the slamming load for axisymmetric geometries during the later entry time based on generalized Wagner theory

ZHAO Jiu-long1,2,YAN Fa-suo1,FAN Lei3
(1.Harbin Engineering University,Harbin 150001,China;2.Shanghai Salvage Co./Ocean C&I,Shanghai 200129,China; 3.Texas A&M University,College Station,Texas 77840,USA)

A numerical code used for calculating the slamming load for axisymmetric geometries based on 3D BEM was mainly introduced.The 3D mechanical model and Wagner’s theory which considering the uplift free surface were both used in the numerical code.Cw which was called wetting factor was introduced to ensure the exact height of the uplift free surface.At the same time,the grid’s movement was also considered and along the interface between the geometry and the free surface one group of grids were modified for a precise wet area.The Nonlinear Bernoulli equation considering the grid’s movement was used for calculating the hydrodynamic pressure and the force could be obtained through integration.Besides,the conception of Virtual surface element was brought here so that the slamming load when the top point of the falling object was under the free surface could be achieved.That is to say,the traditional calculated time process was extended.At the same time,Alaoui slamming test was used to test the code’s validity.

axisymmetric geometries;slamming load;uplift of the free surface; grid’s modification;virtual surface element

TV131.2

A

10.3969/j.issn.1007-7294.2016.07.007

1007-7294（2016）11-1412-08

2016-04-07

趙九龍（1988-），男，工程師，E-mail：tcjiulong@sina.com；閆發(fā)鎖（1977-），男，教授。