复杂外形水下潜伏平台流体动力数值模拟*

王 鹰，袁绪龙，赵文迪，梁丹丹

（1.西北工业大学航海学院，西安 710072；2.青岛西海岸新区海洋发展局，山东 青岛 266000）

0 引言

水下潜伏平台由锚块及系留在水下具有特定功能的平台组成，在水下自主发射及海洋环境监测、水下侦察、通信和作战方面发挥着重要作用。无海流作用时，在正浮力作用下平台呈现垂直状态；有海流存在时，平台在海流绕流的水动力作用下改变初始位置和姿态，产生偏降和倾斜。由于圆柱绕流的卡门涡街效应，会使平台和锚索组成的系统产生周期性的摇摆和晃动。当平台为水雷或者导弹发射筒时，初始位置的改变和姿态变化会影响其初始弹道，降低命中率，当平台用于环境监测等方面会影响数据的准确性。

为了求解水下潜伏平台的运动，必须研究平台自身的流体动力特性以及在海流作用下的卡门涡街效应形成的周期性的流体动力，为系统动力学建模与仿真提供基础数据。

近年来，水下航行器流体动力CFD 计算的研究发展为本文的计算提供了基础，如马烨利用FLUENT 无粘流模型以及动网格技术，数值模拟物体的匀变速运动，并得到其所受的合力，然后根据其受力方程得到物体六自由度运动情况下的附加质量。尤天庆以计算流体动力学方法为基础，结合动网格和动坐标系技术，进行基于雷诺平均N-S 方程的阻尼力计算方法研究，发展针对水下航行体附加力系数和俯仰阻尼力矩系数的数值计算方法。莫慧黠以通用CFD 软件FLUENT 为平台，运用其中的动态网格技术，以三维水下物体为对象，进行了基于3D N-S 方程的全粘流的附加质量计算方法研究。共计算了3 种外形，形成具体的数值计算流程、探讨计算结果的影响因素，并研究外形对流态及附加质量的影响。

由于目前通用CFD 软件的强大功能，包括前处理网格生成软件，可以预见水下航行器流体动力CFD 计算方法将具有更大的从理论到工程实践的可扩展性，比试验方法更易实现，在水下物体航行性能的分析预报等方面将具有广阔的应用前景。

本文针对四联装水下潜伏平台，建立了运动坐标系，介绍了位置力、阻尼力和附加质量计算方法，计算并分析了这种水下潜伏平台的流体动力特性。

1 位置力计算方法与结果

1.1 外形与坐标系定义

本文研究对象为四联装平台，由于实际外形结构尺寸较大，所以采用1∶10 缩比模型进行研究，外形结构如图1 所示。

图1 1∶10 缩比平台外形结构图

平台体坐标系的原点设在平台浮心处，y 轴位于系统对称中心线上，x 轴位于平行与矩形板边的对称面内，垂直于y 轴，z 轴按照右手定则确定。平台运动和海流相对运动形成的攻角定义为相对速度矢量在xoy 平面内的投影与x 轴的夹角，偏下为正。

1.2 计算域选取与网格划分

采用全模三维网格进行计算，计算域为长方体，尺寸为3 500 mm×2 000 mm×2 000 mm，如图2所示，设置来流方向、计算域上下和计算域前后为速度入口，距离平台较远的面为压力出口。

图2 计算域示意图

本文采用NUMECA 公司的HEXPRESS 软件划分网格，利用BOX 标记网格加密区和为特定面指定网格尺寸这两种方法，划分出满足计算需求的网格。整体而言，全局采用非结构的六边形正交各向异性网格填充计算域。对于粘性问题，通过拆分第1层网格的方法，在物面层附近生成各向异性的高质量结构式分布的边界层网格。网格划分时需要注意以下几点：

1）细化平台表面网格

当平台沉浸在海流中，处于全沾湿状态，沾湿部分产生的流体动力为受力主要来源，因此，有必要在平台壳体上布置较密网格，提高流体动力计算精度。因此，需要利用HEXPRESS 的BOX 技术，对整个平台进行网格加密，由于靠近压力出口的部分，速度场和压力场变化较大，所以加大靠近压力出口区域BOX 范围。

2）细化平台连接板、盖子区域的网格

由前文可知，平台是由上、中、下3 块连接板将4 个圆柱筒并且上下各有4 个盖子的复杂四联装平台，与普通航行器不同之处在于不是回转体且3 块连接板的受力较其他位置大得多，所以将平台每部分进行单独网格划分，对连接板网格进行细化。

结合以上考虑，计算域的大部分区域，网格单元都接近长方体，且通过网格的长细比控制，最终生成的网格长细比不大于2。网格数量为305 万左右，网格质量分析显示无负体积、凹面和扭曲网格，且正交性高于80°的网格单元占总数的95%以上。模型壁面网格如图3 所示，XOZ 面上的网格如图4所示。

图3 模型表面网格

图4 XOZ 面上网格

数值求解模型及各参数见表1，各边界条件的参数设置见表2。

表1 模型参数设置

表2 边界条件设置

1.3 计算结果与分析

本文取海流1 kn 的速度进行分析，参考面积为平台横截面面积。

1.3.1 平台阻力特性

平台阻力系数随攻角变化的线性拟合曲线如图5 所示。

图5 阻力系数

1.3.2 平台升力特性

平台升力系数随攻角变化的线性拟合曲线如下页图6 所示。

图6 升力系数

1.3.3 平台俯仰力矩特性

平台俯仰力矩系数随攻角变化的线性拟合曲线如图7 所示。

图7 俯仰力矩系数

可见，在-10°～10°攻角范围内俯仰力矩系数随攻角线性变化，线性拟合得到mz=-0.070 7-0.013 5α。由于平台外形上下不对称，所以在零攻角时，存在一定的升力和俯仰力矩。

2 阻尼力计算方法与结果

2.1 计算原理概述

旋转导数可以采用悬臂水池试验模型获得，也可以利用CFD 软件模拟定常回转运动，获得流体动力系数随回转角速度变化的关系曲线，再线性拟合得到阻尼力系数对无量纲旋转角速度的旋转导数，本文分别计算了潜伏平台绕y 轴和绕z 轴旋转的阻尼力系数旋转导数。

2.2 计算域选取与网格划分

2.2.1 绕y 轴旋转

圆柱形计算域如图8 所示，利用HEXPRESS 进行网格划分，利用自定义函数UDF 指定模型的旋转角速度，采用动计算域进行非定常计算，计算得出不同旋转角速度下的力矩来线性拟合得出旋转导数，计算工况如表3 所示。

图8 圆柱形计算域

表3 计算工况

利用FLUENT 软件提供的旋转参考框架，指定旋转参考框架的旋转角速度，恒定线速度，通过同时改变旋转角速度和旋转半径，来计算平台不同角速度下的流体动力系数，线性拟合得到旋转导数。

计算域如图9 所示，计算工况如表4 所示。

图9 绕z 轴旋转计算域

表4 计算工况表

绕z 轴旋转无量纲角速度与阻力系数之间关系如下页图10 所示。

图10 无量纲角速度与阻力系数关系

绕z 轴旋转无量纲角速度与俯仰力矩系数之间关系如图11 所示。

图11 无量纲角速度与俯仰力矩系数关系

绕y 轴旋转无量纲角速度与偏航力矩系数之间关系如图12 所示。

图12 无量纲角速度与偏航力矩系数关系

3 附加质量计算方法与结果

3.1 计算原理概述

本文利用FLUENT，运用动计算域方法进行非定常计算，利用自定义函数UDF 指定初速度和加速度，监控并保存流体对平台的作用力和力矩；通过受力分析，提取流体动力参数，根据平台外形特点；本文需要计算沿x 轴、y 轴的加速平动和绕x 轴、y轴的加速转动。

3.2 计算域选取与网格划分

沿y 轴加速平动的计算域如图13，沿x 轴加速平动的计算域大小与图13 相同，只是右侧为压力出口，其他5 个面都为压力入口，沿x 轴加速转动的计算域如图14，沿y 轴加速转动的计算域如图15，计算工况如表5 所示。

表5 工况表

图13 沿y 轴加速平动计算域

图14 沿x 轴加速转动计算域

图15 沿y 轴加速转动计算域

3.3 计算结果与分析

4 结论

本文以复杂外形平台为研究对象，运用CFD 计算方法，对平台位置力、阻尼力和附加质量进行了仿真计算。计算了-10°∶10°攻角下的位置力，获得了阻力系数、升力系数和俯仰力矩系数随攻角变化的关系式。根据平台外形特点，计算得出平台主要旋转导数和附加质量。

表6 附加质量计算