0 引　言

带缆水下机器人系统的运动控制与水动力分析一直是研究的热点，建立准确的数学模型对控制律的设计以及水动力变化规律的分析至关重要；同时，水下机器人在水中运动过程中会受到脐带缆张力、螺旋桨推进力及水动力荷载等各种外力的综合作用，受力情况变幻莫测。带缆水下机器人在水中的运动，主要通过脐带缆收放动作和调节螺旋桨转速来实现的，在脐带缆张力与螺旋桨推进力共同作用下，机器人运动过程中又会受到水动力荷载的作用，因而水下机器人控制系统发出的控制力与其所受到的水动力之间存在着复杂的内在耦合关系。探究水动力和控制力之间的内在联系有助于更好地理解带缆水下机器人系统完整的受力状态和水动力变化规律，同时也可以对带缆水下机器人系统向更加精确的运动控制层面提供一定的借鉴。如何对一定控制操纵下水下机器人的水动力进行定量分析和准确预测是现阶段的研究难点之一。

在路径规划和轨迹跟踪问题中，经典控制算法已经得到了广泛的应用，各类新型混合控制算法如雨后春笋般百花齐放，对带缆水下机器人控制系统的设计和应用达到了空前高度与热度。其中较为成熟的有水下机器人定深运动控制、航向和位置控制以及推进器速度控制的PID算法^{[1 − 2]}；轨迹跟踪问题中的混合滑动模态控制^{[3 − 5]}；用于参数调节的自适应控制^[6]；控制系统稳定性分析的鲁棒控制等等^[7]。带缆水下机器人在水中的运动包含3个维度、6个自由度，关于水下机器人的水动力分析中，对系统数学模型又有较高的计算精度要求。例如通过试验确定相关水动力参数^{[8 − 9]}、对水下机器人系统现有的数学模型不断完善与改进^{[10 − 11]}等，相对而言对控制力和水动力之间内在联系的研究并不多见。在不同的控制方法下，机器人运动过程中水动力响应复杂多变，而对模型精度要求越高越高的今天，随着计算机技术高速发展，数值仿真模拟计算的精度和速度远超以往，数值计算凭借可以克服环境干扰、人为操作因素以及难以在实际实验中获得的水动力参数变化规律等优势，一跃成为带缆水下机器人系统水动力计算分析的主流方法，目前仍为诸多研究学者使用。

通过回顾以往文献可知，当前尚未有相关研究将控制力与水动力联系在一起。本文设计了一种带缆水下机器人系统控制力与水动力交互耦合的计算方法，通过Fortran语言、控制程序设计和水动力计算软件Fluent三者之间的混编，在轨迹跟踪仿真模拟计算过程中共享实时计算数据，结合有限差分法和有限体积法实时更新带缆水下机器人系统位移偏差、运动参数和水动力荷载，实现了带缆水下机器人系统在运动控制过程中水动力响应的精准定量分析。为带缆水下机器人系统在一定控制操纵下水动力响应的变化规律起到了一定的预测作用。

1 数学模型

本文研究的带缆水下机器人系统主要由脐带缆、机器人主体和螺旋桨3部分组成，构建数学模型时，需要综合考虑三者之间水动力因素的耦合影响，即对脐带缆、机器人主体和螺旋桨三者运动状态以不同的数学方程描述，再耦合成完整的带缆水下机器人系统水动力数学模型。本文建立的带缆水下机器人系统简图如图1所示。

数学模型的构建从脐带缆局部坐标系(t, n, b)导出的控制方程开始，上末端与固定坐标系(X, Y, Z)中拖船速度耦合，下末端拖点位置在机器人局部坐标系(x, y, z)中与机器人运动速度耦合，即脐带缆上、下末端与拖船、机器人速度耦合关系分别为：

$ \left( {{v_{{t}}},{v_{{n}}},{v_{{b}}}} \right) = \left( {{V_{{X}}},{V_{{Y}}},{V_{{Z}}}} \right){\boldsymbol{D}} ，$

(1)

$ {\boldsymbol{ED}}{{\boldsymbol{V}}_{{T}}} = \left[ {{\boldsymbol{V}} + {\boldsymbol{\Omega }} \times {{\boldsymbol{R}}_{{T}}}} \right]。$

(2)

式中，脐带缆在上末端点三维速度分量(v_t, v_n, v_b)与拖船在该点三维线速度(V_X, V_Y, V_Z)通过矩阵D转换后相一致；脐带缆下末端点（拖点）处三维线速度V_T，经过转换矩阵D和E二次转换后，与拖点处水下机器人三维线速度V=(u, v, w)加上机器人三维角速度Ω=(p, q, r)沿拖点坐标R_T=(x_T, y_T, z_T)方向转换之和相一致。D、E分别为拖船固定坐标系向脐带缆局部坐标系转换矩阵和拖船坐标系向机器人坐标系转换矩阵，具体展开分别为：

$ {\boldsymbol{D}} = \left[ {\begin{array}{*{20}{c}} {\cos \vartheta \cos \varphi }&{ - \cos \vartheta \sin \varphi }&{\sin \vartheta } \\ {\sin \vartheta \cos \varphi }&{ - \sin \vartheta \sin \varphi }&{ - \cos \vartheta } \\ {\sin \varphi }&{\cos \varphi }&0 \end{array}} \right] ，$

(3)

$ {\boldsymbol{E}} = \left[ {\begin{array}{*{20}{c}} {\cos \theta \cos \psi } \\ { - \sin \psi \cos \phi + \sin \phi \sin \theta \cos \psi }\\ {\sin \phi \sin \psi + \cos \phi \cos \psi \sin \theta } \end{array}} \right. \qquad\qquad\qquad$

$\quad \quad\left.{\begin{array}{*{20}{c}} {\cos \theta \sin \psi }&{ - \sin \theta } \\ {\cos \phi \cos \psi + \sin \phi \sin \theta \sin \psi }&{\sin \phi \cos \theta } \\ { - \sin \phi \cos \psi + \cos \phi \sin \theta \sin \psi }&{\cos \phi \cos \theta } \end{array}} \right] 。$

(4)

式中：$ \vartheta $、φ为脐带缆局部坐标系与拖船固定坐标系之间的方向角；θ、ϕ、ψ分别为水下机器人的纵摇角、横摇角和艏摇角。

脐带缆控制方程按照式(5)～式(6)给出^[12]

$ {\boldsymbol{MY'}} = {\boldsymbol{N}}\dot{\boldsymbol{Y}}+ {\boldsymbol{Q}}, $

(5)

$ {\boldsymbol{Y}} = {\left( {T,{v_{\text{t}}},{v_{\text{n}}},{v_{\text{b}}},\vartheta ,\phi } \right)^{\text{T}}},{\text{ }}{\boldsymbol{Y}}' = \frac{{\partial {\boldsymbol{Y}}}}{{\partial s}},\dot{\boldsymbol{Y}}= \frac{{\partial {\boldsymbol{Y}}}}{{\partial t}} 。$

(6)

式中：T为脐带缆张力；s为脐带缆原始长度；t为时间。M、N、Q均为系数矩阵，具体展开式为^[13]

$ {\boldsymbol{M}} = \left[ {\begin{array}{*{20}{c}} 1&0&0&0&0&0 \\ 0&1&0&0&{{v_{{b}}}\cos \varphi }&{ - {v_{{n}}}} \\ 0&0&1&0&{ - {v_{{b}}}\sin \varphi }&{{v_{{t}}}} \\ 0&0&0&1&{ - {v_{{t}}}\cos \varphi + {v_{{n}}}\sin \varphi }&0 \\ 0&0&0&0&{ - T\cos \varphi }&0 \\ 0&0&0&0&0&T \end{array}} \right], $

(7)

$ {\boldsymbol{N}} = \left[ {\begin{array}{*{20}{c}} { - \displaystyle\frac{{em{v_{{t}}}}}{{1 + eT}}} & m & 0 & 0 & {{m_1}{v_{{b}}}\cos \varphi } & { - {m_1}{v_{{n}}}} \\ e & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & {1 + eT} \\ 0 & 0 & 0 & 0 & { - (1 + eT)\cos \varphi } & 0 \\ { - \displaystyle\frac{{e{m_1}{v_{{b}}}}}{{1 + eT}}} & 0 & 0 & {{m_1}} & {{m_1}{v_{{n}}}\sin \varphi - m{v_{{t}}}\cos \varphi } & 0 \\ { - \displaystyle\frac{{e{m_1}{v_{{n}}}}}{{1 + eT}}} & 0 & {{m_1}} & 0 & { - {m_1}{v_{{b}}}\sin \varphi } & {m{v_{{t}}}} \end{array}} \right], $

(8)

$ {\boldsymbol{Q}} = \left[ {\begin{array}{*{20}{c}} { - w\sin \varphi + \displaystyle\frac{1}{2}\rho {\text π} d{C_{{t}}}{v_{{t}}}\left| {{v_{{t}}}} \right|{{\left( {1 + eT} \right)}^{1/2}}} \\ 0 \\ 0 \\ 0 \\ {\displaystyle\frac{1}{2}\rho d{C_{{n}}}{{\left( {v_{{n}}^2 + v_{{b}}^2} \right)}^{1/2}}{v_{{b}}}{{\left( {1 + eT} \right)}^{1/2}}} \\ { - w\cos \varphi + \displaystyle\frac{1}{2}\rho d{C_{{n}}}{{\left( {v_{{n}}^2 + v_{{b}}^2} \right)}^{1/2}}{v_{{n}}}{{\left( {1 + eT} \right)}^{1/2}}} \end{array}} \right]。$

(9)

式中：m为脐带缆单位长度质量；m_l=m +ρA，ρ为水的密度，A为未伸长脐带缆的横截面积；e=1/EA，E为脐带缆的杨氏模量；w为每单位脐带缆长度的淹没重量；d为脐带缆直径。

本文利用有限差分方法将式(5)在空间和时间上进行离散，空间上将脐带缆沿长度方向划分成一系列连续的微段单元，时间上将整个仿真时长分割成连续的时间步，得到差分方程^[14]：

$ \begin{split} &\left[ {{\boldsymbol{M}}_{j + 1}^{i + 1} + {\boldsymbol{M}}_j^{i + 1}} \right]\frac{{{\boldsymbol{Y}}_{j + 1}^{i + 1} - {\boldsymbol{Y}}_j^{i + 1}}}{{\Delta {S_j}}} + \left[ {{\boldsymbol{M}}_{j + 1}^i + {\boldsymbol{M}}_j^i} \right]\frac{{{\boldsymbol{Y}}_{j + 1}^i - {\boldsymbol{Y}}_j^i}}{{\Delta {S_j}}} = \\&\;\; \left[ {{\boldsymbol{N}}_{j + 1}^{i + 1} + {\boldsymbol{N}}_{j + 1}^i} \right]\frac{{{\boldsymbol{Y}}_{j + 1}^{i + 1} - {\boldsymbol{Y}}_{j + 1}^i}}{{\Delta {t^i}}} + \left[ {{\boldsymbol{N}}_j^{i + 1} + {\boldsymbol{N}}_j^i} \right]\frac{{{\boldsymbol{Y}}_j^{i + 1} - {\boldsymbol{Y}}_j^i}}{{\Delta {t^i}}} +\\&\;\; {\boldsymbol{Q}}_{j + 1}^{i + 1} + {\boldsymbol{Q}}_j^{i + 1} + {\boldsymbol{Q}}_{j + 1}^i + {\boldsymbol{Q}}_j^i 。\\[-1pt] \end{split} $

(10)

式中：i为时间步数；j为空间步数；将脐带缆沿长度方向平均分成j段，每一微段长度为△S_j，计算时间步长为△tⁱ，脐带缆节点编号为0～N_d；其中0为上末端点，N_d为下末端点，每一个脐带缆微段节点处包含6个待求未知变量，即Y_j=(T, v_t, v_n, v_b, $ \vartheta $, φ) (j=0～N_d)，因此类似式(10)的算式一共有6N_d个；脐带缆的节点数量一共有(N_d+1)个，每一个节点处都有6个待求未知变量，因此Y_j中共有6(N_d+1)个待求的未知变量。为求解此方程组，引入式(1)和式(2)，此时又有9个多余未知变量V=(u, v, w)、Ω=(p, q, r)以及θ、ϕ、ψ被引入，需要额外补充9个方程才能使方程组闭合可解。机器人六自由度运动方程按照式(11)～(16)展开：

$ {m[\dot u - vr + wq - {x_G}({q^2} + {r^2}) + {y_G}(pq - \dot r) + {z_G}(pr + \dot q)] = {F_{{x}}}, }$

(11)

$ {m[\dot v + ur - wp + {x_G}(pq + \dot r) - {y_G}({p^2} + {r^2}) + {z_G}(qr - \dot p)] = {F_{{y}}}, }$

(12)

$ {m[\dot w - uq + vp + {x_G}(pr - \dot q) + {y_G}(qr + \dot p) - {z_G}({p^2} + {q^2})] = {F_{{z}}}, }$

(13)

$ \begin{split} {I_{{x}}}\dot p +& \left( {{I_{{z}}} - {I_{{y}}}} \right)qr + {I_{{{xy}}}}\left( {pr - \dot q} \right) - {I_{{{yz}}}}\left( {{q^2} - {r^2}} \right) - {I_{{{xz}}}}\left( {pq + \dot r} \right) + \\& m\left[ {{y_{{G}}}\left( {\dot w - uq + vp} \right) - {z_{{G}}}\left( {\dot v + ur - wp} \right)} \right] = {M_{{x}}} ,\\[-3pt] \end{split} $

(14)

$ \begin{split} {I_{{y}}}\dot q +& \left( {{I_{{x}}} - {I_{{z}}}} \right)pr - {I_{{{xy}}}}\left( {qr + \dot p} \right) + {I_{{{yz}}}}\left( {pq - \dot r} \right) + {I_{{{xz}}}}\left( {{p^2} - {r^2}} \right) - \\& m\left[ {{x_{{G}}}\left( {\dot w - uq + vp} \right) - {z_{{G}}}\left( {\dot u - vr - wq} \right)} \right] = {M_{{y}}} ,\\[-3pt] \end{split} $

(15)

$ \begin{split} I_z\dot{r}+ & \left(I_y-I_x\right)pq-I_{xy}\left(p^2-q^2\right)-I_{yz}\left(pr+\dot{q}\right)+I_{xz}\left(qr-\dot{p}\right)+ \\ &m\left[x_G\left(\dot{v}+ur-wp\right)-y_G\left(\dot{u}-vr+wq\right)\right]=M_z。\end{split} $

(16)

在固定惯性坐标系中，水下机器人3个欧拉角的变化速率和在水下机器人局部坐标系中机器人主体角速度的关系为：

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

(17)

式中：(x_G, y_G, z_G)为水下机器人在其局部坐标系中的重心坐标；I_x、I_y、I_z和I_xy、I_yz、I_xz分别为机器人主体质量惯性矩和交叉惯性矩；等式左边表示机器人运动过程中在机器人局部坐标系（非惯性系）所受到的惯性力，等式右边项F₀=(F_x, F_y, F_z)及M₀=(M_x, M_y, M_z)为机器人运动过程中受到的合外力与合外力矩，包含

$ {{\boldsymbol{F}}_0} = {{\boldsymbol{F}}_{{R}}} + {{\boldsymbol{F}}_{{T}}} + {{\boldsymbol{F}}_{{H}}} + {{\boldsymbol{F}}_{{C}}} ，$

(18)

$ \boldsymbol{M}_0=\boldsymbol{M}_R+\boldsymbol{M}_T+\boldsymbol{M}_H+\boldsymbol{M}_C。$

(19)

式中：F_R 为机器人静水回复力；F_H为机器人运动过程中受到的水动力；F_C为螺旋桨发出的控制推进力。脐带缆作为连接工作母船和水下机器人的纽带，担负着传输控制信号与数据的任务，它对机器人运动的影响不可忽略。在机器人局部坐标系中，脐带缆张力F_T和力矩分别为：

$ \mathit{{\boldsymbol{F}}_{T}} = - {\boldsymbol{ED}}T，$

(20)

$ \mathit{\boldsymbol{M}_{T}}=\mathit{\boldsymbol{r}_{t}}\times\mathit{\boldsymbol{F}_{T}}。$

(21)

式中：r_t为脐带缆与机器人主体连接点到机器人主体重心之间的距离。

2 算法设计

式(18)中，当水下机器人重心和浮心不在同一铅垂线上时，静水回复力F_R就会被动产生；F_T通过收放缆动作结合求解方程(10)获得；当水下机器人在水中运动时，会被动受到作用在机器人主体上的水动力荷载F_H，与机器人主体运动的线速度、角速度及螺旋桨旋转产生的流场均有关系；螺旋桨推进力F_C则是通过轨迹跟踪控制算法输出的转速增量在计算流体力学软件Fluent中计算得到。从式(18)和式(19)中可以看出，水下机器人运动过程中主体受到的合外力，既包含了主动调节作用的控制力，又包含了未知被动承受的水动力，而且水动力的变化与主动力的作用方式息息相关。因此，脐带缆张力、机器人水动力与螺旋桨推进力之间存在着复杂的耦合关系，本文通过程序混编计算来探寻三者之间的内在联系。

2.1 基于Fortran语言的耦合方程组求解程序

在Fluent中求解带缆水下机器人系统水动力时，需要预先获得脐带缆、机器人主体以及控制螺旋桨的运动参数，这就需要对脐带缆和机器人主体耦合动力方程组进行编程计算求解。本文以Visual Studio集成编译平台Fortran语言模块对脐带缆和机器人耦合动力方程组进行系统求解，获得脐带缆与机器人实时运动参数，然后再将数据保存并传输至Fluent软件中计算水动力。

Fortran主程序与Fluent水动力计算模块混编数值求解过程中，利用程序执行判断文件“JUDGE.txt”来实现系统方程组求解程序、水动力计算模块以及控制程序三者之间的运行或暂停。当“JUDGE.txt”判断文件中的变量值等于“开”时，Fortran主程序求解脐带缆和水下机器人耦合方程组，获得脐带缆和机器人实时运动参数，然后控制程序执行脐带缆前馈控制律和螺旋桨转速PID控制律，根据目标运动轨迹坐标输出脐带缆收放长度及螺旋桨转速增量，期间Fluent水动力计算模块暂停运行；Fortran主程序与控制程序执行完毕之后，保存计算结果，同时改写“JUDGE.txt”文件中的变量为“关”，计算流程如图2所示。

2.2 控制程序计算

脐带缆收放过程会产生一定的张力作用于机器人主体，限制机器人在深度方向上浮或下潜；螺旋桨转速的改变则会对对机器人运动状态提供合适的推进力，推动机器人沿指定路径运动，二者共同组成水下机器人的整体控制机构。

前馈控制方法高度依赖于被控对象精确的数学模型，调节过程仅与输入的干扰项相关，不涉及反馈环节从而大大降低系统响应延迟，理论上可以实现无差控制。本文设计水下机器人沿指定轨迹在一定时间内进行跟踪运动，控制程序包含调节脐带缆长度的前馈控制方法以及计算各个螺旋桨转速的PID算法。具体过程为，在每一个时间步开始时首先执行前馈控制律调节脐带缆长度，将相邻2个时间步内机器人指定运动轨迹的位移量作为干扰项给到控制律输入端，前馈算法为消除此干扰而执行脐带缆收放动作；随后Fortran主程序求解机器人与脐带缆耦合动力方程组，获得机器人与脐带缆实时运动参数，计算机器人实时位移偏差；最后PID控制律通过位移偏差输出各个螺旋桨转速的增量以减小位移偏差，直至满足精度要求，至此一个时间步内的控制程序执行完毕。前馈控制及PID方法原理如图3和图4所示。

根据位移偏差设计收放缆前馈调节律和控制螺旋桨转速的PID控制律，数学表达式为：

$ {\mathrm{d}}l\left( {nt} \right) = {\mathrm{d}}l\left( {nt - 1} \right) + \delta \left[ {P\left( {nt} \right) - P\left( {nt - 1} \right)} \right]，$

(22)

$ {\Delta r\left( {nt} \right) = {K_P}\Delta e\left( {nt} \right) + {K_I}\Delta e\left( {nt} \right) + {K_D}\left[ {\Delta e\left( {nt} \right) - \Delta e\left( {nt - 1} \right)} \right]。} $

(23)

式中：dl为前馈调节过程中脐带缆微段单元的平均长度；nt为时间步；P为指定运动轨迹的位移坐标；δ为修正系数，无量纲；△r为螺旋桨转速的增量；在不同的运动阶段，K_P、K_I、K_D的取值有所不同。

2.3 FLUENT水动力计算

在每个时间步的计算中，脐带缆和机器人的运动参数通过求解二者的耦合动力方程组(5)、式(11)～式(16)获得；控制程序则是通过位移偏差输出脐带缆收放长度和各个螺旋桨转速的增量；而脐带缆张力、机器人主体受到的水动力荷载以及各个螺旋桨发出的推进力都是由Fluent计算得到的，因此通过方程求解主程序、控制程序和Fluent水动力模块三者混编，实现了控制力与水动力的交互耦合计算。

当Fortran主程序与控制程序计算完毕时，Fluent读取开关判断文件“JUDEGE.txt”变量为“关”，获得Fortran主程序计算的脐带缆和水下机器人运动参数、以及控制程序根据指定位置轨迹输出的各个螺旋桨转速增量，通过用户自定义函数（Users Defined Functions，UDF）计算脐带缆、机器人主体及螺旋桨在运动过程中所受的水动力，以此来实现Fluent、Fortran主程序和控制程序的数据共享和交换，完成带缆水下机器人系统既定的运动控制目标。

Fluent软件自身提供了一套宏函数来帮助定义UDF，也称之为Fluent语言，均以DEFINE_为开头，共有三类：通用的、离散相和多相的。本文为实现Fortran主程序与Fluent之间数据互通和共享一共用到了4个UDF宏函数，均为通用宏，数据传输与共享过程具体展开如下，水动力计算流程如图5所示。

首先，读取程序执行判断文件JUDGE.txt的宏函数为：

1号UDF：DEFINE_EXECUTE_AT_END (name)

本文研究的轨迹跟踪问题在每一时间步内Fluent的计算均为流场复杂多变的瞬态求解过程，在每个时间步计算结束之后执行该宏函数来读取程序执行判断文件“JUDGE.txt”中的变量值，以此实现Fluent软件的运行与暂停功能。

其次，Fortran主程序和控制程序计算得到的水下机器人与脐带缆运动参数以及各个螺旋桨转速，将会被保存并传递到Fluent软件中计算水下机器人主体受到的水动力和各个螺旋桨发出的推进力。Fluent软件用来读取这些参数并将参数值赋予对应的几何模型相应构件的宏函数为：

2号UDF：DEFINE_TRANSIENT_PROFILE (name, current_time)

用于定义瞬态运动参数，将Fortran主程序和控制程序计算得到的水下机器人与脐带缆运动参数以及各个螺旋桨的转速值赋予对应的水下机器人系统几何模型各个对应部件，完成参数的赋值过程。

水下机器人按照指定轨迹运动过程中，形心位置坐标会不断变化。Fluent读取水下机器人形心处三维线速度和角速度、PID控制程序输出的各螺旋桨转速，以及计算机器人主体形心位置水动力和螺旋桨几何中心的推进力，则用到了Fluent的另外2个UDF宏函数：

3号UDF：DEFINE_CG_MOTION (name, dt, vel, omega, time, dtime)

用以指定特定动态区域运动，以及指定每一个时间步内的线速度和角速度，Fluent利用这些速度参数更新动态区域的节点位置。

4号UDF：Compute_Force_And_Moment (domain, thread, x_cg, force, moment, TRUE)

用于捕捉水下机器人主体外壁的力与力矩，其中x_cg是水下机器人的几何中心，force是水下机器人在3个维度受到的水动力，moment是水下机器人运动过程中所受外力对几何中心的力矩，该宏函数实时捕捉水下机器人组合体中心位置x_cg。Fluent实时计算得到水下机器人主体受到的水动力载荷和各个螺旋桨发出的推进力，然后保存并输出，以实现Fluent软件和Fortran主程序以及控制程序之间的数据交换。

3 仿 真 分 析 3.1 模型与算法有效性验证

分别根据已有文献[15 − 16]拖曳体深度控制试验和球形水下机器人水箱位置控制实验，建立对应的几何模型，将几何模型与数学模型关联起来，在与文献[15 − 16]相同的实验条件下进行仿真模拟，对拖曳体的相关参数的监测、对球形水下机器人指定位置信息的控制响应参数分别同实验数据做对比，结果如图6以及表1所示，表1中X轴为水下机器人水平方向前进或后退运动，Z轴为机器人沿纵垂方向的上浮或下潜运动。仿真条件选择标准的k-ε湍流模型、分离求解器SIMPLE算法；设置来流沿X轴负方向、大小为0.2 m/s，时间步长为0.1 s；仿真环境为Inter(R) Xeon(R) CPU E5-2620 v3 at 2.40 GHz。

从图6可知，水下拖曳体深度运动控制过程中，横摇角、纵摇角及深度运动轨迹的仿真模拟结果变化规律与实验结果变化规律相一致，说明通过Fluent水动力软件模拟计算的结果可以较好地反映真实情况；而对于球形水下机器人轨迹控制过程中的关键参数响应，通过对比表1可以发现，系统超调百分比OS%、上升时间T_r、下降时间T_f以及调节时间T_s的数值均比实验值大，是因为仿真计算涉及到螺旋桨旋转流场的变化以及机器人线运动产生周围流场的变化，流场由变化至稳定需要一定的时间积累，这就导致了仿真实验中控制系统响应具有一定的时间延迟性，真是实验中的位置信息则是通过布置在水箱中的相机通过图像处理技术实时监测得到的，不涉及运动流场的变化过程，因此仿真值均超过实验值，但稳态误差e_ss的数值仿真模拟计算结果优于实际实验结果，这是因为数值仿真模拟计算不涉及人为因素，即没有系统误差和随机误差，仿真精度主要以参数设置和计算条件决定，而实际物理实验会存在测量误差与扰动因素等作用的影响，使得稳态误差要高于数值仿真模拟计算。通过2项对比试验可知，通过Fluent软件计算水动力、Fortran主程序求解运动方程组以及控制程序输出的控制量三者混编计算，可以较好地反映真实情况，算法有效性与可靠性得到了验证。

3.2 轨迹跟踪及水动力分析

根据数学模型和控制力与水动力耦合交互算法，对带缆水下机器人系统进行轨迹跟踪仿真模拟，研究带缆水下机器人系统在运动控制过程中脐带缆张力、螺旋桨推进力及机器人主体所受水动力三者之间的关系。

几何模型如图7所示，采用鱼类状主体，两侧对称布置一对导管螺旋桨，中央腔体内布置一个螺旋桨。基本参数如表2所示。

轨迹跟踪给定水下机器人运动路径，仿真模拟中控制操作过程分为3个阶段：

第1阶段（0～3 s）：操纵机器人沿纵垂方向（Z轴）和水平方向（X轴）进行复合运动。其中，以收放缆前馈控制操纵机器人沿Z轴负方向做上浮运动，位移为0.2 m；同时1、2号导管螺旋桨发出的推进力驱动机器人沿X轴负方向以−0.1 m/s的速度匀速前进。

第2阶段（3～9 s）：沿X轴负方向以−0.1 m/s的速度做固定深度的匀速前进。

第3阶段（9～12 s）：通过收放缆操纵机器人沿Z轴正方向向下做下潜运动，位移为0.2 m；同时保持机器人沿X轴负方向以−0.1 m/s的速度匀速前进。

写成数学表达形式：

$\begin{split} &X = 0.8 - 0.1t，{\text{ }}0 < t \leqslant 12 \text{，} \\&Z = \left\{ \begin{gathered} 7.96 - \frac{{0.2}}{3}t \cdot {\alpha _0}，{\text{ }}0 < t \leqslant 3 ，\\ {\text{7}}{\text{.76 }}，3 < t \leqslant 9 ，\\ 7.76 + \frac{{0.2}}{3}t \cdot {\alpha _0}，{\text{ }}9 < t \leqslant 12。\\ \end{gathered} \right.\end{split} $

(24)

仿真运动轨迹及误差如图8所示。可知，机器人水平轨迹的绝对误差始终为正，变化规律为连续3次的先增大后减小，在第2阶段（3～9 s）中的6.5 s附近出现了水平轨迹绝对误差最大值，约为+0.3 m，说明机器人在图8中沿X轴负方向的实时位置坐标始终落后于指定轨迹的位置坐标；机器人沿Z轴方向运动的纵垂轨迹绝对误差同样始终为非负值，说明机器人在上浮运动过程中沿Z轴方向实时位置坐标始终落后于指定轨迹的位置坐标，但是在下潜运动过程中机器人沿Z轴方向实时位置坐标始终领先于指定轨迹的位置坐标。机器人实际运动轨迹在定深运动向下潜运动转变的时刻，由于控制器的超调作用加之机器人的运动惯性，使水平轨迹误差出现了极值，同时实际运动轨迹发生了振荡。

仿真模拟过程中脐带缆长度以及张力随时间变化的关系如图9所示。

从图9(a)可知，脐带缆实际长度随时间的变化大致均分为3个时间段：0～4 s、4～8 s以及8～12 s，与水下机器人实际运动的纵垂轨迹变化规律大体一致。由于PID算法具有一定的延迟性，加之3号螺旋桨的协同推进作用，共同影响了水下机器人脐带缆收放过程，导致脐带缆前馈调节的实际执行时间为4 s，超过了指定时间3 s；8～12 s具有相似的规律出现。

从图9(b)可知，脐带缆张力随时间做周期性振荡变化，周期约为3 s。由于水下机器人轨迹跟踪运动是从流场稳定的静止时刻开始，因此起始阶段脐带缆张力振荡范围较大，约为77～110 kN；一段时间后，水下机器人的运动状态趋于匀速，脐带缆张力振荡范围大大减小，随后由于水下机器人运动轨迹发生变化，导致脐带缆张力振荡更加剧烈，然后又趋于平稳，以此循环。

螺旋桨的转速及其对应的推力如图10所示。

可以看出，螺旋桨转速与其发出的实时推力一一对应，结合图7可知，3个螺旋桨正转时，均会产生沿坐标轴负方向的推进力；螺旋桨反转，推力沿坐标轴正方向，且转速越大，推力越大。

水下机器人轨迹跟踪仿真模拟过程中，实时平动速度如图11所示。

结合图10和图11可以看到：

1）水下机器人实时运动线速度分量与螺旋桨转速的变化有一定的对应关系。由于螺旋桨旋转方向与推力方向时时相反、来流速度设置沿X轴负方向以及水下机器人沿X轴的设定运动速度可知，1、2号导管螺旋桨的推力用来克服来流速度和脐带缆张力沿x轴方向的分量以控制水下机器人按照指定方式运动。

2）水下机器人沿X轴方向指定运动方式为−0.1 m/s的匀速运动，而实际上机器人沿X轴方向平动速度不断变化。在0～3 s的时间段内沿X轴负方向的实时运动速度由零开始逐渐增大，3 s时达到约−0.18 m/s的峰值；由于PID算法具有一定的延迟性与超调性，使得随后的1 s内机器人沿X轴方向实时运动速度减小到约−0.12 m/s；4～9 s，水下机器人沿X轴负方向的运动速度表现为先增大、而后减小再到反向最大、再减小的剧烈变化过程。这是由PID算法调节螺旋桨实时转速的欠灵敏性、以及机器人在变化的流场中运动过程的复杂性决定的。

3）水下机器人沿Z轴方向实时运动速度与3号螺旋桨转速变化规律基本一致，实际运动速度分成4 s一段的三段匀速运动。3号螺旋桨在0～4 s内发出的推力为正，机器人沿纵垂方向做上浮运动；4～8 s螺旋桨转速近似为0，因而无推力，机器人保持定深运动；8～12 s推力为负，机器人做下潜运动；结合图9中脐带缆张力始终大于0可知，3号螺旋桨推力与脐带缆张力始终构成一对相互制衡的平衡力偶，保证水下机器人运动过程中的稳定性。

4）结合脐带缆长度和张力变化特性的结果（见图9）来看，机器人沿X轴方向的运动方式对脐带缆前馈调节的控制效果有时间上的延迟性，而对控制精度影响较小。

在脐带缆收放前馈调节与PID算法控制螺旋桨转速的双重策略下，轨迹跟踪仿真模拟过程中机器人主体受到的水动力载荷如图12所示。

0～4 s阶段，机器人在脐带缆前馈调节与3号螺旋桨协同推进的共同作用下沿Z轴负方向做上浮运动，受到了沿Z轴正方向的水动力载荷；运动初始阶段，因脐带缆前馈调节动作的执行，机器人沿Z轴方向水动力载荷较大，约为20 N，随后迅速降低到10 N左右并缓慢减小，第4 s末降到约7 N；同时在第4 s时刻机器人沿Z轴方向水动力响应出现了振荡，其原因是第4 s末实际收放缆调节动作停止，水下机器人沿Z轴的实际运动速度突然变为0，同时3号螺旋桨转速也减至0，二者共同引起了水下机器人组合体周围流场的突然变化，因此机器人沿Z轴方向水动力载荷也出现了振荡变化；4～8 s时间段，机器人做定深运动，受到的沿Z轴方向的水动力载荷在0附近小范围变化；8～12 s，机器人沿Z轴正方向做下潜运动，受到的水动力载荷为负，同时在第8 s时刻下潜运动初，与第4 s时出现了同样的水动力响应的振荡变化现象，此后稳定在−6 N左右。

关于机器人沿X轴方向的水动力载荷，在0～3 s的时间段内，机器人沿X轴方向指定运动方式为−0.1 m/ s的匀速运动，其中0～2 s机器人从流场稳定条件下开始运动，通过PID控制律1、2号螺旋桨开始旋转输出推力，推动水下机器人沿X轴负方向前进，受到的水动力载荷为正，最大值不超过2 N，相对稳定；2～3 s水动力载荷变大，第3 s时达到16 N左右的峰值；3～4 s由于1、2号螺旋桨转速变小使得机器人沿X轴运动速度减小，因此沿X轴方向的水动力载荷也由16 N下降到6 N左右；4～12 s水下机器人沿X轴方向保持−0.1 m/ s的匀速运动，但由于调节1、2号螺旋桨转速的PID控制律存在一定的时间延迟性以及控制器的超调效应，导致1、2号螺旋桨实际转速在3～9 s出现了剧烈的振荡变化（图10(a)），对应水下机器人沿X轴方向的实时运动速度同样出现了振荡变化；随着螺旋桨转速方向的改变以及机器人运动速度方向的突然改变，沿X轴水动力荷载在第8 s左右由正变负，同时由于机器人组合体周围复合流场是随着机器人运动速度和螺旋桨转速的改变而实时变化的，流场的变化是一个复杂的时间累积过程，因此机器人主体沿X轴方向的水动力载荷在第10 s左右出现了突变：由−10 N突然变化到+18 N然后再减小到+4 N左右。

4 结　语

本文提出一种带缆水下机器人系统水动力与控制力交互耦合计算方法，用于分析带缆水下机器人在一定操作条件下的控制效果和水动力响应特性。为深入探究控制力和水动力的内在规律、以及更加准确地预测水下机器人在控制动作下的动力特性与流场变化分析提供一种思路。

以数学模型、几何模型为基础，通过Fortran语言、控制算法及水动力计算软件Fluent三者混编实现脐带缆、机器人主体以及螺旋桨之间的数据共享和交换，将水下机器人系统运动控制过程中复杂多变的受力状态合理地呈现出来，得到的主要结论有：

1）脐带缆前馈调节过程依赖于水下机器人系统精确的数学模型，机器人沿纵垂方向（Z轴）的运动控制过程较为稳定（见图11），控制精度受PID算法影响不大。收放缆前馈调节的最大相对误差不超过15%；脐带缆张力值随机器人运动状态在一定范围内振荡变化，机器人在下潜运动过程中，脐带缆张力大于上浮运动。

2）螺旋桨实时推力与其转速成正比，转速越大，螺旋桨发出的推进力越大，机器人实时运动线速度与各个螺旋桨的转速变化一一对应。

3）PID控制律调节螺旋桨转速具有一定的延迟性与超调效应，这种影响使脐带缆前馈调节的实际执行时间比指定控制调节时间延长了1 s左右；脐带缆收放过程也影响水下机器人沿X轴方向的运动，使其位移和速度出现了较大的超调现象。

4）水下机器人在轨迹跟踪仿真模拟过程中，主体水动力响应主要受机器人平动速度的变化引起的周围流场的改变、以及3个螺旋桨转速改变引起旋转流场变化二者共同影响。