计算流体力学(computation fluid dynamic, CFD)是水下航行器最重要的总体性能之一, 研究航行器的操纵性对于保证其航行安全、充分发挥水下航行器工作效能, 都有着极其重要的意义^[1]。准确高效的获得水动力系数是开展水下航行器操纵性研究与设计方案改进的重要前提。常用的水动力系数求取方法包括经验公式估算方法、拘束模型试验方法、自航模型试验方法。经验估算方法具有很强的经济性和实用性，但其适用性和精确度难以保证。模型试验方法是当前计算水动力系数较为精确的一种方法，但是由于其试验周期长，耗费成本高，在水下航行器设计初期并不适用^[2]。随着计算机技术和CFD技术的飞速发展，人们逐渐用数值模拟代替真实的模型试验，分析的成本大幅降低，周期大大缩短。大量学者对数值模拟求解水动力系数做了相关研究。Racine^[3]基于重叠网格方法模拟了扁平体水下航行器NNEMO的斜航试验、回转臂试验和平面运动机构试验，求解了11个主要水动力系数。Kim^[4]、Ahmad^[5]及Hiroyoshi^[6]等利用CFD的重叠网格技术模拟平面运动机构试验，通过与真实水池中的试验结果对比证明了数值模拟方法的可靠性。林兆伟等^[7]提出了一种在三维操纵性数值水池中模拟测试水下航行器水动力系数的计算方法，以解决数值计算过程中的精度和效率的问题。文献[8-11]利用CFD软件结合湍流模型、动网格等技术对水下航行器的拘束平面运动进行模拟，进而分析其水动力性能。文献[12-14]为提高旋转导数的预报精度，采用添加动量源项的方法将水下航行器旋臂试验的数值模拟转化为定常问题进行求解，计算结果与试验数据吻合良好。这些研究成果充分说明数值模拟进行水下航行器水动力研究的可靠性和广阔应用前景。

基于平面机构运动模拟、拘束模型试验等方法对求取水下航行器水动力系数的研究已经相当成熟，且应用极为广泛。但是，上述方法的不足之处在于：首先，计算周期较长，因为需要进行多种运动频率或航行速度下的仿真分析，以便求取系数的平均值进而提高结果的准确性；其次，求解过程中每个系数自身都带有一定的误差，因此，分别获取导致动力学模型整体的可靠性不可评估，也难以对建立的动力学模型进行改进或简化。

本文以水下航行器为例，采用3种方法计算其水动力系数，并提出了一种高效的空间拘束运动模拟方法。

1 数值计算方法

为了确定航行器的位置和姿态，研究其运动状态及艇体的水动力特性，本文引入固定坐标系O_E-ξηξ以及运动坐标系O_b-xyz，如图 1所示。固定坐标系的原点可以任意选取计算域中某一点；运动坐标系的原点位于航行器的重心处。基于载体坐标系描述航行器运动时常用参量的符号如表 1所示。

本文水下航行器模型由主艇体、围壳、围壳舵及十字型艉舵4部分组成。围壳舵、十字型艉舵翼断面采用NACA0012翼型。该航行器模型由哈尔滨工程大学设计，参考了德国的U212 A，日本的苍龙号等多个典型的潜艇布局。模型艇长为3.08 m，艇高为0.37 m，直径为0.3 m，重量为220 kg。

计算域设置为一长方体区域，艏部距入口边界2倍艇长；艉部距出口边界5倍艇长；四周距边界2倍艇长。

选取自适应切割体网格技术对计算域进行网格离散，该网格技术能够确保计算域内绝大部分网格为计算性能优良的六面体网格，同时通过边界层网格设置可对艇体近壁面网格进行控制。计算中合理的布置第1层网格的高度，保证20~100，并应用整体动网格技术来模拟水下航行器的六自由度空间运动。在网格生成过程中考虑到计算潜艇结构的复杂性，在对网格相关参数设置时采取在近艇体区域设定较密的网格然后以一定膨胀系数外推。同时适当的加密各附体处的网格以便较准确的捕捉到该处的流动特性。为了避免网格粗糙时引起的计算结果的误差过大，文中在考虑计算资源的基础上，尽可能的对网格进行细化，总的网格数为590万。生成的网格模型如图 2所示。

2 平面运动机构试验模拟

平面运动机构试验是所有拘束模型试验方法中能够获得水动力系数最多的一种方法，通过纯横荡运动、纯摇艏运动、纯升沉运动和纯俯仰运动几个工况的重复模拟，可以获得线/角速度，以及线/角加速度相关的水动力系数。

2.1 纯横荡运动

本文模拟循环水槽中的平面运动机构试验，纯横荡运动中，来流速度沿x方向为U，仅在水平面内的y方向给艇体施加一个按正弦规律变化的位移。其运动参数如下：

$ \left\{ \begin{array}{l} \eta = a\sin \left( {\omega t} \right)\\ v = \dot \eta = a\omega \cos \left( {\omega t} \right) \end{array} \right. $

(1)

式中:η为水下航行器在水平面y方向的位移，a为振幅，ω为振动的圆频率。艇体的振荡频率依次取0.1、0.15、0.2 Hz, 轴向来流速度恒定1.2 m/s，振幅取固定值a为0.04 m。文中计算的小振幅的平面运动，在满足线性假设的前提下水平面内的运动方程为：

$ \left\{ \begin{array}{l} Y = {Y_{\dot v}}\dot v + {Y_v}v + {Y_{\dot r}}\dot r + {Y_r}r\\ N = {N_{\dot v}}\dot v + {N_v}v + {N_{\dot r}}\dot r + {N_r}r \end{array} \right. $

(2)

将横荡运动参数设置带入上述运动方程中并进行无因次化，得到水动力表达式如下：

$ \left\{ \begin{array}{l} Y' = - \frac{{a{\omega ^2}L}}{{{\mathit{\boldsymbol{U}}^2}}}{{Y'}_{{\rm{\dot v}}}}\sin \left( {\omega t} \right) + \frac{{a\omega }}{\mathit{\boldsymbol{U}}}{{Y'}_{\rm{v}}}\cos \left( {\omega t} \right)\\ N' = - \frac{{a{\omega ^2}L}}{{{\mathit{\boldsymbol{U}}^2}}}{{N'}_{{\rm{\dot v}}}}\sin \left( {\omega t} \right) + \frac{{a\omega }}{\mathit{\boldsymbol{U}}}{{N'}_{\rm{v}}}\cos \left( {\omega t} \right) \end{array} \right. $

(3)

分别计算每种工况下艇体所受的侧向力和转首力矩的变化情况。利用最小二乘方法即可求得不同工况下的水动力系数Y′_v、Y′_v、N′_v、N′_v，通过不同工况的重复模拟求解水动力系数的平均值来消除频率对试验结果的影响。

2.2 纯摇艏运动

数值模拟水下航行器纯摇艏运动是在艇模上施加一个艏向角的正弦变化，运动时要求模型迎流方向与中心处轨迹曲线相切。相应的运动参数为：

$ \left\{ \begin{array}{l} \psi = {\psi _0}\sin \left( {\omega t} \right)\\ r = \dot \psi = {\psi _0}\omega \cos \left( {\omega t} \right) \end{array} \right. $

(4)

式中ψ₀为艇体的转艏幅值。

将参数设置代入运动线性方程中得到艇体做纯摇艏运动时的无因次化水动力为

$ \left\{ \begin{array}{l} Y' = - \frac{{{\psi _0}{\omega ^2}{L^2}}}{{{U^2}}}{{Y'}_{{\rm{\dot r}}}}\sin \left( {\omega t} \right) + \frac{{{\psi _0}\omega L}}{U}{{Y'}_{\rm{r}}}\cos \left( {\omega t} \right)\\ N' = - \frac{{{\psi _0}{\omega ^2}{L^2}}}{{{U^2}}}{{N'}_{{\rm{\dot r}}}}\sin \left( {\omega t} \right) + \frac{{{\psi _0}\omega L}}{U}{{N'}_{\rm{r}}}\cos \left( {\omega t} \right) \end{array} \right. $

(5)

与横荡运动相似，通过运动方程的拟合，并将各工况下的计算结果求取平均值得到水动力导数Y′_r、Y′_r、N′_r、N′_r。垂直面内的拘束运动求解方法与水平面的相似，此处不再详述。

3 空间拘束运动模拟

通过模拟平面运动机构试验可以求得动力学模型中的部分水动力系数，但是需要进行多次数值计算。此外，为了保证所得系数的可靠性，计算中需要拟合不同频率的运动工况，通过求取其平均值作为最终的结果。这增加了计算周期，并且限制了水下航行器设计的进度。为解决这一问题，本文提出一种高效的水动力系数计算方法-空间拘束运动模拟方法。利用该方法仅通过一次计算即可求得动力学模型中的全部水动力系数。

3.1 运动方程

水下航行器六自由度空间运动方程以格特勒用于潜艇模拟研究的标准运动方程为基础^[15]，该方程含108个水动力系数。在现实中由于试验的困难很难全部得到^[14]。此外，水下航行器做定速直航等弱机动运动时，由于整体受到的水的作用力变化较小，动力学模型中的非线性项影响很小。因此，本文仅通过数值模拟求取与艇体速度和加速度直接相关的线性水动力系数，在不考虑舵翼影响的前提下，格特勒的标准运动方程可简化为：

$ X = {X_{{\rm{\dot u}}}}\dot u + {X_{{\rm{uu}}}}{u^2} $

(6)

$ Y = {Y_{{\rm{\dot r}}}}\dot r + {Y_{{\rm{\dot v}}}}\dot v + {Y_{\rm{r}}}ur + {Y_{\rm{v}}}uv $

(7)

$ Z = {Z_{{\rm{\dot q}}}}\dot q + {Z_{{\rm{\dot w}}}}\dot w + {Z_{\rm{q}}}uq + {Z_{\rm{w}}}uw + {Z_0} $

(8)

$ K = {K_{{\rm{\dot p}}}}\dot p + {K_{\rm{p}}}up + {K_{{\rm{\dot v}}}}\dot v + {K_{\rm{v}}}\mathit{uv} $

(9)

$ M = {M_{{\rm{\dot q}}}}\dot q + {M_{{\rm{\dot w}}}}\dot w + {M_{\rm{q}}}uq + {M_{\rm{w}}}uw + {M_0} $

(10)

$ N = {N_{{\rm{\dot r}}}}\dot r + {N_{{\rm{\dot v}}}}\dot v + {N_{\rm{r}}}\mathit{ur + }{\mathit{N}_{\rm{v}}}\mathit{uv} $

(11)

尽管线性方程略去了几十个不易计算和测量的耦合项水动力系数，但其仍然有足够的拟合能力对水下航行器小扰动状态下的直航运动进行动力学建模，后续的空间拘束运动模拟结果将进一步证明。

3.2 模拟方法

与平面拘束运动数值模拟方法不同，空间拘束运动模拟是一种在现实环境中较难实现，但在数值环境中很容易完成的模拟方法。在潜艇的6个自由度上人为定义互不相同的6个不规则扰动。使艇体在各个自由度上的速度和加速度都不为0，而且呈现不规则的变化。其定义形式可表示为：

$ \varphi = {\varphi _0} + {\varphi _{\rm{p}}}\sin \left( {2{\rm{ \mathsf{ π} }}{f_{{\varphi _1}}}t} \right)\sin \left( {2{\rm{ \mathsf{ π} }}{f_{{\varphi _2}}}t} \right)\sin \left( {2{\rm{ \mathsf{ π} }}{f_{{\varphi _3}}}t} \right) $

(12)

式中：φ表示艇体的任一运动参数u、v、w、p、q、r；φ₀是一个持续的速度，对于潜艇的x方向来说，φ₀代表其巡航速度，对于其他方向的线速度和角速度φ₀=0；φ_p表示运动中所加扰动的幅值，计算中合理的设置该参数使得运动限制于小幅度的振动。3个正弦函数的乘积是施加于艇体上的一个不规则的扰动，用于调整艇体的运动状态, 使其更接近于潜艇在深海中的运动。合理的调整式中正弦函数的圆频率使得各方向的运动尽可能不相关。为了便于与试验结果的比较，除速度u和角速度p外的运动参数，其设定参考拘束模型试验的振幅和频率。

3.3 相关性分析

在空间拘束运动模拟的基础上采用最小二乘法对水下航行器的水动力系数进行辨识。为保证样本数据的不相关性，在参数设定后选取部分样本点对各速度之间的相关性进行分析。皮尔逊相关系数为：

$ r = \frac{{\sum\limits_{i = 1}^n {\left( {{X_i} - \bar X} \right)\left( {{Y_i} - \bar Y} \right)} }}{{\sqrt {\sum\limits_{i = 1}^n {{{\left( {{X_i} - \bar X} \right)}^2}} } \sqrt {\sum\limits_{i = 1}^n {{{\left( {{Y_i} - \bar Y} \right)}^2}} } }} $

(13)

式中：X为X的均值；Y为Y的均值。

依次计算出各方向速度之间的相关性系数，根据该系数确定各方向运动是否相关进而判断设置的圆频率参数的可行性。该系数的取值[-1, 1]，当r＜0时，对应的2个变量是负相关的；若r＞0, 则二者是正相关的；当r=0时表示各方向速度不存在相关关系。文中设定的圆频率参数得到的速度之间相关系数几乎为0，可见水下航行器各方向的运动相关性很小。

利用CFD仿真分析软件STAR-CCM+模拟水下航行器的空间拘束运动，在计算收敛后监测得到每一时刻艇体各方向的受力和力矩值X、Y、Z、K、M、N，运动时间持续150 s。图 3所示是60~80 s内艇体的线速度以及角速度的变化情况。

3.4 拟合优度检验

回归模型的准确度是判断设定的圆频率参数的否合理性。本文以15~120 s的计算数据为基准，设定了置信区间为95%。运用最小二乘法求解线性运动方程的回归模型，并基于剩余20%的数据对模型的准确性进行交叉验证，确保模型的拟合度良好。引入度量拟合优度的统计量确定系数R^2[16-17]，作为衡量回归模型拟合程度的指标。R²值越接近于1，说明建立的回归模型精度越好。文中艇体各方向的受力及力矩拟合得到的回归模型其系数R²均大于0.95，说明该模型误差较小，能够满足后续计算的需求。

4 计算结果与分析 4.1 水动力试验

为了验证空间拘束运动模拟计算水动力系数方法的可靠性，在循环水槽对潜艇模型开展平面运动机构试验，将试验结果与数值模拟结果进行对比验证。平面运动机构试验包括纯横荡、纯摇艏、纯升沉以及纯俯仰试验，试验中艇体模型的安装情况如图 4所示。

各个试验中流速均为1.2 m/s，振幅为0.04 m；振荡频率区间为0.1 ~ 0.2 Hz。平面运动机构试验和空间拘束运动模拟计算得到的水动力系数对比结果如表 2所示。

4.2 2种数值方法的比较

基于平面拘束运动模拟方法可以求得速度和角速度相关的线性水动力系数，但是为了降低误差需要通过多次计算模拟不同工况下的运动性能，求取系数的平均值，如表 3所示。若使用如下的计算工作站配置：双CPU(Intel Xeon E5-2640 @ 2.50 GHz)且内存储器为48 G。采用该方法完成水动力系数的计算需要进行不少于12个算例的模拟，计算周期超过10 d。而相同的条件下基于空间拘束运动模拟方法仅通过一个算例、一次计算就可求得所需的全套水动力系数，计算时间约为35 h。

比较平面运动机构试验模拟与文中提出的新方法的计算结果，可以看出对于大部分水动力系数，基于2种CFD方法求得的值较为接近。因此，在设计初期，可以用空间拘束运动模拟方法代替现有的常规数值模拟方法，快速求得设计艇型的水动力系数。分析计算值与试验结果的误差，可以看出大部分数值估计的水动力系数值具有较好的精度，能够满足工程实际的需要，进一步说明了空间拘束运动模拟方法的可行性。

解耦方法更注重单个系数的精度，基于该方法计算整合出来的动力学模型整体精度无法评估。本方法更注重动力学模型整体精度，高效计算后可以评估动力学模型的整体精度并进行深入分析。虽然表 2中基于本文方法求取的部分水动力系数误差较平面运动机构试验模拟结果的误差稍大，但是依然能够保证动力学模型的整体精度。

图 5为某时间段内艇体的受力及力矩变化曲线图。基于本文方法求得的水动力系数线性拟合得到的曲线与基于试验的EFD结果进行比较，可以看出2条曲线很接近，充分说明了所建立动力学模型的可靠性。上述回归模型的建立是基于航行器简化后的线性空间运动方程。图中的另一条曲线为CFD数值模拟结果，为数值模拟得到的真实运动受力分析。CFD的线性拟合结果与其相比可以看出线性方程能准确的反应艇体的水动力变化趋势，且其计算数值与艇体的实际模拟值相当接近，说明水下航行器在做小扰动状态下的直航运动时，非线性项的影响较小，可以忽略不计。

4.3 计算误差分析

计算结果中系数Y′_r和N′_v与试验值虽然在同一个量级，但是误差偏大。两种CFD的拟合结果较为接近。造成该误差的原因可能有以下几个方面：1)模型试验中的尺度效应、试验仪器以及PMM运动机构中两支杆运动引起的误差^[18]；2)数值计算中网格、湍流模型等不确定因素引起的误差；3)模型的主艇体横截面尺寸基本满足不大于水槽截面积的3%，但是该模型带有较大的指挥台围壳，当艇体艏部下倾槽底时带来了较大的“堵水”效应，因此引起了试验值测定结果的误差^[19]；4)依据文献[20]中对水动力导数的灵敏度与其辨识精度的关系的分析，对灵敏度较大的水动力导数往往能够得到比较接近真值的估计值，而对模型贡献较小的参数一般很难得到。因此，文中进一步分析了Y方向受力和绕Z轴转矩方程中，几个水动力系数的比重，如图 6所示。图中CFD和EFD曲线分别是依据数值模拟结果和水动力系数实验数据拟合得到的。可见在上述拘束运动中，与其余参数相比，系数Y′_r和N′_v所对应的项在方程中所占的比重微小，图中几乎为一条趋近于0的线。因此，要想准确求得系数Y′_r和N′_v的值，需要分析该系数作用显著的其他运动。

5 结论

1) 该方法计算所得的大部分水动力系数误差小于15%，能够满足工程应用的需求^[21]。与现有数值方法相比，空间拘束运动模拟方法在保证计算精度的同时，极大的缩减了水动力系数计算的周期。

2) 水池试验中由于模型带有较大的指挥台围壳，当艇体艏部下倾槽底时带来了较大的“堵水”效应，故使试验值产生一定的误差。在今后此类带附体试验中，可以考虑尽量增大缩尺比，以减小模型的尺寸。

3) 基于该方法可依据水动力系数的敏感性进一步对水下航行器的动力学模型进行简化。

4) 空间拘束运动方法不局限于小扰动状态时线性水动力系数的求解，同样可应用于各种复杂工况下水下航行器动力学模型的建立。

本研究为水下航行器操纵运动数学模型的建立提供了一种新的途径，为其操纵性预报和仿真系统的设计提供了便利。