旋臂水池试验数值仿真影响因素分析

(深海载人装备国家重点实验室,江苏 无锡 214082)

新型水下潜器外形各异，在设计研发过程中水动力设计无母型参考，潜艇相关规范也不适用,只能依靠水池试验和CFD仿真。在整个研发周期内，尤其是在初期多方案选型时，动稳性评估都是必不可少的关键环节。尽管旋臂水池试验能提供稳定性评估需要的旋转导数，但是，试验周期长，费用高，常常无法及时地用于初期的线型优化设计中。国内对旋臂水池试验的数值仿真方法主要采用多参考系、动量源项法和滑移网格法。前两种都是采用稳态的方法来模拟旋臂水池中物体的运动过程，而滑移网格法采用瞬态模拟，理论上与物理实际的相近程度更高。虽然该方法计算量较大，但是，当前计算能力已经达到了较高的水平，计算能力已经不是突出的问题。因此，考虑采用某商业软件使用滑移网格法对suboff标模进行数值仿真。对旋臂水池的仿真研究[1-3]未见具体的数值计算设置。因此，考虑对计算域大小、边界层厚度、湍流模式、临界雷诺数等因素对仿真的影响进行系统性分析，分析探讨较为准确的旋臂水池数值模拟方法。

1 研究模型及符号定义

采用的计算模型为美国国防高级研究计划署(DARPA)[4-5]提出的SUBOFF AFF-1和AFF-8。SUBOFF主体是水滴形的轴对称体，全附体还包括指挥台围壳和尾翼，见图1。

图1 SUBOFF模型

采用潜艇操纵性规范中规定坐标系、量纲一的量化及表达方法。为方便与试验结果对比，参照文献中的规定，取特征长度L=4.261 m。坐标系原点位于浮心位置处。浮心位于主体轴线上，距艏2.013 m。计算中的工况参数也与试验中的相同，水温18℃，密度ρ=998.55 kg/m3，动力黏度取μ=0.001 053 kg/(m·s)。三向力和力矩系数分别以Cx、Cy、Cz、Cmx、Cmy和Cmz表达，均以特征长度L、旋转线速度V量纲一的量化。无因次角速度r′=rL/V。

2 数值计算方法

2.1 控制方程与湍流模型

不可压缩黏性流体的连续性方程和RANS方程可写成如下形式。

(1)

(2)

式(1)、(2)是不封闭的，需要寻求补充关系(湍流模型)使问题封闭，采用realizablek-ε或SSTk-ω湍流模型进行湍流封闭[6]。

2.2 求解区域与边界条件

为减小网格数量，将整个环形计算区域分割成2个同心的圆环，见图2。大圆环为静止区域，内部小圆环为旋转区域。艇体处于小圆环内，跟随内部圆环围绕轴心旋转，实现艇体的旋转模拟，动静区域通过内外圆环交界面传递数据。旋转半径分别取10、13、15、20和25 m。

图2 计算域及计算域参数示意

2.3 网格建模与生成

采用非结构化网格划分方法对光体和全附体线型进行网格划分，见图3。在边界层附近划分棱柱层网格，并对翼、围壳、艇体周围等区域进行局部加密。全附体计算时，网格总数约200万。

图3 圆环计算域网格划分

2.4 差分格式及算法

采用有限体积法(FVM)离散控制方程和湍流模型。对流项采用二阶迎风格式，扩散项采用中心差分格式。对时间项采用一阶离散格式，每个时间步旋转0.5°，每步内部迭代6次。对压力速度方程采用分离式计算方法。

3 影响因素分析

3.1 外部计算域圆大小的影响

直航仿真中，为了减小远场边界的干扰，通常需要将艇体距离四周边界2倍艇长以上。而旋臂水池仿真中艇体旋转导致流动更加复杂。为了确定合适的边界条件，分别取外部大圆环的半径差Ro等于L、1.5L和2.0L进行计算，结果见表1。

表1 Ro对光体和全附体仿真结果的影响

从表1可以看到，增大Ro对光体和全附体的所有系数的影响非常微弱。可见Ro超过特征长度1倍时已经满足计算要求。

3.2 内部计算域大小的影响

滑移网格法中，艇体周围被小的圆环包裹。为探讨旋转内部计算域对仿真结果的影响，分别取内部圆环半径Ri为2H(H为艇高)和3H对光体和全附体进行计算，结果见图4。

图4 Ri对光体和全附体计算结果的影响

由图4可见，Ri的增加对光体和全附体的Cy和Cmz几乎没有影响，2H已经满足计算需求。尽管Ri对Cx的计算有一定影响，但是程度较小，不足0.4%。此外，Cx不用于求线性旋转导数，是非重点关注量。增大Ri会大幅增加网格数量，计算效率大减。而当Ri从2H继续减小时，会导致艇体与边界距离太短，布置网格数量较少，容易降低网格质量。因此，综合认为旋臂水池仿真中Ri超过2H即可。

3.3 边界层厚度的影响

根据直航阻力计算时的经验，通常第一层边界层厚度对湍流的计算有一定影响。为此，采用realizablek-ε湍流模式，对第一层网格厚度分别取1.2、2.4、3.6和4.8 mm，相对应的y+分别为40、80、120和160。并都取6层棱柱层网格，增长率默认为1.3，其余设置均相同。计算结果见图5。

图5 y+对光体和全附体计算的影响

由图5可见，不同y+下，光体和全附体的力和力矩系数几乎相同，表明其影响均较小。可见，在旋臂水池的模拟中，采用realizablek-ε湍流模式下y+对力和力矩的影响较小，取值在40～160内均可。

3.4 临界雷诺数

旋臂试验中，通常需要先进行临界雷诺数测量试验，而临界雷诺数通常是通过水动力系数的稳定来判断，不同标准可能会导致其数值有一定的变化。

由上节结论可知，采用realizablek-ε时，y+在合适的范围内，对仿真结果影响较小。因此，在不同雷诺数计算时，仅需调整边界层网格厚度保持相当的y+即可。对不同雷诺数下的工况进行仿真时，均采用realizablek-ε湍流模式，网格数量基本相同，仅对边界层网格调整保持y+基本相同即可。从图6计算结果可以看到，光体和全附体下Cx、Cy和Cmz都随着线速度增加逐渐平稳。Cx和Cy在超过4 m/s后变化幅度较小，初步可以认为4 m/s已经到达临界雷诺数，此时雷诺数约17×106，比文献[5]中的14×106稍大。

图6 雷诺数对光体和全附体计算的影响

对临界雷诺数附近不同速度进行线性拟合得到旋转导数见表2。其中，光体的拟合结果见图7。

图7 光体线性拟合求旋转导数

表2 不同线速度下旋转导数计算结果对比

由表2可见，线速度对旋转导数的计算有一定的影响，对光体的影响大于对全附体，影响在5%以内，对全附体的影响则在2%以内。从数值上，随着线速度的增加，光体旋转导数则逐渐减小，全附体则逐渐增大。在超过4 m/s时，线速度对旋转导数的影响可以忽略。

3.5 湍流模式

采用上文研究得到的计算设置，取文献中的线速度6.5 kn，采用不同的旋转半径对光体和全附体进行水平面的旋转仿真。采用2种湍流模式，以对比其影响，见表3。

表3 CFD仿真结果与试验结果对比

从表3中可以看到，realizablek-ε和SSTk-ω模拟得到N′r结果几乎完全相同，而Y′r则有一定的不同,前者与试验结果符合更好。光体模型的计算结果较全附体差，这是由于光体外形光顺，分离点较不固定的原因。分离点的变动对艇体流场和受力产生较大影响，对分离点的准确求取较为困难。

4 结论

以Suboff为研究对象，通过仿真对比，发现计算内部旋转域纵向与艇体距离大于2倍艇高，外部计算域纵向与艇体距离了超过2倍艇长时，计算域大小对计算结果的影响可以忽略。y+在40～160时，对仿真结果的影响不显著。雷诺数对旋转导数的计算有一定影响，且对光体的影响较全附体大，必须超过临界雷诺数以消除其影响。realizablek-ε和SSTk-ω两种湍流模式下得到的N′r几乎完全相同，而Y′r有一定的不同。总体来看前者与试验符合程度更高。全附体N′r模拟结果与试验结果符合非常好，其余系数误差基本在30%以内。尽管以上因素会对计算结果产生一定的影响，但是，波动较小，显现出较强的可靠性，可以采用该方法用于初期的方案优选。