基于CFD方法的船舶运动响应仿真

韩翔希+冯志强+邱昂+符妃+范少涛

摘 要：本研究采用CFD方法建立三维波浪水池，并在所模拟的数值波浪水池中计算Wigley船型顶浪航行在规则波时的运动响应。采用RANS方程求解非定常的不可压缩粘性流体，并采用VOF方法对自由面进行动态模拟。通过编写UDF设置边界入口速度、波高及在波浪水池尾部1～2倍波长区域消波。采用动网格技术结合六自由度求解器模拟船舶的纵摇与垂荡运动。将采用CFD方法模拟得到的Wigley船运动响应与通过边界元法所得到的势流结果进行对比，发现两者吻合良好。本文所探索的方法可有效计算船舶在规则波中的运动响应，便捷高效，易于测量和控制，并能细致描述船体周围的流场情况。同时，也可为模拟分析其他海洋结构物在波浪中的运动响应具有参考价值，丰富与扩展数值波浪水池的应用。

关键词：数值波浪水池；船舶运动响应；VOF；六自由度求解器；动网格

中图分类号：U661.3 文献标识码：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技术的发展，基于CFD方法的数值造波技术逐渐成为研究热点。Harlow等人提出的MAC方法和Hirt等提出的VOF方法，推动了粘性自由表面数值计算技术的快速发展。此后，许多学者在粘性自由液面的数值水池上做了大量工作。Wang采用VOF方法模拟出二维数值波浪水池，并开展波浪对近海平台底部冲击过程的研究；齐鹏和王永学在前人二维波浪水池的基础上，采用VOF的方法，根据造波机的原理建立三维数值波浪水池；Xing-Kaeding等人基于RANS方法模拟计算船舶在波浪中的运动响应，自由面同样采用VOF方法进行处理。

二维数值水槽的建立及分析已经相当成熟，但应用到工程项目中带有局限性，所以建立三维波浪数值模型分析计算波浪与结构物相互作用具有特别重要的意义。另外，目前计算波浪作用在结构物上的载荷大多通过线性势流理论解决，但线性势流理论有明显的局限性，很难考虑粘性和非线性的影响，如果能够模拟高阶非线性波且结果可靠，则可用数值波浪水池模拟各种波况下、特别是极端波况下的波浪运动特性。本研究基于CFD建立三维数值波浪水池，自由液面采用VOF处理，通过编写UDF设置边界入口速度及动量方程中添加源项实现阻尼消波，采用动网格技术结合六自由度求解器模拟船舶的纵摇与垂荡运动，模拟了规则波顶浪中三维船体纵摇与垂荡运动的计算，并与势流理论边界元法得到的结果相比较，两者吻合很好。通过试验，本方法能准确预测规则波顶浪中船舶的纵摇与垂荡运动响应，能细致描述船舶周围的流场。本方法也可为模拟分析其他的海洋结构物在波浪中的运动响应具有参考价值，丰富与扩展数值波浪水池的应用。endprint

2 数学模型及数值方法

2.1 控制方程

本文采用两相不可压的RANS方程求解非定常不可压粘性流体，RANS方程采用 SST k-ω湍流模型封闭，自由液面采用VOF法处理，消波区域采用在动量方程中添加源项的办法进行消波。

连续方程： （1）

RANS方程： （2）

式中：U为速度矢量；流体密度定义为 ，其中体积分数aq表示单元内第q相流体占的体积与该单元的体积之比，本文采用气液两相模拟q=1，2，并且

；μ为相体积分数平均的动力粘性系数。

Menter提出了SST k-ω 两方程模型，它在近壁面保留了原始k-ω的模型，在远离壁面的地方应用了k-ω模型，其k方程和ω方程如下：

（3）

（4）

自由液面采用VOF法处理，aq满足方程：

（5）

2.2 数值造波和消波方法

本文采用在边界设置入口速度的方法造波，需要编写UDF控制入口速度以及波高。

根据线性波理论，小振幅波的速度场为：

（6）

（7）

小振幅波的自由液面高度

（8）

式中：U0为船舶航速；A是波幅；k为波数；ω是波浪圆频率；ωe是遭遇频率；x轴为波浪传播方向。

为了避免入射波与反射波相互干扰，导致波浪出现破碎或不规则现象，采用消波的方法来减弱反射波。本文采用由Baker等人在二维边界条件模式中首次提出的阻尼消波技术，从而可以实现在有限计算区域内原本在无限计算域才能实现的消波效果，这就是所谓的数值阻尼消波原理。本文在计算域的尾部区域设置人工阻尼消波区，根据其他学者的经验，消波区的长度取波长的1-2倍。

在阻尼消波区域内，动量方程为：

（9）

（10）

式中：μ（x）为在阻尼段起点为零的单调递增函数（线性递增、指数递增等）。

3 模型建立

3.1 计算区域

用来计算规则波顶浪中Wigley船模纵摇与垂荡运动的数值波浪水池尺度为：-1.5Lpp

3.2 边界条件及FLUENT设置

数值波浪水池的边界条件设置如下：水池前端为速度入口，水池后端为压力出口，上边界为压力进口，下边界和两个侧面为对称边界。本文利用FLUENT软件、使用VOF方法、结合SST k-ω模型、压力耦合方程组的SIMPLE算法，对流项采用二阶迎风差分格式，求解非定常状态下的紊流问题。应用DEFINE-PROFILE宏设置速度入口和VOF控制二相分布；应用DEFINE SOURCE宏添加动量源项实现数值波浪水池尾部区域的消波功能。采用动网格技术结合六自由度求解器模拟船舶的纵摇与垂荡运动，最终模拟Wigley船模在规则波中顶浪航行时的纵摇与垂荡运动。

4 结果与分析

对航行于波浪中的船舶，一般情况下顶浪时的纵摇和垂荡最严重，因此本文以模拟计算顶浪中船舶纵摇和垂荡运动响应作为所建立数值波浪水池的应用例子。考虑到与试验结果进行对比分析，在计算船模运动响应时，入射波采用线性规则波。设波长λ=1.55 m、4.5 m和6.36 m，波幅A=0.045 m，通過模拟计算分别得到船模的垂荡和纵摇的时历曲线，如图3和图4所示。

从图3和图4可以看出，在短波情况下（λ/L=0.516 7），其垂荡的幅值相比波浪的波幅较小，纵摇角度同样很小；在中等波长情况下（λ/L=1.5），垂荡运动幅值并未达到最大，但纵摇固有周期约等于波浪对船体的扰动力周期，发生谐摇，纵摇幅值达到最大；在最大波长情况下（λ/L=2.12），其垂荡运动幅值比中等波长情况大，纵摇运动幅值比中等波长情况小，此时已经处于随波逐流状态。

可以看出，对于垂荡和纵摇运动而言，规则波对船体扰动力的大小与波长对船长之比有很大关系：λ/L越小，垂荡和纵摇越缓和，故大船在小的波浪上垂荡和纵摇都不会很大；最大纵摇角度出现在中等波长下，此时波长略大于船长，最大波长的纵摇幅值介于小波长和中等波长之间；对于最大垂荡幅值，笔者用势流理论边界元法得到在入射波波长为20～28 m、波幅0.045 m时，船模垂荡运动的幅值约等于波浪的幅值。

基于CFD方法的数值波浪水池可以处理复杂的运动响应问题，其计算结果考虑了粘性及非线性因素等的影响。为了较容易定量分析和比较计算结果，本文仅考虑线性因素，即只考虑船模的运动幅度是线性的，在这种情况下用最小二乘法以正余弦函数对运动幅度时历进行拟合，可以求出船体垂荡和纵摇幅值。

规则波顶浪中船模的纵摇和垂荡最严重，可以将垂荡和纵摇运动表示成如下形式：

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

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

式中：ZG和θ分别为船体的垂荡和纵摇运动随时间变化的函数；Za和θa分别为船体垂荡运动幅值和纵摇运动幅值；ωe为遭遇频率；ε3和ε5为相应的相位。

在波幅A=0.045 m、顶浪海况下零航速时Wigley船模的垂荡和纵摇运动响应幅值曲线，如图5和图6所示，图中同时给出了边界元法得到的势流理论结果。从图5、图6可知，CFD数值计算结果与势流理论解吻合良好，CFD计算考虑了粘性影响，故CFD数值解总体上比势流解要小。

圖7给出了波幅A=0.045 m和λ/L=0.516 7、λ/L=1.5、λ/L=2.12三种情况下，在第11 s时Wigley船模附近数值水池的波形图。从图中可以看出尾部消波区域自由面相对平整，说明本文使用的消波方法效果很好。

在短波海况λ=1.55 m，A=0.045 m和Fr数分别为0.1、0.15、0.2、0.25、0.3时的运动响应曲线，见图8和图9。由图8可知，在Fr=0.1时垂荡运动幅值最大；在Fr>0.1时船模垂荡运动幅值随着航速的增加而减缓。由图9可知，纵摇运动的最大幅值出现在无航速情况下，并且纵摇运动幅值随着航速的增加而减缓。

图10给出了A=0.045 m、λ=1.55 m、Fr=0.2时，在一个周期内（t=0 T， t=0.25 T， t=0.5 T， t=0.75 T） Wigley船模底部的压力分布云图。从图中可以看出各时刻船模底部压力值过渡缓和，说明本文对Wigley船模附近网格划分处理得当。

图11给出了同一个周期内不同时刻自由面波形的变化以及Wigley船模的垂荡和纵摇。从图中可以清晰看到水线面的位置，整个船长范围内波峰与波谷十分清晰且与图10中对应时刻船底压力分布是一致的，即处于波峰的船体部分对应的底部压力值高，波谷部分则相反。

图12给出了A=0.045 m、λ=1.55 m、不同Fr数对应的波形图。从图中可以看出，在航速较小情况下，数值水池波形完整，船尾基本上没有出现兴波；随着航速加大，船行波渐渐明显，由船行波与规则波的相互作用叠加形成的波形和原来规则波在船后尾流区域错分开来，这些现象都与实际中观察到的现象一致。造成这些现象的原因是船模航行时波流共同作用于弯曲的船体，沿船体表面的压力分布不一样，导致船体周围的水面升高或下降，在重力和惯性作用下在船后形成实际的船波，并与规则波相互作用。当Fr=0.25和0.3时，可以明显看到船体尾部形成的凯尔文角，整个船波系基本集中在凯尔文角所限定的扇形面范围内。

综上所述，本文CFD模拟计算结果与势流理论边界元法得到的结果吻合较好。一方面，相对于势流理论的计算，CFD方法更能真实地反映模拟流场，计算结果更加精确；另一方面，CFD方法更加便捷高效，具有易测量和控制及能细致描述船体周围的流场等情况优点，因此采用CFD方法预报波浪中船模的性能比物理模型试验更具有吸引力。

5 结语

本研究基于CFD方法，建立了三维粘性数值波浪水池，并用其模拟了规则波顶浪中以多个航速前进的Wigley船模的纵摇与垂荡运动响应。结果表明：本文的数值模拟方法能准确地反映波浪与船模之间的干扰作用，并能给出较精确的垂荡和纵摇运动幅值，表明了在数值波浪水池中计算船舶运动响应的有效性；与物理模型实验相比，本文的方法更加便捷高效，具有易测量和控制及能细致描述船体周围的流场情况等优点，同时，该方法也为分析其他的海洋结构物在波浪中的运动响应具有参考价值，丰富与扩展了数值波浪水池的应用。

参考文献

[1] 戴世强.认识水波，驾驭水波-评介陶建华的专著《水波的数值模

拟》[J].水动力学研究与进展，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] 齐鹏，王永学.三维数值波浪水池技术与应用 [J]. 大连理工大学学报

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] 孙大鹏，李玉成.数值水槽内的阻尼消波和波浪变形计算 [J].海洋

工程，2000，18 （2）.

[14] 高学平，曾广冬，张亚.不规则波浪数值水槽的造波和阻尼消波 [J].

海洋学报 2002，24 （2）.

[15] 齐鹏，王永学.三维数值波浪水池技术与应用 [J].大连理工大学学

报 2003，43 （6）.endprint