0 引　言

21世纪是海洋的世纪，随着科技进步，为了提升对海洋资源的探索和水下环境的监测能力，进一步认识海洋、保护海洋，水下无人航行器(AUV)逐渐成为科技发达国家关注的热点。根据AUV的排水体积，可将其分为便携式、轻型、重型和矩形4种，其中便携式AUV具有良好的机动性能，广泛应用于军事和民用领域^[1]。AUV在水下作业过程中，其六自由度运动会遭遇各种水下环境扰动，为了能稳定完成水下作业并且自主回收，AUV需要更高的性能保障以及良好的操纵控制能力。在AUV的设计过程中，水动力性能是其设计所要参考的重要依据。要想对水动力性能和实现AUV的操纵控制，就需要获取AUV的水动力系数，水动力系数可以通过AUV不同的运动获取。

目前，获取水动力系数主要近似估算法、船模试验法、理论计算法、数值模拟法和系统参数辨识法5种，其中数值模拟计算水动力系数的方法已经广泛应用于实际工程中。王庆云等^[2]针对Suboff模型，采用流体力学软件Fluent对不同尾操纵面进行不同航速下的阻力分析，并且完成平面运动模拟，最终通过傅离叶级数拟合获取水动力系数。程健^[3]基于Fluent，定义UDF动网格，采用标准k-ω湍流模型，通过建立固定坐标系和运动坐标系，分析了航行器平面运动方程；采用数值模拟获取水下机器人水动力系数，并近似计算部分耦合水动力系数，完成水下机器人的稳定性分析和滑模控制。王雪梅^[4]利用数值模拟计算法，通过Fluent进行数值模拟和水动力系数的近似计算，完成水动力系数的计算，并根据对水下航行器受力分析完成水下航行器的操纵控制。孙铭泽^[5]通过雷诺方程以及离散网格，实现了水下航行器的平面运动虚拟仿真，其水动力系数和实际试验值误差较小，达到工程应用的要求。董苗苗^[6]通过数值模拟获得便携式水下无人航行器水动力系数，完成了AUV的定常运动和非定常运动数值模拟，验证了STAR CCM+数值模拟的有效性。宋男^[7]通过建立水下机器人模型，数值模拟了匀速直线运动和匀加速直线运动，计算出相应的力和力矩获取水动力系数，并对水动力性能分析，验证了数值模拟计算的可靠性。

本文基于数值模拟计算的方法，使用新一代流体力学仿真软件STAR-CCM+^[8]，通过模拟AUV定常运动（匀速直航运动和斜航运动）和非定常运动（纯横荡、纯首摇、纯升沉和纯俯仰）、测出AUV的力和力矩，通过Matlab数据拟合，完成便携式AUV水动力系数的测定，通过对比数值模拟数据误差，分析误差来源并完成初步验证。

1 数值计算模型与设置 1.1 AUV模型

使用以REMUS 100为母艇设计的模型，运用STAR-CCM+对其PMM进行数值模拟计算，航行器裸壳体和声呐换能器的纵切面外形尺寸如图1所示。

十字舵选用流线型NACA0012舵面，其最前端到首部最前端的直线距离为1.189 m。

1.2 数值计算设置

相比于Fluent等传统软件，STAR-CCM+自带平面运动机构（PMM）可以直接模拟AUV的不同工况。本文运用流体力学软件STAR-CCM+和平面运动机构（PMM）对AUV平面运动进行数值模拟，PMM只能对单个平面进行数值计算，一般用于模拟水平面运动，将AUV绕着x轴方向旋转90°，使用PMM完成垂直面运动的数值模拟^[9]。

在计算过程中，湍流模型采用标准 $k - \varepsilon$ 两个方程的模型。为了提高收敛性，将尾翼后端的锐角进行微处理成弧形，防止出现网格单元偏斜角^[10]。选择长方体形状的计算区域，尾流场设为4倍船长，首部和四周为3倍船长，并在AUV周围建立一个球形加密区。流体域的入口面和四周边界都设为速度入口，出口为压力出口。每个步长中进行5次迭代计算；在边界层设置中，棱柱层数为5，棱柱层延伸为1.5，棱柱层总厚度为10 mm，最终壁面Y+值大于30，边界层设置符合要求。逐步渐进地对AUV周围网格进行加密，并对十字控制舵和声呐附近单独加密，进行网格无关性检验之后，最终生成的网格总数约为64万。网格偏斜角是指面网格面积矢量 $a$ （面法向），偏斜角过大会导致求解器收敛性降低，降低计算的精度和稳定性，因此在网格的整个生成过程中，要尽量保证网格单元偏斜角小于85°。

2 定常运动数值模拟

定常运动主要包括匀速直航运动和斜航运动。

2.1 匀速直航运动

匀速直航运动时，AUV的纵向力X可表示为：

$X = {X_{uu}}{u^2} 。$

(1)

本文研究的AUV最大航速为5 kn，故平移速度的范围为0.3～3.0 m/s，记录迭代收敛后的阻力数据,使用Matlab进行数据拟合之后，获得的水动力系数 ${X_{uu}} = - 1.669$ 。

2.2 斜航运动

斜航运动属于平面运动，根据对不同平面划分可分为水平面斜航运动和垂直面斜航运动。数值模拟水平面斜航运动则将AUV绕z轴旋转，使其具有一定的漂角β，数值模拟垂直面斜航运动则将AUV绕y轴旋转，使其具有一定的攻角α。

2.2.1 水平面斜航

设航行器o $\xi$ 轴的前进速度V=1.54 m/s，通过改变漂角的大小产生不同的横向速度v，分别计算漂角β从−10°到10°，每次间隔2°运动下的横向力Y和转首力矩N。水平面斜航运动时，首向角不变，横向力Y和转首力矩N的计算式简化为：

$\left\{ \begin{gathered} Y = {Y_v}v + {Y_{v|v|}}v|v|，\\ N = {N_v}v + {N_{v|v|}}v|v| 。\\ \end{gathered} \right.$

(2)

通过Matlab数据拟合，可测出 ${Y_v}$ ， ${Y_{v{\text{|v|}}}}$ ， ${N}_{v}和{N}_{v\left|v\right|}$ 4个水动力系数，力和力矩的计算结果如图2所示。

拟合之后可得到 ${Y_v} = - 33.2$ ， ${Y_{v{\text{|v|}}}} = - 90.78$ ， ${N_v}= - 22.14$ ， ${N_{v|v|}} = 34.44$ 。

2.2.2 垂直面斜航

通过改变攻角的大小产生不同的垂向速度w，前进速度V=1.54 m/s，分别计算攻角α从−10°到10°，每次间隔2°运动下的垂向力Z和纵倾力矩M。垂直面斜航运动时，纵倾角不变，垂向力Z和纵倾力矩M的计算式为：

$\left\{ \begin{gathered} Z = {Z_w}w + {Z_{w|w|}}w|w|，\\ M = {M_w}w + {M_{w|w|}}w|w|。\\ \end{gathered} \right.$

(3)

通过数据拟合，可测出 ${Z_w}$ ， ${Z_{{\text{w|w|}}}}$ ， ${M}_{w}和{M}_{\text{w|w|}}$ 四个水动力系数，力和力矩的拟合曲线如图3所示。

用Matlab拟合工具箱对垂向力Z和纵摇力矩M进行非线性拟合， ${Z_w} = - 36.89$ ， ${Z_{{\text{w|w|}}}} = - 73.96$ ， ${M_w} = 21.55$ ， ${M_{{\text{w|w|}}}} = - 33.92$ 。

3 非定常运动数值模拟

定常运动无法获得与加速度、角速度相关的水动力系数，相关水动力系数可通过平面运动机构（PMM）模拟AUV非定常运动获取。根据AUV六自由度运动规律，水平面有纯横荡和纯首摇运动，垂直面有纯升沉和纯俯仰运动^[11]。使用动网格技术数值模拟AUV的4个不同工况的运动，在数值模拟中，AUV的航速V=1 m/s，摆动幅值为0.04 m，计算步长设为0.02 s，每个步长迭代5次，分别取0.2 Hz，0.25 Hz，0.312 5 Hz，0.4 Hz，0.5 Hz，0.625 Hz共6个不同的频率进行计算，每个频率都计算6个以上的完整周期。

3.1 水平面非定常数值模拟 3.1.1 纯横荡运动

船模拘束运动的纯横荡运动方程为：

$\left\{ {\begin{array}{*{20}{l}} {\eta = a\sin \omega t} ，\\ {v = \dot \eta = a\omega \cos \omega t}，\\ {\dot v = - a{\omega ^2}\sin \omega t}，\\ {\psi = r = \dot r = 0} 。\end{array}} \right.$

(4)

其中： $\eta$ 为AUV的重心横向位移； $a$ 是横向位移的振幅； $\omega$ 为横荡振荡的角频率。由于 $\eta$ 振幅较小，非线性项可以忽略，因此纯横荡运动下的横向力Y和转首力矩N可表示为：

$\left\{ {\begin{array}{*{20}{l}} {Y = {Y_{\dot v}}\dot v + {Y_v}v}，\\ {N = {N_{\dot v}}\dot v + {N_v}v}。\end{array}} \right.$

(5)

由于单个工况要6个不同频率的数据都是正弦函数，故本文只列举频率为0.625 Hz时的数据，横向力Y和转首力矩N曲线如图4所示。

3.1.2 纯首摇运动

在PMM设置中，在纯横荡的基础上添转艏角，激活纯横摆之后的首向角方程为 $\psi = {\psi _0}\cos (\omega t)$ ， ${\psi _0} = {{a\omega } / V}$ ，所以纯首摇运动方程为：

$\left\{ \begin{gathered} \psi = {\psi _0}\cos \omega t，\\ r = \dot \psi = - {\psi _0}\omega \sin \omega t，\\ \dot r = - {\psi _0}{\omega ^2}\cos \omega t ，\\ v = \dot v = 0 。\end{gathered} \right.$

(6)

纯首摇运动时横向力Y和转首力矩N可表示为：

$\left\{ {\begin{array}{*{20}{l}} {Y = {Y_{\dot r}}\dot r + {Y_r}r}，\\ {N = {N_{\dot r}}\dot r + {N_r}r} 。\end{array}} \right.$

(7)

纯首摇运动时，频率为0.625 Hz时的横向力Y和转首力矩N曲线如图5所示。

3.2 垂直面非定常数值模拟

在STAR-CCM+中将航行器模型绕着x轴方向旋转90°，用平面运动机构（PMM）模拟AUV的垂直面运动。

3.2.1 纯升沉运动

纯横荡运动方程为：

$\left\{ \begin{gathered} \zeta = {a_\zeta }\sin \omega t，\\ w = \dot \zeta = {a_\zeta }\omega \cos \omega t，\\ \dot w = - {a_\zeta }{\omega ^2}\sin \omega t ，\\ \theta = q = \dot q = 0 。\end{gathered} \right.$

(8)

其中： $\zeta$ 为AUV的重心横向位移； ${a_\zeta }$ 为横向位移的振幅。纯升沉运动下的垂向力Z和纵摇力矩M可表示为：

$\left\{ \begin{gathered} Z = {Z_{\dot w}}\dot w + {Z_w}w，\\ M = {M_{\dot w}}\dot w + {M_w}w 。\end{gathered} \right.$

(9)

纯升沉运动时，频率为0.625 Hz时的垂向力Z和纵摇力矩M曲线如图6所示。

3.2.2 纯俯仰运动

纯俯仰运动方程为：

$\left\{ \begin{gathered} \theta = {\theta _0}\cos \omega t，\\ q = \dot \theta = - {\theta _0}\omega \sin \omega t，\\ \dot q = - {\theta _0}{\omega ^2}\cos \omega t ，\\ w = \dot w = 0 。\end{gathered} \right.$

(10)

纯俯仰运动时垂向力Z和纵摇力矩M可表示为：

$\left\{ \begin{gathered} Z = {Z_{\dot q}}\dot q + {Z_{_q}}q，\\ M = {M_{\dot q}}\dot q + {M_q}q 。\end{gathered} \right.$

(11)

纯俯仰运动时，频率为0.625 Hz时的垂向力Z和纵摇力矩M曲线如图7所示。

4 数值计算误差分析

本文非定常运动数值模拟中，使用至强W-2145CPU(主频3.9GHz)的工作站，经过24 h连续模拟计算完成单个工况6个频率的数据计算，无因式水动力系数与REMUS 100的水动力系数汇总如表1所示。

${Y}_{v}{}^{\prime }，{N}_{v}{}^{\prime }，{Z}_{w}{}^{\prime }和{M}_{w}{}^{\prime }$ ，既可通过垂直面斜航运动获得，又可通过纯横荡运动试验得到。对比表中数据，在纯横荡数值模拟中，除了 ${M_w}^\prime$ 误差稍大，其他3个系数的误差都更小。总体来说，纯横荡运动数值模拟除 ${N_v}^\prime$ 以外，其他水动力系数误差都在30%以内，纯升沉运动数值模拟除 ${M_w}^\prime$ 以外，其他水动力系数误差都在都在22%以内，纯首摇运动数值模拟都在一个数量级，纯俯仰运动的数值模拟除 ${M_{\dot q}}^\prime$ 以外，都在一个数量级。数值模拟计算值与REMUS 100的实际值符号都一致，计算数据都在一个数量级上，使用CFD技术理论计算水动力具有一定的可行性。

在纯首摇和纯俯仰运动中，用 ${\psi _0} = {{a\omega } \mathord{\left/ {\vphantom {{a\omega } V}} \right. } V}$ 近似为 ${\psi _0} = \arctan ({{a\omega } \mathord{\left/ {\vphantom {{a\omega } V}} \right. } V})$ ，为了满足近似关系， ${{a\omega } \mathord{\left/ {\vphantom {{a\omega } V}} \right. } V}$ 的数值就得很小，在实际仿真中，振幅和航速一定时，对应的角频率越小，实际误差就越小。

通过对数据误差进行分析，发现在数值模拟计算的过程中，仿真频率、时间步长和边界层设置等都对最终的计算精度有着很大的影响，其中边界层设置要使壁面Y+值大于30；仿真频率和时间步长越小，计算精度越高。在分析仿真频率对计算精度的影响时，对纯横荡运动进行不同频率组的计算（其他设置一样），将频率取小，取0.1 Hz，0.15 Hz，0.2 Hz，0.25 Hz，0.312 5 Hz和0.4 Hz六个频率（低频率组）进行计算，测出的数据和原来6个频率（高频率组）测出的数据进行对比，测出数据如表2所示。

可以看出，纯横荡运动中，频率变低时，水动力系数的误差都变小。试验数据表明，仿真频率的减少提升了数值模拟计算的精度。

5 结　语

本文主要对AUV水动力性能进行研究，通过AUV的定常运动数值模拟，改变AUV的攻角和漂角，分别进行直航运动、斜航运动数值模拟计算；分析非定常运动的性质，通过平面运动机构（PMM）直接模拟航行器的纯横荡、纯首摇运动。在纯升沉和纯俯仰的数值模拟当中，将AUV模型绕着x轴方向旋转90°，利用PMM运动机制，完成垂直面的数值模拟，获取航行器的力和力矩，通过Matlab进行数据拟合，从而获取AUV的水动力系数，并于REMUS 100实际水动力系数进行对比。在对数据误差分析之后，发现在数值模拟计算的过程中，仿真频率、时间步长和边界层设置等都对最终的计算精度有着很大的影响。通过对纯横荡运动不同频率组进行计算，计算数据表明，低频率的数值模拟有利于提升数值计算精度。通过误差分析以及采用低频率组数据的初步验证，为进一步提升数值模拟计算精度奠定了基础。