CN108460206A - 一种波浪滑翔器运动预测方法 - Google Patents

Abstract

本发明公开了一种波浪滑翔器运动预测方法，属于海洋航行器波浪滑翔器的动力学建模领域。本发明包括：建立坐标系，选取自由度；计算刚体质量矩阵与附加力；计算科氏力矩阵；计算惯性水动力矩阵；计算类科氏力矩阵；计算回复力矩阵；计算包含推力、舵力的控制力矩阵；计算包含粘性水动力，风浪流力的阻尼力矩阵；求解波浪滑翔器合力；合力与刚体质量矩阵的逆阵相乘，计算广义加速度；进行数值积分求解广义速度，并进一步求解姿态角以及波浪滑翔器在大地坐标系下的位置。本发明提供的波浪滑翔器运动预测方法，物理意义清晰，分析计算简单，能够有效反映波浪滑翔器的动力学特性，能够有效预测波浪滑翔器的运动状态。

Description

一种波浪滑翔器运动预测方法

技术领域

本发明属于海洋航行器波浪滑翔器的动力学建模领域，具体涉及一种波浪滑翔器运动预测方法。

背景技术

波浪滑翔器是新型的海上观测平台，它通过吸收海洋中的波浪能，将其转化为自身航行的动力，通过太阳能板，将太阳能转化为电能为自身设备供电。它具有续航能力强，运营成本低廉，恶劣海况下适应性强、运行噪声小、无污染等众多优点，为各类海洋信息的观测与研究提供了新途径，现已广泛应用于长时间海洋环境探测、生物迁徙研究、气象信息预报等领域。

然而，波浪滑翔器特殊的刚柔多体结构、特殊的潜体结构使得动力学分析与运动控制相关研究较为困难。

文献(Kraus.Wave Glider Dynamic Modeling，Parameter Identification andSimulation[D].University of Hawaii，2012)中对船舶空间运动数学模型建模方法进行修正以适用于波浪滑翔器，构建了一种波浪滑翔器六自由度操纵性模型。文献(卢旭.波浪滑翔器总体技术研究[D].哈尔滨工程大学硕士学位论文，2015).中基于Fossen动力学模型建立了一种波浪滑翔器六自由度操纵性模型。然而，其中各个矩阵如刚体质量矩阵、科氏向心力矩阵等的物理意义和取值不够明确，与本发明中刚体质量矩阵、科氏向心力、恢复力矩阵等的求解方法有明显不同。除此之外，本发明可以设置浮体垂向运动，便于完成不同海浪环境下的动力学仿真与动力学特性分析。除选取的六个主要自由度外，本发明可间接得到系统重心垂向速度、潜体重心垂向速度等运动量，并将其引入动力学建模过程中，使得本发明中波浪滑翔器动力学建模更加科学严谨，能够更加真实刻画波浪滑翔器动力学特性。

文献(贾丽娟.波浪动力滑翔机双体结构工作机理与动力学行为研究.国家海洋技术中心硕士学位论文.2014)和文献(Qi Zhanfeng.Dynamic Modeling and MotionSimulation for Wave Glider[J].Applied Mechanics and Materials.2013，397：285-290.)基于凯恩方法建立波浪滑翔器动力学模型，而本发明基于Fossen动力学模型建立波浪滑翔器动力学模型，二者理论基础具有本质不同。

发明内容

波浪滑翔器特殊的刚柔混合多体结构使得传统单刚体海洋航行器动力学模型无法直接应用于波浪滑翔器动力学分析，本发明提供一种波浪滑翔器运动预测方法，对单刚体海洋航行器动力学模型在建模假设条件下拓展，使其适用于波浪滑翔器动力学分析，进而完成波浪滑翔器运动状态预测。

本发明是这样实现的：

一种波浪滑翔器运动预测方法，其特征在于，包含以下步骤：

步骤一建立坐标系，选取自由度；

建立大地坐标系，系统坐标系，浮体坐标系，潜体坐标系；大地坐标系原点取为地球表面一固定点，向北为x轴正向，向东为y轴正向，向地心为z轴正向；系统坐标系原点位于系统的重心，x⁰方向垂直于系索，并与潜体中纵剖面位于同一竖直面内指向波浪滑翔器前进方向，z⁰方向取系索方向，由浮体端指向潜体端，y⁰符合右手定则；浮体坐标系原点为系索与浮体连接处，x^F指向浮体艏部方向，y^F指向浮体右舷，z^F指向下；潜体坐标系原点为系索与潜体连接处，x^G指向潜体艏部方向，y^G指向潜体右舷，z^G指向下；选择系统重心在系统坐标系下纵向速度u，横向速度v，系索的横倾角φ，系索的纵倾角θ，浮体的艏向角ψ^F，潜体的纵倾角ψ^G为六个计算的自由度；选取的六个自由度构成广义位置η＝[x y φ θ ψ^F ψ^G]^T，广义速度ν＝[u v p q r^F r^G]^T，广义加速度

步骤二计算刚体质量矩阵M与附加力

刚体质量矩阵M为6*6矩阵，计算过程如下：

考虑系统坐标系原点为系统重心，不考虑浮体潜体自身质量分布，刚体质量矩阵M第1-4行第1-4列元素取Fossen动力学刚体质量矩阵第1-4行第1-4列元素为：

其中，m为系统质量，为浮体质量与潜体质量之和m＝m_F+m_G，I_x＝m_Fd_F ²+m_Gd_G ²I_y＝m_Fd_F ²+m_Gd_G ²，d_F为浮体与系索连接点距离系统重心的长度，d_G为潜体与系索连接点距离系统重心的长度；

在浮体和潜体坐标系下分别根据Fossen动力学求得浮体和潜体各自的转艏运动相关联的惯性力，分别转换至系统坐标系，求和得到二者合力，分解为质量系数矩阵与广义加速度乘积，所述质量系数矩阵为刚体质量矩阵M第1-2行第5-6列元素为：

其中，为浮体重心在浮体坐标系下纵向坐标，为潜体重心在潜体坐标系下纵向坐标；γ，χ为换元表达，γ＝c(θ)s(Δψ)，χ＝s(θ)s(Δψ)s(φ)+c(Δψ)c(φ)，c(·)与s(·)为为简化表达形式采用的求余弦值与正弦值的表达形式；Δψ＝ψ^G-ψ^F，是潜体艏向角与浮体艏向角的差；

在浮体和潜体坐标系下分别根据Fossen动力学求得浮体和潜体各自的转艏运动相关联的惯性力，分别转换至系统坐标系，分别求解对于系统重心的力矩，求和得到二者合力矩，分解为质量系数矩阵与广义加速度乘积，所述质量系数矩阵为刚体质量矩阵M第3-4行第5-6列元素为：

计算浮体在浮体坐标系下的纵向和横向加速度，根据Fossen动力学求得浮体转艏惯性力矩，分解为质量系数矩阵与广义加速度乘积，所述质量系数矩阵为刚体质量矩阵M第5行元素为：

包括一项附加力：

其中，为浮体相对于浮体坐标系Z轴的转动惯量；为换元表达， 为系统重心在系统坐标系下垂向加速度；r为系统坐标系转艏角速度；

计算潜体在潜体坐标系下的纵向和横向加速度，根据Fossen动力学求得潜体转艏惯性力矩，分解为质量系数矩阵与广义加速度乘积，所述质量系数矩阵为刚体质量矩阵M第6行元素为：

包括一项附加力：

其中，为潜体相对于潜体坐标系Z轴的转动惯量；

合并得到刚体质量矩阵：

附加力项为：

其中：

κ，μ为换元表达；κ＝c(θ)s(φ)，μ＝c(θ)c(φ)；为浮体重心在浮体坐标系下垂向加速度；

系统坐标系转艏角速度r根据几何关系按下式计算：

步骤三计算科氏力矩阵Force_C；

科氏力矩阵Force_C为6*1矩阵，计算过程如下：

计算系统重心纵向、横向运动、垂向运动、系索横摇、纵摇运动产生的科氏力；取Fossen动力学科氏力系数矩阵项1-4行1-4列元素，与广义速度项相乘，得科氏力为：

[mwq -mwp mwv -mwu 0 0]^T

其中，w为系统重心在系统坐标系下垂向速度；

计算与浮体和潜体各自的转艏运动相关的系统纵向和垂向运动的科氏力；在浮体和潜体坐标系下分别根据Fossen动力学求得浮体和潜体各自的转艏运动相关联的科氏力，分别转换至系统坐标系，求和得到二者合力为：

其中，u^F、v^F、u^G、v^G分别为浮体和潜体分别在浮体坐标系和潜体坐标系下的纵向和横向速度；α，β为换元表达，α＝c(θ)c(Δψ)，β＝s(θ)c(Δψ)s(φ)-s(Δψ)c(φ)；

在浮体和潜体坐标系下分别根据Fossen动力学求得浮体和潜体各自的转艏运动相关联的科氏力，分别转换至系统坐标系，分别求解对于系统重心的力矩，求和得到二者合力矩为：

其中，Γ为换元表达，Γ＝s(θ)s(φ)；

计算浮体在浮体坐标系下的纵向和横向速度，根据Fossen动力学求得浮体转艏科氏力矩为：

计算潜体在潜体坐标系下的纵向和横向速度，根据Fossen动力学求得潜体转艏科氏力矩为：

合并得到科氏力矩阵Force_C；

其中系统重心在系统坐标系下垂向速度w根据浮体重心在浮体坐标系下垂向速度w^F推算：

浮体和潜体分别在浮体坐标系和潜体坐标系下的纵向和横向速度u^F、v^F、u^G、v^G根据系统重心运动与系索横倾角、纵倾角计算：

u^G＝c(θ)(u+qd_G)+Γ(v-pd_G)+Λw

v^G＝c(φ)(v-pd_G)-s(φ)w

其中，Λ为换元表达，

步骤四计算惯性水动力矩阵Force_M_A；

惯性水动力矩阵Force_M_A为6*1矩阵，计算过程如下：

分别在浮体坐标系和潜体坐标系计算浮体和潜体受到的纵向、横向、垂向惯性水动力，和艏摇惯性水动力矩；其中计算过程中不考虑浮体和潜体的横摇、纵摇运动对上述惯性水动力或惯性水动力矩的影响；

浮体受到x方向的惯性水动力：

浮体受到y方向的惯性水动力：

浮体受到z方向的惯性水动力：

浮体转艏的惯性水动力矩：

潜体受到x方向的惯性水动力：

潜体受到y方向的惯性水动力：

潜体受到z方向的惯性水动力：

潜体转艏的惯性水动力矩：

其中，为浮体惯性水动力系数， 为潜体惯性水动力系数；

计算中用到的浮体重心在浮体坐标系下垂向加速度选为与步骤二中计算刚体质量矩阵M时相同；浮体的纵向、横向加速度，潜体的纵向、横向、垂向加速度如下式：

浮体坐标系和潜体坐标系各自所受力矩作为浮体和潜体分别的转艏力矩；浮体坐标系和潜体坐标系各自所受力，分别转换至系统坐标系下，求合力得到对于系统纵向、横向运动的力，求合力矩得到对于系索横摇、纵摇运动的力矩：

步骤五计算类科氏力矩阵Force_C_A；

Force_C_A为6*1矩阵，计算过程如下：

分别在浮体坐标系和潜体坐标系计算浮体和潜体收到的纵向、横向、垂向类科氏力，和艏摇类科氏力矩；

浮体受到x方向的类科氏力：

浮体受到y方向的类科氏力：

浮体受到z方向的类科氏力：

Z_F_C_A＝0

浮体转艏的类科氏力矩：

潜体受到x方向的类科氏力：

潜体受到y方向的类科氏力：

潜体受到z方向的类科氏力：

Z_G_C_A＝0

潜体转艏的类科氏力矩：

其中，

为浮体惯性水动力系数， 为潜体惯性水动力系数；

计算中用到的浮体重心在浮体坐标系下纵向、横向、垂向速度，潜体重心在潜体坐标系下纵向、横向、垂向速度，选为与步骤三中计算科氏力矩阵Force_C时相同；

浮体坐标系和潜体坐标系各自所受力矩作为浮体和潜体分别的转艏力矩；浮体坐标系和潜体坐标系各自所受力，分别转换至系统坐标系下，求合力得到对于系统纵向、横向运动的力，求合力矩得到对于系索横摇、纵摇运动的力矩：

步骤六计算回复力矩阵G；

回复力矩阵G，只包含对于系索横摇、纵摇运动的回复力矩，计算过程为：浮体在浮体坐标系下受竖直向上与潜体湿重等值的力，潜体在潜体坐标系下受竖直向下与潜体湿重等值的力；分别转换至系统坐标系下后，求合力矩得到对于系索横摇、纵摇运动的回复力矩；对于系统整体横摇、纵摇的恢复力矩分别为和其中为潜体湿重；

步骤七计算包含推力、舵力的控制力矩阵T；

分别在浮体坐标系和潜体坐标系计算浮体和潜体受到的纵向、横向、垂向力，和转艏力矩；浮体在浮体坐标系下和潜体在潜体坐标系下的转艏力矩作为浮体和潜体分别的转艏力矩；浮体在浮体坐标系下和潜体在潜体坐标系下受力，分别转换至系统坐标系下，求合力得到对于系统纵向、横向运动的力，求合力矩得到对于系索横摇、纵摇运动的力矩；

步骤八计算包含粘性水动力，风浪流力的阻尼力矩阵D；

分别在浮体坐标系和潜体坐标系计算浮体和潜体受到的纵向、横向、垂向力，和转艏力矩；浮体在浮体坐标系下和潜体在潜体坐标系下的转艏力矩作为浮体和潜体分别的转艏力矩；浮体在浮体坐标系下和潜体在潜体坐标系下受力，分别转换至系统坐标系下，求合力得到对于系统纵向、横向运动的力，求合力矩得到对于系索横摇、纵摇运动的力矩；

步骤九求解波浪滑翔器受合力

步骤十合力与刚体质量矩阵的逆阵相乘，计算广义加速度A＝M^-1·F；

步骤十一进行数值积分求解广义速度，并进一步求解姿态角以及波浪滑翔器在大地坐标系下的位置。

相比于已有波浪滑翔器动力学模型，本发明的有益效果在于：

1、对单刚体海洋航行器动力学模型在建模假设条件下拓展，使其适用于波浪滑翔器，能够用于波浪滑翔器动力学分析与运动状态预测。

2、可以设置浮体垂向运动，便于完成不同海浪环境下的波浪滑翔器动力学分析与运动状态预测。

3、除选取的六个主要自由度外，本发明可间接得到系统重心垂向速度、潜体重心垂向速度等运动量，并将其引入动力学建模过程中，使得本发明中波浪滑翔器动力学模型更加科学严谨，能够更加真实刻画波浪滑翔器动力学特性。

附图说明

图1为本发明波浪滑翔器运动预测方法流程图；

图2为本发明波浪滑翔器动力学模型坐标系示意图。

具体实施方式

下面结合附图举例对本发明做更详细地描述：

结合图1，本发明提供的波浪滑翔器运动预测方法，包含以下步骤：

(一)建立坐标系，选取自由度；

结合图2，其中坐标系建立的具体方法是：建立大地坐标系，系统坐标系，浮体坐标系，潜体坐标系。大地坐标系原点取为地球表面一固定点，向北为x轴正向，向东为y轴正向，向地心为z轴正向。系统坐标系原点位于系统的重心，x⁰方向垂直于系索，并与潜体中纵剖面位于同一竖直面内指向波浪滑翔器前进方向，z⁰方向取系索方向，由浮体端指向潜体端，y⁰符合右手定则。浮体坐标系原点为系索与浮体连接处，x^F指向浮体艏部方向，y^F指向浮体右舷，z^F指向下。潜体坐标系原点为系索与潜体连接处，x^G指向潜体艏部方向，y^G指向潜体右舷，z^G指向下。选择系统重心在系统坐标系下纵向速度u，横向速度v，系索的横倾角φ，系索的纵倾角θ，浮体的艏向角ψ^F，潜体的纵倾角ψ^G为六个计算的自由度。选取的六个自由度构成广义位置η＝[x y φ θ ψ^F ψ^G]^T，广义速度ν＝[u v p q r^F r^G]^T，广义加速度

波浪滑翔器运动状态预测基于波浪滑翔器动力学模型，所述动力学模型的假设条件具体是：

假设1.大地是理想的平面，即不考虑大地的曲率和地球的自转，那么大地坐标系就成为惯性参考系；

假设2.为了方便利用已知的刚体运动学模型对波浪滑翔器进行分析，考虑到实际工作情况中波浪驱动水面机器人的系索在潜体的拉力下大部分时间是绷直的，则假设系索一直处于绷紧状态；

假设3.考虑到实际航行中波浪滑翔器浮体和潜体的转艏角速度较低，不考虑浮体和潜体各自的转艏运动产生的定轴性效应；

假设4.系统重心在系统坐标系下纵向速度u，横向速度v，系索的横倾角φ，系索的纵倾角θ视为反映波浪滑翔器系统整体运动状态的自由度，浮体的艏向角ψ^F，潜体的纵倾角ψ^G视为波浪滑翔器的浮体和潜体局部的自由度，在计算系统整体运动状态相关受力时，不考虑浮体和潜体的形状。

假设5.假设系统为“缓慢运动”，即系统加速度的变化不剧烈，所述加速度包括各个运动状态量的加速度或角加速度。

(二)计算刚体质量矩阵M与附加力

刚体质量矩阵M为6*6矩阵，计算过程包含如下步骤：

1、考虑到系统坐标系原点为系统重心，不考虑浮体潜体自身质量分布，刚体质量矩阵M第1-4行第1-4列元素取Fossen动力学刚体质量矩阵第1-4行第1-4列元素为

其中，m为系统质量，为浮体质量(包含浮体自身质量与惯性质量)与潜体质量(包含潜体自身质量与惯性质量)之和m＝m_F+m_G，I_x＝m_Fd_F ²+m_Gd_G ²I_y＝m_Fd_F ²+m_Gd_G ²，d_F为浮体与系索连接点距离系统重心的长度，d_G为潜体与系索连接点距离系统重心的长度。

2、在浮体和潜体坐标系下分别根据Fossen动力学求得浮体和潜体各自的转艏运动相关联的惯性力，分别转换至系统坐标系，求和得到二者合力，分解为质量系数矩阵与广义加速度乘积，所述质量系数矩阵为刚体质量矩阵M第1-2行第5-6列元素为

其中，为浮体重心在浮体坐标系下纵向坐标，为潜体重心在潜体坐标系下纵向坐标。γ，χ为换元表达，γ＝c(θ)s(Δψ)，χ＝s(θ)s(Δψ)s(φ)+c(Δψ)c(φ)，c(·)与s(·)为为简化表达形式采用的求余弦值与正弦值的表达形式。Δψ＝ψ^G-ψ^F，是潜体艏向角与浮体艏向角的差。

3、在浮体和潜体坐标系下分别根据Fossen动力学求得浮体和潜体各自的转艏运动相关联的惯性力，分别转换至系统坐标系，分别求解对于系统重心的力矩，求和得到二者合力矩，分解为质量系数矩阵与广义加速度乘积，所述质量系数矩阵为刚体质量矩阵M第3-4行第5-6列元素为

4、计算浮体在浮体坐标系下的纵向和横向加速度，根据Fossen动力学求得浮体转艏惯性力矩，分解为质量系数矩阵与广义加速度乘积，所述质量系数矩阵为刚体质量矩阵M第5行元素为

还包括一项附加力

其中，为浮体相对于浮体坐标系Z轴的转动惯量。为换元表达， 为系统重心在系统坐标系下垂向加速度。r为系统坐标系转艏角速度。

5、计算潜体在潜体坐标系下的纵向和横向加速度，根据Fossen动力学求得潜体转艏惯性力矩，分解为质量系数矩阵与广义加速度乘积，所述质量系数矩阵为刚体质量矩阵M第6行元素为

还包括一项附加力

其中，为潜体相对于潜体坐标系Z轴的转动惯量。

6、合并得到刚体质量矩阵

附加力项为

所述的计算刚体质量矩阵M的第4、5、6步中，其中系统重心在系统坐标系下垂向加速度根据刚体运动规律知识，根据浮体重心在浮体坐标系下垂向加速度与其他已知运动状态量推算，如下式。

其中κ，μ为换元表达。κ＝c(θ)s(φ)，μ＝c(θ)c(φ)。浮体重心在浮体坐标系下垂向加速度其选取可根据不同应用场景与仿真需求根据经验设定，可固定为0，或假设为与波浪表面运动相同，或假设为相对波浪表面运动存在固定滞后时间。

其中系统坐标系转艏角速度r根据几何关系按下式计算：

(三)计算科氏力矩阵Force_C；

科氏力矩阵Force_C为6*1矩阵，计算过程包含如下步骤：

1、计算系统重心纵向、横向运动、垂向运动、系索横摇、纵摇运动产生的科氏力。取Fossen动力学科氏力系数矩阵项1-4行1-4列元素，与广义速度项相乘，得相关的科氏力为：

[mwq -mwp mwv -mwu 0 0]^T

其中，w为系统重心在系统坐标系下垂向速度。

2、计算与浮体和潜体各自的转艏运动相关的系统纵向和垂向运动的科氏力。在浮体和潜体坐标系下分别根据Fossen动力学求得浮体和潜体各自的转艏运动相关联的科氏力，分别转换至系统坐标系，求和得到二者合力为：

其中，u^F、v^F、u^G、v^G分别为浮体和潜体分别在浮体坐标系和潜体坐标系下的纵向和横向速度。α，β为换元表达，α＝c(θ)c(Δψ)，β＝s(θ)c(Δψ)s(φ)-s(Δψ)c(φ)。

3、在浮体和潜体坐标系下分别根据Fossen动力学求得浮体和潜体各自的转艏运动相关联的科氏力，分别转换至系统坐标系，分别求解对于系统重心的力矩，求和得到二者合力矩为：

其中，Γ为换元表达，Γ＝s(θ)s(φ)。

4、计算浮体在浮体坐标系下的纵向和横向速度，根据Fossen动力学求得浮体转艏科氏力矩为：

5、计算潜体在潜体坐标系下的纵向和横向速度，根据Fossen动力学求得潜体转艏科氏力矩为：

6、合并得到科氏力矩阵Force_C。

所述的计算科氏力矩阵Force_C的步骤中，其中系统重心在系统坐标系下垂向速度w根据浮体重心在浮体坐标系下垂向速度w^F推算。

浮体和潜体分别在浮体坐标系和潜体坐标系下的纵向和横向速度u^F、v^F、u^G、v^G根据系统重心运动与系索横倾角、纵倾角计算。

u^G＝c(θ)(u+qd_G)+Γ(v-pd_G)+Λw

v^G＝c(φ)(v-pd_G)-s(φ)w

其中，Λ为换元表达，Λ＝s(θ)c(φ)。

(四)计算惯性水动力矩阵Force_M_A；

惯性水动力矩阵Force_M_A为6*1矩阵，计算过程包含如下步骤：

1、分别在浮体坐标系和潜体坐标系计算浮体和潜体受到的纵向、横向、垂向惯性水动力，和艏摇惯性水动力矩。其中计算过程中不考虑浮体和潜体的横摇、纵摇运动对上述惯性水动力或惯性水动力矩的影响。

浮体受到x方向的惯性水动力

浮体受到y方向的惯性水动力

浮体受到z方向的惯性水动力

浮体转艏的惯性水动力矩

潜体受到x方向的惯性水动力

潜体受到y方向的惯性水动力

潜体受到z方向的惯性水动力

潜体转艏的惯性水动力矩

其中，为浮体惯性水动力系数， 为潜体惯性水动力系数。本领域技术人员可参考相关文献，或利用经验公式，或利用计算流体力学方法，或利用试验方法等选取。

计算中用到的浮体重心在浮体坐标系下垂向加速度选为与所述波浪滑翔器运动预测方法第(二)步中计算刚体质量矩阵M时相同。浮体的纵向、横向加速度，潜体的纵向、横向、垂向加速度根据刚体运动规律知识，根据系统重心纵向、横向、垂向加速度与与其他已知运动状态量推算，如下式。

计算中用到的加速度项均取为上一时刻加速度，所述加速度包括各个运动状态量的加速度或角加速度。在系统“缓慢运动”假设下，系统的加速度变化不剧烈，取较小的计算步长，采用上一时刻加速度计算当前时刻惯性水动力，误差在可控范围内。

2、浮体坐标系和潜体坐标系各自所受力矩作为浮体和潜体分别的转艏力矩；浮体坐标系和潜体坐标系各自所受力，分别转换至系统坐标系下，求合力得到对于系统纵向、横向运动的力，求合力矩得到对于系索横摇、纵摇运动的力矩。

(五)计算类科氏力矩阵Force_C_A。Force_C_A为6*1矩阵，计算过程包含如下步骤：

1、参考Fossen动力学中类科氏力系数矩阵，分别在浮体坐标系和潜体坐标系计算浮体和潜体收到的纵向、横向、垂向类科氏力，和艏摇类科氏力矩。其中计算过程中不考虑浮体横摇、纵摇运动对上述类科氏力或类科氏力矩的影响。

浮体受到x方向的类科氏力

浮体受到y方向的类科氏力

浮体受到z方向的类科氏力

Z_F_C_A＝0

浮体转艏的类科氏力矩

潜体受到x方向的类科氏力

潜体受到y方向的类科氏力

潜体受到z方向的类科氏力

Z_G_C_A＝0

潜体转艏的类科氏力矩

其中，

其中，为浮体惯性水动力系数， 为潜体惯性水动力系数。本领域技术人员可参考相关文献，或利用经验公式，或利用计算流体力学方法，或利用试验方法等选取。

计算中用到的浮体重心在浮体坐标系下纵向、横向、垂向速度，潜体重心在潜体坐标系下纵向、横向、垂向速度，选为与所述波浪滑翔器运动预测方法第(三)步中计算科氏力矩阵Force_C时相同。

2、浮体坐标系和潜体坐标系各自所受力矩作为浮体和潜体分别的转艏力矩；浮体坐标系和潜体坐标系各自所受力，分别转换至系统坐标系下，求合力得到对于系统纵向、横向运动的力，求合力矩得到对于系索横摇、纵摇运动的力矩。

(六)计算回复力矩阵G；回复力矩阵G，只包含对于系索横摇、纵摇运动的回复力矩，计算步骤为：浮体在浮体坐标系下受竖直向上与潜体湿重等值的力，潜体在潜体坐标系下受竖直向下与潜体湿重等值的力；分别转换至系统坐标系下后，求合力矩得到对于系索横摇、纵摇运动的回复力矩。对于系统整体横摇、纵摇的恢复力矩分别为和其中为潜体湿重。

(七)计算包含推力、舵力的控制力矩阵T；

(八)计算包含粘性水动力，风浪流力的阻尼力矩阵D；

所述的波浪滑翔器运动预测方法，其中第(七)、(八)步计算控制力矩阵T与阻尼力矩阵D，计算步骤为：分别在浮体坐标系和潜体坐标系计算推力、舵力、粘性水动力；浮体在浮体坐标系下和潜体在潜体坐标系下的转艏力矩作为浮体和潜体分别的转艏力矩；浮体在浮体坐标系下和潜体在潜体坐标系下受力，分别转换至系统坐标系下，求合力得到对于系统纵向、横向运动的力，求合力矩得到对于系索横摇、纵摇运动的力矩。

其中，在浮体坐标系和潜体坐标系计算推力、舵力、粘性水动力等的计算依赖于流体力学知识，可通过计算流体力学或实验等方法获得。

(九)求解波浪滑翔器受合力

(十)合力与刚体质量矩阵的逆阵相乘，计算广义加速度A＝M^-1·F；

(十一)进行数值积分求解广义速度，并进一步求解姿态角以及波浪滑翔器在大地坐标系下的位置。

本领域技术人员按照本发明提供的波浪滑翔器运动预测方法的步骤，能够较为容易地建立波浪滑翔器动力学模型，并能够在不付出创新性劳动的前提下较为容易地基于MATLAB环境或Visual C++环境等本领域常用技术方法实现波浪滑翔器动力学特性数值仿真，并进一步分析波浪滑翔器的动力学特性。

Claims (1)

1.一种波浪滑翔器运动预测方法，其特征在于，包含以下步骤：

步骤一 建立坐标系，选取自由度；

建立大地坐标系，系统坐标系，浮体坐标系，潜体坐标系；大地坐标系原点取为地球表面一固定点，向北为x轴正向，向东为y轴正向，向地心为z轴正向；系统坐标系原点位于系统的重心，x⁰方向垂直于系索，并与潜体中纵剖面位于同一竖直面内指向波浪滑翔器前进方向，z⁰方向取系索方向，由浮体端指向潜体端，y⁰符合右手定则；浮体坐标系原点为系索与浮体连接处，x^F指向浮体艏部方向，y^F指向浮体右舷，z^F指向下；潜体坐标系原点为系索与潜体连接处，x^G指向潜体艏部方向，y^G指向潜体右舷，z^G指向下；选择系统重心在系统坐标系下纵向速度u，横向速度v，系索的横倾角φ，系索的纵倾角θ，浮体的艏向角ψ^F，潜体的纵倾角ψ^G为六个计算的自由度；选取的六个自由度构成广义位置η＝[x y φ θ ψ^F ψ^G]^T，广义速度ν＝[u v p q r^F r^G]^T，广义加速度

步骤二 计算刚体质量矩阵M与附加力

刚体质量矩阵M为6*6矩阵，计算过程如下：

考虑系统坐标系原点为系统重心，不考虑浮体潜体自身质量分布，刚体质量矩阵M第1-4行第1-4列元素取Fossen动力学刚体质量矩阵第1-4行第1-4列元素为：

其中，m为系统质量，为浮体质量与潜体质量之和m＝m_F+m_G，I_x＝m_Fd_F ²+m_Gd_G ²I_y＝m_Fd_F ²+m_Gd_G ²，d_F为浮体与系索连接点距离系统重心的长度，d_G为潜体与系索连接点距离系统重心的长度；

在浮体和潜体坐标系下分别根据Fossen动力学求得浮体和潜体各自的转艏运动相关联的惯性力，分别转换至系统坐标系，求和得到二者合力，分解为质量系数矩阵与广义加速度乘积，所述质量系数矩阵为刚体质量矩阵M第1-2行第5-6列元素为：

其中，为浮体重心在浮体坐标系下纵向坐标，为潜体重心在潜体坐标系下纵向坐标；γ，χ为换元表达，γ＝c(θ)s(Δψ)，χ＝s(θ)s(Δψ)s(φ)+c(Δψ)c(φ)，c(·)与s(·)为为简化表达形式采用的求余弦值与正弦值的表达形式；Δψ＝ψ^G-ψ^F，是潜体艏向角与浮体艏向角的差；

在浮体和潜体坐标系下分别根据Fossen动力学求得浮体和潜体各自的转艏运动相关联的惯性力，分别转换至系统坐标系，分别求解对于系统重心的力矩，求和得到二者合力矩，分解为质量系数矩阵与广义加速度乘积，所述质量系数矩阵为刚体质量矩阵M第3-4行第5-6列元素为：

计算浮体在浮体坐标系下的纵向和横向加速度，根据Fossen动力学求得浮体转艏惯性力矩，分解为质量系数矩阵与广义加速度乘积，所述质量系数矩阵为刚体质量矩阵M第5行元素为：

包括一项附加力：

其中，为浮体相对于浮体坐标系Z轴的转动惯量；为换元表达， 为系统重心在系统坐标系下垂向加速度；r为系统坐标系转艏角速度；

计算潜体在潜体坐标系下的纵向和横向加速度，根据Fossen动力学求得潜体转艏惯性力矩，分解为质量系数矩阵与广义加速度乘积，所述质量系数矩阵为刚体质量矩阵M第6行元素为：

包括一项附加力：

其中，为潜体相对于潜体坐标系Z轴的转动惯量；

合并得到刚体质量矩阵：

附加力项为：

其中：

κ，μ为换元表达；κ＝c(θ)s(φ)，μ＝c(θ)c(φ)；为浮体重心在浮体坐标系下垂向加速度；

系统坐标系转艏角速度r根据几何关系按下式计算：

步骤三 计算科氏力矩阵Force_C；

科氏力矩阵Force_C为6*1矩阵，计算过程如下：

计算系统重心纵向、横向运动、垂向运动、系索横摇、纵摇运动产生的科氏力；取Fossen动力学科氏力系数矩阵项1-4行1-4列元素，与广义速度项相乘，得科氏力为：

[mwq -mwp mwv -mwu 0 0]^T

其中，w为系统重心在系统坐标系下垂向速度；

计算与浮体和潜体各自的转艏运动相关的系统纵向和垂向运动的科氏力；在浮体和潜体坐标系下分别根据Fossen动力学求得浮体和潜体各自的转艏运动相关联的科氏力，分别转换至系统坐标系，求和得到二者合力为：

其中，u^F、v^F、u^G、v^G分别为浮体和潜体分别在浮体坐标系和潜体坐标系下的纵向和横向速度；α，β为换元表达，α＝c(θ)c(Δψ)，β＝s(θ)c(Δψ)s(φ)-s(Δψ)c(φ)；

在浮体和潜体坐标系下分别根据Fossen动力学求得浮体和潜体各自的转艏运动相关联的科氏力，分别转换至系统坐标系，分别求解对于系统重心的力矩，求和得到二者合力矩为：

其中，Γ为换元表达，Γ＝s(θ)s(φ)；

计算浮体在浮体坐标系下的纵向和横向速度，根据Fossen动力学求得浮体转艏科氏力矩为：

计算潜体在潜体坐标系下的纵向和横向速度，根据Fossen动力学求得潜体转艏科氏力矩为：

合并得到科氏力矩阵Force_C；

其中系统重心在系统坐标系下垂向速度w根据浮体重心在浮体坐标系下垂向速度w^F推算：

浮体和潜体分别在浮体坐标系和潜体坐标系下的纵向和横向速度u^F、v^F、u^G、v^G根据系统重心运动与系索横倾角、纵倾角计算：

u^G＝c(θ)(u+qd_G)+Γ(v-pd_G)+Λw

v^G＝c(φ)(v-pd_G)-s(φ)w

其中，Λ为换元表达，Λ＝s(θ)c(φ)；

步骤四 计算惯性水动力矩阵Force_M_A；

惯性水动力矩阵Force_M_A为6*1矩阵，计算过程如下：

分别在浮体坐标系和潜体坐标系计算浮体和潜体受到的纵向、横向、垂向惯性水动力，和艏摇惯性水动力矩；其中计算过程中不考虑浮体和潜体的横摇、纵摇运动对上述惯性水动力或惯性水动力矩的影响；

浮体受到x方向的惯性水动力：

浮体受到y方向的惯性水动力：

浮体受到z方向的惯性水动力：

浮体转艏的惯性水动力矩：

潜体受到x方向的惯性水动力：

潜体受到y方向的惯性水动力：

潜体受到z方向的惯性水动力：

潜体转艏的惯性水动力矩：

其中，为浮体惯性水动力系数， 为潜体惯性水动力系数；

计算中用到的浮体重心在浮体坐标系下垂向加速度选为与步骤二中计算刚体质量矩阵M时相同；浮体的纵向、横向加速度，潜体的纵向、横向、垂向加速度如下式：

浮体坐标系和潜体坐标系各自所受力矩作为浮体和潜体分别的转艏力矩；浮体坐标系和潜体坐标系各自所受力，分别转换至系统坐标系下，求合力得到对于系统纵向、横向运动的力，求合力矩得到对于系索横摇、纵摇运动的力矩：

步骤五 计算类科氏力矩阵Force_C_A；

Force_C_A为6*1矩阵，计算过程如下：

分别在浮体坐标系和潜体坐标系计算浮体和潜体收到的纵向、横向、垂向类科氏力，和艏摇类科氏力矩；

浮体受到x方向的类科氏力：

浮体受到y方向的类科氏力：

浮体受到z方向的类科氏力：

Z_F_C_A＝0

浮体转艏的类科氏力矩：

潜体受到x方向的类科氏力：

潜体受到y方向的类科氏力：

潜体受到z方向的类科氏力：

Z_G_C_A＝0

潜体转艏的类科氏力矩：

其中，

为浮体惯性水动力系数， 为潜体惯性水动力系数；

计算中用到的浮体重心在浮体坐标系下纵向、横向、垂向速度，潜体重心在潜体坐标系下纵向、横向、垂向速度，选为与步骤三中计算科氏力矩阵Force_C时相同；

浮体坐标系和潜体坐标系各自所受力矩作为浮体和潜体分别的转艏力矩；浮体坐标系和潜体坐标系各自所受力，分别转换至系统坐标系下，求合力得到对于系统纵向、横向运动的力，求合力矩得到对于系索横摇、纵摇运动的力矩：

步骤六 计算回复力矩阵G；

回复力矩阵G，只包含对于系索横摇、纵摇运动的回复力矩，计算过程为：浮体在浮体坐标系下受竖直向上与潜体湿重等值的力，潜体在潜体坐标系下受竖直向下与潜体湿重等值的力；分别转换至系统坐标系下后，求合力矩得到对于系索横摇、纵摇运动的回复力矩；对于系统整体横摇、纵摇的恢复力矩分别为和其中为潜体湿重；

步骤七 计算包含推力、舵力的控制力矩阵T；

分别在浮体坐标系和潜体坐标系计算浮体和潜体受到的纵向、横向、垂向力，和转艏力矩；浮体在浮体坐标系下和潜体在潜体坐标系下的转艏力矩作为浮体和潜体分别的转艏力矩；浮体在浮体坐标系下和潜体在潜体坐标系下受力，分别转换至系统坐标系下，求合力得到对于系统纵向、横向运动的力，求合力矩得到对于系索横摇、纵摇运动的力矩；

步骤八 计算包含粘性水动力，风浪流力的阻尼力矩阵D；

分别在浮体坐标系和潜体坐标系计算浮体和潜体受到的纵向、横向、垂向力，和转艏力矩；浮体在浮体坐标系下和潜体在潜体坐标系下的转艏力矩作为浮体和潜体分别的转艏力矩；浮体在浮体坐标系下和潜体在潜体坐标系下受力，分别转换至系统坐标系下，求合力得到对于系统纵向、横向运动的力，求合力矩得到对于系索横摇、纵摇运动的力矩；

步骤九 求解波浪滑翔器受合力

步骤十 合力与刚体质量矩阵的逆阵相乘，计算广义加速度A＝M^-1·F；

步骤十一 进行数值积分求解广义速度，并进一步求解姿态角以及波浪滑翔器在大地坐标系下的位置。