基于重叠网格的船舶运动计算方法研究

徐伟光，赵发明，何术龙

（中国船舶科学研究中心，江苏 无锡 214082）

基于重叠网格的船舶运动计算方法研究

徐伟光，赵发明，何术龙

（中国船舶科学研究中心，江苏 无锡 214082）

文章应用重叠网格技术和单相Level-set捕捉自由面方法，结合数值造波技术和六自由度运动方程，形成了基于RANS方程的船舶运动响应数值计算方法，并在中国船舶科学研究中心自主研发的船舶专用粘流计算软件OShip上实现。应用该方法计算了DTMB5512船模在迎浪规则波中的运动响应，并通过对计算得到的船模运动时历曲线进行傅里叶变换，提取其运动的幅值和相位，得到了不同航速下船模垂荡和纵摇运动的频响曲线，所得结果与试验结果吻合较好。

重叠网格；六自由度；传递函数和相位

0 引 言

船舶的耐波性作为衡量现代化船舶航行性能的重要指标之一，对船舶的安全性、作战或作业使用性以及适居性都会产生巨大的影响。船舶在波浪中航行的时候，船舶航行的姿态将会严重影响到船舶的阻力性能以及螺旋桨的推进效率，大幅度的运动会致使螺旋桨露出水面，造成飞车现象，同时甲板上浪将影响到船上设备的正常工作，大幅的摇荡还会大大降低对船员和乘客的舒适性。对此，如何准确地预报船舶的耐波性能将是本文研究的重要任务之一。

传统的耐波性预报方法主要基于势流理论，由于其简单、计算效率高的优点得到了广泛的应用，但在处理强非线性问题上还存在一定的局限性，如甲板上浪、自由面破碎等现象，且这些问题都会对船舶的运动产生较大的影响。近些年，随着计算机技术的飞速发展，CFD方法得到了越来越多的应用，由于考虑了流体的粘性效应和非线性因素，CFD方法在准确预报船舶耐波性上体现出了巨大的优势。

国内外已经开展了许多关于船舶运动的粘性数值模拟，Sat等人[2]采用密度函数法对规则波中W-igley船型和S60船型的运动响应进行了数值模拟；Carrica和Wilson等人[3]采用重叠网格技术模拟计算了DTMB 5512船型在中高航速下的大幅度运动响应情况；吴乘胜[4]对基于Rans方程的数值波浪水池技术展开了初步研究，其中对wigley及DTMB5512波浪运动响应进行了预报。

本文将采用Rans方法和重叠网格[5]，对DTMB5512船模在不同航速下，顶浪规则波中的运动响应进行预报，考察船模在较大的波长范围内（0.5L＜λ＜2.5L）的运动特性。本文采用的求解器为中国船舶科学研究中心自主研发的船舶专用粘流计算软件OShip。

1 数学模型及数值方法

1.1 控制方程与湍流模型

在求解不可压流体的非定常运动时，控制方程采用RANS方程，湍流模型为SST型k-ω二方程模型，并采用PISO[6]算法进行RANS方程和连续性方程的求解。

1.2 重叠网格方法

OShip软件采用重叠网格方法来处理运动问题。重叠网格是允许各个网格区域相互重叠、嵌套，并在各个子网格间通过插值来传递流场信息的一种结构网格方法。它将复杂的流动区域分成几何边界比较简单的子区域，极大简化了复杂形状物体网格的生成，提高了结构网格对外形的适应能力。重叠网格生成分为以下三个步骤：（1）寻点：在网格空间对已知点或区域进行搜索，得到与之相关的网格单元。（2）挖洞：在给定的网格中，将落在指定区域中的网格单元从参与流场计算的网格单元集合中剔除。（3）插值：建立各重叠网格间的耦合关系，由相邻贡献单元获得待插值点信息。

图1 重叠网格示意图Fig.1 Demonstration of overset grid

如图1所示，根据贡献单元的位置与插值点的相对位置，可以求得贡献点与待插值点间的插值系数[7]。同时根据这些插值系数，我们可以得到被插值点的任意变量值φ：

式中：ξi，ηi，ζi为插值系数，φi为贡献点上的变量值。

1.3 自由面处理

本文采用单相Level-set方法[8]进行自由面的捕捉，自由面的位置通过Level-set函数的函数值φ确定。单相Level-set方法仅需求解在距离函数φ≤0（表示水相）区域内的流场，空气相中的流场速度则通过速度扩展（velocity extension）的方法来计算的。两相流界面的过渡问题通过显性施加阶跃条件进行处理。此外，在气体中，只需要布置少许网格来满足计算条件，因此相比两相方法，计算资源的消耗大大减小，计算稳定性增加。

1.4 造波方法

本文的计算工况均为线性规则波，可通过定义初始边界条件来实现，在任意时刻的某一点的波浪幅值可以按（2）式表达为：

式中：ξ为自由液面上各点在任意时刻的垂向位置，a为波幅，ω为波浪的自然频率，ωe为遭遇频率，k为波数。在上述波浪初始条件下，速度和压力的表达式为：

式中：U（x，y，z，t），W（x，y，z，t）分别为流体中的点在任意时刻水平方向和垂直方向的速度分量，U0为船模的速度。

1.5 六自由度运动

为了求解运动方程，需要建立两个坐标系，一个是以地面为参考的固定坐标系。另一个是载体坐标系，其原点放在船舶重心处，船尾方向为x轴正方向，右舷方向为y轴正方向，垂直于水线面向上为z轴正方向。

固定坐标系和载体坐标系间的转化关系参见文献[3]。载体坐标系下，船舶六自由度运动方程为：

式中：X，Y，Z，K，M和N分别为纵荡、横荡、垂荡力、横摇、纵摇及艏摇力矩。为线加速度，为欧拉角加速度。六自由度运动求解的流程如图3所示，运动网格重叠结果如图2所示。从图4中可以看出，当船模运动后，重叠网格方法仅需重新计算背景网格即与船模贴体网格的插值关系，而不用对网格重新生成，这样既保证了船体网格的质量保持不变同时也保证了自由面处网格的密度，从而提高了运动求的精度。

图2 固定坐标系和载体坐标系Fig.2 Earth and ship fixed reference systems

图3 六自由度运动求解流程图Fig.3 Solution strategy

图4 动网格重叠示意图Fig.4 Demonstration of dynamic overset grid

1.6 离散方法

控制方程采用体积中心有限差分格式进行离散，时间项采用2阶欧拉向后差分格式、对流项采用2阶混合差分格式，粘流项采用2阶中心差分格式离散。Level-set对流项同样采用2阶混合差分格式。

2 数值计算

2.1 计算对象

本文选用DTMB5512船模作为计算对象。DTMB5512是ITTC推荐的标准模型之一，具有大量的试验数据和计算结果。DTMB5512的缩尺比为1/46.6。

2.2 计算网格

本文的算例中的船体和计算域左右对称，因此在生成网格时只需生成其中的一半。计算区域如图5所示，计算域尺度如下：

图5 计算域网格示意图Fig.5 Computational grid

a.入口—模型首部前0.6 Lpp处；b.出口—模型为尾部1.5 Lpp处；c.侧边界—模型侧方1.0 Lpp处；d.上边界—水线以上0.25 Lpp处；e.下边界—水线以下1.0 Lpp处。

2.3 网格生成

为了能够准确计算船舶在波浪中的运动，在划分背景网格时，保证一个船长内至少含有60个网格。整个计算区域的网格总数为1 224 000，计算网格如图6所示。

图6 船首、尾网格Fig.6 Local grid

2.4 计算工况

数值模拟中的各计算工况的参数均列于表1中。

表1 DTMB5512耐波性计算工况Tab.1 Summary of simulation conditions

3 计算结果及分析

3.1 数据处理方法

船模在顶浪规则波中作垂荡和纵摇的耦合运动，且船模的运动和受力具有周期性，因此将船模的运动用傅里叶级数展开可以写为：

式中：x3n，x5n和γ3n，γ5n分别为船模垂荡和纵摇运动的n阶幅值和相位，而垂荡和纵摇运动响应的传递函数可表达为（1阶运动）：

式中：γ31，γ51为传递函数的相位角。

船模试验中规定：t=0时规则波的波峰正好经过船模重心的纵向位置LCG，而本文数值计算中t=0时，LCG处未必对应着波峰，因此在计算相位时还需进行修正：

式中：γζLCG为t=0时LCG处的波浪相位，γζ0为t=0时入口处的波浪相位，s为入口至LCG的距离，γ3，γ5则为修正后的相位角。

3.2 传递函数及运动响应相位

在幅值和相位的处理中，采用船模运动稳定后有效的数据进行分析，采用快速傅里叶变换进行离散数据的处理从而获得船模运动的1阶幅值和对应的相位，将由FFT获得的垂荡和纵摇运动的1阶幅值代入（7）式和（8）式中，得到垂荡和纵摇运动的传递函数及运动响应相位，图7-14给出了本文计算值与试验值和切片法计算值（SMP）的比较结果。

图7 Fr=0.28垂荡运动传递函数Fig.7 Heave transfer functions at Fr=0.28

图8 Fr=0.28纵摇运动传递函数Fig.8 Pitch transfer functions at Fr=0.28

图9 Fr=0.28垂荡运动响应相位Fig.9 Heave phases at Fr=0.28

图10 Fr=0.28纵摇运动响应相位Fig.10 Pitch phases at Fr=0.28

图11 Fr=0.41垂荡运动传递函数Fig.11 Heave transfer functions at Fr=0.41

图12 Fr=0.41纵摇运动传递函数Fig.12 Pitch transfer functions at Fr=0.41

图13 Fr=0.41垂荡运动响应相位Fig.13 Heave phases at Fr=0.41

图14 Fr=0.41纵摇运动响应相位Fig.14 Pitch phases at Fr=0.41

图15 垂荡运动传递函数随遭遇频率变化Fig.15 Heave transfer functions for all Fr versus fe

图16 纵摇运动传递函数随遭遇频率变化Fig.16 Pitch transfer functions for all Fr versus fe

由图7-10的结果可得出如下结论：

（1）当Fr=0.28时，本文的计算值与试验值符合较好，计算所得船模运动随波长变化的趋势与试验值一致，本文计算结果优于切片法（SMP）的计算值。

（2）当λ增大时，垂荡运动的传递函数值也随之增大，在λ/L=1.3附近垂荡的传递函数值到达峰值点。随着波长进一步增加，传递函数值有所下降，而当λ/L＞1.5时，随着波长的增加，传递函数值再次平稳地上升。同时，纵摇运动的传递函数值在λ/L=1.5附近达到峰值点，之后随着波长的增加而减小最终趋于平缓。另一方面，从计算结果来看，峰值处传递函数计算值要略小于试验值，但峰值点出现位置捕捉得较为准确，相比之下切片法所计算的峰值位置向右偏移较为严重，且峰值也存在较大误差。

（3）船模的垂荡和纵摇运动存在着一定的相位差，γx3要超前γx5，图9、10中，随着波长的增加，γx3、γx5也随之减小，并分别在λ/L=0.85和0.7附近达到最小值，之后随着波长的增加而增加，最后保持不变。最终γx3=0表明船模的运动与波浪同步。计算结果表明，在峰值点附近本文的相位计算值要稍大于试验值，即相位要超前一些，但从整体趋势上来看二者符合得较好，而切片法的计算值在峰值附近要偏小很多。

图11-14为Fr=0.41时的计算结果，可以看出：

（1）垂荡和纵摇运动传递函数分别在λ/L=1.4和λ/L=1.6处达到峰值点，之后传递函数值随着波长的增加而减小。本文的计算值与试验值符合得较好，准确地捕捉到了峰值点的位置，同时曲线的趋势与试验值一致，但在峰值附近本文的计算值要偏小一些。相比之下切片法的计算值误差则较大，且其计算值的峰值点右移。

（2）垂荡和纵摇运动的相位γx3和γx5的谷值点分别出现在λ/L=1.0以及λ/L=0.8附近。本文相位计算值在波长较大的区域内与试验值符合较好，在短波范围内则稍差，且在1.0＜λ/L＜1.5范围内计算值稍微偏大，但整体趋势与试验值符合得较好。

图15和16给出了运动计算结果：当Fr=0.28时，垂荡峰值点出现在λ/L=1.3附近，此时fe=1.014；当Fr=0.41时，λ/L=1.4，此时fe=1.130，两个峰值处的遭遇频率并不相同，这表明对于DTMB 5512船模来说，某一航速下的峰值点不一定位于固有频率处，即发生共振并非对应着最大运动幅值。

根据Lewis，Journee，Fonseca和Soares等人[9-10]的结论，对于绝大多数情况，某一航速下，垂荡和纵摇运动的峰值发生在L/λ=0.75（λ/L=1.33）附近，此时波浪扰动力最大，进一步当遭遇频率与船模的固有频率相等即fe=fn，两个条件同时满足时，船舶将产生最大运动响应。fn计算如（12）式：

式中：ω3由Lloyd（1989）按照（13）式给出：

式中：C33为垂荡恢复力，A33为垂荡附加质量。由于Fr=0.41时的遭遇频率fe=1.130更加接近船模的固有频率fn，因此该航速下的传递函数曲线更加陡峭，幅值也更大，本文的计算结果也很好地验证了这一点。

4 结 论

本文基于重叠网格技术，对波浪中DTMB 5512船模的运动响应进行了数值模拟，并将计算结果与试验值和切片法（SMP）的计算值进行了比较和验证，得到以下结论：

（1）本文采用的基于重叠网格的数值方法在预报船模运动响应时得到了稳定、收敛的计算结果，证明了该方法的稳定性。

（2）本文计算了两种航速Fr=0.28和Fr=0.41的运动响应，得到的传递函数和相位与试验值都符合得较好，能够较为准确地捕捉到峰值点位置，但计算值在峰值点附近偏小，而切片法（SMP）的计算值误差则较大。

（3）对于DTMB 5512船模来说，某一航速下的峰值点不一定位于固有频率处。对于绝大多数情况，某一航速下，垂荡和纵摇运动的峰值发生在L/λ=0.75（λ/L=1.33）附近，并且当遭遇频率与船模的固有频率相等fe=fn时，将产生最大的运动响应。本文的计算结果在一定程度上验证了这一点，不足的地方将在今后的工作中继续完善。

综上，本文的计算结果证明了本文采用的基于重叠网格技术的数值方法在单体船耐波性预报上是适用、稳定以及有效的。该方法在多体船及多体复合船耐波性预报上的推广将在未来的工作中完成。

[1]周秀红,赵发明,王丽艳.船舶粘流计算软件“Oship”开发[J].中国造船,2014,55(1):90-103. Zhou Xiuhong,Zhao Faming,Wang Liyan.Development of software‘Oship’for open/overlap ship viscous caculation[J]. Shipbulding of China,2014,55(1):90-103.

[2]Miyata S H,Sato T.CFD simulation of 3-dimensional motion of a ship in waves:Application to an advancing ship in regular heading waves[J].Journal of Marine Science and Technology,1999,4(3):108-116.

[3]Carrica P M,Wilson R V,Noack R W,et al.Ship motions using single-phase Level-set with dynamic overset grids[J]. Computers&Fluids,2007,36(9):1415-1433.

[4]吴乘胜.基于Rans方程的数值水池研发及其应用研究[D].无锡:中国船舶科学研究中心,2014. Wu Chengsheng.Development of a numerical wave tank based on RANS equations and its application[D].Wuxi:China Ship Science Research Center,2014.

[5]赵发明,高成君,夏 琼.重叠网格在船舶CFD中的应用研究[J].船舶力学,2011,15(4):332-341. Zhao Faming,Gao Chengjun,Xia Qiong.Overlap grid research on the application of ship CFD[J].Journal of Ship Mechanics,2011,15(4):332-341.

[6]Issa R I.Solution of the implicitly discretised fluid flow equations by operator-splitting[J].Journal of Computational Physics, 1986,62(1):40-65.

[7]Wilson R,Carrica P,Hyma M,Stern F.A single-phase level set method with application to breaking waves and forward speed diffraction problem[C]//25th Symposium on Naval Hydrodynamics St.Canada,2004.

[8]Lewis E V.Principles of naval architecture[J].The Society of Naval Architects and Marine Engineers,New York,1999.

[9]Journee J M J.Experiments and calculations on four wigley hullforms[R].Delft University of Technology,Ship Hydromechanics Lab,Report 1999,No.909.

[10]Fonseac N,Soares C G.Experimental investigation of the onliner effects on vertical motions and loads of a containership in regular wave[J].Journal of Ship Research,2004b,48:118-147.

Research on calculation method for ship motion based on dynamic overset approach

XU Wei-guang,ZHAO Fa-ming,HE Shu-long
(China Ship Scientific Research Center,Wuxi 214082,China)

This research applies dynamic overset grid and single phase Level-set method and accommodates numerical wave making technique and 6 DOF motion equations to form the numerical method for calculating ship motion response,and this method is realized by using self-developed calculation software-Oship developed by China Ship Scientific Research Center.Motion responses of DTMB 5512 advancing in regular head waves are computed using this method.Through the fast Fourier transform of time history of ship motions, amplitudes and phases of ship motions are acquired to determine the heave and pitch transfer functions and phases of DTMB 5512 at different speed.The simulation results show fair agreements with experiments.

overset grid;6 DOF;transfer function and phase

U661.3

：Adoi:10.3969/j.issn.1007-7294.2016.07.005

1007-7294（2016）07-0824-09

2015-08-14

徐伟光（1989-），男，硕士研究生，E-mail：654050813@qq.com；赵发明（1979-），男，高级工程师。