0 引　言

水下作业机器人（ROV）在海洋科学、资源勘探和海洋工程等领域发挥着越来越重要的作用。水下作业机器人的设计和操作受到水动力效应的显著影响，水动力系数是描述水下作业机器人在水中运动时所受到的阻力、升力和力矩等水动力效应的重要参数。而水下作业机器人的机箱、浮体、螺旋桨以及所有作业工具直接安装于机器人框架上，导致其水动力特性复杂。因此，机器人水动力特性的准确计算对机器人控制系统设计以及其运动性能预报至关重要，如何有效精确获取ROV的水动力系数是水下作业机器人运动性能研究的重要方向。

目前获取水动力系数的方法有经验公式计算方法、操纵性模型试验方法和计算流体力学（CFD）数值模拟方法。经验公式计算实用性和经济性强，但是难以保证精度和适用性；操纵性模型试验方法准确性高，但是试验周期长、成本高，不适用于水下作业机器人的初期设计；随着计算机和流体力学的发展，以计算机为基础的CFD软件为水下作业机器人计算水动力提供了一种新的方法。CFD数值模拟方法是目前为止使用广泛、成本低、周期短、比较可靠的一种获取水动力系数的方法。Guo等^[1]采用STAR-CCM+对四自由度船体进行直航运动、斜航运动和回转运动开展数值模拟，将拟合得到相应的水动力系数运用到“Z”字形运动仿真中。Liu等^[2]以KCS集装箱船为研究对象，对其进行了拘束模试验模拟，得到了船模所有线性和非线性水动力系数，并对其操纵性进行预测。刘冬雨等^[3]采用STAR-CCM+软件对无人双体船进行PMM（平面运动机构）运动仿真试验，通过对比不同拟合方式下的运动仿真轨迹的吻合度，选取出最符合实际的拟合方法。刘晨飞等^[4]采用CFD技术对KVLCC2船模进行拘束模试验，获得精确的水动力系数，利用龙格-库塔法对其进行操纵性运动仿真。黎欣昌等^[5]利用STAR-CCM+软件模拟AUV的小振幅PMM运动，将数值计算的水动力系数与REMUS100实际系数进行对比。结果表明，低频振荡的数值计算有助于提高水动力系数的精度。

本文聚焦于水下作业机器人的水动力系数和运动预报。首先，采用软件STAR-CCM+对复杂构型水下作业机器人开展直航运动、斜航运动和平面运动机构运动模拟，得到水下作业机器人简单运动下各自由度上的水动力和力矩；然后，针对不同的水动力导数型式，采用最小二乘法进行水动力导数辨识；最后，基于CFD软件计算的水下作业机器人二维回转运动的水动力计算结果，探讨不同水动力导数型式下复杂构型水下作业机器人水动力模型的精确性。

1 数值方法与模型建立 1.1 CFD数值模拟方法

模拟水下作业机器人在无界流中运动时，将周围流体认为是不可压缩的粘性流，通过有限体积法求解流体控制方程获得流域解，CFD中的连续性方程和雷诺平均N-S方程^[6]如下：

$\frac{{\partial {u_i}}}{{\partial {x_i}}} = 0 ，$

(1)

$\rho \frac{{\partial {u_i}}}{{\partial t}} + \rho {u_j}\frac{{\partial {u_i}}}{{\partial {x_j}}} = - \frac{{\partial p}}{{\partial {x_i}}} + \mu \frac{\partial }{{\partial {x_j}}}\left( {\frac{{\partial {u_i}}}{{\partial {x_j}}}} \right) + {S_i} 。$

(2)

式中：u_i和u_j为速度分量；x_i和x_j为位移分量；ρ为流体密度；p为流体微元上的正压力；μ为动力粘性系数；t为时间；S_i为广义源项。

1.2 计算对象

本文研究对象是一款采用扁平化设计的复杂构型水下作业机器人，其总体外观如图1所示，扁平化设计使之具有更灵巧的操纵性。8个螺旋桨推进器实现ROV水下六自由度运动，水平方向的4个推进器对称分布，竖直方向的4个推进器均匀分布在顶部4个角。浮力材料可为水下作业机器人提供100 kg的可变载荷，调整载荷大小，使浮力和重力尽可能的相等。ROV的基本尺寸如表1所示。

1.3 数值计算域

为了便于计算，对上述ROV原型进行简化，基于STAR-CCM+软件建立三维数值水池模型。坐标原点位于ROV质心位置。流体域尺寸为：−3L<x<4L，−3L<y<3L，−3L<z<3L（其中L为ROV的长度）。如图2所示，将ROV的前进方向面设为速度入口，其对立面设为压力出口，数值水池上下表面设为无滑移壁面，其左右表面设为对称平面，ROV表面设为无滑移壁面。

采用可实现k-ε湍流模型求解连续性方程和动量方程。采用重叠网格模拟ROV的全局运动，采用切割体网格对ROV进行网格划分，采用棱柱层网格模拟边界层，为了精确捕捉流场变化，对ROV周围进行局部网格加密，划分后的网格数量约216.7万，网格划分后的ROV模型如图3所示。在模拟单自由度运动时，采用DFBI（Dynamic Fluid Body Interaction）模块计算ROV所受载荷，以对应的来流速度模拟ROV在该方向上的运动速度；在模拟PMM运动和ROV复杂多维运动时，采用动网格技术，对相关边界施加运动条件，从而实现精确模拟相应的运动，以保证模拟的准确性和计算效率。时间求解器采用隐式非稳态，与显式方法相比，隐式方法更稳定，具有较好的耗散性和数值稳定性。采用分离流方法处理边界层产生的流动分离现象。

2 ROV数学模型

为了准确描述ROV的空间位置，本文建立了以水下作业机器人的质心为原点的运动坐标系Oxyz和以地心为原点的惯性坐标系Eξηζ，ROV的运动参数定义如表2所示。

在描述ROV的空间位置时，可以先在惯性坐标系下建立运动方程，然后利用坐标系的转换矩阵完成惯性坐标系Eξηζ向运动坐标系Oxyz下的运动方程转换，获得水下作业机器人的空间运动方程，其转换公式为：

$\left[\begin{gathered}\xi \\ \eta \\ \zeta \\ \end{gathered}\right]=T\left[\begin{gathered}x \\ y \\ z \\ \end{gathered}\right]。$

(3)

式中：T为转换矩阵。

同理，ROV的姿态变换也可以由下式将惯性坐标系转变为运动坐标系：

$\left[ \begin{gathered} {\dot \varphi } \\ {\dot \theta } \\ {\dot \psi } \\ \end{gathered} \right] = \left[ {\begin{array}{*{20}{c}} 1&{\sin\varphi \tan\theta }&{\cos\varphi \tan\theta } \\ 0&{\cos\varphi }&{ - \sin\varphi } \\ 0&{{\hbox{${\sin\varphi }$} \mathord{\left/ {\vphantom {{\sin\varphi } {\cos\theta }}}\right.} \hbox{${\cos\theta }$}}}&{{\hbox{${\cos\varphi }$} \mathord{\left/ {\vphantom {{\cos\varphi } {\cos\theta }}}\right.} \hbox{${\cos\theta }$}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} p \\ q \\ r \end{array}} \right]。$

(4)

式中：$\dot{\phi }、\dot{\theta }、\dot{\psi }$为ROV分别绕3个坐标轴旋转的角度对时间t的导数，即角速度。

水下作业机器人是一个非线性的复杂系统，其在水中受到的力包括：水动力、静力、推进器力、脐带缆作用力以及外部环境的干扰力，ROV的六自由度动力学方程可以表示为：

$\boldsymbol{M\dot{v}}+\boldsymbol{C}\left(\boldsymbol{v}\right)\boldsymbol{v}+\boldsymbol{D}\left(\boldsymbol{v}\right)\boldsymbol{v}+\boldsymbol{g}\left(\eta\right)=f。$

(5)

式中：M为ROV的质量矩阵；v为ROV六自由度上的速度和角速度矩阵；$\dot v$为ROV六自由度上的加速度和角加速度矩阵；C(v)为附加科氏向心力矩阵；D(v)为ROV的阻尼矩阵；g(η)为静力矩阵；f为ROV在水下受到的所有力的总和。

水下作业机器人在做二维回转运动时，只考虑纯纵荡、纯横荡和纯首摇3种模态下的一阶水动力导数，以此建立机器人的水动力方程^[7]：

${\boldsymbol{X}}={\boldsymbol{X}}_0+X_{\dot{u}}\dot{u}+X_uu+X_{\dot{v}}\dot{v}+X_vv+X_{\dot{r}}\dot{r}+X_rr，$

(6)

${\boldsymbol{Y}}={\boldsymbol{Y}}_0+Y_{\dot{u}}\dot{u}+Y_uu+Y_{\dot{v}}\dot{v}+Y_vv+Y_{\dot{r}}\dot{r}+Y_rr，$

(7)

${\boldsymbol{N}}={\boldsymbol{N}}_0+N_{\dot{u}}\dot{u}+N_uu+N_{\dot{v}}\dot{v}+N_vv+N_{\dot{r}}\dot{r}+N_rr。$

(8)

式中：${X}_{\dot{u}}、{X}_{\dot{v}}、{X}_{\dot{r}}、{Y}_{\dot{u}}、{Y}_{\dot{v}}、{Y}_{\dot{r}}、{N}_{\dot{u}}、{N}_{\dot{v}}、{N}_{\dot{r}}$为不同运动下的加速度系数；${X}_{u}、{X}_{v}、{X}_{r}、{Y}_{u}、{Y}_{v}、{Y}_{r}、{N}_{u}、{N}_{v}、{N}_{r}$为不同运动下的速度系数；$\boldsymbol{X}_0、\boldsymbol{Y}_0、\boldsymbol{Z}_0$为ROV在静水状态下的受力；$\dot u$为纵向加速度；$\dot v$为横向加速度；$\dot r$为绕z轴旋转的角加速度。

然而，对于复杂构型的水下作业机器人，仅考虑一阶水动力系数所建立的水动力模型往往不够精确，一般需要综合考虑非线性的二阶水动力系数或者其他运动下的水动力系数。

3 数值模拟计算结果

水动力系数包括与加速度、角加速度相关的惯性类水动力系数和与速度、角速度有关的粘性类水动力系数。粘性类水动力系数通过直航、斜航运动获取，惯性类水动力系数通过PMM运动获取^[8]。

3.1 直航运动

直航阻力试验是根据试验模型在不同流速下受到的纵向阻力来确定直航阻力系数。工况设置速度入口处来流速度范围为0.5～1.5 m/s，每隔0.25 m/s取一个工况。计算得到直航运动阻力，阻力的计算公式如下：

$C_d=\frac{\boldsymbol{\boldsymbol{F}}}{\dfrac{1}{2}\rho V^2L^2}。$

(9)

式中：F为直航阻力；ρ为液体密度，在这里取1030 kg/m³，L为机器人的特征长度，即1.8 m；V为来流速度；${C_d}$为阻力系数。

数值模拟结果如表3所示，可以看出ROV的直航阻力大小与来流速度呈正相关。根据表3给出的结果，利用Matlab基于最小二乘法进行二次曲线拟合，拟合曲线如图4所示。所得拟合曲线方程取二阶导数，再将其无量纲化处理得到水平面直航阻力系数$X_{uu}^{'} = 0.58$。

3.2 斜航运动

水平面斜航试验数值模拟就是将水下作业机器人在流体计算域中绕z轴旋转一定的角度，即漂角β。给定来流速度V，此时，在运动坐标系下，水下作业机器人分别有纵向和横向的速度分量^[9]：

$\left\{ {\begin{array}{*{20}{l}} {u = V\cos \beta } ，\\ {v = V\sin \beta } 。\end{array}} \right.$

(10)

式中：u、v分别为ROV在纵向和横向上的速度分量；V为来流速度；$\beta$为漂角。

改变漂角β的大小，相当于机器人以相同的水平速度V和不同的横向速度v运动，得到机器人在y方向上的横向力Y和绕z轴的横向力矩N。将来流速度V分别设为0.5、0.75、1、1.25、1.5 m/s，每个速度对应的漂角β区间在[−10°,10°]之间，增量为2°，共55个工况。

对STAR-CCM+中监测得到的横向力和横向力矩数据用Matlab基于最小二乘法进行三次曲线拟合，得出曲线方程，求出曲线在v=0处的斜率，得到相对应的水动力系数。通过回归分析，得到水平面斜航运动无因次化水动力系数，如表4所示。

3.3 PMM运动

平面运动机构试验是通过给定来流速度，在机器人的6个自由度上增加正弦运动，改变不同的振荡频率，模拟水下作业机器人做纯升沉、纯首摇、纯横荡、纯横摇、纯俯仰和纯纵荡运动，获取其对应的水动力系数。

3.3.1 纯横荡运动

在PMM运动中，纯横荡运动是在x方向给定一个来流速度V=0.8 m/s，同时在y方向上增加一个正弦运动，模拟振荡频率f=0.2、0.4、0.6、0.8、1.0时的5个周期的升沉运动。纯横荡运动的运动参数方程如下：

$\left\{ {\begin{array}{*{20}{l}} {y = {y_0}\sin \omega t} ，\\ {v = \dot y = {y_0}\omega \cos \omega t} ，\\ {\dot v = \ddot y = - {y_0}{\omega ^2}\sin \omega t}，\\ {p = \dot p = 0} 。\end{array}} \right.$

(11)

式中：y₀为振荡的最大振幅，取0.1；v为横向速度；p为横摇角；$\dot p$为横摇角速度；ω为圆频率；t为时间。

将ROV纯横荡运动的方程进行简化，只考虑横向力和力矩，其水动力方程为：

$\left\{\begin{array}{*{20}{l}}\boldsymbol{Y}=\boldsymbol{Y}_0+Y_{\dot{v}}\dot{v}+Y_vv，\\ \boldsymbol{N}=\boldsymbol{N}_0+N_{\dot{v}}\dot{v}+N_vv。\end{array}\right.$

(12)

将参数方程代入水动力方程，并将其进行无因次化处理，得到纯横荡运动无因次化水动力方程：

$\left\{ {\begin{array}{*{20}{l}} {{\boldsymbol Y^{'}} = \boldsymbol Y_0^{'} + {Y_a}\sin \omega t + {Y_b}\cos \omega t} ，\\ {{\boldsymbol N^{'}} = \boldsymbol N_0^{'} + {N_a}\sin \omega t + {N_b}\cos \omega t} 。\end{array}} \right.$

(13)

式中，${Y_a} = - \displaystyle\frac{{{y_0}{\omega ^2}{L^2}}}{{{V^2}}}Y_{\dot v}^{'}$；${Y_b} = \displaystyle\frac{{{y_0}\omega L}}{V}Y_v^{'}$；${N_a} = - \displaystyle\frac{{{y_0}{\omega ^2}{L^2}}}{{{V^2}}}N_{\dot v}^{'}$；${N_b} = \displaystyle\frac{{{y_0}\omega L}}{V}N_v^{'}$$\boldsymbol{Y}^{'}、\boldsymbol{N}^{'}$分别为无因次化后的横向力和回转力矩；$\boldsymbol{Y}_0^{'}、\boldsymbol{N}_0^{'}$分别为静水状态下无因次化后的横向力和回转力矩。

将在STAR-CCM+中数值模拟监测得到的横向力和力矩进行无量纲化，其无量纲化数据随时间的周期变化如图5所示。

使用Matlab中的工具箱cftool对不同频率下模拟得到的力和力矩进行傅里叶级数拟合，拟合曲线如图6所示，拟合得到${Y}_{a}、{Y}_{b}、{N}_{a}、{N}_{b}$如表5所示。

傅里叶级数拟合得到的一系列${Y}_{a}、{Y}_{b}、{N}_{a}、{N}_{b}$，再利用Matlab基于最小二乘法进行三次曲线拟合，得到的曲线方程的一阶导数即为线性水动力系数，拟合曲线如图7所示，表6为纯横荡运动的水动力系数。

3.3.2 纯艏摇运动

给定入口速度V=0.8 m/s，同时施加一个绕z轴旋转的正弦运动，开展纯首摇运动模拟，与纯横荡运动类似，其运动方程参考文献[7]的表述。

图8为ROV纯艏摇运动下振荡频率f=1时一个周期内的速度分布情况，整体上看ROV在各个时段的最大速度值相近，最小速度都出现在了ROV的顶部，但是在一个周期内的速度还是出现了明显的周期性变化。

将在STAR-CCM+中采集到的数据利用Matlab基于最小二乘法进行三次曲线拟合，得到的曲线方程的一阶导数即为线性水动力系数，拟合曲线如图9所示。根据拟合曲线得到纯首摇运动无因次水动力系数，如表7所示。

4 复杂运动结果分析与验证 4.1 CFD网格收敛性分析

数值计算网格数量和网格质量都会直接影响数值模拟计算结果。本文采用ITTC^[10]推荐的方法对ROV复杂运动条件下的计算域进行网格收敛性分析。在动网格技术下，以恒定角速度r=0.1 rad/s做二维回转运动，回转半径为5 m^[11]。

首先选择一个基础网格尺寸，然后再用基础网格的$\sqrt 2$倍进行网格细化，最终划分了162万、315万、442万等3套粗、中、细网格用于收敛性分析。图10为使用3套不同网格模拟ROV的纵向力和横向力的时程曲线。与一般船体结构不同，本文研究的水下作业机器人结构复杂，容易在众多区域产生微小漩涡，从而导致了复杂运动中受力的不稳定。如图10所示，不同网格的受力曲线仍存在不同程度的波动。为了更直观地比较不同网格对受力的影响，对所得的稳定状态下的纵向力和横向力取平均值，即图10中的水平线。由图可知，中网格和细网格的数值相近，为了计算的精确度，本文采用细网格所得结果进行验证。

4.2 二维回转运动受力预报

根据4.1节的分析，本文采用422万的细网格数量进行模拟ROV在二维水平面上的回转运动，计算结果取平均值，纵向力X为−182.5 N，横向力Y为−17.7 N，回转力矩N为−24.4 N。为探讨不同水动力系数型式所建立的水动力模型的精确性，将在STAR-CCM+中用简单的运动模拟得到的水动力系数代入到ROV的运动学方程中，使用式(6)～式(8)在Matlab中计算复杂的二维回转运动中水下作业机器人受到的力和力矩。表8为在计算水动力时考虑不同水动力系数型式的3种方案。

图11为不同方案下Matlab水动力系数计算结果及其误差变化，为了方便展示与比较，将纵向力数据减小为原来的0.1。

结合表9和图11，在使用方案1计算水动力时，纵向力和回转力矩的计算只考虑由自身速度和回转角速度引起的一阶和二阶水动力系数，计算结果误差在13%左右。横向力需同时考虑纵向速度和角速度，但由于其较小，对水动力系数的变动较敏感。

现有的研究多集中于结构对称的水下作业机器人，通常忽略了$\boldsymbol{X}_0、\boldsymbol{Y}_0、\boldsymbol{N}_0$的影响。因此，方案2是在方案1的基础上，进一步考虑了水下作业机器人在无扰动情况下的受力$\boldsymbol{X}_0、\boldsymbol{Y}_0、\boldsymbol{N}_0$。研究发现，尽管$\boldsymbol{X}_0、\boldsymbol{Y}_0、\boldsymbol{N}_0$的数值较小，但是影响不可忽视。将纵向力X增加了$\boldsymbol{X}_0$项后，误差降低了3%；而横向力和回转力矩在增加$\boldsymbol{Y}_0、\boldsymbol{N}_0$项后，计算精度也有所提高。

方案3是在方案2的基础上增加二阶非线性项。在计算纵向力和横向力时，分别增加二阶水动力系数${X_{rr}}$和${Y_{uu}}$可减少约1%的误差。在计算回转力矩时，增加由纵向速度引起的一阶回转水动力系数${N_u}$项后，可显著提高计算精度，这是由于系数${N_u}$对纵向流体载荷更敏感。在本计算中，加入一阶项${N_u}$后计算结果已经足够准确，且速度的平方项影响较小，可以不用考虑二阶项${N_{uu}}$的影响。

5 结　语

本文基于CFD软件STAR-CCM+，针对结构复杂的水下作业机器人，开展直航运动、斜航运动、平面运动机构运动和二维回转运动模拟，并通过Matlab搭建四自由度ROV运动仿真平台，所得结论如下：

1）采用CFD软件模拟低速、低频下的直航、斜航、纯横荡、纯艏摇运动，能够获得相应的水动力和力矩，进而拟合得到复杂构型水下作业机器人的惯性类水动力系数和粘性类水动力系数。使用得到的水动力系数能够实现四自由度ROV运动仿真平台的搭建，通过比较，证明所得水动力系数的可靠性。

2）搭建复杂构型的水下作业机器人运动仿真平台时，需要考虑无扰动时的力和力矩。计算纵向力时，需要考虑由回转角速度引起的二阶纵向力系数${X_{rr}}$项。在计算横向力时，应综合考虑由纵向速度和回转角速度引起的水动力系数${Y_u}$和${Y_r}$。在计算回转力矩时，无需考虑二阶速度系数的影响，仅考虑一阶纵向运动速度系数${N_u}$可以将计算误差降低到4.49%左右。

本文为复杂构型水下作业机器人获取水动力系数提供了一种有效的方法，为复杂构型水下作业机器人的运动仿真及操纵性预报奠定了基础。