文 | 语君辞
编辑 | 语君辞
前言：
目前，模型规模的测试和模拟用于预测和优化船舶的设计，以确保最佳的动力性能。然而，在转换雷诺数差异的模型试验数据的过程中，尺度效应会导致很大的误差。
尽管有经验公式可以用来纠正这一错误，但尾流自适应螺旋桨设计和节能装置的节能效果分析存在严重问题。 随着流体力学计算机软件、硬件技术和数值模拟方法的最新发展，基于计算流体动力学（CFD）的全尺寸船舶动力性能正在积极研究中。
全尺寸数值模拟有助于预测节能设备的性能；并揭示了它们的节能机理。全尺寸船舶的动力性能计算已经逐步发展，然而，目前还没有对全尺寸船舶动力性能的仿真策略进行相关研究，尤其是对全尺寸仿真的关键问题，如时间步长和壁Y+值。
因此，我们将普通货船REGAL作为本研究的研究对象，通过对各种网格数、时间步长和壁Y+值的计算分析，构建了全尺寸船舶阻力模拟策略。
在此基础上，模拟了不同时间步长的自推进特性，并通过比较海试结果进行了验证。最后，分析了不同时间步长下的非定常螺旋桨推力和转矩，总结了不同要求下的全尺寸船舶自推进仿真策略。
一、数值模型与方法
1.普通货船REGAL
图1 普通货船REGAL的试验与仿真几何模型
2015年10月，LR利用从伊斯坦布尔到瓦尔纳的普通货船航程获得了一套全面的试验数据。2016年，LR船舶性能集团组织了一次全面的流体动力学研讨会，并公布了REGAL的完整试验数据。
图1显示了REGAL的测试和模拟模型，这是一艘单螺旋桨推进船，甲板上方有桅杆和上部结构。根据研讨会的建议，在本研究的模拟中考虑了上部结构的影响，省略了桅杆。通用货船REGAL的参数和海上试验条件如表1所示。
表1
2.计算域和边界条件
如图2所示，通过考虑船体的对称性，对半艘船进行了阻力模拟。船舶阻力计算的计算域为矩形，其长度、宽度和高度为4.5LPP×1。分别为5LPP×3LPP。从船到域的入口和顶部的距离为1，从侧面到1.5 LPP，从底部到2 LPP，和从出口到2.5 LPP。
计算域的入口和出口分别设置为速度入口和压力出口，侧面设置为对称平面。为了与研讨会的结果进行比较，船体和螺旋桨表面被设置为光滑的墙壁。
图2 阻力模拟域（左）和螺旋桨域（右）
在自推进模拟中，对整艘船进行了建模，计算域的宽度是阻力模拟域的两倍。如图2所示，螺旋桨旋转是通过滑动网格实现的，旋转域的直径为1.15D。
3.网格生成
STAR-CCM+的网格生成器用于划分计算域的网格。为了提高数值模拟的准确性，有必要细化计算域的特定区域。图3显示了仿真领域关键部分的网格结构。可以看出，对自由曲面的网格进行了细化，以更好地模拟自由曲面。
图3 网格分布：（a）开尔文波系，（b）船头，（c）边界层，（d）螺旋桨
同时，对船头、船尾、螺旋桨等表面曲率变化较大的区域也进行了细化。阻力（网格G1）和自推进计算中使用的网格参数列于表2中。
二、电阻模拟
基于REGAL进行了全尺寸船舶阻力仿真。 仿真的重点是网格数、时间步长和墙Y+值的影响，以获得仿真策略 。模拟的船速为14节，相应的弗劳德数为0.196。
表2
在预测多种工况的阻力之前，对纵倾和下沉进行了预测，并在计算中允许船舶绕Y轴的垂直运动和旋转。在随后的多条件阻力计算中，船舶被预先调整到上述姿态预测结果，然后进行固定姿态的阻力计算。
一方面，可以避免不同工况下可能出现的姿态差异，进而影响Y+和时间步长的单变量分析。另一方面，可以实现加速收敛和减少预测时间。预测的倾斜度和下沉度分别为0.0344°和0.194 m，与车间结果一致。电阻计算条件如表3所示。
表3
1.网格不确定性分析
采用笛卡尔剖切单元法生成由棱柱层网格和修剪网格组成的计算网格。为了消除Y+值的影响，在进行网格不确定性分析时，有必要保持棱柱层网格的拓扑形式不变。网格细化比率r被设置为1.2。表4列出了三组网格参数和CPU模拟次数。
表4
在劳氏船级社中，大多数给定的时间步长分布在0.0025–0.01 LPP/v的范围内。Y+值的范围统计如表5所示。基于此，在网格不确定性研究中选择了相对较小的时间步长和Y+值，分别为0.0025 LPP/v和100。
表5
基于上述三组网格，预测了REGAL在14节时的阻力，表6总结了数据。从表中可以看出，三组网格可以很好地预测船舶阻力，精度为0.5%。这表明这种网格拓扑结构可以很好的模拟全尺寸船舶阻力。
表6
表7总结了REGAL船舶阻力的网格收敛验证结果。从表7中可以看出，网格收敛比小于1， 这表明网格拓扑是单调收敛的 。网格不确定性为1.99%S1（其中S1是基于网格G1的模拟结果），这表明了令人满意的网格收敛行为。
表7
图4显示了基于不同网格模拟的自由表面波形的结果，并比较了Y/LPP=0.1时的纵波切割。WaterElevation是无量纲波高，定义为WaterElevation=h/LPP，h是波高。
可以看出， 基于三组网格的波形模拟结果非常相似，单个波峰和波谷之间只有细微的差异；此外，船头和船尾波浪系统发展良好，可以传输到远离船尾的位置2 LPP。
然而，具有最小网格数的G3波形模拟结果具有受网格拓扑影响的迹线。此外，随着网格数量的减少，船尾波浪系统的波高略有增加；从而改变船舶阻力。
图4 在不同网格下Y/LPP=0.1处切割的自由表面波形和纵波的比较
图5将基于本研究中使用的网格G1的Y/LPP=0.1时的波形曲线与劳氏船级社的结果进行了比较。可以看出，本研究中模拟的纵波切口与劳氏船级社中相同工况（即代码2、代码9、代码10.2、代码10.3、代码16、代码17和代码19）下的计算结果非常接近。
上述工况计算结果之间的波形峰值误差均在15%以内。这也表明了阻力模拟策略在本研究中应用的可行性。后续关于时间步长和墙Y+值的选择的研究将基于网格G1。
图5 本研究与Lloyd’s研讨会在Y/LPP=0.1时纵波切割的比较
2.关于设置时间步长的研究
为了研究适用于全尺寸电阻模拟的时间步长，选择了四种不同的时间步长：0.024、0.048、0.096和0.192s，分别对应于0.00125、0.0025、0.005和0.01LPP/v。在设置时间步长的研究中，Y+值也被设置为100。
表8
表8列出了四种不同时间步长的电阻模拟结果。在本次海试中，REGAL的阻力没有使用LR进行测量。因此，使用LR（2016）中相同工况（代码2、代码9、代码10.2、代码10.3、代码16、代码17、代码19）的平均阻力值来验证本研究中的模拟结果。
表8总结了LR车间四个时间步长的电阻模拟结果和结果的总体误差在5%以内。与其他两个小的时间步长相比，0.005和0.01LPP/v时间步长期间的船舶阻力显著增加。
图6 REGAL在不同时间步长下的船舶阻力分量
图6显示了船舶阻力分量随时间步长的变化。从图中可以看出，船舶摩擦阻力几乎不随时间步长的变化，剩余阻力受时间步长的影响很大。当t=0.0025 LPP/v增加到0.005 LPP/v时，剩余阻力突然变化，这导致在0.005和0.01 LPP/v的时间步长下船舶阻力的误差更大。
为了分析船舶剩余阻力变化的原因，分析了不同时间步长下的波形；结果如图7和图8所示。从图中可以看出，时间步长对波形的演变和捕获具有更大的影响。
图7 REGAL在不同时间步长下的自由表面波形
当时间步长为0.00125 LPP/v和0.0025 LPP/v时，总体波形和纵波切割模拟结果相对接近。然而，当时间步长为0.005和0.01LPP/v时，会出现自由表面捕获问题，并且波形不能正确和充分地发展，这会影响造波电阻。
图8 不同时间步长之间Y/LPP=0.1处纵波切割的比较
CFL数是评估数值计算的一个重要条件，其在自由表面上的分布如图9所示。可以看出，自由表面上的CFL数与该位置的网格大小密切相关（图3（a））。最小网格处的CFL数最大，在0.0025 LPP/v和0.005 LPP/v的时间步长下分别约为0.8和1.6。
图9 不同时间步长下自由表面上的CFL数
图10显示了REGAL在不同时间步长下的标称尾流场（无轮毂）。模拟结果清楚地捕捉到了该船型钩形速度等值线的存在。
然而， 当时间步长增大时，尾流场的速度分布之间会出现一些差异 。提取了螺旋桨圆盘上不同螺旋桨半径下的轴向速度；结果如图11所示。图11中横坐标θ的定义如图10（a）所示。
图10 REGAL在不同时间步长下的标称尾流场
从图11可以看出，当时间步长为0.00125和0.0025 LPP/v时，速度分布几乎相同，但当它们为0.005 LPP/v或0.01 LPP/v，r/r=0.3和r/r=0.5时的速度值显著增加。螺旋桨盘上的平均无量纲轴向速度和不同时间步长下的不同半径如图12所示。
当时间步长大小为0.0025、0.005 LPP/v或0.01 LPP/v时，螺旋桨盘上的平均轴向速度与当时间步长尺寸为0.00125 LPP/v的平均轴向速率之差分别为0.11%、2.87%或3.24%。
图11 不同半径下轴向速度的周向分布
基于对不同时间步长的船舶阻力分量、波形和尾流场的分析，全尺寸船舶（低速船型）阻力模拟的时间步长应确保自由表面的CFL数不大于1，建议值为0.0025 LPP/v。
3.选择墙Y+值
当外推船舶阻力时，摩擦阻力的预测通常基于ITTC1957线。随着全尺寸仿真的发展，许多研究人员确定了全尺寸船舶摩擦阻力仿真值与ITTC1957线之间存在一定误差。
图12 不同时间步长下螺旋桨盘（左）和不同半径（右）上的平均轴向速度
为了分析适用于全尺寸模拟的船面壁Y+值，参考LR车间中选择的壁Y+的范围（表5），选择了六个不同的目标壁Y+数值（棱镜层厚度不变，同时调整数量（10-14）和拉伸（1.05-1.2）边界层网格）：60、100、200、300、500和800。
图13 REGAL表面上不同的墙Y+值
在壁Y+值收敛研究中，模拟基于网格G1，时间步长为0.048 s，对应于0.0025LPP/v。REGAL表面上不同的壁Y+分布如图13所示。可以看出，船体表面上的Y+值总体上略小于目标Y+值。Y+值在水线以下船体表面的分布范围如表9所示。
图14显示了不同Y+条件下半船阻力的时间历程。可以看出，电阻倾向于在100 s左右收敛，每个工况的模拟物理时间为200 s，电阻取100–200 s间隔的平均值。
表9
对不同壁Y+值的船舶阻力分量进行了监测，剩余阻力和摩擦阻力随壁Y+的变化如图15所示。可以看出，当壁Y+值为60、100或200时，残余电阻几乎稳定；当壁Y+值大于300时发生轻微的增加。
当其低于300时，模拟摩擦阻力基本保持不变，这与Katsui线和Grigson线都很好地一致。当Y+为500–800时，船舶的摩擦阻力在一定程度上增加。
图14 不同Y+值下的半船阻力的时间历程
当壁Y+值为100、200、300、500或800时，船舶的总阻力与壁Y+为60时的总阻力之差分别为0.15%、0.23%、1.07%、1.43%或1.86%。
图15 REGAL具有不同壁Y+值的船舶阻力组件
图16表示不同壁Y+值下REGAL上的表面摩擦系数。可以看出，船头和船尾的表面摩擦系数随着壁面Y+值的增加而增加，特别是当Y+大于200时。这与不同壁Y+下边界层流动模拟的差异有关。
图16 REGAL在不同壁面Y+值下的表面摩擦系数
图17显示了不同壁Y+值下船尾压力系数的分布（Song等人，引文2018）。可以看出，当Y+为300、500和800时，船尾的压力在一定程度上降低，这增加了船头和船尾之间的压力差，从而 增加了粘性压力阻力 。
图17 REGAL在不同壁Y+值下的压力系数
图18和图19分别表示在不同壁Y+条件下，Y/LPP=0.1处切割的自由表面波形和纵波。可以看出，在不同Y+值条件下，波形的计算结果非常相似，只有在弓形波形上有一些差异，这也与边界层流动的演变有关。不同工况下表面摩擦系数、压差和波形的变化解释了图15中船舶阻力的差异。
图18 REGAL在不同壁Y+值下的自由表面波形
图20显示了REGAL在不同壁Y+值下的标称尾流场。从图中可以看出，当壁Y+值大于300时，涡流结构的形状显示出一定的差异。提取了螺旋桨圆盘上不同螺旋桨半径下的速度；结果如图21所示。
图19 不同壁Y+条件下Y/LPP=0.1处纵波切割的比较
可以看出，随着r/r的降低，不同条件下的速度差异更大。换句话说，越靠近螺旋桨盘的中心，Y+值对速度的影响就越大。Y+为500和800时的速度分布与其他工况明显不同。
为了定量分析壁面Y+值的影响，获得了不同半径（不安装假毂）的螺旋桨盘的平均无量纲轴向速度；结果如图22所示。当壁Y+小于200时，螺旋桨盘表面和不同半径的平均速度保持稳定。
图20 REGAL在不同壁面Y+值下的标称尾流场
当壁面Y+大于200时，r/r=0.3时的平均速度增加，r/r=0.5和r/r=0.7时的平均速率减小，螺旋桨盘上的平均速率在一定程度上减小。
图21 不同半径下轴向速度的周向分布
当壁面Y+值为100、200、300、500或800时，螺旋桨盘上的平均轴向速度与当壁面Y+值为60时的轴向速度之差分别为0.08%、-0.03%、0.61%、1.23%或1.84%。
图22 不同壁Y+值下螺旋桨盘（左）和不同半径（右）上的平均轴向速度
根据本节阻力、桨盘速度和流场的模拟结果，进行全尺寸船舶模拟时，Y+值在200以内是合适的。
三、自推进模拟
在阻力仿真策略的基础上，对自推进策略进行了研究。船体Y+值和网格拓扑与阻力模拟中保持相同。自推进模拟的网格数为3022万。根据劳氏船级社给出的海试结果，当螺旋桨转速为106.4转/分时，船舶速度为13.05节（对应的弗劳德数为0.182）。
图23 不同时间步长的模拟螺旋桨推力和扭矩
在本研究中，选择的五个时间步长分别为0.000783 s（0.000038 LPP/v）、0.001566 s（0.000076 LPP/v。本节中的图23-25显示了1°/步的所有模拟结果。
图24 不同时间步长下的模拟推力推导分数
表10和图23比较了推力和扭矩的计算结果和试验结果。从表11可以明显看出，推力和扭矩对时间步长的敏感度相对较低。对于所有五个时间步长，推力和扭矩的误差都小于4%。
通过分析，推力和扭矩随着时间步长的增加而减小。在0.5–2°/步的情况下，推力和扭矩的变化相对较小，但在4°/步时会突然下降。
图25 不同时间步长下滑动网格界面上的CFL数
需要注意的是，不同螺旋桨转速下的阻力为370–375 kN，不同时间步长下的阻力和推力之差在−2.48–3.94%之间。这两个极限相对接近，这表明该船处于自推进状态。
图24将本研究中推力推导分数的模拟结果与劳氏船级社在不同时间步长下的结果进行了比较。可以看出，对于相同的时间步长，本研究的结果与劳氏船级社的结果一致。
表10
Code8.1使用虚拟磁盘模型来模拟自推进，并在相对较长的时间步长内获得了可接受的推力推导分数；这就是虚拟磁盘模型的优点。通过模拟基于螺旋桨旋转的自推进，当时间步长小于2°/步（0.000152 LPP/v）时，推力推导分数更可靠。
表11
图25显示了滑动网格界面上的CFL数量分布。可以看出，2°/步和4°/步的CFL数分别约为0.7和1.4。由于两种时间步长下的推力和扭矩之间存在相当大的差异， 应确保滑动网格界面上的CFL数不大于1。
图26 阻力（a）和自推进（b）条件下REGAL上的压力分布
图26和27分别显示了不同条件下REGAL上压力系数和表面摩擦系数的分布。可以看出，阻力和自推进条件之间的两种分布的差异主要在船尾。这是因为螺旋桨的运行导致其前方的水速度增加，从而降低了压力系数并增加了船尾的摩擦系数。
图27 阻力（a）和自推进（b）条件下REGAL上的表面摩擦系数
图28显示了阻力和自推进条件下螺旋桨盘下游横截面X/D=0.2315处的速度分布。在阻力模拟中，该型船的钩形速度剖面被清晰地捕捉到。自推进模拟的速度等值线图清楚地表明了尾流的不对称性。此外，轮毂涡流引起的低速区域也在本节中得到了很好的显示。
图28 阻力（左）和自推进（右）条件下螺旋桨盘下游横截面X/D=0.2315的速度等值线图
图29描述了螺旋桨涡流结构，其中等表面由Q=0.5（||Ω||2−||S||2）=50表示，其中S和Ω分别表示应变率和涡度率大小。
从图中可以清楚地观察到螺旋桨尾流的非定常特性。因此，在自推进模拟中，许多研究人员除了关注自推进因素外，还关注螺旋桨的非稳态特性。
螺旋桨激振力是船尾结构振动的主要来源，也是导致船体结构损坏的激振力引起的振动。基于快速傅立叶变换计算了不同时间步长下的非定常推力和转矩；结果如图30所示。
图29 螺旋桨和方向舵周围的旋涡结构
图30显示，REGAL船后面的螺旋桨在轴频率（SF，1.773 Hz）和叶片通过频率（BPF，7.093 Hz）下具有波动振幅。然而，对于每个时间步长大小，在频域结果中存在一定的差距。SF激振力主要由机械不平衡引起，这是由于螺旋桨的安装过程或船尾不对称造成的。
图31和图32分别比较了本研究中模拟的SF和BPF以及劳氏船级社的波动幅度。可以看出，在某些工况下，本研究中的模拟结果与车间结果具有一定程度的一致性。
劳氏船级社中代码3、代码9和代码16的振幅似乎不取决于时间步长，这可能是由于数值模拟策略之间的差异。
图30 不同时间步长的推力（a）和转矩（b）的频域结果
对于本研究中的模拟结果，激振力的振幅显示出时间步长的影响，尤其是当8°/步长的计算结果不同时，推力的SF幅度和转矩的BPF幅度随着时间步长的大小而增加。同时，推力的BPF幅度和SF幅度随着时间步长的增加而减小。
此外， 推力和转矩在0.5°/步和1°/步时的波动幅度具有很高的相似性 。因此，在求解螺旋桨诱导激振力时，建议时间步长不应大于1°/步（0.000076 LPP/v）。
图31 SF不同时间步长下的推力（a）和转矩（b）的幅度
结论
从这项研究中可以得出以下结论，Katsui和Grigson线与模拟的摩擦阻力吻合良好，可用于验证摩擦阻力模拟值。
壁Y+值将影响边界层流动的捕获、船舶压力场的演变和波形，从而影响船舶阻力。根据REGAL的模拟结果，当模拟全尺寸船舶模拟的动力性能时，Y+值在200以内是合适的。
图32 不同时间步长下推力（a）和扭矩（b）的BPF振幅
全尺寸船舶阻力模拟的波形演变很大程度上受时间步长的影响。电阻模拟的时间步长应满足自由表面CFL数不大于1的要求；推荐值为0.0025 LPP/v。
目前，用于全尺寸船舶海上试验的开放测试数据集很少。REGAL公布的劳氏船级社试验数据可以为全尺寸船舶的动力性能模拟提供验证数据。此外，它还可以增加研究人员对全尺寸船舶CFD的信心。
参考文献：
1.阻力、动力、耐波性和5415机动的模型和全尺寸URANS模拟，船舶研究杂志，斯特恩，2009年
2.KCS KRISO集装箱船采用离散螺旋桨的全尺寸自推进计算，计算机和流体，卡斯特罗，2011年
3.船舶阻力部件和形状因子的尺度效应，海洋工程，德米雷尔，2020年
4.用于性能预测的精确光滑摩擦线，皇家造船工程师学会汇刊，格里森，1993年
5.全尺寸计算流体动力学模拟对船舶的重要性，国际造船与海洋工程杂志，霍奇基奇，2013年