无人水翼航行器运动的非线性建模
无人水翼航行器运动的非线性建模
无人水翼航行器的结构和运动特点
为了提高航行器快速性，应尽量提高推进器的推力和减小航行器航行的阻力。从 节能的角度考虑，减少航行器航行时的阻力是提高快速性的好办法。航行器阻力包括水阻力和空气阻力。由于空气阻力比航行器水阻力小得多，因此主要考虑航行器的水阻力，其中包括摩擦阻力、涡流阻力(形状阻力)和兴波阻力三个部分。涡流阻力和兴波阻力可以通过改善航行器的线型而得以改进。
水翼附着了一个可以通过控制系统调节的襟翼。在不同航速和海面状况的条 件下，通过调整襟翼可以实现对升力的控制。通过对水翼实施控制，从而产生与波浪干扰力相反的力，以抵消或减轻波浪的干扰所引起的水翼船垂荡和纵摇运动， 保证水翼船在风浪中平稳地航行。虽然各种水翼船实际控制系统的具体构成各有不同，但是大都具有与下图相类似的结构原理。
无人水翼航行器的非线性运动数学模型
高速水翼航行器的实际运动异常复杂，在一般情况下具有6个自由度，其中三个线性运动和三个转动运动 ，分别为：横荡 、 纵荡 、 垂荡和横摇 、纵摇 、艏 摇 。其中横荡 、纵荡 、垂荡分别为沿X，Y，Z轴的平直运动 ，横摇 、纵摇和艏摇分别为绕X，Y，Z轴的转动运动 ，运动名称定义见下表所示：
| 运动 | 坐标轴 | 力与力矩 | 速度/角速度 | 位置/姿态 |
|---|---|---|---|---|
| 横荡（进退） | 沿 x 轴方向运动 | X | u | x |
| 纵荡（横移） | 沿 y 轴方向运动 | Y | v | y |
| 升沉（垂荡） | 沿 z 轴方向运动 | Z | w | z |
| 横摇 | 旋转轴： x 轴 | K | p | $\varphi $ |
| 纵揺 | 旋转轴： y 轴 | M | q | $\theta $ |
| 艏摇 | 旋转轴： z 轴 | N | r | $ \psi$ |
（1） 附体坐标系
附体坐标系$O_{b} $ $X_{b} $$Y_{b} $$Z_{b} $ (或称为“随船运动坐标系”):为了描述航行器航行过程中的姿态，需要定义一个与航行器固连，即随着航行器一 起运动的附体坐标系 $O_{b} $ $X_{b} $$Y_{b} $$Z_{b} $ 。 它的原点$O_{b} $位于航行器重心，$O_{b} $ $X_{b} $轴在中纵剖面内平行于船体基线，指向船首； $O_{b} $ $Y_{b} $ 轴垂直于中纵剖面，指向右； $O_{b} $ $Z_{b} $ 轴在中纵剖面内垂直于船体基线， 指向龙骨。
（2） 空间坐标系
空间固定坐标系O’X’Y’Z’: 随航行器运动坐标系$O_{b} $ $X_{b} $$Y_{b} $$Z_{b} $是建立船舶质心运动模型的基本坐标系，不能反映航行器运动欧拉角的变化。为此需要建立一个固定在地面的固定坐标系O‘X’Y‘Z’, 其 中O‘Z’轴垂直向下，O’X ‘轴平行于水平面，坐标原点O ‘在t=0 时刻与附体坐标系的原点$O_{b} $ 重合。
$$
m(\ddot{\xi} + v\dot{\theta}) = -\sum_{i=1}^{2} (F_{fi} + F_{fpi}) \cos\theta + mg \cos\theta
$$
$$
J_{yy} \ddot{\theta} = Fd - M_f - \sum_{i=1}^{2} (F_{fi} + F_{fpi}) d_i \cos\theta + \sum_{i=1}^{2} (F_{fi} + F_{fpi}) d_i \sin\theta
$$
这个方程表示的是在$Y_{b} $ 轴上，航行器受到的力矩与它的角加速度之间的关系。$F_{fi} $ 为前后水翼( i = 1 为前翼，i = 2 为后翼) 所产生的力;$F_{fpi} $ 为襟翼所产生的力,后翼无襟翼，固定角度; $F_{Di} $ 为水翼所产生的阻力,水翼所产生的阻力及其力矩可以忽略; f 为连接杆在水中运动所产生的阻力; F 为电机推力。θ为纵摇角，ξ 为无人水翼航行器垂直于水平面升沉量。其中：
- $J_{yy}$是航行器绕$Y_{b} $ 轴的转动惯量;
- $\ddot{\theta }$ 是航行器的纵摇加速度;
- Fd是电机推力产生的力矩;
- $M_{f}$是连接杆在水中运动产生的阻力产生的力矩,为fd;
- $F_{fi}$是水翼产生的力;
- $F_{fpi}$是襟翼产生的力;
- $F_{Di}$是水翼产生的阻力，可以忽略;
- d和$d_{i}$是相应力臂的长度。
加入随机海浪，产生波动，$Z_{W}$ 和 $M_{W}$ 是随机海浪的力和力矩, $Z_{c}$ 和 $M_{\mathrm{c}}$ 是控制力和力矩。将部分参数展开表示后如下:
$$
m(\ddot{\xi} + V \dot{\theta}) = -\sum_{i=1}^{2} (F_{fi} + F_{fpi}) \cos\theta + mg \cos\theta + Z_w
$$
$$
J_{yy} \ddot{\theta} = Fd - M_f - \sum_{i=1}^{2} (F_{fi} + F_{fpi}) d_i \cos\theta + \sum_{i=1}^{2} (F_{fi} + F_{fpi}) d_i \sin\theta + M_w
$$
将航行器的非线性数学模型改写成如下形式：
$$
Z(\ddot{\xi}, \dot{\xi}, \xi, \ddot{\theta}, \dot{\theta}, \theta) = m(\ddot{\xi} + V \dot{\theta}) + \sum_{i=1}^{2} (F_{fi} + F_{fpi}) \cos\theta - mg \cos\theta - Z_w = 0
$$
$$
M(\ddot{\xi}, \dot{\xi}, \xi, \ddot{\theta}, \dot{\theta}, \theta) = J_{yy} \ddot{\theta} - M_w + \sum_{i=1}^{2} (F_{fi} + F_{fpi}) d_i \cos\theta - \sum_{i=1}^{2} (F_{fi} + F_{fpi}) d_i \sin\theta - Fd + M_f = 0
$$
线性化：
水翼所产生的作用力
在航行器的非线性数学模型中与水翼有关的作用力有两个，分别是水翼所产 生的升力和惯性力，即$F_{fi}$=$L_{fi}$+$F_{ai}$。
水翼升力的计算
水翼所产生的阻力及其力矩可以忽略，而仅考虑水翼所产生的升力$L_{fi}$, 故水翼的升力可表示为:
$$
L_{fi}=1/2\rho S_{i}V^{2} C_{Li}
$$
其中, $\rho$ 为流体的密度, $S_{i}$ 是水翼投影的面积, 可由 $S_{i}=l_{i} \cdot b_{i}$ 求得, $l_{i}$ 是水翼的翼展, $b_{i}$ 是水翼的翼弦。
水翼的升力系数 $C_{L i}$ 有时会受到水翼浸于流体中的深度的影响, 也称作 “浅浸效应”。 由于全浸式水翼艇的水翼接近液面时, 水和大气交界会使水翼的特性, 即升力系数的斜 率 $d C_{L i} / d \alpha_{i}$ 变小。
水翼升力系数的求法可以采用 $K . L$. Wadlin 提出的方法进行求解:
$$
C_{L i}=\frac{2 K_{2 i} \pi \lambda_{i} \alpha_{i}}{\lambda_{i}+2 K_{2 i}+1}+\frac{8}{3}\left(1-\frac{\lambda_{i}}{10}\right) \sin ^{2} \alpha_{i} \cos \alpha_{i}
$$
式中, $\lambda_{i}$ 是展弦比, $K_{2 i}$ 是二维深度修正因子, $\alpha_{i}$ 是水翼的攻角。二维深度修正因 子 $K_{2 i}$ 的表达式如下:
$$
K_{2 i}=\frac{\left(4 f_{i}\right)^{2}+8 f_{i} \sin \alpha_{i}+1}{\left(4 f_{i}\right)^{2}+8 f_{i} \sin \alpha_{i}+2}
$$
式中, $f_{i}$ 是水翼在 $1 / 4$ 翼弦处水翼的浸没深度, $f_{i}$ 根据水翼处于自由液面时具有的 临界条件可以用以下方式计算:
$$
f_{i}=\frac{d_{i}}{b_{i}}+\frac{\frac{d_{i}}{4} \sin \alpha_{i}}{0.05+\frac{d_{i}}{b_{i}}}
$$
式中, $d_{i}$ 为水翼的导边到自由液面的长度, 计算方法如下:
$$
d_{i}=d_{0 i}+\xi-\left(x_{f i}-x_{G}\right) \sin \theta
$$
攻角 $\alpha_{i}$ 的计算方法如下:
$$
\alpha_i = \alpha_{si} - \alpha_{0i} + \frac{\dot{\xi} - (x_{fi} - x_g) \dot{\theta} - \dot{\zeta}_i}{V} + \theta
$$
其中, $\alpha_{s i}$ 为水翼的安装角, $\alpha_{0 i}$ 是水翼的零升力攻角, $\dot{\zeta}_{i}$ 是波浪在水翼处的运动速度垂直于水面方向的分量。
在参考论文中，对于全浸式水翼艇的水翼升力系数由于计算过于繁杂, 本文根据经验将其定为常 数, 取值 0.33^[1]^，其襟翼计算的升力系数参考如下：
水翼的惯性力
航行器在水中运动时可以看作是一个刚体，受到的惯性力除了航行器本身的以外还有周围流体的作用。航行器运动时周围的流体质点就会受到扰动，就会将流体反作用力施加在航行器上。流体作用力大小正比于运动加速度，方向与运动加速度方向相反， 也称作惯性力，大小如下：
$$
{F_{a i}=-m_{f i}\left(\ddot{\xi}+V \dot{\theta}-d_{i} \ddot{\theta}-\ddot{\zeta}_{i}\right)}
$$
$m_{f i}$ 为水翼的附加质量，它是水翼的惯性阻力与对应的水翼垂向运动加速度的比例系数。由于水的密度较大，所以水翼的附加质量不可忽略。其中, ${m_{f i}=\frac{\pi}{2} \cdot \frac{b_{i}^{2}}{4} \cdot l_{i} \cdot \rho, l_{i}}$ 是水翼的翼展, $b_{i}$ 是水翼的翼弦, V为水翼航行器的航速, $d_{i}$ 是水翼到船体的重心的距离, $\ddot{\zeta}_{i}$ 是波浪在水翼处的运动加速度垂直于水面方向的分量。
水翼航行器转动惯量的估算
一个长方体关于其对角线的转动惯量可以通过以下公式计算：
$$
I = (M/12) * (h² + w² + l²)
$$
其中：
- I 是转动惯量
- M 是质量
- h、w、l 分别是高度、宽度和长度
假设我们有一个质量为1.5千克的长方体，其长度为1米，宽度和高度都是0.3米，我们可以将这些数值代入公式：
$$
I = (1.5/12) * (0.3² + 0.3² + 1²)
$$
$$
I = 0.125 kg*m²
$$
这就是这个长方体关于其对角线的转动惯量。注意，这是在假设质量均匀分布的情况下的转动惯量。如果质量分布不均匀，转动惯量会有所不同。
刚体转动惯量是惯性张量在特定轴线上的取值，而惯性张量则是一个包含了刚体在不同轴线上转动惯量的矩阵。即Iyy取0.14.
对于我们所设计的水翼航行器航，行器结构参数如下表所示。
潜行速度10km/h（2.8m/s）,滑行速度20km/h(5.6m/s),即V=5.6m/s。
| 名称 | 数值 |
|---|---|
| 航行器质量 m/kg | 1.5 |
| 重力加速度 g /( m·$s ^{-2} $) | 9.8 |
| 航行器转动惯量 Jyy /( kg·$m^{-2} $ ) | 0.14（估计） |
| 连接杆长 d/m | 0.09 |
| 连接杆1/前翼距航行器重心距离 d1 /m | 0.3（估计） |
| 连接杆2/后翼距航行器重心距离 d2 /m | 0.3（估计） |
| 电机距航行器重心距离 d2 /m | 0.3（估计） |
| 前翼展长 L1 /m | 0.135 |
| 前翼弦长 b1 /m | 0.0324 |
| 后翼展长 L2 /m | 0.1016 |
| 后翼弦长 b2/m | 0.0168 |
线性化后的方程
$$
Z_{\ddot{\xi}} = m + \sum_{i=1}^{2} m_{fi}
$$
$$
Z_{\dot{\xi}} = \frac{1}{2} \rho V^2 \sum_{i=1}^{2} \left( S_i \frac{\partial C_{Li}}{\partial \dot{\xi}} \right)
$$
$$
Z_\xi = \frac{1}{2} \rho V^2 \sum_{i=1}^{2} \left( S_i \frac{\partial C_{Li}}{\partial \dot{\xi}} \right)
$$
$$
Z_{\ddot{\theta}} = -\sum_{i=1}^{2} (m_{fi} d_i)
$$
$$
Z_{\dot{\theta}} = \frac{1}{2} \rho V^2 \sum_{i=1}^{2} \left( S_i \frac{\partial C_{Li}}{\partial \dot{\theta}} \right)
$$
$$
Z_\theta = \frac{1}{2} \rho V^2 \sum_{i=1}^{2} \left( S_i \frac{\partial C_{Li}}{\partial \dot{\theta}} \right)
$$
$$
Z_{fp-\alpha i} = \frac{1}{2} \rho V^2 S_{fpi} \frac{\partial C_{fp-Li}(\alpha_{fpi})}{\partial \alpha_{fpi}}
$$
$$
M_{\ddot{\xi}} = -\sum_{i=1}^{2} (m_{fi} d_i)
$$
$$
M_{\dot{\xi}} = -\frac{1}{2} \rho V^2 \sum_{i=1}^{2} \left( S_i \frac{\partial C_{Li}}{\partial \dot{\xi}} d_i \right)
$$
$$
M_\xi = -\frac{1}{2} \rho V^2 \sum_{i=1}^{2} \left( S_i \frac{\partial C_{Li}}{\partial \dot{\xi}} d_i \right)
$$
$$
M_{\ddot{\theta}} = J_{yy} + \sum_{i=1}^{2} (m_{fi} d_i^2)
$$
$$
M_{\dot{\theta}} = \frac{1}{2} \rho V^2 \sum_{i=1}^{2} \left( S_i \frac{\partial C_{Li}}{\partial \dot{\xi}} d_i \right) - \sum_{i=1}^{2} (m_{fi} V d_i)
$$
$$
M_\theta = -\frac{1}{2} \rho V^2 \sum_{i=1}^{2} \left( S_i \frac{\partial C_{Li}}{\partial \dot{\xi}} d_i \right)
$$
$$
M_{fp-\alpha i} = -\frac{1}{2} \rho V^2 S_{fpi} \frac{\partial C_{fp-Li}(\alpha_{fpi})}{\partial \alpha_{fpi}} d_i
$$
通过MATLAB进行计算得到：
- 前后水翼的质量 m_f1 和 m_f2: 0.055652 0.011261
- 水翼的面积 S_1 和 S_2: 0.004374 0.0017069
- 襟翼的面积 S_fp: 0.0008748
- $Z_{\ddot{\xi}}$: 1.5669
- $Z_{\dot{\xi}}$: 31.4649
- $Z_{\xi}$: 31.4649
- $Z_{\ddot{\theta}}$: -0.020074
- $Z_{\dot{\theta}}$: 31.4649
- $Z_{\theta}$: 31.4649
- $Z_{f p-\alpha i}$: 1.0973
- $M_{\ddot{\xi}}$: -0.020074
- $M_{\dot{\xi}}$: -9.4395
- $M_{\xi}$: -9.4395
- $M_{\ddot{\theta}}$: 0.14602
- $M_{\dot{\theta}}$: 9.3271
- $M_{\theta}$: -9.4395
- $M_{f p-\alpha i}$: -0.65841
状态空间方程
状态向量 $x$ 和控制输入向量 $u$ 已经给出，同时还给出了扰动向量 $W$。
将得到的方程整理成 $\dot{x}=Ax+Bu+DW$ 的形式，其中 $A$ 是系统矩阵，$B$ 是输入矩阵，$D$ 是扰动矩阵。把每一个方程中的线性项和非线性项分开,确定矩阵 $A$, $B$ 和 $D$。
由于方程中有常数项，也就是 $(x_{1}-0.045)$ 和 $(x_{3}-0.0349)$。可以通过在每个方程中加上这些常数项的系数，然后从右侧减去相同的项来消除这些项。 所以，状态方程可以写成以下的形式：
[1]杨洋. 全浸式水翼艇纵向运动鲁棒滑模重复控制设计[D].哈尔滨工程大学,2020.
wechat- alipay
