韓翔希+馮志強(qiáng)+邱昂+符妃+范少濤

摘 要：本研究采用CFD方法建立三維波浪水池，并在所模擬的數(shù)值波浪水池中計(jì)算Wigley船型頂浪航行在規(guī)則波時的運(yùn)動響應(yīng)。采用RANS方程求解非定常的不可壓縮粘性流體，并采用VOF方法對自由面進(jìn)行動態(tài)模擬。通過編寫UDF設(shè)置邊界入口速度、波高及在波浪水池尾部1～2倍波長區(qū)域消波。采用動網(wǎng)格技術(shù)結(jié)合六自由度求解器模擬船舶的縱搖與垂蕩運(yùn)動。將采用CFD方法模擬得到的Wigley船運(yùn)動響應(yīng)與通過邊界元法所得到的勢流結(jié)果進(jìn)行對比，發(fā)現(xiàn)兩者吻合良好。本文所探索的方法可有效計(jì)算船舶在規(guī)則波中的運(yùn)動響應(yīng)，便捷高效，易于測量和控制，并能細(xì)致描述船體周圍的流場情況。同時，也可為模擬分析其他海洋結(jié)構(gòu)物在波浪中的運(yùn)動響應(yīng)具有參考價值，豐富與擴(kuò)展數(shù)值波浪水池的應(yīng)用。

關(guān)鍵詞：數(shù)值波浪水池；船舶運(yùn)動響應(yīng)；VOF；六自由度求解器；動網(wǎng)格

中圖分類號：U661.3 文獻(xiàn)標(biāo)識碼：A

Simulation of Ship Motion Response Based on CFD

HAN Xiangxi1，2， FENG Zhiqiang1，2， QIU Ang3， FU Fei1，2 FAN Shaotao4

（ 1. Guangxi Engineering Technology Research Center of Marine Digital Design and Advanced Manufacturing， Qinzhou 535099；

2. Qinzhou University， Qinzhou 535099； 3. Guangdong Shipping Science Research Institute， Guangzhou 510000；

4. China Energy Engineering Group Guangdong Electric Power Design Institute Co.， Ltd. Guangzhou 510663 ）

Abstract： This paper establishes a three-dimensional （3D） numerical wave tank based on the theory of viscous flow to simulate the unsteady motion response of a Wigley ship navigating in regular heading waves. The governing equations， Reynolds Averaged Navier-Stokes and continuity equations are discretized by finite volume method， A six-degrees-of-freedom （6 DOF） solver is employed to predict the motions of the ship， and volume of fluid （VOF） method is adopted to capture the nonlinear free surface by writing user-defined functions （UDF） to set inlet boundary. The outgoing waves are dissipated inside an artificial damping zone located at the rear part （about 1-2 wave lengths） of the wave tank. The numerical simulation results show good agreement with the theoretical results. This research can be used to further analyze and predict hydrodynamic performance of the ship and marine floating structures in waves and help to extend the applications of the numerical wave tank.

Key words： Numerical wave tank； Ship motion response； Volume of fluid； Six DOF Solver； Dynamic mesh

1 引言

近年來，由于CFD技術(shù)的發(fā)展，基于CFD方法的數(shù)值造波技術(shù)逐漸成為研究熱點(diǎn)。Harlow等人提出的MAC方法和Hirt等提出的VOF方法，推動了粘性自由表面數(shù)值計(jì)算技術(shù)的快速發(fā)展。此后，許多學(xué)者在粘性自由液面的數(shù)值水池上做了大量工作。Wang采用VOF方法模擬出二維數(shù)值波浪水池，并開展波浪對近海平臺底部沖擊過程的研究；齊鵬和王永學(xué)在前人二維波浪水池的基礎(chǔ)上，采用VOF的方法，根據(jù)造波機(jī)的原理建立三維數(shù)值波浪水池；Xing-Kaeding等人基于RANS方法模擬計(jì)算船舶在波浪中的運(yùn)動響應(yīng)，自由面同樣采用VOF方法進(jìn)行處理。

二維數(shù)值水槽的建立及分析已經(jīng)相當(dāng)成熟，但應(yīng)用到工程項(xiàng)目中帶有局限性，所以建立三維波浪數(shù)值模型分析計(jì)算波浪與結(jié)構(gòu)物相互作用具有特別重要的意義。另外，目前計(jì)算波浪作用在結(jié)構(gòu)物上的載荷大多通過線性勢流理論解決，但線性勢流理論有明顯的局限性，很難考慮粘性和非線性的影響，如果能夠模擬高階非線性波且結(jié)果可靠，則可用數(shù)值波浪水池模擬各種波況下、特別是極端波況下的波浪運(yùn)動特性。本研究基于CFD建立三維數(shù)值波浪水池，自由液面采用VOF處理，通過編寫UDF設(shè)置邊界入口速度及動量方程中添加源項(xiàng)實(shí)現(xiàn)阻尼消波，采用動網(wǎng)格技術(shù)結(jié)合六自由度求解器模擬船舶的縱搖與垂蕩運(yùn)動，模擬了規(guī)則波頂浪中三維船體縱搖與垂蕩運(yùn)動的計(jì)算，并與勢流理論邊界元法得到的結(jié)果相比較，兩者吻合很好。通過試驗(yàn)，本方法能準(zhǔn)確預(yù)測規(guī)則波頂浪中船舶的縱搖與垂蕩運(yùn)動響應(yīng)，能細(xì)致描述船舶周圍的流場。本方法也可為模擬分析其他的海洋結(jié)構(gòu)物在波浪中的運(yùn)動響應(yīng)具有參考價值，豐富與擴(kuò)展數(shù)值波浪水池的應(yīng)用。endprint

2 數(shù)學(xué)模型及數(shù)值方法

2.1 控制方程

本文采用兩相不可壓的RANS方程求解非定常不可壓粘性流體，RANS方程采用 SST k-ω湍流模型封閉，自由液面采用VOF法處理，消波區(qū)域采用在動量方程中添加源項(xiàng)的辦法進(jìn)行消波。

連續(xù)方程： （1）

RANS方程： （2）

式中：U為速度矢量；流體密度定義為 ，其中體積分?jǐn)?shù)aq表示單元內(nèi)第q相流體占的體積與該單元的體積之比，本文采用氣液兩相模擬q=1，2，并且

；μ為相體積分?jǐn)?shù)平均的動力粘性系數(shù)。

Menter提出了SST k-ω 兩方程模型，它在近壁面保留了原始k-ω的模型，在遠(yuǎn)離壁面的地方應(yīng)用了k-ω模型，其k方程和ω方程如下：

（3）

（4）

自由液面采用VOF法處理，aq滿足方程：

（5）

2.2 數(shù)值造波和消波方法

本文采用在邊界設(shè)置入口速度的方法造波，需要編寫UDF控制入口速度以及波高。

根據(jù)線性波理論，小振幅波的速度場為：

（6）

（7）

小振幅波的自由液面高度

（8）

式中：U0為船舶航速；A是波幅；k為波數(shù)；ω是波浪圓頻率；ωe是遭遇頻率；x軸為波浪傳播方向。

為了避免入射波與反射波相互干擾，導(dǎo)致波浪出現(xiàn)破碎或不規(guī)則現(xiàn)象，采用消波的方法來減弱反射波。本文采用由Baker等人在二維邊界條件模式中首次提出的阻尼消波技術(shù)，從而可以實(shí)現(xiàn)在有限計(jì)算區(qū)域內(nèi)原本在無限計(jì)算域才能實(shí)現(xiàn)的消波效果，這就是所謂的數(shù)值阻尼消波原理。本文在計(jì)算域的尾部區(qū)域設(shè)置人工阻尼消波區(qū)，根據(jù)其他學(xué)者的經(jīng)驗(yàn)，消波區(qū)的長度取波長的1-2倍。

在阻尼消波區(qū)域內(nèi)，動量方程為：

（9）

（10）

式中：μ（x）為在阻尼段起點(diǎn)為零的單調(diào)遞增函數(shù)（線性遞增、指數(shù)遞增等）。

3 模型建立

3.1 計(jì)算區(qū)域

用來計(jì)算規(guī)則波頂浪中Wigley船?？v搖與垂蕩運(yùn)動的數(shù)值波浪水池尺度為：-1.5Lpp

3.2 邊界條件及FLUENT設(shè)置

數(shù)值波浪水池的邊界條件設(shè)置如下：水池前端為速度入口，水池后端為壓力出口，上邊界為壓力進(jìn)口，下邊界和兩個側(cè)面為對稱邊界。本文利用FLUENT軟件、使用VOF方法、結(jié)合SST k-ω模型、壓力耦合方程組的SIMPLE算法，對流項(xiàng)采用二階迎風(fēng)差分格式，求解非定常狀態(tài)下的紊流問題。應(yīng)用DEFINE-PROFILE宏設(shè)置速度入口和VOF控制二相分布；應(yīng)用DEFINE SOURCE宏添加動量源項(xiàng)實(shí)現(xiàn)數(shù)值波浪水池尾部區(qū)域的消波功能。采用動網(wǎng)格技術(shù)結(jié)合六自由度求解器模擬船舶的縱搖與垂蕩運(yùn)動，最終模擬Wigley船模在規(guī)則波中頂浪航行時的縱搖與垂蕩運(yùn)動。

4 結(jié)果與分析

對航行于波浪中的船舶，一般情況下頂浪時的縱搖和垂蕩最嚴(yán)重，因此本文以模擬計(jì)算頂浪中船舶縱搖和垂蕩運(yùn)動響應(yīng)作為所建立數(shù)值波浪水池的應(yīng)用例子?？紤]到與試驗(yàn)結(jié)果進(jìn)行對比分析，在計(jì)算船模運(yùn)動響應(yīng)時，入射波采用線性規(guī)則波。設(shè)波長λ=1.55 m、4.5 m和6.36 m，波幅A=0.045 m，通過模擬計(jì)算分別得到船模的垂蕩和縱搖的時歷曲線，如圖3和圖4所示。

從圖3和圖4可以看出，在短波情況下（λ/L=0.516 7），其垂蕩的幅值相比波浪的波幅較小，縱搖角度同樣很?。辉谥械炔ㄩL情況下（λ/L=1.5），垂蕩運(yùn)動幅值并未達(dá)到最大，但縱搖固有周期約等于波浪對船體的擾動力周期，發(fā)生諧搖，縱搖幅值達(dá)到最大；在最大波長情況下（λ/L=2.12），其垂蕩運(yùn)動幅值比中等波長情況大，縱搖運(yùn)動幅值比中等波長情況小，此時已經(jīng)處于隨波逐流狀態(tài)。

可以看出，對于垂蕩和縱搖運(yùn)動而言，規(guī)則波對船體擾動力的大小與波長對船長之比有很大關(guān)系：λ/L越小，垂蕩和縱搖越緩和，故大船在小的波浪上垂蕩和縱搖都不會很大；最大縱搖角度出現(xiàn)在中等波長下，此時波長略大于船長，最大波長的縱搖幅值介于小波長和中等波長之間；對于最大垂蕩幅值，筆者用勢流理論邊界元法得到在入射波波長為20～28 m、波幅0.045 m時，船模垂蕩運(yùn)動的幅值約等于波浪的幅值。

基于CFD方法的數(shù)值波浪水池可以處理復(fù)雜的運(yùn)動響應(yīng)問題，其計(jì)算結(jié)果考慮了粘性及非線性因素等的影響。為了較容易定量分析和比較計(jì)算結(jié)果，本文僅考慮線性因素，即只考慮船模的運(yùn)動幅度是線性的，在這種情況下用最小二乘法以正余弦函數(shù)對運(yùn)動幅度時歷進(jìn)行擬合，可以求出船體垂蕩和縱搖幅值。

規(guī)則波頂浪中船模的縱搖和垂蕩最嚴(yán)重，可以將垂蕩和縱搖運(yùn)動表示成如下形式：

zG= za cos（ωet=ε3） （11）

θ= θa cos（ωet=ε3） （12）

式中：ZG和θ分別為船體的垂蕩和縱搖運(yùn)動隨時間變化的函數(shù)；Za和θa分別為船體垂蕩運(yùn)動幅值和縱搖運(yùn)動幅值；ωe為遭遇頻率；ε3和ε5為相應(yīng)的相位。

在波幅A=0.045 m、頂浪海況下零航速時Wigley船模的垂蕩和縱搖運(yùn)動響應(yīng)幅值曲線，如圖5和圖6所示，圖中同時給出了邊界元法得到的勢流理論結(jié)果。從圖5、圖6可知，CFD數(shù)值計(jì)算結(jié)果與勢流理論解吻合良好，CFD計(jì)算考慮了粘性影響，故CFD數(shù)值解總體上比勢流解要小。

圖7給出了波幅A=0.045 m和λ/L=0.516 7、λ/L=1.5、λ/L=2.12三種情況下，在第11 s時Wigley船模附近數(shù)值水池的波形圖。從圖中可以看出尾部消波區(qū)域自由面相對平整，說明本文使用的消波方法效果很好。

在短波海況λ=1.55 m，A=0.045 m和Fr數(shù)分別為0.1、0.15、0.2、0.25、0.3時的運(yùn)動響應(yīng)曲線，見圖8和圖9。由圖8可知，在Fr=0.1時垂蕩運(yùn)動幅值最大；在Fr>0.1時船模垂蕩運(yùn)動幅值隨著航速的增加而減緩。由圖9可知，縱搖運(yùn)動的最大幅值出現(xiàn)在無航速情況下，并且縱搖運(yùn)動幅值隨著航速的增加而減緩。

圖10給出了A=0.045 m、λ=1.55 m、Fr=0.2時，在一個周期內(nèi)（t=0 T， t=0.25 T， t=0.5 T， t=0.75 T） Wigley船模底部的壓力分布云圖。從圖中可以看出各時刻船模底部壓力值過渡緩和，說明本文對Wigley船模附近網(wǎng)格劃分處理得當(dāng)。

圖11給出了同一個周期內(nèi)不同時刻自由面波形的變化以及Wigley船模的垂蕩和縱搖。從圖中可以清晰看到水線面的位置，整個船長范圍內(nèi)波峰與波谷十分清晰且與圖10中對應(yīng)時刻船底壓力分布是一致的，即處于波峰的船體部分對應(yīng)的底部壓力值高，波谷部分則相反。

圖12給出了A=0.045 m、λ=1.55 m、不同F(xiàn)r數(shù)對應(yīng)的波形圖。從圖中可以看出，在航速較小情況下，數(shù)值水池波形完整，船尾基本上沒有出現(xiàn)興波；隨著航速加大，船行波漸漸明顯，由船行波與規(guī)則波的相互作用疊加形成的波形和原來規(guī)則波在船后尾流區(qū)域錯分開來，這些現(xiàn)象都與實(shí)際中觀察到的現(xiàn)象一致。造成這些現(xiàn)象的原因是船模航行時波流共同作用于彎曲的船體，沿船體表面的壓力分布不一樣，導(dǎo)致船體周圍的水面升高或下降，在重力和慣性作用下在船后形成實(shí)際的船波，并與規(guī)則波相互作用。當(dāng)Fr=0.25和0.3時，可以明顯看到船體尾部形成的凱爾文角，整個船波系基本集中在凱爾文角所限定的扇形面范圍內(nèi)。

綜上所述，本文CFD模擬計(jì)算結(jié)果與勢流理論邊界元法得到的結(jié)果吻合較好。一方面，相對于勢流理論的計(jì)算，CFD方法更能真實(shí)地反映模擬流場，計(jì)算結(jié)果更加精確；另一方面，CFD方法更加便捷高效，具有易測量和控制及能細(xì)致描述船體周圍的流場等情況優(yōu)點(diǎn)，因此采用CFD方法預(yù)報波浪中船模的性能比物理模型試驗(yàn)更具有吸引力。

5 結(jié)語

本研究基于CFD方法，建立了三維粘性數(shù)值波浪水池，并用其模擬了規(guī)則波頂浪中以多個航速前進(jìn)的Wigley船模的縱搖與垂蕩運(yùn)動響應(yīng)。結(jié)果表明：本文的數(shù)值模擬方法能準(zhǔn)確地反映波浪與船模之間的干擾作用，并能給出較精確的垂蕩和縱搖運(yùn)動幅值，表明了在數(shù)值波浪水池中計(jì)算船舶運(yùn)動響應(yīng)的有效性；與物理模型實(shí)驗(yàn)相比，本文的方法更加便捷高效，具有易測量和控制及能細(xì)致描述船體周圍的流場情況等優(yōu)點(diǎn)，同時，該方法也為分析其他的海洋結(jié)構(gòu)物在波浪中的運(yùn)動響應(yīng)具有參考價值，豐富與擴(kuò)展了數(shù)值波浪水池的應(yīng)用。

參考文獻(xiàn)

[1] 戴世強(qiáng).認(rèn)識水波，駕馭水波-評介陶建華的專著《水波的數(shù)值模

擬》[J].水動力學(xué)研究與進(jìn)展，2007，22 （2）.

[2] Longuet-Higgins MS， Cokelet CD. The deformation of steep surface waves on

water， a numerical method of computation [C]. Proceedings of the Royal Society

of London 1976， A350.

[3] Baudic SF， Williams AN， Kareem A. A two-dimensional numerical wave

flume-Part 1： Nonlinear wave generation propagation and absorption [J]. Journal

of Offshore Mechanics and Arctic Engineering 2001， 123.

[4] Ryu S， Kim MH， Lynett PJ. Fully nonlinear wave-current interactions and

kinematics by a BEM-based numerical wave tank [J]. Computational Mechanics

2003， 32.

[5] Harlow FH， Welch JE. Numerical calculation of time-dependent viscous

incompressible flow of fluid with free surface [J]. Physics of Fluids 1965， 8.

[6] Hirt CW， Nichols BD. Volume of fluid （VOF） method for the dynamics of free

boundaries [J]. Journal of Computational Physics 1981， 3.

[7] Wang YX. Numerical wave channel with absorbed wave maker [J]. China

Ocean Engineering 1995， 9（2）.

[8] Wang YX. Wave slamming on offshore structure [A]. Proceedings of the 6th

International Off shore and Polar Engineering Conference： Vol 3 [C]. Los

Angeles， US， 1996. 1922196.

[9] 齊鵬，王永學(xué).三維數(shù)值波浪水池技術(shù)與應(yīng)用 [J]. 大連理工大學(xué)學(xué)報

2003， 43（6）.

[10] Xing-Kaeding Y， Jensen G， Hadzic I， Peric M. Simulation of flow-induced

ship motions in waves using a RANSE method [J]. Ship Technology Research

2004， 51.

[11] Menter FR. Zonal two equation k-w turbulence models for aerodynamic flows

[R]. AIAA-93-2906， 1993.

[12] Baker GR， Meiron DI， Orszag SA. Application of generalized vortex method to

nonlinear free-surface flows [C]. Proceedings of the 3rd International

Conference on Numerical Ship Hydrodynamics. Paris， France， 1981.

[13] 孫大鵬，李玉成.數(shù)值水槽內(nèi)的阻尼消波和波浪變形計(jì)算 [J].海洋

工程，2000，18 （2）.

[14] 高學(xué)平，曾廣冬，張亞.不規(guī)則波浪數(shù)值水槽的造波和阻尼消波 [J].

海洋學(xué)報 2002，24 （2）.

[15] 齊鵬，王永學(xué).三維數(shù)值波浪水池技術(shù)與應(yīng)用 [J].大連理工大學(xué)學(xué)

報 2003，43 （6）.endprint