段文洋，陳紀(jì)康，趙彬彬

（哈爾濱工程大學(xué)船舶工程學(xué)院，黑龍江哈爾濱150001）

基于泰勒展開(kāi)邊界元法的深水浮體二階平均漂移力計(jì)算

段文洋，陳紀(jì)康，趙彬彬

（哈爾濱工程大學(xué)船舶工程學(xué)院，黑龍江哈爾濱150001）

為了消除傳統(tǒng)低階面元法計(jì)算非光滑邊界單元處的切向誘導(dǎo)速度精度低的缺點(diǎn)，提出了一階泰勒展開(kāi)邊界元方法。對(duì)邊界積分方程中偶極強(qiáng)度在單元內(nèi)進(jìn)行泰勒展開(kāi)，保留一階導(dǎo)數(shù)項(xiàng)，并對(duì)邊界場(chǎng)點(diǎn)求2個(gè)正交方向的切向?qū)?shù)來(lái)封閉方程組，從而可同時(shí)求解單元中的速度勢(shì)及其切向速度。數(shù)值計(jì)算表明，該方法可提高浮體邊界拐角處切向速度和水線處波面爬高的計(jì)算精度，進(jìn)而保證浮體二階平均波浪漂移力的計(jì)算精度和收斂性。

二階波浪力；泰勒展開(kāi)；邊界元方法；深水浮體；柯欽函數(shù)；近場(chǎng)公式；遠(yuǎn)場(chǎng)公式

Hess等［1］將邊界元方法拓展到求解浮體水動(dòng)力學(xué)問(wèn)題中。從邊界積分方程數(shù)值離散看，邊界元法可分為低階邊界元方法（常值面元法）和高階邊界元方法。

對(duì)于低階邊界元方法，通常用2種方法求解浮體水動(dòng)力問(wèn)題，即：利用源-偶混合分布方法求解速度勢(shì)，分布源方法求解誘導(dǎo)速度［2］，或者直接利用分布源方法求解速度勢(shì)和誘導(dǎo)速度［3］。對(duì)于任意邊界，只要面元數(shù)足夠，這2種方法求解速度勢(shì)都能得到令人滿意的精度，但是第1種方法的收斂性更好［4］。而對(duì)于非光滑邊界，2種方法計(jì)算得到的非光滑邊界面元處的誘導(dǎo)速度均不精確。這導(dǎo)致由速度平方項(xiàng)和水線波浪爬高產(chǎn)生的二階定常波浪載荷隨著單元數(shù)量增加不收斂［2］。

在相同面元數(shù)下高階邊界元方法極大改進(jìn)了計(jì)算精度。Choi等［5?8］表明高階邊界元方法計(jì)算收斂較快。但是對(duì)于邊界具有拐角且需要求解速度勢(shì)導(dǎo)數(shù)的情況，由于某些邊值問(wèn)題在角點(diǎn)存在奇異性，并且高階邊界元法需要處理復(fù)雜的高階奇異積分，往往并不能給出滿意的計(jì)算結(jié)果［8］。

針對(duì)上述問(wèn)題，段文洋［9］提出了泰勒展開(kāi)邊界元方法（Taylor expansion boundary element method，TE?BEM）的思想，即對(duì)源-偶混合分布方法中偶極強(qiáng)度（即速度勢(shì)）進(jìn)行泰勒展開(kāi)，并保留一階導(dǎo)數(shù)項(xiàng)（即切向速度），未知數(shù)是速度勢(shì)及其切向?qū)?shù)，如果離散面元數(shù)為N，則對(duì)于三維問(wèn)題，TEBEM未知數(shù)的個(gè)數(shù)為3N。所以需要再引入2N個(gè)代數(shù)方程來(lái)封閉TEBEM的方程組。對(duì)此，本文提出對(duì)格林第三公式中場(chǎng)點(diǎn)，考察其沿著法向趨于邊界時(shí)垂直于法線平面內(nèi)相互正交的2個(gè)方向上的一階導(dǎo)數(shù)特性，發(fā)現(xiàn)其積分擁有同偶極分布類似的特性，即主值系數(shù)為2π。這樣就提供了2組共2N個(gè)代數(shù)方程。因?yàn)樵谏鲜鎏├照归_(kāi)過(guò)程中僅僅展開(kāi)到一階，所以本文方法稱為一階TEBEM。

對(duì)于海洋大型浮體波浪載荷計(jì)算，和解析解計(jì)算結(jié)果的比較表明，本文提出的方法無(wú)論對(duì)于光滑邊界還是邊界突變處的速度勢(shì)和切向速度都具有極高的計(jì)算精度，浮體二階波浪定常力收斂很快。

1 深水浮體二階平均波漂力的定解問(wèn)題

假設(shè)流體無(wú)粘，不可壓縮且流動(dòng)無(wú)旋。本文討論深水浮體水動(dòng)力的頻域解，右手直角坐標(biāo)系原點(diǎn)位于靜水面，z軸垂直向上，流場(chǎng)速度勢(shì)可寫(xiě)為Φ( x，y，z，t )＝Re[ φ( x，y，z)e-iωt]，φ滿足拉普拉斯方程，在靜水面上滿足線性自由面條件。

從總的速度勢(shì)中分出入射波勢(shì)即：φ＝φ0+φP。已知入射波速度勢(shì)表達(dá)式為

式中：ω、a、k0、β分別為入射波浪的角頻率、波幅、波數(shù)及入射波浪向角，g為重力加速度。

根據(jù)疊加原理［4］，擾動(dòng)勢(shì)可分解為7個(gè)部分：

式中：vj為浮體6個(gè)自由度的運(yùn)動(dòng)速度，φj（j＝1，2，…，6）為物體單位速度運(yùn)動(dòng)產(chǎn)生的輻射勢(shì)，φ7為繞射勢(shì)。

擾動(dòng)速度勢(shì)φj（j＝1，2，…，7）定解問(wèn)題為

一是按照“下管一級(jí)”的原則，部屬各級(jí)預(yù)算單位負(fù)有對(duì)本級(jí)和下級(jí)預(yù)算單位動(dòng)態(tài)監(jiān)控的職責(zé)。但各級(jí)預(yù)算單位動(dòng)態(tài)監(jiān)控開(kāi)展工作不夠均衡，部分單位未開(kāi)展動(dòng)態(tài)監(jiān)控工作。二是對(duì)在監(jiān)控中下發(fā)的疑點(diǎn)問(wèn)題，有關(guān)單位的反饋、整改不夠及時(shí)，影響了動(dòng)態(tài)監(jiān)控應(yīng)取得的成效。

一旦求得φ及?φ之后，作用在浮體上的一階和二階水動(dòng)力可根據(jù)式（4）、（5）求得。下面給出浮體的一階水動(dòng)力和二階平均波漂力表述為平均濕表面積分的近場(chǎng)公式：

式中：浮體平動(dòng)位移為T(mén)＝Re( ηTe-iωt)；ηT＝(η1，η2，η3)；角位移為Ω＝Re( ηRe-iωt)；ηR＝ (η，η，η)；總位移為X＝Re(χe-iωt)，χ＝

456(χ1，χ2，χ3)＝ηT+ηR×r；其中波面升高表達(dá)式為：ζ3＝iωφ／g；Awp為水線面面積；ρ表示流體密度；xf是水線面漂心縱向位置坐標(biāo)；r為位置矢量；n表示物面法向，其正方向指向流域外部；wl為浮體平均位置的水線；上標(biāo)＂?＂表示共軛函數(shù)。

而遠(yuǎn)場(chǎng)公式［4］的收斂性更高。其x、y方向表達(dá)式為

式中：υ＝ω2／g。K為柯欽函數(shù)，基于分布源方法的柯欽函數(shù)表達(dá)式可參見(jiàn)文獻(xiàn)［4］。而基于源-偶混合分布方法的柯欽函數(shù)表達(dá)式在下節(jié)給出。

2 一階泰勒展開(kāi)邊界元方法

在流場(chǎng)邊界上運(yùn)用格林第三公式，利用脈動(dòng)源興波格林函數(shù)可以獲得在浮體平均濕表面上表示的流場(chǎng)速度勢(shì)邊界積分方程［4］：

當(dāng)場(chǎng)點(diǎn)沿著法線方向趨近物面時(shí)，可得到基于源-偶混合分布的邊界積分方程，即

為了補(bǔ)充一階TEBEM方法的方程組，對(duì)格林第三公式中場(chǎng)點(diǎn)，考察其沿著法向趨于邊界時(shí)垂直于法線平面內(nèi)相互正交的2個(gè)方向上的一階導(dǎo)數(shù)?？傻玫饺缦路匠探M：

注意方程（9）中的求導(dǎo)為關(guān)于場(chǎng)點(diǎn)處2相互正交方向的切向?qū)?shù)。式中x-和y-為場(chǎng)點(diǎn)所在面元局部坐標(biāo)下的坐標(biāo)，如圖1所示。

將積分邊界S劃分為N個(gè)面元，與傳統(tǒng)的常值面元法不同，對(duì)方程（8）、（9）中的偶極強(qiáng)度在積分邊界劃分后的面元局部坐標(biāo)系下（見(jiàn)圖1）進(jìn)行泰勒展開(kāi)，并保留一階導(dǎo)數(shù)項(xiàng)，即

上式各矩陣中元素表達(dá)式：k＝1，2，3，m＝1，2，3

式中：上標(biāo)i和j表示面元編號(hào)。

一旦利用源-偶混合分布方法求解得到浮體濕表面的速度勢(shì)場(chǎng)后，可借助于柯欽函數(shù)，利用近場(chǎng)速度勢(shì)來(lái)表示無(wú)窮遠(yuǎn)處流場(chǎng)質(zhì)點(diǎn)的速度勢(shì)。再利用遠(yuǎn)場(chǎng)公式（6）可計(jì)算浮體平均波漂載荷。遠(yuǎn)場(chǎng)擾動(dòng)速度勢(shì)表達(dá)式為

式中：r為柱坐標(biāo)系下的徑向坐標(biāo)。

圖1 局部坐標(biāo)系Fig.1 The schematic of the local coordinate

3 用解析解對(duì)TEBEM精度的驗(yàn)證

為了驗(yàn)證TEBEM對(duì)于帶有拐角浮體水動(dòng)力的計(jì)算精度。選擇海洋工程常用的有限吃水圓柱型浮體。對(duì)于無(wú)限水深情況，沒(méi)有找到有限吃水圓柱水動(dòng)力問(wèn)題的解析解。為此，考察有限水深情況下有限吃水圓柱繞射問(wèn)題的解析解［10］，分析當(dāng)吃水增大，直至計(jì)算結(jié)果不受水深變化影響的速度勢(shì)和其導(dǎo)數(shù)穩(wěn)定值代替無(wú)限水深圓柱繞射勢(shì)及其導(dǎo)數(shù)。截?cái)鄨A柱半徑為R，吃水為T(mén)＝1.25R。

選取距離柱面頂部d＝T／50處網(wǎng)格中心點(diǎn)作為研究對(duì)象。利用一階TEBEM方法和分布源方法計(jì)算無(wú)限水深圓柱繞射勢(shì)及其誘導(dǎo)速度。圖2和圖3給出了各點(diǎn)在k0R＝3.26時(shí)z方向無(wú)因次繞射誘導(dǎo)速度實(shí)部φDzc和虛部φDzs。利用入射波相速度做無(wú)因次化處理。從圖中可以看出，一階TEBEM方法求解的z方向的切向誘導(dǎo)速度精度明顯比分布源方法高。

圖2 k0R＝3.26時(shí)，選取點(diǎn)誘導(dǎo)速度實(shí)部Fig.2 Real part of the diffraction induced velocity at selected points，k0R＝3.26

圖3 k0R＝3.26時(shí)，選取點(diǎn)誘導(dǎo)速度虛部Fig.3 Imaginary part of the diffraction induced veloci?ty at selected points，k0R＝3.26

圖4 圓柱的平均漂移力Fig.4 The diffraction mean drift force of cylinder

圖4 給出了圓柱繞射平均漂移力的解析解，分布源方法和一階泰勒展開(kāi)邊界元方法基于近場(chǎng)公式的計(jì)算值，及基于源偶混合分布和分布源方法的遠(yuǎn)場(chǎng)積分公式計(jì)算值。為了方便比較，平均漂移力結(jié)果利用ρga2R作無(wú)因次化處理。圖5給出了4種方法的計(jì)算結(jié)果的相對(duì)誤差。從圖中可以看出，一階TEBEM方法的相對(duì)誤差在整個(gè)波長(zhǎng)范圍內(nèi)控制在4%以內(nèi)。在高頻段，可以看出一階TEBEM方法的相對(duì)誤差明顯小于分布源方法。而2種遠(yuǎn)場(chǎng)公式的計(jì)算結(jié)果相同。與解析解吻合，精度較高。在頻域方法中，會(huì)存在不規(guī)則頻率的現(xiàn)象，本文利用的是擴(kuò)展邊界積分方程的方法消除不規(guī)則頻率［11］，因此計(jì)算結(jié)果中沒(méi)有看到不規(guī)則頻率現(xiàn)象。

圖5 圓柱平均漂移力相對(duì)誤差Fig.5 Relative error of the diffraction mean drift force

圖6 水線積分項(xiàng)貢獻(xiàn)Fig.6 Contribution of the waterline integral

圖7 速度平方項(xiàng)貢獻(xiàn)Fig.7 Contribution of the square velocity item integral

4 實(shí)船二階定常波浪力結(jié)果與分析

本文以KVLCC2船為模型進(jìn)行實(shí)船驗(yàn)證。入射波浪向角為艏斜浪β＝165°。圖9、10給出了KVLCC2在艏斜浪方向入射波浪作用下，一階TEBEM方法和分布源方法計(jì)算得到的縱向平均漂移力結(jié)果。為了方便比較，利用ρga2B2／L作無(wú)因次化處理。圖中L為船長(zhǎng)，B為船寬，λ為波長(zhǎng)。從圖9、10中可以看出，分布源方法與遠(yuǎn)場(chǎng)公式在高頻段收斂于不同的值。2種方法存在一定的誤差。相比較一階TEBEM方法與遠(yuǎn)場(chǎng)公式結(jié)果更為吻合。圖11、12給出了相同工況下，2種方法計(jì)算的KVLCC2橫向平均漂移力結(jié)果。從圖中可看出，分布源方法的計(jì)算結(jié)果在船長(zhǎng)波長(zhǎng)比大于2后發(fā)散，而一階TEBEM方法很快收斂。

圖9 基于一階TEBEM方法計(jì)算的縱向平均漂移力（β＝165°）Fig.9 The surge mean drift force calculated by the 1storder TEBEM method at β＝165°

圖10 基于分布源方法計(jì)算的縱向平均漂移力（β＝165°）Fig.10 The surge mean drift force calculated by the source distribution method at β＝165°

圖11 基于一階TEBEM方法計(jì)算的橫向平均漂移力（β＝165°）Fig.11 The transverse mean drift force calculated by the 1storder TEBEM method at β＝165°

圖12 基于分布源方法計(jì)算的橫向平均漂移力（β＝165°）Fig.12 The transverse mean drift force calculated by the source distribution method at β＝165°

5 結(jié)論

本文提出了一種邊界積分方程新的求解方法，稱之為泰勒展開(kāi)邊界元方法。利用該方法可精確求解邊界突變處的誘導(dǎo)速度，并解決了傳統(tǒng)高階邊界元方法處理奇異積分困難的問(wèn)題，以無(wú)限水深圓柱波浪繞射及6自由運(yùn)動(dòng)的KVLCC2肥大船舶來(lái)驗(yàn)證一階TE?BEM方法的精度，結(jié)論如下：

1）通過(guò)與解析解的比較可知，一階TEBEM方法可解決物面突變處誘導(dǎo)速度精度低的問(wèn)題。與傳統(tǒng)分布源方法相比，即使在相同的計(jì)算量下，一階TEBEM方法的計(jì)算精度仍比分布源方法要高。

2）由于高精度的誘導(dǎo)速度求解，一階TEBEM方法計(jì)算的浮體二階平均波浪漂移力精度更高，尤其體現(xiàn)在高頻段。

本文只是對(duì)偶極強(qiáng)度進(jìn)行一階泰勒展開(kāi)，但是在解決實(shí)際工程問(wèn)題時(shí)，比如船舶波浪增阻問(wèn)題，需要精確求解速度勢(shì)的二階導(dǎo)數(shù)，因此將一階TEBEM方法拓展至二階是未來(lái)需進(jìn)行的工作。

［1］HESS J L，SMITH A M O.Calculation of non?lifting poten?tial flow about arbitrary three?dimensional bodies［J］.Jour?nal of Ship Research，1964，8：22?44.

［2］NEWMAN J N，LEE C H.First and second order wave effects on a submerged spheroid［J］.Journal of Ship Re?search，1991，35（3）：183?190.

［3］CHEN X B.Middle?field formulation for the computation of wave?drift loads［J］.Journal of Engineering Mathematics，2007（59）：61?82.

［4］戴遺山，段文洋.船舶在波浪中運(yùn)動(dòng)的勢(shì)流理論［M］.北京：國(guó)防工業(yè)出版社，2008：28?35.DAI Yishan，DUAN Wenyang.Potential flow theory of ship motions in waves［M］.Beijing：National Defense Industry Press，2008：28?35.

［5］CHOI Y R，HONG S Y，CHOI H S.An analysis of second?order wave forces on floating bodies by using a higher?order boundary element method［J］.Ocean Engineering，2000，28：117?138.

［6］TENG B.Second order wave action on 3D floating body［J］.Journal of Hydrodynamics A，1995，10（3）：316?327.

［7］柏威.非線性波浪與任意三維物體的相互作用［D］.大連：大連理工大學(xué)，2001：83?96.BAI Wei.Non?linear wave interaction with arbitrary 3?D bodies［D］.Dalian：Dalian University of Technology，2001：83?96.

［8］SHAO Y L，ALTINSON O M.Use of body fixed coordinate system in analysis of weakly nonlinear wave body problems［J］.Applied Ocean Research，2010，32：20?33.

［9］DUAN W Y.Taylor expansion boundary element method for floating body hydrodynamics［C］／／Proc 27th Intl Workshop on Water Waves and Floating Bodies.Copenhagen，Den?mark，2012.

［10］CHEN X B.Analytical evaluation of the second?order loads upon a tension leg platform column［Z］.IFP report No.38 739.Rueil?Malmaison，F(xiàn)rance，1991.

［11］ZHU X M.Irregular frequency removal from the boundary integral equation for the wave?body problem［D］.Boston：Massachusetts Institute of Technology，1994：31?35.

Calculation of second?order mean drift loads for the deepwater floating body based on the Taylor expansion boundary element method

DUAN Wenyang，CHEN Jikang，ZHAO Binbin
（College of Shipbuilding Engineering，Harbin Engineering University，Harbin 150001，China）

Taylor expansion boundary element method（1storder TEBEM）is developed for improving the solution accuracy of the tangential induced velocity at the unsmooth corners of the floating body in the traditional low order panel method.This method is based on the framework of the constant panel method to discrete the boundary integra?tion equation（BIE），which mainly applies the Taylor expansion formula to the dipole strength of the BIE on each element and reserves the first order derivatives which are exactly the tangential velocity on the element.Finally it u?tilizes the corresponding tangential derivatives of two tangential directions with respect to the field points on the boundary to close the equations.Numerical results showed that the 1storder TEBEM can accurately solve the tan?gential velocity at the corner of floating body and the height of wave surface of waterline.Compared with constant panel methods，it is found that TEBEM can get higher accuracy and convergence results of the second order wave drift loads.

second order mean wave drift loads；Taylor expansion；boundary element method；deepwater floating body；Kochin function；near?field formulation；far?field formulation

10.3969／j.issn.1006?7043.201406006

http：／／www.cnki.net／kcms／detail／23.1390.U.20150109.1453.001.html

O352

A

1006?7043（2015）03?0302?05

2014?06?05.網(wǎng)絡(luò)出版時(shí)間：2015?01?09.基金項(xiàng)目：國(guó)家自然科學(xué)基金資助項(xiàng)目（11272097）.

段文洋（1967?），男，教授，博士生導(dǎo)師.

段文洋，E?mail：duanwenyang＠hrbeu.edu.cn.