0 引 言

潜艇作为海军重要的作战武器，对于维护国家海权、协同水面舰船作战、进行战略威慑具有重要意义^[1]。获取潜艇水动力系数是评价其操纵性的重要途径。计算流体力学（CFD）作为重要的途径，可以节省经费、缩短研究周期，且由于获得的水动力系数越来越准确，故越来越受到认可。

文献^[2-5]以 Suboff 模型为研究对象对潜艇定常 运动进行了计算，分析了网格和湍流模型对计算精度的影响，并获得了部分线型和非线性水动力系数。文献^[6]通过MRF方法将回转运动计算转化为定常计算，并进行了数值模拟，分析了回转运动的原理。文献^[7-9]通过自编程序对回转运动进行了研究，分析了湍流模型、网格结构、边界层等对计算精度的影响。

以上文献涉及到的数值计算均是基于固联于艇体的运动坐标系，并将非定常问题通过 UDF 功能转化为定常计算。本文以潜艇 Suboff 模型为研究对象，基于大地坐标系，借助刚性网格方法和 Mesh Motion 功能，通过流体计算软件 Fluent，对潜艇直航阻力运动、水平面和垂直面变漂角运动、悬臂回转运动进行数值计算，并将计算得到的水动力系数与试验值进行对比，验证方法的有效性和准确性。

1 基本理论与控制方程 1.1 ALE 方程

在很多运动物体模拟中，几何体的运动可以分解为整个计算域的刚体运动和局部位置的相对运动的叠加。本文用 ALE 方程来描述这种计算域的刚性运动。其在直角坐标系下的积分形式为：

$\frac{\partial }{\partial t}\iiint\limits_{\Omega }{\vec{Q}\text{d}V}+\iint_{\partial \text{ }\Omega }{\vec{F}(\vec{Q})}\cdot \vec{n}\text{d}S=0,$

(1)

式中的解向量为：

$\vec{Q}={{(\rho ,\rho u,\rho v,\rho w,\rho e)}^{\text{T}}}$

(2)

对流通量矢量项 $\overrightarrow F (\overrightarrow Q )$ 在 x、y 和 z 轴的分量为：

$\begin{align} & {{{\vec{F}}}_{x}}(\vec{Q})=\left\{ \begin{matrix} \rho {{U}_{n}} \\ \rho u{{U}_{n}}+p \\ \rho v{{U}_{n}} \\ \rho w{{U}_{n}} \\ (e+p){{U}_{n}}+{{x}_{t}}p \\ \end{matrix} \right\},{{{\vec{F}}}_{y}}(\vec{Q})=\left\{ \begin{matrix} \rho {{V}_{n}} \\ \rho u{{V}_{n}} \\ \rho v{{V}_{n}}+p \\ \rho w{{V}_{n}} \\ (e+p){{V}_{n}}+{{y}_{t}}p \\ \end{matrix} \right\}, \\ & {{{\vec{F}}}_{z}}(\vec{Q})=\left\{ \begin{matrix} \rho {{W}_{n}} \\ \rho u{{W}_{n}} \\ \rho v{{W}_{n}} \\ \rho w{{W}_{n}}+p \\ (e+p){{W}_{n}}+{{z}_{t}}p \\ \end{matrix} \right\}\text{.} \\ \end{align}$

(3)

式中：Ω 为控制体，∂Ω 为控制体单元的边界；dV 为体积微元；$\overrightarrow n $ 为控制体边界外法向单元向量；dS 为面积微元；ρ，u，v，w，e 分别为流体的密度，x，y，z 轴的流体速度和单位体积流体的总内能。

1.2 刚性网格法

对一些特殊的运动形式可以将物体作为刚性物体处理，此时，计算域网格作为刚性网格与物体固连，所以不用在每个时间步重新生成新的网格，基于 ALE 方程便可描述计算域的运动形式。刚性运动网格计算方法与固定计算域相比只是在通量计算时考虑了网格的运动速度。网格运动速度和网格在每个时刻点的位置由如下方法确定。

设物体的转动轴的方程为：

$\frac{{x - {x_0}}}{A} = \frac{{y - {y_0}}}{B} = \frac{{z - {z_0}}}{C} = t\,$

(4)

设网格 i 中心点的坐标为 $m({x_i},{y_i},{z_i})$，常量 $t = \frac{{A({x_i} - {x_0}) + B({y_i} - {y_0}) + C({z_i} - {z_0})}}{{{A^2} + {B^2} + {C^2}}}$，则 m 在转动轴上的投影点坐标 $n({x_n},{y_n},{z_n})$ 为 $n(At + {x_0},Bt + {y_0},Ct + {z_0})$，则网格中心点关于时间 t 的关系式为：

$x = ({x_i} - {x_n})\cos (a(t)) + ({z_i} - {z_n})\sin (a(t)) + {x_n}\,$

(5)

$y = {y_i}\,$

(6)

$z=-({{x}_{i}}-{{x}_{n}})\sin (a(t))+({{z}_{i}}-{{z}_{n}})\cos (a(t))+{{z}_{n}}.$

(7)

2 水动力系数数值计算实现

水动力系数数值计算通过以下数值试验实现：直航阻力数值试验、垂直面变漂角数值试验、水平面变漂角数值试验。通过 Suboff 模型，利用 Fluent 软件实现各种运动。

2.1 数值计算的准备

通过专业网格划分软件 ICEM CFD，采用全局映射式六面体网格离散方法建立网格。计算区域的大小是：距离潜艇前部有 1 L 距离，尾部有 2 L 距离，上下左右为 1 L 距离。根据 RNG k-ε 湍流模型对于近壁面网格处理方法的要求，Y+ 主要集中在（0，10）,网格数目为 360 万。

2.2 直航阻力数值试验

根据纵向直航阻力试验的特点设定计算域运动速度 u₀ = 0.5 m/s，1 m/s，2 m/s，2.5 m/s，3.045 m/s，5.144 m/s，6.091 m/s，7.160 m/s，8.230 m/s，9.254m/s，Turbulent Intensity 和 Turbulent Viscosity Ratio 均设定为 1%。根据直航阻力运动特点，u = u₀，v = 0，w = 0，p = 0，q = 0，r = 0，则六自由度空间运动方程^[10]可简化为：

$X\text{=}\frac{1}{2}\rho {{L}^{2}}{{{{X}'}}_{uu}}{{u}^{2}},\text{Y=}\frac{1}{2}\rho {{L}^{2}}{{{{Y}'}}_{0}}{{u}^{2}},Z\text{=}\frac{1}{2}\rho {{L}^{2}}{{{{Z}'}}_{0}}{{u}^{2}},M\text{=}\frac{1}{2}\rho .{{L}^{3}}{{{{M}'}}_{0}}{{u}^{2}}.$

(8)

其中 u，v，w，p，q，r 分别为计算域横向速度、纵向速度、垂向速度、横倾角速度、纵倾角速度及偏航角速度。通过 Matlab 对计算结果进行二次样条线拟合如图 2 所示。图 3 是流场的速度剖面图。

2.3 垂直面变漂角数值试验

通过计算不同攻角条件下潜艇的受力和力矩，可推算相关水动力系数。假定潜艇的速度为 U = 3.045 m/s，攻角 α 为 ± 2°，± 4°，± 6°，± 8°，± 10°，± 12°，± 15°，± 20°，± 25°，则计算域速度 $u{\rm{ = }}U\cos \alpha $，$w{\rm{ = }}U\sin \alpha $。计算结果如图 3 所示。

由于 Suboff 模型对称性，故 Y₀′ = 0，K₀′ = 0，N₀′ = 0。而 M_ww′ 数值较小可约等于 0。垂向力 Z 和纵倾力矩 M 关于垂向速度 w 的函数为：

$\left\{ \begin{array}{*{35}{l}} Z(w)\text{=}{{Z}_{0}}\text{+}{{Z}_{w}}w\text{+}{{Z}_{w\text{ }\!\!|\!\!\text{ }w\text{ }\!\!|\!\!\text{ }}}w\text{ }\!\!|\!\!\text{ }w\text{ }\!\!|\!\!\text{ ,} \\ M(w)\text{=}{{M}_{0}}\text{+}{{M}_{w}}w\text{+}{{M}_{w\text{ }\!\!|\!\!\text{ }w\text{ }\!\!|\!\!\text{ }}}w\text{ }\!\!|\!\!\text{ }w\text{ }\!\!|\!\!\text{ }\text{.} \\ \end{array} \right.$

(9)

Z、M 关于不同速度 w 可得如下方程组：

$\left\{ \begin{align} & \sum\limits_{i}{({{Z}_{0}}\text{+}{{Z}_{w}}}w\text{+}{{Z}_{w\text{ }\!\!|\!\!\text{ }w\text{ }\!\!|\!\!\text{ }}}w\text{ }\!\!|\!\!\text{ }w\text{ }\!\!|\!\!\text{ -}{{\text{Z}}_{i}})\text{=0}, \\ & \sum\limits_{i}{({{Z}_{0}}\text{+}{{Z}_{w}}w\text{+}{{Z}_{w\text{ }\!\!|\!\!\text{ }w\text{ }\!\!|\!\!\text{ }}}w\text{ }\!\!|\!\!\text{ }w\text{ }\!\!|\!\!\text{ -}{{\text{Z}}_{i}}){{\text{w}}_{i}}\text{=0}}, \\ & \sum\limits_{i}{({{Z}_{0}}\text{+}{{Z}_{w}}w\text{+}{{Z}_{w\text{ }\!\!|\!\!\text{ }w\text{ }\!\!|\!\!\text{ }}}w\text{ }\!\!|\!\!\text{ }w\text{ }\!\!|\!\!\text{ -}{{\text{Z}}_{i}}){{\text{w}}_{i}}\text{ }\!\!|\!\!\text{ }{{\text{w}}_{i}}\text{ }\!\!|\!\!\text{ =0}}, \\ & \sum\limits_{i}{({{\text{M}}_{0}}\text{+}{{\text{M}}_{w}}}w\text{+}{{\text{M}}_{w\text{ }\!\!|\!\!\text{ }w\text{ }\!\!|\!\!\text{ }}}w\text{ }\!\!|\!\!\text{ }w\text{ }\!\!|\!\!\text{ -}{{\text{M}}_{i}})\text{=0}, \\ & \sum\limits_{i}{({{\text{M}}_{0}}\text{+}{{\text{M}}_{w}}}w\text{+}{{\text{M}}_{w\text{ }\!\!|\!\!\text{ }w\text{ }\!\!|\!\!\text{ }}}w\text{ }\!\!|\!\!\text{ }w\text{ }\!\!|\!\!\text{ -}{{\text{M}}_{i}}){{\text{w}}_{i}}\text{=0}, \\ & \sum\limits_{i}{({{\text{M}}_{0}}\text{+}{{\text{M}}_{w}}}w\text{+}{{\text{M}}_{w\text{ }\!\!|\!\!\text{ }w\text{ }\!\!|\!\!\text{ }}}w\text{ }\!\!|\!\!\text{ }w\text{ }\!\!|\!\!\text{ -}{{\text{M}}_{i}}){{\text{w}}_{i}}\text{ }\!\!|\!\!\text{ }{{\text{w}}_{i}}\text{ }\!\!|\!\!\text{ =0}\text{.} \\ \end{align} \right.$

(10)

解方程组并进行无因次化，可得水动力系数。

2.4 水平面变漂角数值试验

假定潜艇的速度为 U = 3.045 m/s，漂角 β 为 ± 2°，± 4°，± 6°，± 8°，± 10°，± 12°，± 14°，± 16°，由垂直面变漂角 β 试验特点可知，u ≠ 0，v ≠ 0，w = 0，p = 0，q = 0，r = 0。计算域运动速度 ${{u = }}U \cdot \cos \beta $，$v{{ = }}U \cdot \sin \beta $。

由于 Suboff 模型对称性，故 Y₀′ = 0，N₀′ = 0。横向力 Y 和偏航力矩 N 关于横向速度 v 的函数为：

$\left\{ \begin{array}{*{35}{l}} Y(v)={{Y}_{v}}v+{{Y}_{v\text{ }\!\!|\!\!\text{ }v\text{ }\!\!|\!\!\text{ }}}v\text{ }\!\!|\!\!\text{ }v\text{ }\!\!|\!\!\text{ }, \\ N(v)={{N}_{v}}v+{{N}_{v\text{ }\!\!|\!\!\text{ }v\text{ }\!\!|\!\!\text{ }}}v\text{ }\!\!|\!\!\text{ }v\text{ }\!\!|\!\!\text{ }\text{.} \\ \end{array} \right.$

(11)

Y，N 关于不同速度 v 可得如下方程组：

$\left\{ \begin{align} & \sum\limits_{i}{({{Y}_{v}}v}+{{Y}_{v\text{ }\!\!|\!\!\text{ }v\text{ }\!\!|\!\!\text{ }}}v\text{ }\!\!|\!\!\text{ }v\text{ }\!\!|\!\!\text{ -}{{Y}_{i}})\text{=0}, \\ & \sum\limits_{i}{({{Y}_{v}}v\text{+}{{Y}_{v\text{ }\!\!|\!\!\text{ }v\text{ }\!\!|\!\!\text{ }}}v\text{ }\!\!|\!\!\text{ }v\text{ }\!\!|\!\!\text{ -}{{Y}_{i}}){{\text{v}}_{i}}\text{=0}}, \\ & \sum\limits_{i}{({{Y}_{v}}v\text{+}{{Y}_{v\text{ }\!\!|\!\!\text{ }v\text{ }\!\!|\!\!\text{ }}}v\text{ }\!\!|\!\!\text{ }v\text{ }\!\!|\!\!\text{ -}{{Y}_{i}}){{\text{v}}_{i}}\text{ }\!\!|\!\!\text{ }{{\text{v}}_{i}}\text{ }\!\!|\!\!\text{ =0}}, \\ & \sum\limits_{i}{({{N}_{v}}v}\text{+}{{N}_{v|v|}}v|v|-{{N}_{i}})\text{=0}, \\ & \sum\limits_{i}{({{N}_{v}}}v+{{N}_{v|v|}}v|v|-{{N}_{i}}){{v}_{i}}=0, \\ & \sum\limits_{i}{({{N}_{v}}}v+{{N}_{v|v|}}v|v|-{{N}_{i}}){{v}_{i}}|{{v}_{i}}\text{ }\!\!|\!\!\text{ =0}\text{.} \\ \end{align} \right.$

(12)

解方程组并进行无因次化，可得水动力系数。

2.5 计算结果校验

通过直航阻力数值试验及水平面变漂角数值试验获得的水动力导数计算结果误差为 3.7%，满足计算工程需求。通过垂直面变漂角数值试验获得的水动力导数计算结果误差为 7%，这是由于当攻角较大时，通过湍流模型计算不能准确捕捉到流动分离；而且全附体潜艇大攻角运动时，潜艇指挥台围壳导致背流区分离，而且有主附体干扰的存在，同时指挥台围壳对流场的干扰也会对翼产生影响，影响舵效及伴流场。

3 结 语

本文根据每种运动的不同，对潜艇六自由度空间运动方程进行化简，获得了直航阻力运动、垂直面变漂角运动、水平面变漂角运动的空间运动方程。利用 Suboff 模型，基于刚性网格法，通过设定计算域的运动形式，实现了潜艇直航阻力运动、垂直面变漂角运动、水平面变漂角运动的数值模拟，通过对计算结果进行拟合得到了线性水动力系数和非线性水动力系数，通过与相关试验结果进行对比可知，计算误差在 7% 左右。