深远海养殖工船最小推进功率研究

2021-12-18陈作钢

海洋工程 2021年6期

吕磊，陈作钢, 3，代燚, 3

(1. 上海交通大学海洋工程国家重点实验室，上海 200240; 2. 上海交通大学船舶海洋与建筑工程学院，上海 200240; 3. 高新船舶与深海开发装备协同创新中心(船海协创中心)，上海 200240)

深远海养殖近年来获得了广泛的关注。在深远海海域开展水产养殖，现代化、工业化、机械化的养殖装备是关系到养殖成败的关键因素[1]。深远海养殖工船作为深远海养殖的一个发展方向，集成了繁育、养殖、加工、冷冻冷藏等鱼货物供给的一条龙功能，有效地推进海洋渔业的转型发展[2]。养殖工船通常系泊工作于水质和风浪条件合适的深远海海域，并汲取海水至养殖液舱[3]。当收到恶劣海况预报时，养殖工船可依靠其自航能力提前撤离危险区域，达到防风避台的目的。为了在设计初期确保养殖工船具有恶劣海况下的生存能力，最小推进功率的校核是必不可少的关键设计步骤。

2009年国际海事组织(IMO)海洋环境保护委员会(MEPC)提出船舶能效设计指数(EEDI)[4]，旨在提高船舶能效水平，减少船舶业的温室气体排放。船舶的设计建造者为了满足EEDI的要求，较为有效的方法就是降低船舶装机功率和航速。这种方法同时也使得船舶在恶劣海况下可能会出现因功率不足而影响操纵性能的情况，容易导致航行安全问题。针对这一问题，MEPC 第64届会议通过了“恶劣海况下船舶维持操纵性的最小推进功率临时导则”[5]，并于第65届会议对其进行了完善，即MEPC.232(65)决议(以下简称为决议)。决议明确要求船舶主机推进功率应不小于在恶劣海况下保持船舶操纵性所需要的最小推进功率。

根据决议的说明，最小推进功率有两种评估方法:等级1的线评估法和等级2的简化评估法。线评估方法较为严苛，对于肥大型船舶来说一般难以满足[6]。简化评估方法则通过计算恶劣海况下船舶在迎风迎浪时以最小前进速度航行时的阻力，进而间接得到所需的最小推进功率。这种方法将阻力分为4部分:静水阻力、空气阻力、附体阻力和波浪增阻，并给出了恶劣海况的风浪参数。其中，前3项阻力分量都可由决议推荐的经验公式计算得到，波浪增阻的计算则需要通过模型试验或其他可靠的方法得到。综合来看，波浪增阻的准确预报是最小推进功率校核的关键。

除了船模试验外，波浪增阻的研究方法有理论计算法和经验公式法。理论计算法在势流理论框架下分为远场法和近场法。远场法基于船体绕辐射波能在远场总动量变化得到定常二阶力，近场法是通过直接积分船体湿表面二阶压力得到[7]。经验公式方法从理论推导出发，通过试验数据对公式局部修正。这些方法对船型和航速都有一定的适用范围，且无法考虑各种复杂流动现象导致的强非线性因素。计算流体力学方法(computational fluid dynamics, 简称CFD)不仅可以充分考虑流体的黏性作用，还可以将流体的非线性因素计算在内。CFD方法与试验流体力学(experimental fluid dynamics，简称EFD)方法相比，成本更低，又可以给出船体周围复杂流场的信息，近年来成为求解船舶与海洋工程水动力问题的一种重要手段。

沈志荣等[8]利用naoe-FOAM-SJTU求解器计算分析了Wigley III型船模在迎浪中的运动响应及波浪增阻，并与试验数据和势流计算结果进行了对比。方昭昭等[9]基于FLUENT求解器建立数值波浪水池，就规则波顶浪中航行船舶的运动与波浪增阻进行数值模拟计算，计算结果和试验值吻合较好。曹阳等[10]使用重叠网格方法对KVLCC2船型迎浪航行进行了数值模拟研究，证明了CFD方法能够在全波长范围内更准确地预报船舶在波浪上增阻。这些研究使用不同的数值模拟方法计算了各类船型在模型尺度下的波浪增阻问题，表明使用CFD方法计算波浪增阻问题是切实可行的。但针对养殖工船最小推进功率校核的实际问题，模型尺度的计算将会带来尺度效应的影响，若进行尺度修正又会带来误差，更合适的方法是直接进行波浪增阻的实尺度数值模拟计算。

以十万吨级深远海养殖工船为研究对象，使用STAR-CCM+商业求解器，以实尺度CFD方法得到了其在规定海况下的波浪增阻计算结果，并对其最小推进功率进行了计算校核。为验证计算方法的可靠性，进行了该船型在迎浪航行状态下的EFD船模试验，并与几组模型尺度下的CFD数值计算结果进行了对比。最终结果表明实尺度CFD数值模拟方法可以用于养殖工船最小推进功率的校核，该船型使用简化评估方法校核最小推进功率更容易满足规范要求。

1 EFD试验准备和CFD计算的前处理

1.1 模型参数

研究对象为十万吨级深远海养殖工船，主要船型参数见表1。

表1 主要船型参数Tab. 1 Main ship parameters

EFD模型试验及与其对照的CFD模型尺度计算所用的船模缩尺比为1∶80，船体模型严格依据型线图按上述缩尺比加工，模型几何尺寸误差不超过2 mm，吃水误差不超过1 mm。直接用于最小推进功率校核的CFD实尺度算例则使用实际船型参数。图1为试验模型及CFD计算模型。

图1 试验模型和计算模型Fig. 1 Test model and calculation model

1.2 模型试验和模拟计算海况

根据MEPC.232(65)决议要求，最小推进功率等级2的简化评估方法基于这样的原则:如果船舶具有足够的装机功率在迎风迎浪中以一定的前进速度移动，则船舶也能够在任何其他方向的波浪和风中保持航向，即满足最小前进速度要求也就是满足航向保持要求[11]。

首先需要确定船舶在迎风迎浪时保证航向稳定性的最小航速。最小前进速度定义为:

VS=max{2.058 m/s,Vref-5.144(AR%-0.9)}

(1)

(2)

ALS=LppTm[1+25(BWL/Lpp)2]

(3)

式中:VS和Vref分别为最小前进速度和最小基准速度，Vref由表2线性插值得到。AR%为舵面积AR占船体侧面浸水面积ALS的百分比，又称为舵面积系数。

表2中AF为船体纵向受风面积，AL为船体横向受风面积。经计算，此养殖工船最小前进速度VS为2.953 m/s，对应的弗劳德数为0.061。

表2 最小基准航向保持速度定义Tab. 2 Definition of minimum reference course keeping speed

设置6组不同波长的波浪参数进行EFD模型试验，并选取其中3组参数进行CFD模型尺度计算并与之进行对比，以验证CFD实尺度计算方法的可靠性。为试验的方便，采取等波高的方法设计波浪参数。再设置5组实尺度CFD计算海况进行数值计算，采用等波陡的方法设计波浪参数，得到的波浪增阻计算结果直接用于最小推进功率的校核。所有的EFD模型试验和CFD模拟计算波浪参数如表3所示。

表3 EFD模型试验和CFD模拟计算的波浪参数Tab. 3 Wave parameters of EFD model test and CFD simulation

1.3 EFD模型试验设置

EFD模型试验在风洞循环水槽实验室进行，图2为风洞循环水槽试验装置总体布局及循环水槽全景。试验中使用的主要试验设备有:循环水槽、随动式4自由度适航仪、日章电机LMC-1541单分量天平、惯量调节架等。其中循环水槽试验段长8.0 m，宽3.0 m，水深1.6 m，最高流速3 m/s；配备4座叶轮泵，可进行分层流试验；具有造波、消波能力，可生成波高0.1 m的规则波、不规则波和内波；试验段上游装有表面流加速装置，弥补了边界层速度亏缺；配备了水自动过滤装置和水位自动调节装置。该试验设施可完全满足试验要求。

图2 主要试验装置全景Fig. 2 Panorama of the main test device

试验原理需保证模型与实船在几何学上保持几何相似，运动学上保持弗劳德数相似。试验进行之前，先使用惯量调节架对试验船模的重心位置和惯量进行调节，使其与设定参数相一致。试验时将模型置于水槽中线处，与来流方向一致。水槽测试段前端造波机运动，形成浪流联合条件。随动式4自由度适航仪自动控制调节电机扭矩大小，使得模型在波浪中随动，得到模型垂荡和纵摇的运动幅值时历变化。

1.4 CFD模拟计算设置

1.4.1 控制方程和物理模型

模拟计算的控制方程为连续性方程和雷诺平均Navier-Stokes方程(RANS):

(4)

(5)

(6)

(7)

(8)

其中，体积分数aq表示单元内第q相流体占总体积的比例分数。流体体积的输运方程为:

(9)

为了求解雷诺应力项，引入可实现K-Epsilon模型来封闭控制方程。造波方法采用边界输入法，根据Stokes五阶波的解析解，在计算域的边界直接设定水质点速度和波面瞬时高度。定义纵向坐标x对应于波浪传播方向，垂向坐标z为铅垂向上的方向。根据Stokes五阶波的解析解，在计算域的边界分别设定水质点x方向速度ux、z方向速度uz和波面瞬时高度ηw为:

(10)

(11)

(12)

式中:ω、d和k分别为圆频率、水深和波数。其余各项系数可参考文献[12]。在出口边界设置宽度为两倍波长的消波区，使用阻尼消波法消除反射波的影响。此外，使用HRIC(high-resolution interface capturing)高分辨率接触面捕捉方法提高自由液面的捕捉精度。

使用有限体积法(FVM)将控制方程离散化。所有数值插分、微分和积分都是基于二阶近似。所得的耦合方程组被线性化，并由隐式非定常求解器求解。计算时使用SIMPLE方法以获取压力值并校正速度。

1.4.2 计算域和网格设置

为减少计算量，根据船舶在迎浪状态下流场特征和运动响应的对称性，只取船中纵剖面一侧的流域进行研究。坐标原点设在船体重心处，计算域大小及各边界条件如图3所示。

图3 计算域设置Fig. 3 Computing domain settings

使用重叠网格方法及DFBI (displaying dynamic fluid body interaction)模型实现对船体迎浪时的垂荡与纵倾运动的模拟。重叠区域和背景区域的网格采用切割体非结构网格形式，基础尺寸设为1/10LPP；船体边界层采用棱柱层网格进行捕捉，棱柱层增长率设为1.2。模型尺度计算时需调整近壁面第一层网格的高度以使船体表面Y+在30至100之间，确保壁面函数方法可以较好地求解出近壁面区域的流动情况。实尺度模拟计算的雷诺数相较模型尺度来说要大715倍，为避免第一层网格高度过小，放宽Y+到300至500范围内，并适当增加棱柱层数。在背景网格区域设置多层嵌套的加密区，使网格尺寸平滑过渡。在船艏和船艉桨轴等流场变化较快的区域设置适当的加密区，以精确捕捉流场特征。此外，为保证波浪的多周期稳定传播，还需对自由液面处网格进行加密。在波浪起伏范围内，Z方向网格尺寸设为0.05H，H为波高，保证20个网格的解析度；X方向的网格尺寸设为0.4H；Y方向网格尺寸设为4H，并设置多层加密区进行过渡。由于自由液面处加密区以及棱柱层设置的不同，模型尺度算例的网格数在9×106至12×106之间，实船尺度算例的网格在11×106至17×106之间。图4为实船尺度算例中纵对称面船体附近的网格示意，船体表面Y+的值也显示在图中。

图4 实尺度CFD计算网格设置和船体表面Y+分布Fig. 4 Full scale CFD calculation grid setting and hull surface Y+ distribution

2 EFD试验结果和CFD计算的后处理

为了定量分析和比较计算结果，采用傅里叶级数展开的方式对数据进行分析处理。某一物理量的时域φ(t)可表示为:

(13)

式中:ωe为船模遭遇频率，γi为第i阶系数对应的初始相位。将船舶运动中垂荡和纵摇时历曲线依据公式(13)进行傅里叶级数展开，可以得到不同阶数的幅值以及相位。取傅里叶级数中的一阶响应系数，可以将船舶的垂荡运动z和纵摇运动θ表示成如下形式:

z=zacos(ωet+εz)

(14)

θ=θacos(ωet+εθ)

(15)

(16)

(17)

(18)

其中，a是入射波波幅；k=2π/λ是波数；Rwave为规则波中的波浪增阻，Rwave=FX,wave+FX,calm，FX,wave和FX,calm为船体在波浪和静水中阻力的平均值。其中垂荡和纵摇以船体重心为参考点。

(19)

式中:FX为沿X方向总阻力，S为船体在静水中的湿表面面积。此外，在中纵剖面船艏2LPP处设置数值探针，以监测入射波波幅随时间的变化情况，确保波浪参数的准确。

2.1 EFD阻力试验结果和模型尺度CFD计算结果

图5 EFD与CFD结果对比Fig. 5 Comparison of EFD and CFD results

3组模型尺度CFD计算结果与EFD试验结果相比，在λ=0.2LPP的入射波激励下3个参数的计算值略大于试验值。随着入射波波长的增大，CFD计算结果均小于试验值。除了λ=0.2LPP入射波对应的波浪增阻系数以外，其余数据点计算值和试验值的误差均在10%以内，符合工程精度要求，说明使用CFD方法计算船在波浪中的运动响应和阻力是可行的。

2.2 实尺度CFD计算结果

图6 入射波波长0.5LPP算例监测结果Fig. 6 Monitoring results of the incident wave wavelength 0.5LPP case

从监测结果可知，当计算稳定后，各参数时历值呈周期性变化。其中入射波最大波幅的平均值为1.061 m，与斯托克斯五阶波理论值1.085 m相比小2.30%，衰减幅度在可接受范围内。升沉和纵摇运动传递函数变化趋势接近线性，总阻力系数时历曲线则有明显的非线性特征。图7为自由液面波形图，由于船型较为肥大且航速较低，在船艏可看出明显的堆水现象，显示出非线性特征。同时在重叠网格和背景网格的交界处，可以看出流场有一定的不连续，这是由于重叠网格的插值和网格大小的变化造成的，说明在网格的设置方面还有改进空间。

图7 自由液面波形Fig. 7 Free surface wave pattern

将所有算例的波浪增阻系数σaw进行汇总，如图8所示。

图8 实尺度CFD波浪增阻系数计算结果 Fig. 8 Full-scale CFD wave resistance increase coefficient calculation results

由结果分析可知，波浪增阻峰值对应的入射波波长在0.5Lpp附近，增加此波长附近工况有利于提高波浪增阻的计算精度。CFD实尺度计算结果符合客观规律，较为可靠。

2.3 最小推进功率校核

MEPC.232(65)决议中的等级1线评估法给出的最小功率线值PLV可由式(20)计算:

PLV=a×(DWT)+b

(20)

式中：DWT为载重吨。对于养殖工船，参数a和b分别取0.068 9和3 253.0，计算得到PLV为11 610 kW。

等级2的简化评估法将船舶在恶劣海况下螺旋桨的推力T分为4部分，分别为静水阻力Rcw、空气阻力Rair、附体阻力Rapp和波浪增阻Raw，并考虑推力减额分数t:

T=(Rcw+Rapp+Rair+Raw)/(1-t)

(21)

其中，静水中的裸露船体阻力Rcw和空气动力学阻力Rair均可由推荐的经验公式计算得到；由于附体阻力Rapp一般占比较小，可由简单等效方法计算。波浪增阻Raw需要由水池试验或等效的可靠方法获得。研究采用实尺度CFD方法计算得到波浪增阻系数进行谱分析得到。

2.3.1 确定恶劣海况参数

为计算十万吨级深远海养殖工船在极端情况下的生存能力，需要对恶劣海况进行规定，并计算在此海况下养殖工船的阻力和推力性能。最小推进功率导则所规定的恶劣海况参数需根据表4来确定。

表4 恶劣海况参数Tab. 4 Severe sea conditions parameters

根据此深远海养殖工船船长线性插值，得到对应的恶劣海况有效波高HS为5.230 m，谱峰周期TP在7 s至15 s之间，平均风速Vw为18.406 m/s。

2.3.2 计算静水阻力Rcw

静水阻力的计算根据相当平板理论，采用1957年ITTC公式来计算静水阻力Rcw，计算公式如下:

(22)

(23)

(24)

式中:k为船体形状因子，可由经验公式(24)得到；S为船体湿表面积；CF为摩擦阻力系数；CB为方形系数。

2.3.3 计算附体阻力Rapp

由于附体阻力占总阻力比例较小，一般可以使用简单等效方法计算。将所有附体的湿表面积Sapp加入公式(25)船体的湿表面积中，用与静水阻力相同的计算方法得到附体阻力的大小，即:

(25)

2.3.4 计算空气阻力Rair

空气阻力Rair可以采用式(26)计算:

(26)

式中:AF为船体正面受风面积；Vrw为船与风的相对速度，Vrw=VS+VW；Cair为空气阻力系数，可以由试验或者经验给出。这里采用经验公式得到Cair为:

(27)

式中:AL为船体横向受风面积；C为船舶横向受风面积的形心距船中的位置，靠近船艏为正。

2.3.5 计算波浪增阻Raw

规定恶劣海况不规则波浪环境下的波浪增阻Raw可根据不同频率下规则波中波浪增阻Rwawe通过谱分析方法得到，而规则波中波浪增阻Rwawe可根据公式(18)由波浪增阻系数σaw得到，所以Raw可表示为:

(28)

式中:a为入射波波幅；S(ω)为波浪谱密度函数。最小推进功率导则中推荐波浪谱密度函数采用1978年第15届ITTC提出的Jonswap平均波能谱，其定义和相关系数表达式为:

(29)

(30)

ω0=4.85/T1

(31)

T1=4.85TP/(2π)

(32)

(33)

其中，HS为有效波高，T1为平均周期，TP为谱峰周期。有效波高HS和谱峰周期TP由上述恶劣海况定义表得到。

2.3.6 确定最小推进功率

当螺旋桨所需推力T确定后，可由式(34)通过螺旋桨敞水特征曲线插值得到进速系数J，具体为:

(34)

式中：KT(J)为推力分数，DP为螺旋桨直径，进速ua=VS(1-w)，w为伴流分数。研究螺旋桨敞水特性曲线参考MAU4-40桨标准图谱，进速系数J根据曲线插值得到。伴流分数参考最小推进功率导则的推荐值选为0.35，推力减额系数保守估计为t=0.7w。

推力系数KT、转矩系数KQ和进速系数J分别定义为：

(35)

(36)

(37)

有了进速系数J，螺旋桨转速n、收到功率PD以及主机功率PS可由式(38)到式(40)计算得到：

n=ua/(JDP)

(38)

(39)

PS=PD/ηs

(40)

式中：ηs为传输效率，对于尾机型船可取0.98。由于柴油机转矩转速特征曲线对可用功率的限制，还需结合曲线检验转矩是否超过最大转矩，以最终确定所需的最小装机功率值。计算结果见表5。

表5 最小推进功率简化评估法计算结果Tab. 5 Calculation results of the minimum push power simplified evaluation method

由计算结果可知，对于此深远海养殖工船来说，最小推进功率校核时所规定的恶劣海况下波浪增阻占总阻力比例较大，在28%至57%之间。谱峰周期TP在8 s至10 s之间所对应的海况较为危险，此时总阻力和所需功率较大，转矩达到最大转矩的63.9%。等级2简化评估法得到的最小推进功率最大值为6 832.870 kW，为等级1线评估法给出的最小功率线值PLV的58.9%，在装机功率不变的情况下，使用简化评估法进行校核更容易达到规范要求。