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

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

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

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

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

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

軸對(duì)稱體；砰擊載荷；液面升高；網(wǎng)格截?cái)?；虛擬面元

0 引 言

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

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

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

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

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

1 數(shù)值算法

1.1 改進(jìn)Wagner方法在砰擊載荷計(jì)算中的應(yīng)用原理

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

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

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

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

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

通過求解方程組（3）便可以求解出各面元中心點(diǎn)處的分布源密度，進(jìn)一步可得到個(gè)面元中心點(diǎn)的速度勢(shì)：

緊接著可以求出各面元中心點(diǎn)處的誘導(dǎo)速度：

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

計(jì)算得到各時(shí)刻對(duì)應(yīng)的繞流壓力。

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

1.2 自由液面處理方法

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

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

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

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

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

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

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

1.3 面元網(wǎng)格的運(yùn)動(dòng)及重構(gòu)

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

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

圖3 自由液面交界處的網(wǎng)格截?cái)嗉爸貥?gòu)方法Fig.3 The gridding truncation method nearby the boundary of the free surface and the object

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

1.4 砰擊后時(shí)段算法適用性研究

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

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

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

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

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

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

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

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

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

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

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

3 結(jié) 語

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

（2）算法中考慮了網(wǎng)格運(yùn)動(dòng)，在自由液面附近對(duì)網(wǎng)格進(jìn)行截?cái)嘀貥?gòu)處理，精準(zhǔn)了每一時(shí)間步對(duì)應(yīng)的水下濕面積。

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

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

圖8 半球阻力系數(shù)時(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]戴遺山,段文洋.船舶在波浪中運(yùn)動(dòng)的勢(shì)流理論[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-），男，教授。