0 引　言

水下无人航行器（UUV）是海洋资源开发、海洋环境监测及海洋生态保护等的关键装备之一。近年来，随着人类对深海资源开发的不断深入，UUV由于在军事与科研方面有重要应用，引起了广泛关注，研究价值和意义也日益凸显^[1-2]。

操纵性能是UUV总体性能的主要指标之一，直接影响着航行器执行作业任务的能力，在航行器初步设计和运动控制参数选取时，操纵性能的预报尤为重要^[3]。建立合理的空间运动数学模型是研究UUV操纵性及控制系统设计的基础，王波^[4]、赵金鑫^[5]等参考潜艇操纵性研究的方法搭建了非线性的UUV空间运动数学模型进行了仿真预报，段斐等^[6]针对REMUS模型中的推力（力矩）和舵力（力矩）难于获取等问题，提出了一种修正的REMUS模型并完成了航行器的运动仿真预报。戴君锐等^[7]采用六自由度模型完成了UUV操纵运动仿真，并与K-T模型的仿真结果进行了对比。徐得志等^[8]参考舰船操纵性理论研究方法，采用四阶龙格-库塔算法对UUV垂向操纵运动进行了研究。聂为彪等^[9]建立了UUV水平面内的动力学方程，应用Matlab编程对研究对象的水平面操纵运动进行了预报。

本文基于水下航行力学基本原理，参考鱼雷操纵性基本理论^[10]，在小冲角、小侧滑角、小机动运动条件下导出了UUV六自由度空间运动的线性数学模型，并搭建了操纵运动仿真预报平台，对某型UUV单平面典型操纵运动特性进行仿真预报，为UUV水动力布局和控制系统的设计提供一定的技术指导与理论支撑。

1 UUV空间运动数学模型 1.1 基本假设

1）UUV是刚体，并关于 $ xoy $ 平面对称；

2）流体是无粘不可压缩的；

3）流体为无界流场，在航行体运动之前是静止的；

4）坐标原点与UUV浮心重合；

5）不考虑地球的自转和地球的曲率，近似认为地面坐标系为惯性坐标系。

1.2 坐标系选择与运动学参数

采用如图1所示坐标系，包括地面坐标系 $ E{X_e}{Y_e}{Z_e} $ 和航行体坐标系 $ Oxyz $ 。地面坐标系与大地相连，航行体坐标系与UUV相连，坐标原点与UUV浮心重合。

地面坐标系 $ E{X_e}{Y_e}{Z_e} $ 中的坐标（ $ {X_e},{Y_e},{Z_e} $ ）和航行体坐标系 $ Oxyz $ 中的坐标（ $ x,y,z $ ）满足式（1）所示，其中 ${\boldsymbol{C}}_{\mathbf{E}}^{\mathbf{O}}$ 为航行体坐标系到地面坐标系的转换矩阵。

$ \left[ \begin{gathered} {X_e} \\ {Y_e} \\ {Z_e} \\ \end{gathered} \right] = {\boldsymbol{C}}_{\mathbf{E}}^{\mathbf{O}} \left[ \begin{gathered} x \\ y \\ z \\ \end{gathered} \right]，$

(1)

$ \begin{split}{\boldsymbol{C}}_{\mathbf{E}}^{\mathbf{O}} =& \left[ \begin{array}{cc} \cos \theta \cos \psi & \sin \psi \sin \varphi - \sin \theta \cos \psi \cos \varphi \\ \sin \theta &\cos \varphi \cos \theta \\ - \cos \theta \sin \psi &\cos \psi \sin \varphi + \sin \theta \sin \psi \cos \varphi \end{array} \right. \\ & \left. \begin{array}{c} \sin \psi \cos \varphi + \sin \theta \cos \psi \sin \varphi \\ - \sin \varphi \cos \theta \\ \cos \psi \cos \varphi - \sin \theta \sin \psi \sin \varphi \end{array} \right] 。\end{split}$

(2)

运动参数包括：

1）UUV在地面坐标系中的位置（ $ {X_e},{Y_e},{Z_e} $ ）和姿态角（ $ \varphi $ ， $ \psi $ ， $ \theta $ ）；

2）航行体坐标系下的UUV速度（ $ {v_x},{v_y},{v_z} $ ）和角速度（ $ {\omega _x} $ , $ {\omega _y} $ , $ {\omega _z} $ ）；

3）UUV速度 $ v $ 和速度的流体动力角 $ \alpha $ ， $\ \beta $ 关系如下式：

$ \left\{ \begin{gathered} v = \sqrt {{v_x}^2 + {v_y}^2 + {v_z}^2} ，\\ \alpha = \arctan ( - {v_y}/{v_x})，\\ \beta = \arctan ({v_z}/\sqrt {{v_x}^2 + } {v_y}^2)。\\ \end{gathered} \right. $

(3)

1.3 运动学方程

根据UUV速度在地面坐标系和航行体坐标系下的转换关系，可得平动运动学方程：

$ {\left[ {{X_e}\;{Y_e}\;{Z_e}} \right]^{{\rm{T}}} } = {\boldsymbol{C}}_{\mathbf{E}}^{\mathbf{O}} {\left[ {{v_x}\;{v_y}\;{v_z}} \right]^{{\rm{T}}} }。$

(4)

式中： $ {\boldsymbol{C}}_{\mathbf{E}}^{\mathbf{O}} $ 为航行体坐标系到地面坐标系的转换矩阵，如式（2）所示。

由UUV旋转角在地面坐标系和航行体坐标系下的关系，可得转动运动学方程：

$ \left[ \begin{gathered} {\dot \varphi } \\ {\dot \psi } \\ {\dot \theta } \\ \end{gathered} \right] = \left[ \begin{gathered} 1\;\;\;\; - \cos \varphi \tan \theta \;\;\;\;\;\;\sin \varphi \tan \theta \\ 0\;\;\;\;\;\cos \varphi /\cos \theta \;\;\; - \sin \varphi /\cos \theta \\ 0\;\;\;\;\;\;\;\;\;\;\;\sin \varphi \;\;\;\;\;\;\;\;\;\;\;\cos \varphi \\ \end{gathered} \right] \left[ \begin{gathered} {\omega _x} \\ {\omega _y} \\ {\omega _z} \\ \end{gathered} \right]。$

(5)

1.4 动力学方程

根据动量定理和动量矩定理^[11]，可以得到UUV在无流界中的动力学运动方程，并通过系列的推导，可表示为矩阵形式：

$ \begin{split}\left[ \begin{array}{*{20}{c}} {{\dot v}_x} \\ {{\dot v}_y} \\ {{\dot v}_z} \\ {{\dot \omega }_x} \\ {{\dot \omega }_y} \\ {{\dot \omega }_z} \\ \end{array} \right] = {\left[ \begin{array}{*{20}{c}} m+{\lambda _{11}} & 0 & 0 & 0 & m{z_c} & - m{y_c} \\ 0 & m + {\lambda _{22}} & 0 & - m{z_c} & 0 & m{x_c} + {\lambda _{26}} \\ 0 & 0 & m + {\lambda _{33}} & m{y_c} & {\lambda _{35}} - m{x_c} & 0 \\ 0 & - m{z_c} & m{y_c} & {J_{xx}} + {\lambda _{44}} & 0 & 0 \\ m{z_c} & 0 & {\lambda _{35}} - m{x_c} & 0 & {J_{yy}} + {\lambda _{55}} & 0 \\ - m{y_c} & m{x_c} + {\lambda _{26}} & 0 & 0 & 0 & {J_{zz}} + {\lambda _{66}} \\ \end{array} \right]^{ - 1}}\times \left[ \begin{array}{c} {F_{sx}} + {F_{gx}} + {F_{tx}} + T \\ {F_{sy}} + {F_{gy}} + {F_{ty}} \\ {F_{sz}} + {F_{gz}} + {F_{tz}} \\ {M_{sx}} + {M_{gx}} + {M_{tx}} + \Delta {M_{xp}} \\ {M_{sy}} + {M_{gy}} + {M_{ty}} \\ {M_{sz}} + {M_{gz}} + {M_{tz}} \\ \end{array} \right]。\end{split} $

(6)

式中：

1） $ {\left[ {T\;\;\Delta {M_{xp}}} \right]^{\text{T}}} $ 为螺旋桨的推力与扭矩，T为螺旋桨在航行体坐标系中的有效推力； $ \Delta {M_{xp}} $ 为单桨或不平衡对转桨产生的扭矩。

2） $ {\left[ {{F_{gx}}\;\;{F_{gy}}\;\;{F_{gz}}\;\;{M_{gx}}\;\;{M_{gy}}\;{M_{gz}}} \right]^{{\rm{T}}} } $ 为由静水力引起的力与力矩，该静水力与力矩主要由重力G和浮力B产生，如下式：

$ \left[ \begin{gathered} {F_{gx}} \\ {F_{gy}} \\ {F_{gz}} \\ {M_{gx}} \\ {M_{gy}} \\ {M_{gz}} \\ \end{gathered} \right] = \left[ \begin{gathered} {\text{(}}B - G)sin\theta \\ {\text{(}}B - G)\cos \theta \cos \varphi \\ - {\text{(}}B - G)\cos \theta \sin \varphi \\ G\cos \theta ({y_c}\sin \varphi + {z_c}\cos \varphi ) \\ - G({x_c}\cos \theta \sin \varphi + {z_c}\sin \theta ) \\ G({y_c}\sin \theta - {x_c}\cos \theta \cos \varphi ) \\ \end{gathered} \right]。$

(7)

式中： $ {x_c},{y_c},{z_c} $ 为UUV重心在航行体坐标系下的位置。

3） $ {\left[ {{F_{tx}}\;\;{F_{ty}}\;\;{F_{tz}}\;\;{M_{tx}}\;\;{M_{ty}}\;{M_{tz}}} \right]^{{\rm{T}}} } $ 表示除其他力以外的流体动力，如下式：

$ \begin{split}&\left[ \begin{gathered} {F_{tx}}\;\; \\ {F_{ty}}\;\; \\ {F_{tz}}\;\; \\ {M_{tx}} \\ {M_{ty}}\; \\ {M_{tz}} \\ \end{gathered} \right] = \\&\left[ \begin{gathered} - m\left[ {{v_z}{\omega _y} - {v_y}{\omega _z} + {y_c}{\omega _x}{\omega _y} + {z_c}{\omega _x}{\omega _z} - {x_c}({\omega _y}^2 + {\omega _z}^2)} \right]\;\; \\ - m\left[ {{v_x}{\omega _z} - {v_z}{\omega _x} + {x_c}{\omega _x}{\omega _y} + {z_c}{\omega _y}{\omega _z} - {y_c}({\omega _x}^2 + {\omega _z}^2)} \right]\; \\ - m\left[ {{v_y}{\omega _x} - {v_x}{\omega _y} + {x_c}{\omega _x}{\omega _z} + {y_c}{\omega _y}{\omega _z} - {z_c}({\omega _x}^2 + {\omega _y}^2)} \right]\;\; \\ - \left[ {m{y_c}({v_y}{\omega _x} - {v_x}{\omega _y}) + m{z_c}({v_z}{\omega _x} - {v_x}{\omega _z}) + ({J_{zz}} - {J_{yy}}){\omega _y}{\omega _z}} \right] \\ - \left[ {m{z_c}({v_z}{\omega _y} - {v_y}{\omega _z}) + m{x_c}({v_x}{\omega _y} - {v_y}{\omega _x}) + ({J_{xx}} - {J_{zz}}){\omega _z}{\omega _x}} \right]\; \\ - \left[ {m{x_c}({v_x}{\omega _z} - {v_z}{\omega _x}) + m{y_c}({v_y}{\omega _z} - {v_z}{\omega _y}) + ({J_{yy}} - {J_{xx}}){\omega _x}{\omega _y}} \right] \\ \end{gathered} \right]。\end{split}$

(8)

4） $ {\left[ {{F_{sx}}\;\;{F_{sy}}\;\;{F_{sz}}\;\;{M_{sx}}\;\;{M_{sy}}\;{M_{sz}}} \right]^{{\rm{T}}} } $ 表示黏性流体动力，UUV在小冲角、小侧滑角、小机动运动条件下^[10]，可用线性导数来表示，如下式：

$ \left[ \begin{gathered} {F_{sx}}\; \\ {F_{sy}}\; \\ {F_{sz}}\; \\ {M_{sx}}\; \\ {M_{sy}}\; \\ {M_{sz}} \\ \end{gathered} \right] = \frac{1}{2}\rho {v^2}S\left[ \begin{gathered} {C_x}(0) \\ {C_y}(0) + C_y^\alpha \alpha + C_y^{{{\bar \omega }_z}}{{\bar \omega }_z} + C_y^{{\delta _e}}{\delta _e} \\ C_z^\beta \beta + C_z^{{{\bar \omega }_x}}{{\bar \omega }_x} + C_z^{{{\bar \omega }_y}}{{\bar \omega }_y} + C_z^{{\delta _r}}{\delta _r}\; \\ L\left[ {m_x^\beta \beta + m_x^{{{\bar \omega }_x}}{{\bar \omega }_x} + m_x^{{{\bar \omega }_y}}{{\bar \omega }_y} + m_x^{{\delta _r}}{\delta _r}\;} \right]\; \\ L\;\left[ {m_y^\beta \beta + m_y^{{{\bar \omega }_x}}{{\bar \omega }_x} + m_y^{{{\bar \omega }_y}}{{\bar \omega }_y} + m_y^{{\delta _r}}{\delta _r}} \right] \\ L\left[ {{m_z}(0) + m_z^\alpha \alpha + m_z^{{{\bar \omega }_z}}{{\bar \omega }_z} + m_z^{{\delta _e}}{\delta _e}} \right] \\ \end{gathered} \right]。$

(9)

式中： ${\bar \omega _x} = \dfrac{{{\omega _x}L}}{v}；{\bar \omega _y} = \dfrac{{{\omega _y}L}}{v}；{\bar \omega _z} = \dfrac{{{\omega _z}L}}{v}$ ；L为UUV长度； $ \rho $ 为流体密度； $ S $ 为UUV横截面积。

2 UUV操纵运动仿真平台搭建

搭建操纵性仿真预报平台^[12-13]，该平台将航行体空间运动数学模型分为姿态角模型、速度转化模型、粘性流体动力模型、静水力模型、其他力模型、坐标系转换模型等。

仿真预报过程中采用变步长四阶-五阶Runge-Kutta算法，输入UUV基本参数及水动力参数，并赋予初始速度、推力和舵角，运行程序后便可输出对应时刻UUV的姿态角、位置、速度等，生成对应曲线图，即可对UUV的操纵性能进行预报。

3 UUV操纵运动仿真预报 3.1 水平面运动仿真 3.1.1 回转运动预报

给定推力值T＝95.96 N，对应设计航速4 kn，初始速度 $ v = 2\;{\rm{m}}/{\rm{s}} $ ，并分别操垂直舵 $ {\delta _r} = $ 10°，15°，20°，其浮心的运动轨迹详如图2所示，结果如表2所示。

可以看出，在设计航速为4 kn下，舵角越大，回转直径、回转周期和纵距越小。在操大舵角20°时，回转直径为56.3 m，约为航行体长（7 m）的8倍，可以看出该UUV具有较高的水平面回转性能。

3.1.2 水平面Z形操舵预报

给定推力值T，初始速度，对应航速分别为4 kn和8 kn下，对于 $ {\delta _r}/\psi = {15^ \circ }/{15^ \circ } $ 进行Z型操舵运动，反复操舵4次，其偏航角 $ \psi $ 和舵角 $ {\delta _r} $ 随时间的变化轨迹如图3所示。由此得到UUV的初转期 $ {t_a} $ 、超越时间 $ {t_{0v}} $ 、超越偏航角 $ {\psi _{0v}} $ 、全周期 $ {T_z} $ 等特征参数，便可评估UUV水平面航向改变性能。

从表3可以看出，航速8 kn时，UUV的初转期 $ {t_a} $ 、超越时间 $ {t_{0v}} $ 、超越偏航角 $ {\psi _{0v}} $ 、全周期 $ {T_z} $ 均比低航速4 kn时小。此外，还可以看出，该UUV具备良好的应舵性和操纵性，响应时间快。

3.2 垂直面运动仿真 3.2.1 垂直面T形操舵预报

给定推力值T，初始速度，对应航速分别为4 kn和8 kn下，操升降舵 $ {\delta _e} = {15^ \circ } $ ，UUV的俯仰角 $ \theta $ 和深度 $ {Y_e} $ 等都在改变，当俯仰角达到设定俯仰角10°时，立即回舵到初始状态，当俯仰角速度为0时， $ \theta $ 达到稳定值，UUV稳定在一个新的航向上。升降舵 $ {\delta _e} $ 、俯仰角 $ \theta $ 、航行深度 $ {Y_e} $ 随时间的变化轨迹如图4所示。

从表4两个响应航速结果可以看出，航速越大，UUV的初转期 $ {t_a} $ 越小，说明UUV的应舵更快，转首性好，下潜快，但同时超越俯仰角 $ {\theta _{0v}} $ 和越深度 $ {\xi _{0v}} $ 也越大。此外，从图4可以看出，航速越大，航行器稳定在新的航向上的时间也越长。

3.2.2 垂直面纯下潜预报

给定推力值T，初始速度，在设计航速4 kn下，分别操横舵 $ {\delta _e} $ ＝15°，10°进行下潜，其俯仰角 $ \theta $ 和深度 $ {Y_e} $ 随时间的变化轨迹如图5所示。

从图5可以看出，在同一航速为4 kn下，操横舵角越大，UUV转首性越好，下潜得也越快，但对应稳定后的俯仰角也越大，符合UUV的实际操纵运动特性。

4 结　语

本文基于水下航行力学基本原理，参考鱼雷操纵性基本理论，导出了UUV六自由度空间运动线性数学模型，并搭建了操纵运动仿真预报平台，采用变步长四阶-五阶龙格-库塔（Runge-Kutta）算法，对某型UUV的水平面回转、Z形操舵运动以及垂直面T形操舵、纯下潜等运动进行了仿真预报。仿真结果表明，该UUV具备良好的应舵性和操纵性能，响应时间较快，能真实地反映研究对象的操纵运动特性，其结果可为UUV水动力布局和控制系统的设计提供一定的技术指导与理论支撑。