船舶浅水航行下沉量和纵倾的数值计算

吴明 王骁 应荣镕 侯建军 石爱国

（海军大连舰艇学院航海系 大连 116018）

船舶浅水航行下沉量和纵倾的数值计算

吴明 王骁 应荣镕 侯建军 石爱国

（海军大连舰艇学院航海系 大连 116018）

基于CFD方法对船舶浅水航行航态变化问题进行研究。通过同时求解RANS方程和刚体运动学方程对船舶浅水航行流场进行数值模拟，并根据浅水流场空间受限的特点，采用三种动网格相结合的方法，解决船舶在浅水中航行时航态变化的网格更新问题，在数值模拟中还考虑了螺旋桨旋转的影响。文中以S60船模（Cb=0.6）为对象，计算了其浅水航行的下沉量和纵倾值，计算结果与船模水池试验数据对比，吻合良好。

计算流体力学；螺旋桨；浅水；下沉；纵倾

0 引言

船舶在浅水航行时，由于船体周围流场的压力和速度分布不均匀，会产生比较明显的下沉现象，尤其是大型船舶在浅水高速航行时，船体下沉现象更加明显，这已经成为影响船舶浅水航行安全的重要因素之一。为了确保船舶航行安全，许多学者对船舶浅水航行的下沉现象进行了理论研究和试验，最具代表性的是1967年，Tuck利用细长体理论提出的计算公式［1］。后期在Tuck公式的基础上，学者们又发展了许多解析式及半经验的计算公式。第23届ITTC操纵性委员会专门对浅水及受限水域船舶的操纵性问题进行研究总结，并推荐一些实际广泛应用的计算公式［2-5］，对保证船舶浅水航行安全起到一定的指导作用。

虽然针对船舶浅水域航行下沉量的研究已取得不少成果，但这些研究多基于势流理论方法，忽略了流体的粘性影响，而船舶浅水流场与船体及水底的壁面效应密切相关；因此，采用势流方法研究船舶浅水流场有一定的局限性。近年来，计算机技术的飞速发展为采用粘性流方法研究船舶水动力学问题提供了有利的条件，粘性流方法在船舶阻力、耐波性、操纵性等方面取得了不少成果。然而，目前关于考虑浮态变化船舶浅水航行数值模拟研究比较少，主要原因是在浅水流场中，船底与水底之间的空间狭小，采用简单的动网格方式，难以实现船体浮态变化的模拟。此外，现有的船舶浅水粘性流场数值模拟尚未考虑螺旋桨的影响，而受螺旋桨的排出流影响，船艉部流场会发生变化，这将会影响船舶在浅水中的航态，因此也是需要深入研究的问题。

针对上述问题，本文借助CFD商用软件Fluent，采用同时求解描述船舶流场的RANS方程和描述船舶自由运动的刚体运动学方程的方法，结合软件提供的三种动态网格方法，实现带螺旋桨船模浅水中航行流场的数值模拟，给出船模浅水航行时的下沉量和纵倾值，并与船模水池试验数据进行对比验证。

1 CFD数学模型

船舶浅水流场可视为不可压缩粘性流场，采用雷诺平均方法对船舶浅水域流场进行数值模拟，流场连续性方程和动量方程为：

式中：ui为流体时均速度分量，m/s；

p为流体压强，Pa；

fi为流体体积力分量，N；

ρ为流体密度，kg/m3；

ν为流体的运动粘性系数，m2/s；

ui′为相对于时均流速的湍流脉动速度分量，m/s；为雷诺应力。

采用船舶计算流体力学中应用较为广泛的RNG k-ε二方程［6］湍流模型封闭上述方程组。该模型是由重正化群理论推导得出，其湍流动能和湍流耗散率方程为：

式（3）和式（4）中各符号的物理意义详见文献［7］。

自由面的模拟采用VOF方法［8］来追踪。设第q项流体的体积分数为Cq，其输运方程为：

每一控制体的密度ρ为体积分数的平均值：

式中：ρq为第q项流体的密度，单位kg/m3。

2 船舶运动数学模型

2.1 船舶运动控制方程

为了计算船舶在波浪中的升沉和纵倾运动，需定义三个坐标系：第一个是固定坐标系，其空间位置不变，流体的速度、压力等在该坐标系求解；第二个是平移坐标系，其随船体作平移运动；第三个是随船坐标系，其随船体作旋转运动，船体的摇荡运动在该坐标系求解。通常将船体坐标系的原点取在船的重心上，把坐标轴取作与船的中心惯性主轴相重合，这样描述船舶摇荡运动简单方便。

设Λ（φθψ）为船舶旋转的欧拉角，对应的分别为船舶的横摇角、纵摇角和首摇角。如果只考虑角位置关系，向量在固定坐标系与随船坐标系之间的转换矩阵为：

设Ω=（p q r）为角速度在运动坐标系的投影，则角速度在运动坐标系的投影与欧拉角的导数之间的转换关系为：

船体可视为刚体，船舶摇荡运动遵循刚体运动的动量定理和动量矩定理。当重心位于原点时，在运动坐标系，船舶空间运动一般方程为：

式中：F=（X Y Z）为船体所受外力，N；L=（K M N）为船体所受外力合力对原点的力矩，N·m。

2.2 螺旋桨运动控制方程

在考虑螺旋桨的船舶摇荡运动数值模拟中，螺旋桨的运动比较复杂，可分解为两种运动：一是螺旋桨随船体作摇荡运动，二是螺旋桨绕转轴做旋转运动。由于螺旋桨与船体的相对位置固定，螺旋桨与船体的摇荡运动是同步的，所以前者可以通过求解船体运动控制方程获得。而螺旋桨绕转轴的旋转时，螺旋桨的轴向和旋转中心是变量，因此需要给出瞬态的控制方程。

在全局坐标系，设初始时刻螺旋桨旋转中心的坐标为A，船体绕转心旋转后，新的螺旋桨旋转中心坐标为A1=AS。

另取螺旋桨轴线上一点坐标B，船体绕转心旋转后，该点坐标为B1=BS，则新的螺旋桨轴向向量为：

2.3 数值求解步骤

船舶升沉和纵倾运动模拟的关键是如何实现流场与船舶运动的耦合求解，本文数值模拟的步骤如下：

（1）将船模按初始浮态固定，生成计算域和网格，对计算域进行初始化；

（2）开始迭代计算，待流场基本稳定时，加载船体和螺旋桨运动的方程；

（3）在当前时间步长Δt内，根据船体所受垂向力Z和纵摇力矩M，通过求解船体运动控制方程，得到船体垂荡位移的变化量ΔxCG和纵摇角的变化量ΔθCG；

（4）根据螺旋桨自转角速度更新螺旋桨绕轴的旋转角度，同时根据ΔxCG和ΔθCG更新船体和螺旋桨的位置和姿态；

（5）返回第（3）步，进行下一时间段的步长计算，求解航态改变后的流场，并根据计算结果再次更新螺旋桨绕轴的旋转角度，改变船体和螺旋桨的位置和姿态，直到设定的时间步数计算完毕。

通过上述步骤可实现流场与船体运动的耦合求解，实现带螺旋桨船模浅水航行升沉和纵倾运动的数值模拟。

3 数值模拟及结果分析

3.1 算例描述

选用S60船模进行数值模拟，文献［9］对该船模进行浅水试验，给出阻力、纵倾和升沉测试结果。表1为S60船模的主尺度数据，图1为本文采用GAMBIT软件生成的三维船体模型。

表1 S60船模船型参数

图1 S60船模

为了比较螺旋桨对船舶在浅水中航态的影响，本文分别对不带螺旋桨和带螺旋桨S60船模浅水域航行进行数值模拟，不带螺旋桨的数值试验7组，船模前进速度U=1.107 4~1.701 8m/s（对应的浅水傅汝德数Fh=0.5~0.8）；带螺旋桨的数值试验，选用前苏联克雷洛夫船舶研究院的CK3-110型螺旋桨，船模前进速度为U=1.328 8（对应的浅水傅汝德数Fh=0.6），试验中参考物理水池船模自航试验方法，根据不同航速下的船模自航点，确定该航速下螺旋桨的转速。

3.2 计算域及网格划分

根据水池试验的条件设定计算域的尺寸，深吃水比H/T=2，水池宽度为2.09L，计算域的长、宽、深分别为：4L×2.09L×2T，其中入口距船艏1L，出口距船艉2L，自由面距上、下边界分别约0.25L和2T。

船舶摇荡运动的模拟需要依靠动网格技术来实现。Fluent软件提供了三种基本的网格更新算法：弹簧近似光滑法（Spring Based Smoothing）、动态分层法（Dynamic Layering）和局部重构法（Local Remeshing），这三种算法通过拉伸、压缩网格或者增加、减少网格以及局部重新生成网格来解决边界运动的问题。

受限水域船舶下沉量数值计算的难点在于：船舶的升沉和纵倾运动，需要通过动网格技术来实现；而限制水域船体底部与航道底面的空间区域狭小，使网格弹性形变和重构的余地较小，在数值计算中极容易引起网格质量的下降，导致计算发散。因此如何减小或者避免数值计算中的网格质量下降，是受限水域船舶下沉量数值计算的关键。

本文采用三种动网格方案相结合的方法解决上述问题。如图2所示，跟据浅水流场及船模纵倾和升沉运动的特点，采用分区划分网格的方法，生成了结构/非结构混合网格。在船体附近Z1区域生成非结构网格，在该区域采用弹性形变和局部重构法实现船体的纵摇运动；在Z2区域生成结构网格，在该区域采用动态分层法实现船体的升沉运动，同时Z1区域整体也随Z2区域作升沉运动。采用这样的动网格方案可以减小网格弹性形变和重构，从而避免数值计算过程中因网格质量下降而引起计算发散。

图2 网格划分

3.3 离散格式及求解算法

采用有限体积法对微分方程进行离散，压力插值采用PRESTO格式，动量、湍流动能和湍流耗散率均采用二阶迎风格式离散，体积分数采用Modified HRIC格式离散。采用SIMPLE算法对流场进行求解。

3.4 数值模拟结果及分析

图3为S60裸船模浅水中不同浅水傅汝德数下的升沉量计算值与船模试验值的比较图，从图中可以看出，随着航速不断增加，船体浅水航行的下沉量不断增大。本文对浅水下沉量的计算与船模水池试验值相比吻合较好。

图4为S60裸船模浅水中不同浅水傅汝德数下的纵倾计算值与船模试验值的比较图。定义正的纵倾值表示船体产生向艉纵倾。从图中可以看出，在Fh＜0.75时，船体有微小的向艏纵倾；Fh＞0.75时，船体开始产生向艉纵倾，通过比较可以看出，本文方法对纵倾的计算值与水池试验值吻合得较好。

图3 浅水中的升沉

图4 浅水中的纵倾

图5 升沉时历（Fh=0.6）

图6 纵倾时历（Fh=0.6）

图5和图6中虚线为S60裸船模的浅水航行下沉量和纵倾的水池试验值，从图中可以看出，与S60裸船模浅水水池试验数据相比，螺旋桨的旋转使船舶浅水航行中的船体下沉量和艉倾增大。

图7和图8为浅水傅汝德数为0.6，船体航态固定时，S60裸船模与带螺旋桨船模艉部船底的压力分布对比图。从图中可以看出，尽管带桨船体在浅水中运动的流场十分复杂，但螺旋桨的排出流使紧贴船艉湿表面的作用力降低。由于船艉的压力是表面积分的结果，因此会导致船体向艉纵倾力矩增大，进而使船体下沉量和艉倾增大。浅水流场中船体的压力分布变化是引起船体的下沉量和纵倾变化的重要原因之一，本文的数值模拟可以捕捉到浅水流场变化的细节，为研究船舶浅水流场提供了一种有效的方法。

图7 裸船模艉底部表面压力等值线

图8 带螺旋桨船模艉底部表面压力等值线

4 结论

本文基于粘性流理论，实现带螺旋桨船模浅水域航行升沉和纵倾运动的数值模拟。数值计算结果表明：

（1）借助CFD软件，采用三种动网格相结合的方法，可以实现考虑船体浮态变化的船舶浅水航行的数值模拟；计算得到的S60裸船模浅水航行下沉量和纵倾与船模水池结果吻合良好，可为浅水域航行下沉量计算提供参考。

（2）螺旋桨的旋转会使船艉附近流场发生变化，引起船艉底部压力下降，这是船体艉纵倾及下沉量增大的重要原因之一。

本文对带螺旋桨船舶浅水域航行下沉量问题进行了初步研究。模拟的船舶航行工况较少，在今后的工作中将进一步研究螺旋桨转速同船舶航行下沉量及纵倾之间的关系。

［1］TUCK E O.Sinkage and trim in shallow water of finite width［J］.Schiffstechnik，1967，14：92-94.

［2］HOOFT JP.The be haviour of a ship in head waves at restricted water depth［J］.International Shipbuilding Progress，1974，21：367.

［3］BARRASSC B.A unified approach to squat calculations for ships［J］.PIANC Bulletin，1979，32：3-10.

［4］DAND IW，Ferguson A.The squat of full ships in shallow water［J］.The Naval Architect，1973（4）：237-255.

［5］ANKUDINOV V，DAGGETT L，HUVALL C，et al.Squat predictions for maneuvering applications［C］//Proceedings ofMARSIM’96，Copenhagen，［s.n.］，1996：467-495.

［6］YAKHOT V，ORSZAG SA，THANGAM S，et al.Development of turbulence models for shear flows by a double expansion technique［J］.Physics of Fluids A，1992（4）：1510-1520.

［7］王福军.计算流体动力学分析［M］.北京：清华大学出版社，2004：124.

［8］HIRTCW，NICHOLSBD.Volume of fluid（VOF）method for the dynamics of free boundaries［J］.Journal of Computational Physics，1981，39：201-225.

［9］JIANG T.Investigation of waves generated by ships in shallow water［C］//Proceeding of the 22nd Symposium on Naval Hydrodynamics，Washington DC：［s.n.］，1998：601-612.

Numerical calculation of ship sink age and trim in shallow water

WU Ming WANG Xiao YING Rong-rong HOU Jian-jun SHIAi-guo
（Dalian Naval Academy，Department of Navigation，Dalian 116018，China）

A study on the navigation behavior for ships in shallow water has been carried out based on CFD.This paper numerically simulates flow field around ship in shallow water by simultaneously solving equations of Reynolds Averaged Navier-Stokes（RANS）and kinematic equations of a rigid body.To settle the mesh update issue for the ships in shallow water when the navigation behavior is changed，this paper adopts the method of coupling three dynamic mesh according to the characteristics of limited space of shallow water.It also considers the influence of the propeller rotation.The numerical simulation is carried out on a S60model（Cb=0.6）to calculate the ship sinkage and trim travelling in shallow water.The numerical results are in good agreement with the experimental data.

CFD；ship；propeller；shallow water；sinkage；trim

U661.3

A

1001-9855（2013）03-0007-05

辽宁省博士科研启动基金项目（20111037），辽宁省优秀青年人才基金项目（2012J21DW027）。

2013-01-10；

2013-03-25

吴明（1984-），男，博士，研究方向：船舶耐波性。

王骁（1980-），男，博士，讲师，研究方向：舰船操纵性。

应荣镕（1984-），男，硕士，讲师，研究方向：舰船操纵性。

侯建军（1973-），男，博士，高级工程师，研究方向：舰船耐波性。

石爱国（1956-），男，教授，研究方向：舰船操纵性。