CN106484969B - 一种大包线强机动飞行器动力学高精度仿真方法 - Google Patents

Abstract

一种大包线强机动飞行器动力学高精度仿真方法，属于飞行器动力学与控制领域。该方法首先建立了大包线强机动飞行器的刚体弹性耦合动力学高阶模型，建模过程仅考虑小幅线性振动假设，充分考虑了飞行器大包线强机动飞行过程中显著且快时变的气动力和力矩、发动机推力和力矩、重力、姿态强机动、气动舵快速运动与结构弹性振动之间的相互耦合影响，模型包含了全面的高阶非线性项，刚体运动与弹性振动之间的耦合影响项、姿态机动和气动舵快速运动对弹性振动和姿态的影响项。因此，该高阶模型能够真实反映大包线强机动飞行器真实状态，利用该高阶模型进行仿真分析，可以用于验证飞行器设计的合理性，验证结果较传统模型更准确可靠。

Description

一种大包线强机动飞行器动力学高精度仿真方法

技术领域

本发明涉及一种大包线强机动飞行器动力学高精度仿真方法，属于飞行器动力学与控制领域。

背景技术

利用飞行器的刚体弹性耦合动力学模型进行仿真，得到飞行器位置、速度、姿态、角速度以及弹性振动模态坐标随时间变化的关系曲线，是从动力学角度验证飞行器设计合理性的一个关键手段。如果姿态或弹性振动曲线发散，则说明设计不合理。

随着航天技术的快速发展，新型高性能飞行器越来越成为各航天强国和大国的关注点。大包线、强机动、高超声速飞行成为新型高性能飞行器的典型特点，与传统弹箭和航空器差异显著。

大包线、高超声速飞行导致新型飞行器受到的气动力矩、环境干扰力矩大，且变化明显。强机动飞行导致飞行器本体姿态角速度和角加速度、气动舵角速度和角加速度均很显著。并且，新型飞行器也多为面对称外形，横航向动力学耦合严重。在气动力矩、环境干扰力矩、飞行器本体和气动舵剧烈运动、横航向耦合严重的情况下，飞行器刚体运动和弹性振动之间耦合特性复杂、强烈，传统的刚体弹性耦合动力学模型由于未全面考虑上述因素的影响，不再完全适用于新型的大包线强机动飞行器特点和设计需求。利用传统的刚体弹性耦合动力学模型进行仿真，得到的结果不能反映飞行器真实状态，难以准确评估飞行器设计的合理性。

发明内容

本发明要解决的技术问题是：克服现有技术的不足，提供一种大包线强机动飞行器动力学高精度仿真方法，能够真实反映大包线强机动飞行器状态，准确评估飞行器设计的合理性。

本发明的技术解决方案是：一种大包线强机动飞行器动力学高精度仿真方法，包括如下步骤：

(1)建立大包线强机动飞行器的刚体弹性耦合动力学高阶模型，所述刚体弹性耦合动力学高阶模型包括地面坐标系下的飞行器平动动力学方程：

飞行器本体坐标系下的飞行器转动动力学方程：

飞行器本体坐标系下的飞行器弹性振动动力学方程：

A_bg为地面坐标系到本体坐标系的坐标转换矩阵；

A_bi为惯性坐标系到本体坐标系的坐标转换矩阵；

A_brk为气动舵k本体坐标系到飞行器本体坐标系的坐标转换矩阵；

A_bv为速度坐标系到本体坐标系的坐标转换矩阵；

A_ei为惯性坐标系到地球坐标系的坐标转换矩阵；

A_gb为本体坐标系到地面坐标系的坐标转换矩阵；

A_ge为地球坐标系到地面坐标系的坐标转换矩阵；

A_grk为气动舵k本体坐标系到地面坐标系的坐标转换矩阵；

A_og为地面坐标系到弹道坐标系的坐标转换矩阵；

A_rkb为飞行器本体坐标系到气动舵k本体坐标系的坐标转换矩阵；

A_vo为弹道坐标系到速度坐标系的坐标转换矩阵；

d_bg为b系原点相对g系原点的矢径在b系中的表示；

d_dg为b系原点相对g系原点的矢径在g系中的表示；

d_gi为g系原点相对i系原点的矢径在i系中的表示；

d_i为b系原点相对i系原点的矢径在i系中的表示；

d_rk为气动舵k本体系原点在b系中的未变形矢径；

F_bg为飞行器弹性振动对飞行器转动的耦合系数矩阵；

F_frk为飞行器弹性振动对气动舵k转动的耦合系数矩阵；

g为重力加速度；

J_bb为飞行器相对b系的转动惯量矩阵；

J_rk为气动舵k相对本体坐标系的转动惯量矩阵；

M为飞行器的总质量；

M_rk为气动舵k的质量；

为飞行器弹性振动的等效质量矩阵；

P_bg为飞行器转动对飞行器平动的耦合系数矩阵；

P_f为飞行器弹性振动对飞行器平动的耦合系数矩阵；

P_rk为气动舵k转动对飞行器平动的耦合系数矩阵；

Q_b为作用于飞行器的气动力矩和推力产生的力矩；

Q_bc为发动机推力产生的力矩在飞行器本体坐标系下的表示；

Q_bq为作用于飞行器上的气动力矩在飞行器本体坐标系下的表示；

Q_dc为发动机推力在飞行器本体坐标系下的表示；

Q_dq为气动力在飞行器本体坐标系下的表示；

Q_rk为气动舵的驱动力矩和作用于气动舵上的干扰力矩对弹性振动产生的广义力；

为气动分布力对飞行器弹性振动模态坐标的广义力；

q_dq为气动分布力在飞行器体坐标系上的表示；

R_{bg_rk}为气动舵k转动对飞行器转动的耦合系数矩阵；

r_b为飞行器机身上任一质点dm相对i系的矢径；

r_rk为气动舵k上任一质点dm相对i系的矢径；

为飞行器本体的模态振型函数矩阵；

为在飞行器本体与发动机连接点处的值；

为在飞行器本体与气动舵k连接点处的值；

η_b为飞行器弹性振动的模态坐标列向量；

ρ_b为飞行器机身上任一质点dm相对b系坐标原点的未变形矢径；

ρ_rk为气动舵k上任一质点dm相对本体系的矢径；

ρ_rk0为气动舵k的质心在其本体系中的表示；

ω_{bg_b}为b系相对g系的角速度在b系中的表示；

ω_{bi_b}＝ω_{bg_b}+A_biω_{gi_i}为b系相对i系的角速度在b系中的表示；

ω_{gi_i}为g系相对i系的角速度在i系中的表示；

ω_{rkb_rk}为r_k系相对b系的角速度在r_k系中的表示；

ω_{rkg_rk}＝ω_{rkb_rk}+A_rkbω_{bg_b}为r_k系相对g系的角速度在r_k系中的表示；

ω_{rki_rk}＝ω_{rkb_rk}+A_rkb(ω_{bg_b}+A_biω_{gi_i})为r_k系相对i系的角速度在r_k系中的表示；

Λ_b为飞行器弹性振动的特征频率对角阵；

ξ_b为飞行器弹性振动的粘性阻尼对角阵；

i系为惯性坐标系，b系为飞行器本体坐标系，r_k系为气动舵k本体坐标系，g系为地面坐标系。

(2)利用步骤(1)得到的高阶模型，根据飞行器位置、速度、姿态、角速度、弹性振动模态坐标以及模态坐标速度的初始值以及除这些量以外的所有其他参数的设计值进行飞行器刚体弹性耦合动力学仿真，得到飞行器位置、速度、姿态、角速度以及弹性振动模态坐标随时间变化的关系曲线，用于从动力学角度验证飞行器设计的合理性。

本发明与现有技术相比的优点在于：

(1)本发明提出的动力学高阶模型是基于普适的力学原理推导得出，除了弹性振动作小幅线性振动假设以外，未做其它简化处理，包含了完整的高阶非线性耦合影响项，能够真实表征大包线强机动飞行器的动力学特征。

(2)本发明提出的动力学高阶模型包含飞行器平动与弹性振动相互耦合影响项、飞行器平动与转动相互耦合影响项(飞行器本体坐标系原点在系统质心时此耦合影响消失)、飞行器转动与弹性振动相互耦合影响项、气动舵运动对飞行器平动的影响项、气动舵运动对飞行器转动的影响项、气动舵运动对弹性振动的影响项、弹性振动对飞行器转动惯量的时变影响项，影响因素全面；模型中的广义力为包括重力和气动力对飞行器平动的合力、气动力和发动机对飞行器转动的合力矩、气动力和发动机对弹性振动的广义力，全面、充分考虑了重力、气动、发动机等各种外部激励对刚体和弹性动力学的影响。基于该模型进行动力学仿真，能够真实反映大包线强机动飞行器状态，准确评估飞行器设计的合理性。

(3)本发明提出的动力学高阶模型还可作为其他简化方程(如控制系统设计采用的线性化方程)推导的依据，相较于传统模型，推导结果更准确可靠。

附图说明

图1为飞行器结构及相关坐标系示意图；

图2为惯性坐标系和地球坐标系示意图；

图3为地球坐标系和地面坐标系示意图；

图4为地面坐标系和弹道坐标系示意图；

图5为速度坐标系和飞行器本体坐标系示意图。

具体实施方式

定义自变量为任意实数a的矩阵函数A_x(a)、A_y(a)、A_z(a)如下：

其中实数为自变量。

定义自变量为任意3维列向量的叉乘矩阵算子如下：

其中为任意。

飞行器包含1个刚体弹性耦合的机身和N个刚性气动舵，结构图如图1所示。为了建立耦合动力学模型，引入如下坐标系：

(1)惯性坐标系o_ix_iy_iz_i(简称i系)

原点o_i在地心；o_ix_i轴在赤道面内指向春分点；o_iy_i轴在赤道面内与o_ix_i轴垂直指向东；o_iz_i轴按右手法则确定。单位向量为i。

(2)飞行器本体坐标系o_bx_by_bz_b(简称b系)

o_bx_by_bz_b为直角坐标系并与飞行器本体固连。原点o_b取在飞行器的质心上；o_bx_b轴与飞行器纵轴重合，指向头部为正；o_by_b轴位于飞行器纵向对称面内与o_bx_b轴垂直，指向上为正；o_bz_b轴按右手法则确定。单位向量为b。

(3)气动舵k本体坐标系o_rkx_rky_rkz_rk(简称rk系)

o_rkx_rky_rkz_rk为直角坐标系并与气动舵k本体固连。原点o_rk取在气动舵k与飞行器连接处的中心上；o_rky_rk轴与气动舵k旋转轴重合，平行于o_by_b轴；o_rkz_rk轴垂直于舵面，指向右为正；o_rkx_rk轴按右手法则确定。单位向量为r_k。

(4)地球坐标系o_ex_ey_ez_e(或称ECEF，简称e系)

o_ex_ey_ez_e为直角坐标系并与地球固连。原点o_e位于地心；o_ex_e轴在赤道面内指向本初子午线；o_ey_e轴在赤道面内与x_e轴垂直，指向东为正；o_ez_e轴按右手法则确定。单位向量为e。

(5)地面坐标系o_gx_gy_gz_g(或称北天东地理坐标系，简称g系)

o_gx_gy_gz_g为直角坐标系并与地球表面固连。原点o_g取在飞行器起飞点上(严格地说，应取在起飞瞬时飞行器的质心上)；o_gx_g轴在当地水平面内沿当地经线的切线方向，指向北为正；o_gy_g轴沿当地地理垂线的方向，向上为正；o_gz_g轴在当地水平面内沿当地纬线的切线方向，指向东为正。单位向量为g。

(6)速度坐标系o_vx_vy_vz_v(简称v系)

o_vx_vy_vz_v为直角坐标系并与飞行器速度矢量固连。原点o_v取在飞行器的质心上；o_vx_v轴与飞行器速度矢量重合；o_vy_v轴位于飞行器纵向对称面内与o_vx_v轴垂直，指向上为正；o_vz_v轴按右手法则确定。

(7)弹道坐标系o_ox_oy_oz_o(简称o系)

o_ox_oy_oz_o为直角坐标系并与飞行器速度矢量固连。原点o_o取在飞行器的质心上；o_ox_o轴与飞行器速度矢量重合；o_oy_o轴位于包含速度矢量的铅垂面内并与o_ox_o轴垂直，指向上为正；o_oz_o轴按右手法则确定。

如图2所示，惯性坐标系和地球坐标系之间的相对方位可由本初子午线的恒星时角α_s0确定，则由惯性坐标系到地球坐标系的坐标转换矩阵为

A_ei＝A_z(α_s0)

如图3所示，地球坐标系和地面坐标系之间的相对方位可由飞行器起飞点的地理经度λ和纬度Φ确定。由地球坐标系到地面坐标系的转换矩阵可表示为

A_ge＝A_x(-90°)A_y(-90°-Φ)A_z(λ)

如图4所示，地面坐标系和弹道坐标系之间的相对方位可由两个角度确定，分别定义如下：

1)弹道倾角θ_v：飞行器质心的速度矢量V(即o_vx_v轴)与水平面之间的夹角。速度矢量指向水平面上方，θ_v角为正；反之为负。

2)弹道偏角ψ_v：飞行器质心的速度矢量V在水平面内投影与地面坐标系的o_gx_g轴正向间的夹角。速度矢量偏东时，ψ_v为正；反之为负。

由地面坐标系到弹道坐标系的坐标转换矩阵为

A_og＝A_z(θ_v)A_y(-ψ_v)

弹道坐标系和速度坐标系之间的相对方位由倾斜角(或倾侧角)γ_v确定。倾侧角：飞行器纵向对称面与铅垂面的夹角。从飞行器尾部向前看，若纵向对称面向右倾斜，则γ_v为正；反之为负。由弹道坐标系到速度坐标系的坐标转换矩阵为

A_vo＝A_x(γ_v)

如图5所示，速度坐标系和飞行器本体坐标系之间的相对方位可由两个角度确定，分别定义如下：

1)攻角α：飞行器质心的速度矢量V(即o_vx_v轴)在飞行器纵向对称面o_bx_by_b上的投影与o_bx_b轴之间的夹角。若o_bx_b轴位于V的投影线的上方(即产生正升力)时，攻角α为正；反之为负。

2)侧滑角β：速度矢量V与飞行器纵向对称面之间的夹角。沿飞行方向视察，若来流从右侧流向飞行器(即产生负侧向力)，则所对应的侧滑角β为正；反之为负。

由速度坐标系到飞行器本体坐标系的坐标转换矩阵为

A_bv＝A_z(α)A_y(β)

从理论力学或分析力学中知道，对于n个自由度的非保守系统的动力学方程可写成如下的拉格朗日方程的形式

或

其中L＝T-U称为拉格朗日函数，Q_j为对应广义坐标q_j的广义力，T和U分别是系统的动能和势能。在微振动理论中，动能T与广义坐标无关(因质量是常量)，即则上述方程可写为

如果考虑阻尼的情况，定义耗散力：当仅考虑粘性阻尼的情况，作用在质点上的阻力是线性非保守力，此力大小与速度的一次方成正比，方向相反。因这种力使机械能耗散，故又称耗散力。

经过推导，如果令Q_Rj为对应于广义坐标q_j的广义耗散力，则

其中D称为耗散能，C是广义阻尼矩阵，q＝[q₁ q₂ … q_n]^T是由q_j(j＝1,2,…,n)构成的广义坐标列阵。

如果单独把耗散能提出来，对上述的拉格朗日方程进行修改

此时上述方程中的Q_j是不包括阻尼力的广义力。

采用混合坐标法对柔性结构进行模化，并依据拉格朗日方程推导刚弹耦合动力学模型，推导过程不做任何其它假设、简化和省略，可以得到大包线强机动飞行器刚体弹性耦合动力学模型。

动力学模型的建立过程：

1、飞行器上各质点速度的导出：

1)飞行器机身

令飞行器机身上任一质点dm相对b系坐标原点的未变形矢径为ρ_b，变形位移为u_b，则该质点相对i系的矢径r_b可表示为

r_b＝d_i+A_ib(ρ_b+u_b)＝A_igd_bg+d_gi+A_ib(ρ_b+u_b)

其中d_i为b系原点相对i系原点的矢径，d_i＝A_igd_bg+d_gi。d_bg为b系原点相对g系原点的矢径，d_gi为g系原点相对i系原点的矢径。

该质点相对于i系的速度υ_b为：

其中

进而有

令

g_b＝ρ_b+u_b，g_bi＝d_i+A_ib(ρ_b+u_b)，

进一步有：

2)气动舵k(k＝1,2,3,…,N)

假设气动舵k上任一质点dm相对本体系的矢径为ρ_rk，则该质点相对i系的矢径r_rk可表示为

r_rk＝A_igd_bg+d_gi+A_ibd_rk+A_ibu_rk+A_irkρ_rk

其中d_rk和u_rk为气动舵k本体系原点在b系中的未变形矢径和弹性变形位移。

进一步求导可得气动舵k上任一质点的速度υ_rk为

其中g_{rk_i}＝d_i+A_ib(u_rk+d_rk)+A_irkρ_rk，g_{rk_b}＝u_rk+d_rk+A_brkρ_rk，

在小幅线性振动的前提下，

2、动能

1)飞行器机身的动能T_b

其中

2)气动舵k的动能T_rk

其中

3)总的动能T

其中

考虑到飞行器本体坐标系的原点位于质心，即可以得到

3、势能U

总势能U为重力势能和弹性势能之和，即

其中Λ_b飞行器弹性振动的特征频率对角阵。

4、耗散能D

飞行器的耗散能D为

其中C_b＝2ξ_bΛ_b，ξ_b为线性粘性阻尼对角阵。

5、动力学方程

选择飞行器本体相对地面坐标系的位置和姿态、气动舵k转角、以及弹性振动的模态坐标为系统的广义坐标，将动能T、势能U、耗散能D的表达式带入Lagrange方程可得：

其中Q_d、Q_b、Q_η和Q_rk分别为对应于位置、姿态角、弹性振动的模态坐标和气动舵k转角的广义力。

将上面三式展开，可得最终的平动、转动、弹性振动方程如下：

1)飞行器平动动力学方程

在地面坐标系下，飞行器平动动力学方程为

2)飞行器转动动力学方程

在飞行器本体坐标系下，飞行器转动动力学方程为

3)飞行器弹性振动动力学方程

在飞行器本体坐标系下，飞行器弹性振动动力学方程为

符号说明：

A_bg为地面坐标系到本体坐标系的坐标转换矩阵；

A_bi为惯性坐标系到本体坐标系的坐标转换矩阵；

A_brk为气动舵k本体坐标系到飞行器本体坐标系的坐标转换矩阵；

A_bv为速度坐标系到本体坐标系的坐标转换矩阵；

A_ei为惯性坐标系到地球坐标系的坐标转换矩阵；

A_gb为本体坐标系到地面坐标系的坐标转换矩阵；

A_ge为地球坐标系到地面坐标系的坐标转换矩阵；

A_grk为气动舵k本体坐标系到地面坐标系的坐标转换矩阵；

A_ib为本体坐标系到惯性坐标系的坐标转换矩阵；

A_ig为地面坐标系到惯性坐标系的坐标转换矩阵；

A_irk为气动舵k本体坐标系到惯性坐标系的坐标转换矩阵；

A_og为地面坐标系到弹道坐标系的坐标转换矩阵；

A_rkb为飞行器本体坐标系到气动舵k本体坐标系的坐标转换矩阵；

A_vo为弹道坐标系到速度坐标系的坐标转换矩阵；

D为飞行器的耗散能；

d_bg为b系原点相对g系原点的矢径在b系中的表示；

d_dg为b系原点相对g系原点的矢径在g系中的表示；

d_gi为g系原点相对i系原点的矢径在i系中的表示；

d_i为b系原点相对i系原点的矢径在i系中的表示；

d_rk为气动舵k本体系原点在b系中的未变形矢径；

F_bg为飞行器弹性振动对飞行器转动的耦合系数矩阵；

F_frk为飞行器弹性振动对气动舵k转动的耦合系数矩阵；

g为重力加速度；

J_bb为飞行器相对b系的转动惯量矩阵；

J_rk为气动舵k相对本体坐标系的转动惯量矩阵；

M为飞行器的总质量；

M_rk为气动舵k的质量；

为飞行器弹性振动的等效质量矩阵；

N为气动舵的个数；

P_bg为飞行器转动对飞行器平动的耦合系数矩阵；

P_f为飞行器弹性振动对飞行器平动的耦合系数矩阵；

P_rk为气动舵k转动对飞行器平动的耦合系数矩阵；

Q_b为作用于飞行器的气动力矩和推力产生的力矩；

Q_bc为发动机推力产生的力矩在飞行器本体坐标系下的表示；

Q_bq为作用于飞行器上的气动力矩在飞行器本体坐标系下的表示；

Q_dc为发动机推力在飞行器本体坐标系下的表示；

Q_dq为气动力在飞行器本体坐标系下的表示；

Q_rk为气动舵的驱动力矩和作用于气动舵上的干扰力矩对弹性振动产生的广义力；

为气动分布力对飞行器弹性振动模态坐标的广义力；

q_dq为气动分布力在飞行器体坐标系上的表示；

R_{bg_rk}为气动舵k转动对飞行器转动的耦合系数矩阵；

r_b为飞行器机身上任一质点dm相对i系的矢径；

r_rk为气动舵k上任一质点dm相对i系的矢径；

T为飞行器的总动能；

T_b为飞行器机身的动能；

T_rk为气动舵k的动能；

t为时间；

U为飞行器的势能；

u_b为飞行器机身上任一质点dm的变形位移；

u_rk为气动舵k上任一质点dm的弹性变形位移；

v_b为飞行器机身上任一质点dm相对i系的速度；

v_rk为气动舵k上任一质点dm相对i系的速度；

α为攻角；

α_s0为本初子午线的恒星时角；

β为侧滑角；

为飞行器本体的模态振型函数矩阵；

为在飞行器本体与发动机连接点处的值；

为在飞行器本体与气动舵k连接点处的值；

γ_v为倾斜角(或倾侧角)；

η_b为飞行器弹性振动的模态坐标列向量；

λ为地理经度；

Φ为地理纬度；

ψ_v为弹道偏角；

θ_r为弹道倾角；

θ_rk为气动舵k的转角；

ρ_b为飞行器机身上任一质点dm相对b系坐标原点的未变形矢径；

ρ_rk为气动舵k上任一质点dm相对本体系的矢径；

ρ_rk0为气动舵k的质心在其本体系中的表示；

ω_{bg_b}为b系相对g系的角速度在b系中的表示；

ω_{bi_b}＝ω_{bg_b}+A_biω_{gi_i}为b系相对i系的角速度在b系中的表示；

ω_{gi_i}为g系相对i系的角速度在i系中的表示；

ω_{rkb_rk}为r_k系相对b系的角速度在r_k系中的表示；

ω_{rkg_rk}＝ω_{rkb_rk}+A_rkbω_{bg_b}为r_k系相对g系的角速度在r_k系中的表示；

ω_{rki_rk}＝ω_{rkb_rk}+A_rkb(ω_{bg_b}+A_biω_{gi_i})为r_k系相对i系的角速度在r_k系中的表示；

Λ_b为飞行器弹性振动的特征频率对角阵；

ξ_b为飞行器弹性振动的粘性阻尼对角阵。

利用得到的高阶模型，根据飞行器位置、速度、姿态、角速度、弹性振动模态坐标以及模态坐标速度的初始值以及这些量以外的所有其他参数的设计值进行飞行器刚体弹性耦合动力学仿真，得到飞行器位置、速度、姿态、角速度以及弹性振动模态坐标随时间变化的关系曲线，从而从动力学角度验证飞行器设计的合理性。

从上面推导高阶模型的过程可以看出本发明的动力学高阶模型具有如下特点：

(1)建模过程仅考虑小幅线性振动假设，不做其它假设。

(2)充分考虑了飞行器大包线强机动飞行过程中显著且快时变的气动力和力矩、发动机推力和力矩、重力、姿态强机动、气动舵快速运动与结构弹性振动之间的相互耦合影响，除弹性振动(小幅线性振动)外其他因素均不做小量假设。模型包含了全面的高阶非线性项，刚体运动与弹性振动之间的耦合影响项、姿态机动和气动舵快速运动对弹性振动和姿态的影响项。

(3)从普适的力学一般原理出发推导耦合模型，推导过程中保留所有高阶项，不做忽略。

因此，该高阶模型能够真实反映大包线强机动飞行器真实状态，利用该高阶模型进行仿真分析，可以用于验证飞行器设计的合理性，验证结果较传统模型更准确可靠，同时该高阶模型也可用于大包线强机动飞行器的各种简化模型推导、分析。

本发明未详细说明部分属于本领域技术人员公知技术。

Claims (1)

1.一种大包线强机动飞行器动力学高精度仿真方法，其特征在于包括如下步骤：

(1)建立大包线强机动飞行器的刚体弹性耦合动力学高阶模型，所述刚体弹性耦合动力学高阶模型包括地面坐标系下的飞行器平动动力学方程：

飞行器本体坐标系下的飞行器转动动力学方程：

飞行器本体坐标系下的飞行器弹性振动动力学方程：

A_bg为地面坐标系到本体坐标系的坐标转换矩阵；

A_bi为惯性坐标系到本体坐标系的坐标转换矩阵；

A_brk为气动舵k本体坐标系到飞行器本体坐标系的坐标转换矩阵；

A_bv为速度坐标系到本体坐标系的坐标转换矩阵；

A_ei为惯性坐标系到地球坐标系的坐标转换矩阵；

A_gb为本体坐标系到地面坐标系的坐标转换矩阵；

A_ge为地球坐标系到地面坐标系的坐标转换矩阵；

A_grk为气动舵k本体坐标系到地面坐标系的坐标转换矩阵；

A_og为地面坐标系到弹道坐标系的坐标转换矩阵；

A_rkb为飞行器本体坐标系到气动舵k本体坐标系的坐标转换矩阵；

A_vo为弹道坐标系到速度坐标系的坐标转换矩阵；

d_bg为b系原点相对g系原点的矢径在b系中的表示；

d_dg为b系原点相对g系原点的矢径在g系中的表示；

d_gi为g系原点相对i系原点的矢径在i系中的表示；

d_i为b系原点相对i系原点的矢径在i系中的表示；

d_rk为气动舵k本体系原点在b系中的未变形矢径；

F_bg为飞行器弹性振动对飞行器转动的耦合系数矩阵；

F_frk为飞行器弹性振动对气动舵k转动的耦合系数矩阵；

g为重力加速度；

J_bb为飞行器相对b系的转动惯量矩阵；

J_rk为气动舵k相对本体坐标系的转动惯量矩阵；

M为飞行器的总质量；

M_rk为气动舵k的质量；

为飞行器弹性振动的等效质量矩阵；

P_bg为飞行器转动对飞行器平动的耦合系数矩阵；

P_f为飞行器弹性振动对飞行器平动的耦合系数矩阵；

P_rk为气动舵k转动对飞行器平动的耦合系数矩阵；

Q_b为作用于飞行器的气动力矩和推力产生的力矩；

Q_bc为发动机推力产生的力矩在飞行器本体坐标系下的表示；

Q_bq为作用于飞行器上的气动力矩在飞行器本体坐标系下的表示；

Q_dc为发动机推力在飞行器本体坐标系下的表示；

Q_dq为气动力在飞行器本体坐标系下的表示；

Q_rk为气动舵的驱动力矩和作用于气动舵上的干扰力矩对弹性振动产生的广义力；

为气动分布力对飞行器弹性振动模态坐标的广义力；

q_dq为气动分布力在飞行器体坐标系上的表示；

R_{bg_rk}为气动舵k转动对飞行器转动的耦合系数矩阵；

r_b为飞行器机身上任一质点dm相对i系的矢径；

r_rk为气动舵k上任一质点dm相对i系的矢径；

为飞行器本体的模态振型函数矩阵；

为在飞行器本体与发动机连接点处的值；

为在飞行器本体与气动舵k连接点处的值；

η_b为飞行器弹性振动的模态坐标列向量；

ρ_b为飞行器机身上任一质点dm相对b系坐标原点的未变形矢径；

ρ_rk为气动舵k上任一质点dm相对本体系的矢径；

ρ_rk0为气动舵k的质心在其本体系中的表示；

ω_{bg_b}为b系相对g系的角速度在b系中的表示；

ω_{bi_b}＝ω_{bg_b}+A_biω_{gi_i}为b系相对i系的角速度在b系中的表示；

ω_{gi_i}为g系相对i系的角速度在i系中的表示；

ω_{rkb_rk}为r_k系相对b系的角速度在r_k系中的表示；

ω_{rkg_rk}＝ω_{rkb_rk}+A_rkbω_{bg_b}为r_k系相对g系的角速度在r_k系中的表示；

ω_{rki_rk}＝ω_{rkb_rk}+A_rkb(ω_{bg_b}+A_biω_{gi_i})为r_k系相对i系的角速度在r_k系中的表示；

Λ_b为飞行器弹性振动的特征频率对角阵；

ξ_b为飞行器弹性振动的粘性阻尼对角阵；

i系为惯性坐标系，b系为飞行器本体坐标系，r_k系为气动舵k本体坐标系，g系为地面坐标系；

(2)利用步骤(1)得到的高阶模型，根据飞行器位置、速度、姿态、角速度、弹性振动模态坐标以及模态坐标速度的初始值以及除这些量以外的所有其他参数的设计值进行飞行器刚体弹性耦合动力学仿真，得到飞行器位置、速度、姿态、角速度以及弹性振动模态坐标随时间变化的关系曲线，用于从动力学角度验证飞行器设计的合理性。