摘要：商船会因装载不同而吃水大幅增减，舰艇会因战损而出现严重纵横倾，如何在浮态急剧变化时估算船舶的快速性，是一项具有普遍意义的课题。为了解决上述问题，采用计算流体力学方法对某型船模进行了大量计及自由表面的不同浮态下的阻力数值计算，在深入分析船体浮态改变时影响船舶阻力诸因素的基础上，对数值计算结果进行了研究，推导出船舶吃水、纵倾和横倾变化时的阻力估算公式。据此进行的数值计算结果与水池试验结果对比，吻合良好，从而验证了提出的计算模式的可行性，为解决船舶任意浮态下的阻力计算问题提供了新的思路。
关键词：船舶、舰船工程；阻力；浮态变化；计算流体动力学
1问题的提出
舰船在运行期间往往不是处于正常浮态：对于民船常常因装载或压载不同而产生吃水、纵倾的变化，对于军舰，还会因弹药物资消耗、战损等原因产生横倾的变化，这些浮态的变化会影响到舰船的快速性，而舰船在吃水、纵倾、横倾变化时的阻力规律，是阻力研究领域的空白。
文章采用船舶计算流体动力学(Computational Fluid Dynamics，CFD)与理论研究相结合的方法，探讨计及自由面的舰船吃水、纵倾、横倾改变时的阻力变化规律，以求低耗、省时、精确地获得舰船几种浮态下的阻力，为解决舰船任意浮态下的阻力问题提供一种新的思路。
2船舶CFD方法
以INSEAN 2340船模为研究对象，用Pro/Engineer软件生成该船体的实体模型，然后在Gambit软件中分别生成正常浮态、不同吃水、纵倾和横倾的计算域网格，最后通过Fluent软件对不同浮态、不同航速时的船模受力进行计算。
2.1计算模型
INSEAN 2340船模为美国海军驱逐舰DDG51模型，是国际船模试验水池会议(International Towing Tank Conference，ITTC)推荐的用于船舶水动力数值模拟研究的基准船模之一，其型值数据见表1。
表1 INSEAN2340船模主要参数
|
缩尺比 |
24.83 |
|
垂线间长/m |
5.72 |
|
船宽/m |
0.76 |
|
吃水/m |
0.248 |
|
排水量/t |
0.549 |
|
方形系数 |
0.506 |
按照与水池试验等比例船模建立计算模型，坐标系的定义为：坐标原点位于船模重心，X轴指向下游，Z轴垂直水平面向上，Y轴指向右舷。计算域设置成立方体形状，各边界的位置为：
入口距船艏1倍船长，出口距船艉3倍船长，顶部边界距水线0.5倍船长，底部边界距水线1倍船长，左、右边界距船中纵剖面1倍船长。
2.2控制方程
对于不可压缩粘性流流动，以张量法表示的控制方程如下：
连续性方程： (1)
RANS方程： (2)
式(1)和式(2)中，ui为流体时均速度分量；p为流体压强；fi为流体体积力分量；ρ为流体密度；v为流体的运动粘性系数；u′i为相对于时均流速的湍流脉动速度分量；为雷诺应力。
2.3湍流模型
经过对多种湍流模型的比较研究，文章在低速、中速和高速数值模拟时分别选用剪切应力输运模型(Shear Stress Transport，SST)SST k-ω、Realizable k-ε模型和重整化群湍流模型(Renormalization Group，RNG)RNG k-ε湍流模型来封闭RANS方程，各模型方程详见文献。
2.4网格划分
采用混和网格技术划分计算域[7]，在船体表面划分边界层，边界层为三棱柱形的半结构网格；在船体附近区域划分四面体形状的非结构网格；远离船体的区域采用结构网格；自由面附近采用尺度较小的结构网格。计算域的网格分布如图1所示。
2.5边界条件
1)入流平面：设置为压力入口条件；
2)出流平面：设置为压力出口条件；
3)物面：船体表面设置为无滑移壁面；
4)外边界：设置为自由滑移壁面。
图1计算域网格划分
2.6自由面的模拟
采用流体何种法(Volume of Fluid，VOF)模拟自由面。VOF方法引进了体积分数的概念，用体积分数表示网格单元内的流体状态(充满、没有或存在自由液面)。
定义第q种流体的体积分数为Cq，则存在3种情况：
1)Cq=0，表示网格单元内不含第q种流体；
2)Cq=1，表示网格单元内充满第q种流体；
3)0<Cq<1，表示网格单元内存在自由面。
体积分数Cq同样需要通过对其控制方程的求解获得。
3吃水变化时的阻力变化规律
3.1数值模拟结果
计算了船模5种吃水、5种不同弗劳德数下的阻力，用弗劳德换算法将船模吃水变化阻力换算成实船阻力并拟合成曲线，如图2所示。
图2吃水变化弗劳德数2阻力曲线
图2中5根曲线，分别为设计吃水，以及吃水ΔT增加0.225m、0.55m、0.825m和1.1m的弗劳德数-阻力曲线。
3.2试验分析
将阻力随吃水变化的规律表述为：
(3)
则： (4)
式(4)中，Δ为排水量。
y可视为弗劳德数和方形系数的函数，有：
(5)
式(5)中，是同一弗劳德数下y的平均值，它是Fr的函数，I是考虑CB随吃水变化的修正因子，由于CB的变化主要是因吃水T变化引起的，因此，I可简化为：
(6)
对…y进行弗劳德数Fr的2阶拟合，可得：
(7)
拟合结果如图3所示。
图3…y随Fr变化曲线
对I进行ΔT/T的1阶拟合，可得：
I=-1.3182×ΔT/T+0.1483 (8)
拟合结果如图4所示。
至此，得到一定条件下的计算吃水变化时阻力的公式：
(9)
式(9)的适用条件为：
1)瘦型带球鼻艏船，CB为0.5左右；
2)Fr为0.138～0.410之间；
3)吃水变化范围为正常吃水T～1.18T。
图4 I随ΔT/T变化曲线
4纵倾变化时的阻力变化规律
变纵倾试验分为艏倾和艉倾两个部分。分别是绕OY轴旋转0.5°、1.0°、1.5°、2.0°、2.5°、-0.5°、-1.0°、-1.5°、-2.0°、-2.5°，其中艉倾为“+”，艏倾为“-”。
4.1数值模拟结果
计算了船模11种纵倾、5种速度下的阻力，用弗劳德换算法将船模纵倾变化阻力换算成实船阻力，并拟合成曲线，如图5、图6所示。
4.2试验分析
使用纵倾引起的相对于舰艇正浮时的阻力增加率，来描述其阻力变化规律。阻力增加率公式为：
(10)
则纵倾变化时阻力的计算公式为：
(11)
式(1)中，R0为正浮态的阻力值；η为纵倾阻力
增加率。系数a0、a1、a2、a3随纵倾角变化值(见表2)。纵倾的中间值的阻力增加率，可通过系数的线性内差获得。
图5艏倾变化弗劳德数-阻力曲线
图6艉倾变化弗劳德数-阻力曲线
表2 纵倾阻力增加率
|
纵倾角/(°) |
阻力增加率公式 |
|
-2.5 |
12.0126Fr3-6.8006Fr2+1.1656Fr+0.0750 |
|
-2.0 |
8.4315Fr3-4.4297Fr2+0.7258Fr+0.0412 |
|
-1.5 |
11.3996Fr3-9.5653Fr2+1.4459Fr-0.0573 |
|
-1.0 |
14.6160Fr3-9.5159Fr2+1.9497Fr-0.1319 |
|
-0.5 |
10.2630Fr3-6.4733Fr2+1.2630Fr-0.1048 |
|
0.5 |
0.4826Fr2-0.2346Fr+0.0209 |
|
1.0 |
0.7046Fr2-0.2315Fr+0.0401 |
|
1.5 |
1.2059Fr2-0.4054Fr+0.0760 |
|
2.0 |
1.5648Fr2-0.4767Fr+0.1085 |
|
2.5 |
1.7637Fr2-0.4960Fr+0.1476 |
5横倾变化时的阻力变化规律
横倾试验分为5组，分别为绕OX轴旋转5°、10°、15°、20°、25°。用弗劳德换算法将船模横倾变化阻力换算成实船阻力，并拟合成曲线(如图7所示)。
图7横倾变化弗劳德数-阻力曲线图
采用和计算纵倾阻力相似的方法，得到计算横倾变化时的阻力：
(12)
式(12)中，R0为正浮态的阻力值；η为横倾阻力增加率。系数b0、b1、b2随横倾角变化值(如表3所示)。
表3 横倾阻力增加率
|
横倾角/(°) |
阻力增加率公式 |
|
5 |
-0.0679Fr2+0.0226Fr+0.0216 |
|
10 |
0.5585Fr2-0.2941Fr+0.0962 |
|
15 |
0.5255Fr2-0.2959Fr+0.1109 |
|
20 |
0.4941Fr2-0.3268Fr+0.1551 |
|
25 |
0.6684Fr2-0.3821Fr+0.1891 |
6验证
为了验证文章数值计算的准确性，对船模进行几组不同航速下的数值计算，将阻力的数值计算结果与ITTC提供的船模水池试验结果进行对比验证(如表4所示)，Fr=0.28时自由面上的压力分布(如图8所示)。
从表4可以看出文章采用数值方法计算出的船舶阻力与船模水池试验结果吻合较好，相对误差在2%以内，可以满足工程需要。
表4 阻力的数值计算值和水池试验值比较
|
弗劳德数 |
数值模拟/N |
水池试验/N |
相对误差/% |
|
0.20 |
21.84 |
21.58 |
1.20 |
|
0.25 |
35.31 |
34.68 |
1.82 |
|
0.28 |
45.85 |
45.18 |
1.48 |
|
0.35 |
81.04 |
80.66 |
0.47 |
|
0.41 |
150.73 |
152.70 |
-1.29 |
图8自由面上的压力分布
为了验证增阻规律的准确性，采用工程上应用比较广泛的基尔斯法对吃水变化阻力公式进行验证(见表5)，该方法的基本思路是：如果新船的船型较母型船有所改变，则它们在相同弗劳德数时的剩余阻力系数分别为CR2和CR1，则：
CR2=J·CR1 (13)
表5 吃水变化阻力计算的验证
|
弗劳德数 |
增深值/m |
公式计算/105N |
基尔斯法/105N |
相对误差/% |
|
0.21 |
0.0155 |
4.1189 |
4.0991 |
0.48 |
|
0.0310 |
4.4283 |
4.4649 |
-0.82 | |
|
0.28 |
0.0155 |
7.0556 |
7.4456 |
-5.24 |
|
0.0310 |
7.6388 |
8.0270 |
-4.84 | |
|
0.41 |
0.0155 |
23.6340 |
22.9600 |
2.93 |
|
0.0310 |
25.7390 |
25.6940 |
0.18 |
式(13)中，J为阻力的修正系数。至于摩擦阻力系数可根据平板摩擦阻力公式计算。
从表5可以看出，这里选择Fr=0.21、Fr=0.28和Fr=0.413种弗劳德数情况进行验证文章公式计算结果与基尔斯法计算结果相对误差最大在-5%左右，从而证明了计算方法的可行性。
作者：吴明，王骁，杨波，石爱国，杨宝璋 来源：中国航海
