常數(shù)分布Rankine源法與二階繞射問題精度研究

徐 剛，段文洋

(1.江蘇科技大學(xué) 船舶與海洋工程學(xué)院，江蘇 鎮(zhèn)江 212003;2.哈爾濱工程大學(xué) 船舶工程學(xué)院，黑龍江 哈爾濱 150001)

雖然很多學(xué)者致力于發(fā)展完全非線性的數(shù)值技術(shù)，但是由于深海平臺的工作環(huán)境及其自身特點(diǎn)，波浪對結(jié)構(gòu)物的二階作用不可忽視，發(fā)展一種有效的處理二階問題的數(shù)值方法是非常有必要的.

在水動力學(xué)領(lǐng)域求解勢流問題有2種常用的方法，即有限元法(finite elementmethod，F(xiàn)EM)和邊界元法(boundary elementmethod，BEM).但是網(wǎng)格生成困難是有限元法的最大缺點(diǎn)，高質(zhì)量的網(wǎng)格劃分通常需要依據(jù)物體和波浪的運(yùn)動［1］.因此，本文采用邊界元法.

在應(yīng)用邊界元法求解線性時域問題時，通常也有2種方法:1)滿足自由面條件的時域格林函數(shù)法，但其含有時間記憶項，在長時歷的數(shù)值模擬時，每一步的計算所需要的時間越來越長，導(dǎo)致計算效率非常低;2)Rankine源法，這里的格林函數(shù)不滿足自由面條件，因此在自由面上也需要分布源匯.Rankine源法的有以下3個優(yōu)點(diǎn):沒有記憶項［2］、可以直接得到自由面上速度勢分布、容易拓展到非線性問題.因此，本文采用基于Rankine源的邊界元法在時域內(nèi)求解非線性的波物相互作用問題.

基于Rankine源的邊界元法需要在距離物體有限遠(yuǎn)的距離處截斷流域，然而在截斷面(人工邊界)上還沒有準(zhǔn)確的完全無反射的輻射邊界條件存在.目前，Orlanski法［3］和阻尼區(qū)法［4］已經(jīng)廣泛應(yīng)用于模擬開路邊界問題.而在開路邊界附近網(wǎng)格不能足夠密時，Orlanski法可能會產(chǎn)生不正確的相速度.Clément［5］提出了一種耦合方法(piston-beach hybrid absorber)來消除波的反射;Boo［6］也提出了一種耦合的方法(an absorbing beach and the stretching technique［7］)來處理不規(guī)則波問題;Wang 等［2］通過結(jié)合阻尼區(qū)法和Sommerfeld-Orlanski法作為輻射條件，研究了多個柱體Trappingmode問題.但阻尼區(qū)法的消波效率取決于阻尼區(qū)寬度與被吸收波長之比，阻尼區(qū)越寬消波效率也就越高，這會導(dǎo)致自由面上網(wǎng)格數(shù)量的急劇增多，降低計算效率，尤其會表現(xiàn)在非線性的不規(guī)則波問題上.

為了尋找一種普適的輻射條件來處理非線性不規(guī)則波問題，本文將采用多次透射公式［8］(multitransmitting formula，MTF)作為遠(yuǎn)方輻射條件.在MTF法中，需引入人工波速來取代波的相速度，通常情況下，沒有必要取人工波速等于物理波速.同時，在自由面上采用了積分格式的自由面條件［9］(integral formof free surface boundary condition，IFBC)，通過此積分格式可以得到自由面上的速度勢分布，它與傳統(tǒng)的有限差分法相比數(shù)值上更穩(wěn)定，目前該方法已經(jīng)成功應(yīng)用于求解線性輻射［10］和不規(guī)則波繞射問題［11］.

為了檢驗(yàn)本文方法在二階問題中的適用性，同時也為了研究文中基于常值分布簡單格林函數(shù)法的數(shù)值精度，文中將對圓球繞流和直立圓柱樁的二階繞射問題進(jìn)行了研究，分析了本文結(jié)果與相關(guān)參考文獻(xiàn)中結(jié)果存在差異的主要原因.

1 控制方程

建立笛卡爾坐標(biāo)系作為參考坐標(biāo)系，將xoy定義在靜水面上，z軸垂直于靜水面向上，如圖1所示.物面用SH表示，其單位法向量n指向物體內(nèi)部.人工輻射邊界SC將流場分為內(nèi)域和外域2個部分.內(nèi)域的邊界S=SF+SH+SB+SC.

圖1 坐標(biāo)系和計算域定義圖Fig.1 Definition sketch

假定浮體是剛性不透水的，海底在z=-h的水平平面上，而且也不透水.假定流體是無粘、無旋不可壓縮的理想流體，則流體的運(yùn)動可以通過速度勢φ來表示，該速度勢在流體域Ωf內(nèi)滿足拉普拉斯方程.而且速度勢在海底SB和物面SH上滿足剛性不透水邊界條件;在z=η處的瞬時自由面SIF上滿足運(yùn)動學(xué)和動力學(xué)邊界條件，其中η表示自由面波高;在距離物體有限遠(yuǎn)處的人工邊界SC上滿足合適的人工輻射邊界條件.

考慮到自由面條件到二階，對于k階水波繞射問題［12］(k依次為1和2)，各階繞射勢φ(k)D在流體域Ωf內(nèi)滿足定解條件:

在海底SB、物面SH和靜水面SF上的滿足的邊界條件分別為

式中:下標(biāo)I和D分別表示入射勢和繞射勢，其中在

人工邊界SC上滿足二階多次透射邊界條件，其具體形式將在后文討論.

2 多次透射邊界條件

廖振鵬［8］應(yīng)用時空外推法給出了外傳波局部無反射邊界條件的一般表達(dá)式MTF，該方法是從波的傳播特性出發(fā)，將人工輻射邊界SC上的某點(diǎn)在t時刻的速度勢，用鄰近流體域Ωf內(nèi)的內(nèi)點(diǎn)在t時刻以前的速度勢來表示.

圖2 多次透射邊界條件Fig.2 Radiation condition on the artificial boundary

多次透射邊界條件如圖2.x軸垂直于人工輻射邊界SC并指向人工邊界的外部區(qū)域，將x軸與人工邊界的交點(diǎn)0設(shè)為人工邊界SC上的外傳點(diǎn)，j點(diǎn)是由0點(diǎn)沿其法線方向指向邊界內(nèi)部區(qū)域的點(diǎn)，j點(diǎn)與0點(diǎn)之間的距離為jCaΔt，這里Ca為人工波速(可不等于物理波速Cx，具體取值范圍可參考文獻(xiàn)［10-11］)，由此可以得到N階MTF的形式:

式中:整數(shù)p表示時間，N表示MTF階數(shù)，其中二項式系數(shù)CNj具有如下形式:

通過計算實(shí)踐［10-11］表明，二階MTF已能解決問題，因此式(7)有如下形式:

為了防止數(shù)值實(shí)現(xiàn)過程中出現(xiàn)低頻失穩(wěn)的現(xiàn)象，需在MTF中引入附加因子γ2，則式(9)可改為

式中:速度勢φ的下標(biāo)0表示人工輻射邊界上的0點(diǎn)，1和2分別表示在流體域內(nèi)部沿0點(diǎn)的法線反向分別前進(jìn)CaΔt和2CaΔt的內(nèi)部點(diǎn)1和2.

3 自由面條件積分格式

通常自由面上的速度勢φ都采用有限差分法來求解，但這種處理方法相對比較復(fù)雜，同時還會出現(xiàn)數(shù)值誤差的累積，易導(dǎo)致計算結(jié)果發(fā)散，數(shù)值穩(wěn)定性比較差，尤其表現(xiàn)在高階導(dǎo)數(shù)問題上.在文中，將采用一種積分格式的自由面條件IFBC［9］，通過此積分格式可以得到自由面上一階和二階繞射勢的表達(dá)式，而且這種積分格式的自由面條件穩(wěn)定性相當(dāng)高:

4 邊界積分方程及離散

利用基于三維格林公式的邊界積分方程可以得到控制域表面上各階繞射勢:

式中:p(x，y，z)為場點(diǎn)，q(ξ，η，ζ)為源點(diǎn)，G(p，q)為格林函數(shù).由于海底是水平的，可以利用海底鏡像.對于后續(xù)討論的圓柱繞射問題格林函數(shù)G的形式如下:

式中:r表示場點(diǎn)與源點(diǎn)之間的距離，點(diǎn)q＇是點(diǎn)q關(guān)于海底的鏡像，因此距離r可用下面2個關(guān)系式來表示:

本文求解積分方程式(14)的數(shù)值方法是低階面元法，在各控制面上劃分單元，且在相應(yīng)單元上速度勢及其法向?qū)?shù)均看成常數(shù)，邊界條件在單元的中心點(diǎn)滿足.式(14)的離散形式如下

其中，人工輻射邊界上的單元速度勢可以用人工邊界條件MTF獲得，具體計算關(guān)系式如下:

式中:nSH、nSF和nSC分別表示控制面SH、SF和SC上的單元數(shù);φDjm表示在時間t=mΔt時第j塊面元上的繞射勢是其法向?qū)?shù);φDj'，m-1和 Dj″，m-2表示在人工邊界處的內(nèi)點(diǎn)j'和j″分別在t=(m-1)Δt和t=(m-2)Δt時的繞射勢;方程中 Sij和Dij為面元之間的影響系數(shù).通過式(18)就可以得到整個控制面上的法向速度，再由式(13)、(14)就可以得到每個單元上的速度分布.

5 水動力計算公式

作用在物體上的水動力F可以直接通過對瞬時物面上的壓力積分獲得，對于二階繞射問題可以寫為［12］

水動力又可分解成如下形式:

式中:F(2)=F(21)+F(22)，w0表示物面與靜水面的交線，F(xiàn)(1)表示一階力，F(xiàn)(2)表示二階倍頻力，F(xiàn)(21)表示由一階勢引起的二階倍頻力，F(xiàn)(22)表示由二階勢引起的二階倍頻力，F(xiàn)0表示二階定常力，符號“﹤﹥”表示時間平均.

6 數(shù)值結(jié)果及討論

由于現(xiàn)有計算二階水動力的方法所得到的結(jié)果之間存在很大的離散度［12-16］，同時為了研究本文方法在二階問題中的適用性，本文以圓球繞流和直立圓柱樁的二階繞射問題為例來檢驗(yàn)本文方法的數(shù)值精度，并通過對數(shù)值結(jié)果的比較和分析來說明現(xiàn)有結(jié)果存在較大離散度的原因.

6.1 圓球繞流求解兩種分布源模型

首先以三維圓球繞流為例，研究不同分布源法［17］的數(shù)值精度，分布源模型如式(13)、(14).對于圓球繞流問題，需要將式(13)、(14)中的速度勢下標(biāo)D去掉，同時不考慮時間因數(shù):

式中:源點(diǎn)(ξ，η，-ζ)是源點(diǎn)(ξ，η，ζ)關(guān)于對稱面xoy的鏡像點(diǎn).

圖3 計算域和網(wǎng)格模型Fig.3 Control domain and mesh of themodel

圖3為圓球繞流模型，圖3(a)為模型1，沒有考慮圓球關(guān)于xoy平面的對稱性(在對稱面上也需要劃分網(wǎng)格，即半球與對稱面交界處網(wǎng)格法向量存在突變，格林函數(shù)G=1/r1);圖3(b)為模型2，利用圓球關(guān)于xoy平面的對稱性(對稱面上不需要劃分網(wǎng)格，即在半球上不存在法向量突變，格林函數(shù)G=1/r1+1/r3)，這里的格林函數(shù)滿足式(28)，其中r1和 r3由式(16)、(29)確定.利用式(13)、(14)可以得到面元中心點(diǎn)上的速度勢以及x、y和z方向上的速度 φx，φy和 φz.

圖4～6為圓球不同位置處(在xoz平面上，x＞0)的速度勢、速度以及其與解析解［18］之間的相對誤差曲線.圖4為xoz面上圓球速度勢隨坐標(biāo)z的變化關(guān)系(x＞0);圖5為xoz面上圓球x方向速度隨坐標(biāo)z的變化關(guān)系(x＞0);圖6為xoz面上圓球z方向速度隨坐標(biāo)z的變化關(guān)系(x＞0).從圖中可以看出，使用模型2的結(jié)果明顯優(yōu)于模型1，得到速度勢和速度相對誤差隨水深z變化不大，而模型1得到的x和z方向的速度隨水深變化非常劇烈，尤其表現(xiàn)在對稱面和圓球交界處(即法向突變處)，這就說明分布源法在控制形狀突變處得到的面元切向速度是不準(zhǔn)確的，存在較大的誤差.

圖4 速度勢及其相對誤差Fig.4 The potential and its relative error

圖5 x方向速度及其相對誤差Fig.5 The velocity in x direction and its relative error

圖6 z方向速度及其相對誤差Fig.6 The velocity in z direction and its relative error

6.2 直立圓柱樁繞射

為了進(jìn)一步研究常數(shù)分布簡單格林函數(shù)法的數(shù)值精度，本文在時域內(nèi)模擬了直立圓柱樁的二階繞射問題.由于海底鏡像，這里的水線處就相當(dāng)于模型1(面元存在法向突變);圓柱與水底交界處就相當(dāng)于模型2.本文的結(jié)果將與文獻(xiàn)［17］中的頻域解析解進(jìn)行比較.設(shè)入射波為周期性的Stokes波［19］，水深為h，柱半徑為a.計算中取人工波速Ca等于一階入射波的相速度Cx，其滿足穩(wěn)定實(shí)現(xiàn)MTF的基本條件，即滿足 Ca的取值范圍［10-11］.

為了避免初始效應(yīng)對數(shù)值結(jié)果的影響，在物面條件式(3)中使用平滑函數(shù)M(t):

初始時刻對應(yīng)于全流場繞射勢和繞射波高為零，平滑時間Tm取3倍的一階入射波周期T.

首先本文通過源-耦混合分布模型求得整個控制面上繞射勢φ及其法向?qū)?shù)φn的值(如無特殊說明后續(xù)圖中關(guān)于一階繞射勢的相關(guān)量均略去了上標(biāo)(1)和下標(biāo)D)，并將物面和自由面上的結(jié)果與頻域解析解進(jìn)行比較.為了便于圖形顯示，本文一般只給出了第10個周期的時歷曲線，同時在圖中標(biāo)明了該點(diǎn)所處的位置.從圖7、8可以看出:無論是速度勢還是其法向?qū)?shù)都與解析解吻合的非常好，這就說明MTF能夠有效地將一階繞射波傳出人工輻射邊界.除有特殊說明外，以下圖中的量都是采用國際單位，t/T表示無量綱的時間;計算模型:水深h=1m，柱半徑a=1 m;入射波波幅 A=1m，波數(shù) k=1.0m-1.

圖7 圓柱上點(diǎn)(1.0 m，0，0.025 m)和(1.0 m，0，0.475 m)的速度勢及其法向?qū)?shù)Fig.7 The potential and its derivative on the cylinder(1.0 m，0，0.025 m)，(1.0 m，0，0.475 m)

圖8 自由面上點(diǎn)(1.052 m，0，0)和(2.622 m，0，0)的速度勢及其法向?qū)?shù)Fig.8 The potential and its derivative on the free surface(1.052 m，0，0)，(2.622 m，0，0)

其次，通過分布源模型將由混合分布模型所得到的法向速度n用來求解每個面元上的速度φx，φy和φz.圖9分別表示圓柱側(cè)面上水線處的點(diǎn)(1.0 m，0，0.025 m)和水底處的點(diǎn)(1.0 m，0，0.047 5 m)的速度曲線;圖10分別表示自由面上水線處的點(diǎn)(1.052 m，0，0)和物面與遠(yuǎn)方人工邊界之間的一點(diǎn)(2.622 m，0，0)的速度曲線.從圖中可以看出:分布源法得到的物面和自由面上面元的法向速度與解析解吻合的比較好，但是面元的切向速度(圓柱z方向的速度，自由面水平方向的速度)與解析解差別較大，尤其表現(xiàn)在控制域形狀突變處(如:水線以及遠(yuǎn)方控制面與自由面交界處).由于物面和自由面上y=0處有φy=0，所以圖中沒有顯示φy的值.導(dǎo)致面元切向速度誤差比較大的原因在于:分布源法在面元法向突變處的精度較低，這與三維圓球繞流問題的結(jié)論是一致的.

圖9 圓柱上點(diǎn)(1.0 m，0，0.025m)和(1.0 m，0，0.475 m)的速度Fig.9 The velocity on the cylinder(1.0 m，0，0.025 m)，(1.0 m，0，0.475 m)

圖10 自由面上點(diǎn)(1.052 m，0，0)和(2.622 m，0，0)的速度Fig.10 The velocity on the free surface(1.052 m，0，0)，(2.622 m，0，0)

本文進(jìn)一步研究網(wǎng)格密度對分布源精度的影響，圖11是物面上網(wǎng)格加密3倍之后，將物面上同一點(diǎn)的切向速度與先前結(jié)果的比較，除了物面上緊挨著自由面一點(diǎn)切向速度僅相位趨于解析解之外，其余各點(diǎn)的幅值和相位都與解析解更接近，這就說明加密網(wǎng)格之后數(shù)值結(jié)果的總體誤差變小，精度提高，但這并不能從根本上解決分布源法在面元法向突變處精度較低的問題.此外，從計算效率的角度考慮，又不能無限增加網(wǎng)格的數(shù)量，這會導(dǎo)致效率急劇降低，不利于文中方法應(yīng)用到實(shí)際工程問題中去.由于二階繞射勢和水動力的求解與物面和自由面上一階速度密切相關(guān)，因此切向速度存在較大誤差必然會導(dǎo)致二階勢也存在較大誤差.

圖11 圓柱上點(diǎn)(1.0m，0，0.025m)和(1.0m，0，0.475m)z方向的速度Fig.11 The velocity on the cylinder in z direction(1.0 m，0，0.025 m)，(1.0 m，0，0.475 m)

本文為了研究切向速度精度對二階水動力的影響，將文獻(xiàn)［17］中一階繞射勢的解析解代入到二階非齊次自由面條件中，如式(6)所示，得到的二階水動力及其相對誤差如圖12所示，從圖中可以看出隨著頻率的增加相對誤差逐漸變大.但圖13中的時歷曲線卻能說明MTF可以作為遠(yuǎn)方邊界條件求解二階水波問題［20-21］.

圖12 x方向的二階力及其相對誤差Fig.12 Second order forces in x direction and its relative error

圖13 二階水動力的時歷曲線Fig.13 The time histories of second order hydrodynamic force

7 結(jié)論

本文對圓球繞流和直立圓柱樁的二階繞射問題進(jìn)行了研究，得到如下結(jié)論:

1)10個周期的數(shù)值模擬沒有出現(xiàn)發(fā)散現(xiàn)象，結(jié)果穩(wěn)定性非常高，說明MTF可以用于時域二階水波問題的數(shù)值模擬;

2)通過將數(shù)值結(jié)果與解析解進(jìn)行比較，可知常數(shù)分布源法在控制形狀突變處求得的面元切向速度精度比較低;

3)通過將二階繞射結(jié)果與參考文獻(xiàn)中的數(shù)值結(jié)果進(jìn)行比較，得出本文方法計算得到的二階水動力與半解析解存在較大差異的原因.

［1］WU G X，TAYLOR R E.The coupled finite element and boundary element analysis of nonlinear interactions between waves and bodies［J］.Ocean Engineering，2003，30:387-400.

［2］WANG C Z，WU G X.Time domain analysis of second-order wave diffraction by an array of vertical cylinders［J］.Journal of Fluids and Structures，2007，23(4):605-631.

［3］ORLANSKI I.A simple boundary condition for unbounded hyperbolic flows［J］.Journal of Computational Physics，1976，21:251-269.

［4］ISRAELI M，ORSZAG S A.Approximation of radiation boundary conditions［J］.Journal of Computational Physics，1981，41:115-135.

［5］CLéMENT A.Coupling of two absorbing boundary conditions for 2-D time-domain simulations of free surface gravity waves［J］.Journal of Computational Physics，1996，126:139-151.

［6］BOO SY.Linear and nonlinear irregular waves and forces in a numerical wave tank［J］.Ocean Engineering，2002，29:475-493.

［7］FORRISTALL G Z.Irregular wave kinematics froma kinematics boundary condition fit(KBCF)［J］.Applied Ocean Research，1985，7:202-212.

［8］廖振鵬.工程波動理論導(dǎo)論［M］.北京:科學(xué)出版社，2004:141-189.

［9］戴遺山，段文洋.船舶在波浪中運(yùn)動的勢流理論［M］.北京:國防工業(yè)出版社，2008:214-217.

［10］XU Gang，DUANWenyang.Time domain simulation of irregular wave diffraction［C］//Proceedings of the 8th International Conference on Hydrodynamics.Nantes，2008.

［11］XU Gang，DUAN Wenyang.Time domain simulation for waterwave radiation by floating structures(Part A.)［J］.Journal of Marine Science and Application，2008，7:226-235.

［12］ISAACSON M，CHEUNG K F.Time domain second-order wave diffraction in three dimensions［J］.Journal ofWaterway，Port，Coastal，and Ocean Engineering，1992，118(5):496-516.

［13］ISAACSON M，NG JY T，CHEUNG K F.Second-order wave radiation of three-dimensional bodies by time-domain method［J］.International Journal of Offshore and Polar Engineering，1993，3(4):264-272.

［14］G?REN?.On the second-order wave radiation of an oscillating vertical circular cylinder in finite-depth water［J］.Journal of ShipResearch，1996，40(3):224-234.

［15］滕斌，柏威.三維浮體二階輻射問題的時域模擬［J］.中國造船，2002，43:289-298.TENG Bin，BAIWei.Simulation of second-order radiation of 3D bodies in time domain［J］.Shipbuilding of China，2002，43:289-298.

［16］TAYLOR R E，HUNG S M.Second order diffraction forces on a vertical cylinder［J］.Applied Ocean Research，1987，9(1):19-30.

［17］戴遺山.艦船在波浪中運(yùn)動的頻域與時域勢流理論［M］.北京:國防工業(yè)出版社，1998:44-54，116-117.

［18］吳望一.流體力學(xué)［M］.北京:北京大學(xué)出版社，1983:124-127.

［19］鄒志利.水波理論及其應(yīng)用［M］.北京:科學(xué)出版社，2005:25-34.

［20］徐剛，段文洋.半潛柱體垂蕩運(yùn)動二階水動力分析［J］.哈爾濱工程大學(xué)學(xué)報，2010，31(4):414-420.XU Gang，DUAN Wenyang.Second-order hydrodynamic analysis of surface-piercing circular cylinders undergoing heavemotion［J］.Journal of Harbin Engineering University，2010，31(4):414-420.

［21］徐剛.不規(guī)則波中浮體二階水動力時域數(shù)值模擬［D］.哈爾濱:哈爾濱工程大學(xué)，2010:117-135.XU Gang.Time-domain simulation of second-order hydrodynamic force on floating bodies in irregular waves［D］.Harbin:Harbin Engineering University，2010:117-135.