水下滑翔机是一种新型的水下航行器，通过改变自身净浮力和质心位置来实现水下滑翔运动。与采用外部驱动装置的传统AUV相比，水下滑翔机具有成本低、能耗低、续航力强等特点，更适合大范围、长航程的工作任务，在海洋资源探测和海洋活动监测等领域具有广阔的应用前景^[1-2]。为降低内部机构调节的能耗，在大多数工作时间内，水下滑翔机处于稳态运动，因此，稳态运动分析通常可为其运动性能优化和控制系统设计提供指导。基于垂直平面内的直线稳态运动方程，王延辉等^[3]分析了滑翔机运动性能与各参数间的关系。同样利用直线稳态运动方程，Leonard等^[4]建立了水下滑翔机的反馈控制方法。将滑翔机的转向速度作为摄动参数，Mahmoudian等^[5]推导了滑翔机稳态转向时转向速度的近似解析表达式，Yang等^[6]推广此方法分析了装有尾舵的滑翔机转向性能。Zhang等^[7-8]建立了一种螺旋稳态运动的快速求解方法，并与Seawing滑翔机的试验结果进行对比验证。赵宝强等^[9-10]对水下滑翔机定常运动进行建模，并采用李雅普诺夫方法进行了稳定性分析。Liu等^[11]通过理论分析和样机试验发现，水下滑翔机自身的水动力特点可以使其产生一种新型的螺旋稳态运动，文中称为反向螺旋运动。

可以看出，水下滑翔机的稳态运动形式主要有直线稳态运动和螺旋稳态运动。相比于直线稳态运动，螺旋稳态运动的力学机理更为复杂。对于螺旋稳态运动，尚未看到其运动速度变化规律的系统分析，而且，现有多数研究将两种稳态运动形式单独讨论，未分析二者间的联系。本文建立水下滑翔机的稳态动力学模型，并给出能够反映两种稳态运动特征的运动速度描述方法。通过系统地改变控制参数，建立了稳态运动速度变化规律的研究方法，并以沈阳自动化研究所研制的Seawing水下滑翔机^{[8, 12]}为例开展研究。

1 坐标系与滑翔机的系统组成

如图 1所示，引入惯性坐标系E-ξηζ和体坐标系o-xyz来描述滑翔机的空间运动。其中，坐标轴E-ζ与重力方向一致，坐标系o-xyz的原点位于滑翔机的浮心，o-x轴沿着滑翔机的机身轴线指向首部方向，o-y轴指向右侧机翼的方向。水下滑翔机可以视为由固定结构、移动滑块、净浮力质量组成的多体系统。其中，固定结构包括机身及与机身固结在一起的所有结构，质量表示为m_rb；移动滑块为可实现平移和旋转运动的内部结构，质量表示为m_p；净浮力质量为浮力调节系统引起的可变化质量，表示为m_b。本文将滑翔机系统视为质量可变，浮力不变的结构，当净浮力质量为零时，滑翔机系统的总质量恰好等于其所受的浮力。向量r_rb、r_p、r_b分别对应m_rb、m_p、m_b的质心在体坐标系中的位置。r_rx为m_p的质心沿着o-x轴的位置，定义m_p沿着o-x轴的转角为γ。其中，γ=0对应m_p的质心处在体坐标的o-xz平面内，此时滑翔机将在垂直平面内做直线运动；γ≠0将导致滑翔机所受的升力偏离垂直平面，升力在水平面中的分量将作为向心力使其产生转向运动。

2 稳态运动的动力学模型

将水下滑翔机在惯性坐标系下的位置和姿态分别定义为b=[ξ_B  η_B  ζ_B]^T和Ω=[φ  θ  ψ]^T，并将其在体坐标系下的线速度和角速度分别定义为υ=[u  v  w]^T和ω=[p  q  r]^T，则水下滑翔机的运动学方程为^[13]

$ \begin{array}{l} \mathit{\boldsymbol{\dot b}} = {\mathit{\boldsymbol{R}}_{{\rm{EB}}}}\mathit{\boldsymbol{\upsilon }}, \\ \mathit{\boldsymbol{ \boldsymbol{\dot \varOmega} }} = {\mathit{\boldsymbol{R}}_{\Omega {\rm{B}}}}\mathit{\boldsymbol{\omega }} \end{array} $

(1)

其中，R_EB和R_BΩ可用水下滑翔机的姿态角来表示。采用简化形式c·=cos(·)，s·=sin(·)，t·=tan(·)表示三角函数，则有

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{R}}_{{\rm{EB}}}} = }\\ {\left( {\begin{array}{*{20}{c}} {{\rm{c}}\theta {\rm{c}}\psi }&{{\rm{s}}\varphi {\rm{s}}\theta {\rm{c}}\psi - {\rm{c}}\varphi {\rm{s}}\psi }&{{\rm{c}}\varphi {\rm{s}}\theta {\rm{c}}\psi + {\rm{s}}\varphi {\rm{s}}\psi }\\ {{\rm{c}}\theta {\rm{s}}\psi }&{{\rm{c}}\varphi {\rm{c}}\psi + {\rm{s}}\varphi {\rm{s}}\theta {\rm{s}}\psi }&{ - {\rm{s}}\varphi {\rm{c}}\psi + {\rm{c}}\varphi {\rm{s}}\theta {\rm{s}}\psi }\\ { - {\rm{s}}\theta }&{{\rm{s}}\varphi {\rm{c}}\theta }&{{\rm{c}}\varphi {\rm{c}}\theta } \end{array}} \right)}\\ {{\mathit{\boldsymbol{R}}_{{\rm{B\Omega }}}} = \left( {\begin{array}{*{20}{c}} 1&0&{ - {\rm{s}}\theta }\\ 0&{{\rm{c}}\varphi }&{{\rm{s}}\varphi {\rm{c}}\theta }\\ 0&{ - {\rm{s}}\varphi }&{{\rm{c}}\varphi {\rm{c}}\theta } \end{array}} \right)} \end{array} $

(2)

利用动量定理和坐标变换，水下滑翔机在体坐标系下的线动量P和角动量Π变化率为

$ \mathit{\boldsymbol{\dot P}} = {\rm{ }}\mathit{\boldsymbol{P\omega }} + \mathit{\boldsymbol{R}}_{{\rm{EB}}}^{\rm{T}}\left( {\mathit{\boldsymbol{G}} + \mathit{\boldsymbol{B}}} \right) + {\mathit{\boldsymbol{F}}_{\rm{B}}} $

(3)

$ \mathit{\boldsymbol{ \boldsymbol{\varPi} }} = \mathit{\boldsymbol{ \boldsymbol{\varPi} \omega }} + \mathit{\boldsymbol{P\upsilon }} + {\mathit{\boldsymbol{r}}_{{\rm{CG}}}} \times \mathit{\boldsymbol{R}}_{{\rm{EB}}}^{\rm{T}}\mathit{\boldsymbol{G}} + {\mathit{\boldsymbol{M}}_{\rm{B}}} $

(4)

其中，G和B分别为水下滑翔机的重力和浮力，r_CG为重心在体坐标系下的位置向量，F_B和M_B分别为体坐标系下水下滑翔机所受的粘性水动力和粘性水动力矩。对于F_B和M_B，可结合滑翔机运动速度和水动力系数来计算^{[8, 13]}。

通过分析水下滑翔机及附连水的总动能，可将线动量P和角动量Π表达为

$ \left( \begin{array}{l} \mathit{\boldsymbol{P}}\\ \mathit{\boldsymbol{ \boldsymbol{\varPi} }} \end{array} \right) = \mathit{\boldsymbol{M}}\left( \begin{array}{l} \mathit{\boldsymbol{\upsilon }}\\ \mathit{\boldsymbol{\omega }} \end{array} \right) = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{M}}_t}}&{{\mathit{\boldsymbol{C}}_t}}\\ {\mathit{\boldsymbol{C}}_t^T}&{{\mathit{\boldsymbol{I}}_t}} \end{array}} \right]\left( \begin{array}{l} \mathit{\boldsymbol{\upsilon }}\\ \mathit{\boldsymbol{\omega }} \end{array} \right) $

(5)

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{M}}_t} = {\mathit{\boldsymbol{M}}_A} + ({m_{{\rm{rb}}}} + {m_{\rm{p}}} + {m_{\rm{b}}})\mathit{\boldsymbol{I}}}\\ {{\mathit{\boldsymbol{C}}_t} = {C_A} - {m_{{\rm{rb}}}}{{\mathit{\boldsymbol{\hat r}}}_{{\rm{rb}}}} - {m_{\rm{p}}}{{\mathit{\boldsymbol{\hat r}}}_{\rm{p}}} - {m_{\rm{b}}}{{\mathit{\boldsymbol{\hat r}}}_{\rm{b}}}}\\ {{I_t} = {J_A} + {J_{{\rm{rb}}}} - {m_{\rm{p}}}{{\mathit{\boldsymbol{\hat r}}}_{\rm{p}}}{{\mathit{\boldsymbol{\hat r}}}_{\rm{p}}} - {m_{\rm{b}}}{{\mathit{\boldsymbol{\hat r}}}_{\rm{b}}}{{\mathit{\boldsymbol{\hat r}}}_{\rm{b}}}} \end{array} $

(6)

式中：M为整个滑翔机及附连水系统的广义质量矩阵，M_A、J_A和C_A分别为附连水在体坐标系下的质量矩阵，转动惯量矩阵和耦合项矩阵，J_rb为固定质量在体坐标系的转动惯量矩阵。

当水下滑翔机稳态运动时，广义质量矩阵M保持不变，对式(5)两端求导并结合式(3)和(4)，可有：

$ \left( \begin{array}{l} {\mathit{\boldsymbol{\dot \upsilon }}}\\ {\mathit{\boldsymbol{\dot \omega }}} \end{array} \right) = {\mathit{\boldsymbol{M}}^{ - 1}}\left( {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{P}} \times \mathit{\boldsymbol{\omega }} + \mathit{\boldsymbol{R}}_{{\rm{EB}}}^{\rm{T}}\left( {\mathit{\boldsymbol{G}} + \mathit{\boldsymbol{B}}} \right) + {\rm{ }}{F_{\rm{B}}}}\\ {\mathit{\boldsymbol{ \boldsymbol{\varPi} }} \times \mathit{\boldsymbol{\omega }} + \mathit{\boldsymbol{P}} \times \mathit{\boldsymbol{\upsilon }} + {\rm{ }}{\mathit{\boldsymbol{r}}_{{\rm{CG}}}} \times \mathit{\boldsymbol{R}}_{{\rm{EB}}}^{\rm{T}}\mathit{\boldsymbol{G}} + {\mathit{\boldsymbol{M}}_{\rm{B}}}} \end{array}} \right) $

(7)

基于式(1)和(7)，给定b，Ω，υ，ω的初值后，即可实现水下滑翔机的稳态运动仿真。

为进一步定量地验证本文动力学模型的正确性，将本文回转半径的分析结果与文献[8]中的结果进行对比，如图 2所示，其中r_rx=0.421 6 m, m_b=0.3 kg。可以看出，两种方法所得结果的变化趋势是一致的，但量值上存在较小的差别，原因为文献[8]中采用递归方法实现稳态螺旋运动的快速分析，在推导过程中忽略了质量矩阵中影响较小的耦合项，而本文方法从滑翔机的运动微分方程入手，更为详细地计入了质量矩阵。总体分析，本文的动力学模型可用于水下滑翔机稳态运动速度的分析。

3 稳态运动速度的描述

当水下滑翔机达到稳态运动时，在体坐标系中的速度分量将保持不变，即${\mathit{\boldsymbol{\dot \upsilon }}}$=0，${\mathit{\boldsymbol{\dot \omega }}}$=0。而体坐标中运动速度与滑翔机空间运动方向的对应关系复杂，不便于理解滑翔机空间运动的物理含义。尽管惯性坐标系中的运动速度描述与滑翔机的空间运动方向相对应，但当滑翔机处于螺旋稳态运动时，沿E-ξ和E-η轴的线速度分量并不会随时间达到稳定。为此，本文从螺旋稳态运动入手，采用线速度v_t、v_ζ和角速度ω_ζ来描述水下滑翔机的稳态运动速度。如图 3所示，Seawing水下滑翔机螺旋运动轨迹，水下滑翔机螺旋运动在E-ξη平面的投影轨迹应为匀速圆周运动，对应轨迹上的某一位置，定义v_t为圆周上的切向线速度，其与E-ξη平面相平行，v_ζ为沿E-ζ方向的垂向线速度，ω_ζ为滑翔机绕中心轴的转向角速度。

将b，Ω，υ，ω的初值取为零向量，并将控制参数取值为γ=90°，r_rx=0.42 m, m_b=0.3 kg，得到滑翔机各速度分量随时间的变化规律，如图 4所示。可以看出沿E-ξ和E-η轴的速度分量v_ξ和v_η处于周期性振荡，而速度分量v_t、v_ζ、ω_ζ将逐渐达到稳态。对滑翔机在垂直平面内的直线稳态运动而言，可将其理解为ω_ζ=0的螺旋稳态运动，此时v_t等价于垂直平面内的水平方向运动速度，v_ζ等价于垂直平面内的垂直方向运动速度。因此，v_t、v_ζ和ω_ζ可实现两种稳态运动中运动速度分量的统一描述，便于后续进行统一分析。

4 稳态运动速度的分析

利用已给出的稳态速度统一描述方法和稳态运动仿真技术，分析水下滑翔机控制输入与运动速度间的变化规律。水下滑翔机的姿态控制输入量包括：移动滑块质心沿o-x轴的位置r_rx、移动滑块绕o-x轴的转角γ和净浮力质量m_b。

根据滑翔机重力与浮力的关系，可将滑翔机的稳态运动分类为上浮稳态运动和下潜稳态运动。当重力大于浮力，即m_b＞0时，滑翔机将进行下潜稳态运动，为保证滑翔机沿着机头的方向滑行，需调节r_rx使滑翔机的重心在浮心的前方。而当重力小于浮力，即m_b＜0时，滑翔机将进行上浮稳态运动，此时需调节r_rx使滑翔机的重心在浮心的后方。以Seawing水下滑翔机为例进行研究，通过系统地改变典型工况中的单一控制输入量，设计一系列稳态运动计算工况，如表 1所示，进而分析滑翔机稳态运动速度分量随控制参数的变化规律。其中，γ=0°对应垂直面内的直线运动；γ=90°对应螺旋运动中回转速度最大的滑块转角；在垂直面内运动时，r_rx=0.42 m和r_rx=0.38 m分别对应下潜和上浮中水平速度的最大值。

通过对组合1的分析，可得到各运动速度分量随转角γ的变化规律，如图 5所示。由于滑翔机外形和质量系统的左右对称性，导致切向速度v_t和垂向速度v_ζ关于γ=0°位置左右对称；转向速度ω_ζ关于(0, 0)点中心对称。随着转角γ的绝对值的增大，转向速度ω_ζ增大，垂向速度v_ζ的绝对值也增大，切向速度v_t先增大后减小。这是由于转角γ的增大将导致滑翔机升力在水平方向的向心力分量增大，在垂直方向分量减小。向心力增大导致滑翔机的转速增大。而升力在垂直方向分量的减小，将导致垂直方向的受力平衡需要更多的阻力，因而垂向速度v_ζ的绝对值增大。

通过对组合2的分析，可得到各运动速度分量随r_rx的变化规律，如图 6所示。可以看出，在垂直面内运动时(γ=0°)，切向速度v_t，垂向速度v_ζ的变化规律与其在螺旋运动时(γ=90°)是一致的；v_t与v_ζ在下潜阶段的变化规律与其在上浮阶段也是一致的；在下潜和上浮时v_t均存在最大值；当r_rx向r_rx=0.35 m和r_rx=0.45 m的两端变化时，v_ζ的绝对值逐渐变大并趋于平缓。转向速度ω_ζ在下潜和上浮阶段的变化规律存在明显的不同，在下潜阶段ω_ζ存在最小值，对应的r_rx恰好在v_t最大值对应的r_rx附近；在上浮阶段，当r_rx向r_rx=0.35 m变化时，ω_ζ逐渐增大并趋于平缓。当水下滑翔机在垂直平面内滑行时，r_rx通常取值在v_t最大值对应的位置，若此时需要调整γ进行转向，结合图 5可以看出，可以考虑在上浮阶段调整γ来获得更大的转向速度。

通过对组合3的分析，可得到各运动速度分量随m_b的变化规律，如图 7所示。可以看出，在多数情况下，速度分量v_t与v_ζ的绝对值随m_b绝对值的增大而增大；但在r_rx=0.38 m, γ=90°的螺旋上浮工况中，v_ζ存在一个转折点，在m_b＜-0.35 kg后，v_ζ的绝对值反而随着m_b绝对值的增大而减小，相应的切向速度v_t绝对值增大的幅度加快；转向速度ω_ζ在下潜和上浮阶段的变化规律存在明显不同；在下潜阶段ω_ζ随m_b绝对值的增大而减小；在上浮阶段ω_ζ随m_b绝对值的增大而增大。

5 结论

1) 采用切向线速度v_t，垂向线速度v_ζ和转向角速度ω_ζ的运动速度描述方法，可实现两种稳态运动形式中运动速度的统一描述。

2) 水下滑翔机的转向角速度ω_ζ在上浮和下潜运动中的变化趋势差别很大，而线速度(v_t、v_ζ)在上浮和下潜运动中的变化趋势相近。

3) 运动速度变化规律可为水下滑翔机的运动控制提供指导。如，当需要螺旋转向时，可以考虑增大转角γ来提高转向速度，同时，在相同转角下，也可考虑在上浮阶段进行转向来提高转向速度。

本文的稳态运动速度分析是以Seawing水下滑翔机为例开展的。而对于不同机构设置的水下滑翔机，本文的研究思路也同样适用。