CN112364571A - 大型复杂耦合航天器动力学模型建模方法 - Google Patents

Abstract

本发明涉及自动控制技术、力学建模领域，为提出一种大型复杂耦合航天器动力学模型建模方法，使建立的数学模型更加贴合实际卫星模型。本发明，大型复杂耦合航天器动力学模型建模方法，步骤如下：S1、建立航天器太阳帆板简化模型；S2、建立液体晃动动力学模型；S3、建立刚液航天器耦合模型；S4、建立的液体晃动模型进行补充，以获得更贴近设计晃动规律的液体晃动数学模型；S5、需要同时考虑S1中柔性附件振动和S2、S3中液体晃动对刚体的耦合影响，建立大型柔性充液航天器动力学模型；S6、最后根据刚柔航天器、刚液航天器模型分析，建立大型柔性充液航天器模型。本发明主要应用于航天器动力学模型建模。

Description

大型复杂耦合航天器动力学模型建模方法

技术领域

本发明涉及自动控制技术、力学建模领域。具体涉及一种大型复杂耦合航天器动力学模型建模方法。

背景技术

随着空间技术的不断发展与航天需求的不断增长，航天器主要沿着两个大的方向发展：一是微小化、集群化，多个微小卫星协同执行某一任务；二是尺寸大、结构复杂，让卫星在太空中长时间执行任务。针对于大型卫星，结构越来越复杂，尺寸也越来越大，但由于发射成本及运载能力限制的原因，在减轻航天器重量的同时，结构的挠性也越来越大。在航天器完成诸如机动、转向和空中对接等动作时，很容易激起柔性结构的振动，这种振动会增加星体姿态调整时间、影响指向精度及航天器精密仪器的正常工作。另外，随着航天器运载能力、机动能力及长寿命等性能的提高，液体燃料质量占航天器总质量的比值不断增大，这种不利影响也随之增大，液体晃动问题在充液航天器的姿态控制系统设计中无法避免，燃料质量的增加容易导致航天器储液箱液体发生强烈晃动，严重影响星体姿态控制的精度与稳定性。因此，对于大型复杂耦合航天器动力学模型的建模主要考虑刚体、柔性附件以及贮液箱中液体晃动的力与力矩的作用。

模型作为控制的基础，对其进行描述与分析有利于航天器控制器的设计。建立大型复杂耦合航天器动力学模型主要分为三个部分：1、分析柔性附件的振动；2、分析液体晃动对刚体运动的影响；3、根据动量矩守恒原理建立大型柔性充液航天器动力学模型。

发明内容

为克服现有技术的不足，本发明旨在提出一种大型复杂耦合航天器动力学模型建模方法，使建立的数学模型更加贴合实际卫星模型，并为航天器控制器的设计提供基础。为此，本发明采取的技术方案是，大型复杂耦合航天器动力学模型建模方法，步骤如下：

S1、先根据太阳能帆板的结构尺寸，忽略对航天器太阳能帆板的动力学影响微弱的小的杆件结构，建立航天器太阳帆板简化模型；

S2、应用动量矩守恒定理建立航天器系统的姿态动力学方程，利用弹簧质量等效模型的前两阶模态对液体晃动部分进行等效，建立液体晃动动力学模型；

S3、采用CFD与等效力学结合的方法，建立刚液航天器耦合模型，具体是采用CFD软件的流体体积函数方法进行三维液体晃动仿真计算，得到椭球体贮箱内流动参数的变化规律，并进行深入分析，从而完善S2采用等效力学模型所得的液体晃动数学模型；

S4、综合考虑S2、S3两种液体建模方法，依据CFD分析软件获得的液体晃动内部参数结果对等效力学方法建立的液体晃动模型进行补充，以获得更贴近设计晃动规律的液体晃动数学模型；

S5、需要同时考虑S1中柔性附件振动和S2、S3中液体晃动对刚体的耦合影响，建立大型柔性充液航天器动力学模型；

S6、最后根据刚柔航天器、刚液航天器模型分析，建立大型柔性充液航天器模型。

步骤S1详细步骤如下：

考虑到太阳能帆板是一个典型的悬臂梁结构，为使分析更加直观，首先忽略刚体运动，只考虑刚体和柔性附件连接处力的作用，单独分析悬臂梁结构，采用假设模态法，将太阳能帆板简化为欧拉-伯努利悬臂梁,令P(x,t)是单位长度悬臂梁的横向外力分布，m(x)是悬臂梁的质量分布，EI(x)是悬臂梁的刚度分布，E是弹性模量，I(x)为悬臂梁在x处的惯性矩阵，w(x,t)是距离悬臂梁原点x处的截面在时间t时刻的纵向位移,M为每一段微元所受扭转力矩，F_s为微元所受的剪切力,对微元dx进行受力与力矩分析，得到如下力和力矩平衡方程：

式中，m(x)＝ρ_wingW_wingH_wing为单位长度质量,考虑到dx的二阶小量影响很小，故将其省略，得到：

梁的弯曲位移与扭转力矩M的关系得到：

将式(3)、(4)代入式(1)，整理为：

方程(5)即为悬臂梁的运动方程，基于该方程，后续将进行梁的固有振动特性分析，进而求得柔性帆板各阶模态的固有频率以及振型函数：

首先，考虑航天器运行在微重力环境下，受到重力很小，且几乎不受到其他外力的影响，因此，悬臂梁被视为自由运动模态，即悬臂梁的横向外力分布P(x,t)＝0，则基于式(5)得：

对于式(6)中的w(x,t)，采用假设模态分析法表示为：

式中，φ_n(x)为模态函数，χ_n(t)为广义坐标；

将式(7)代入式(6)，得到：

对于(8)式，左边对时间为t常数，右边对坐标x为常数，因此为保证(8)成立，其必须等于一个常数，记为Ω²，如式(9)所示，式中Ω为梁的固有振动频率：

利用分离变量方法，将(9)写成两个独立的常微分方程，为：

式中：

从而得到：

式(10)与式(11)是求解自由梁的标准方程；为了获得自由梁的固有频率Ω以及结构振型，需求解方程(10)，(11)，(11)式的通解为：

χ(t)＝A₁sinωt+A₂cosωt (12)

式(10)的通解为：

φ(x)＝De^rx (13)

将式(13)代入式(10)得：

r⁴-β⁴＝0 (14)

解得：

r_1,2＝±β r_3,4＝±iβ (15)

因此式(10)的通解表示为：

φ(x)＝D₁e^βx+D₂e^-βx+D₃e^iβx+D₄e^-iβx (16)

将式(16)转换为三角函数的形式：

φ(x)＝a_n[sin(βx)-sinh(βx)-α_n(cos(βx)-cosh(βx))] (17)

式中，

对式(17)进行归一化，从而求得系数a_n：

考虑到边界条件，由于悬臂梁一段自由，一段固定在航天器主体上，得到悬臂梁的边界条件为：

w(0,t)＝0，w′(0,t)＝0，w″(0,t)＝0，w″′(0,t)＝0

初始条件为：

将上述边界条件、初始条件代入式(7)，得模态函数的边界条件及初始条件：

φ(x)|_x＝0＝0，φ′(x)|_x＝0＝0

φ″(x)|_{x＝L_wing}＝0，φ″′(x)|_{x＝L_wing}＝0

将模态函数的边界条件、初始条件代入(17)得：

cos(βL_wing)·cosh(βL_wing)+1＝0 (19)

式(19)为超越方程，因此无法得到其精确的解，故在此用MATLAB程序来求解此方程，得到较为精确的数值解，从而获得柔性附件的固有频率及振型。

步骤S2详细步骤如下：

假设充液航天器贮液箱模型为椭圆形，令航天器质心坐标为O，储液箱液体分为两个部分，第一部分为不参与晃动的液体，其质量为m_l0，其质心与航天器质心距离为r_l0；第二部分为参与晃动的液体，将其等效为二阶弹簧质量模型并考虑其模态振动，各阶晃动模型参数为：晃动质量m_li，弹簧刚度k_li，阻尼c_li，平衡状态时各阶质量块距离航天器质心的距离为r_zi，晃动状态时各阶质量块距离航天器质心的距离为

η_i＝[η_i1 η_i2]^T，式中，η_i1为晃动质量沿着OX轴的晃动位移，η_i2为晃动质量沿着OY轴的晃动位移，i＝1,2为液体等效过程中取前两阶晃动位移；

等效弹簧质量的动量表示为：

式中，m_l为液体燃料的总质量，r_l为其到质心的距离；

在旋转坐标系下，等效液体晃动的弹簧质量动量守恒，得液体晃动动力学方程为：

步骤S3详细步骤如下：

采用基于有限体积法的VOF多相流模型模拟贮液箱内液体晃动现象，VOF模型中，对第i相流体的体积分数记为α_i，当α_i＝0时，贮液箱内不包含第i相流体；当α_i＝1时，贮液箱中充满第i相流体；当0<α_i<1则控制体处于相界面的位置，控制体中所有的体积分数之和等于1，对于航天器贮液箱内液体分布情况来说，为气液两项流：

α_l+α_g＝1 (22)

式中，α_l为液相体积分数，α_g为气相体积分数；

VOF方法的控制方程组包括：

a.连续性方程

所谓连续性方程，就是对贮液箱中任一点来说，流向该节点的流量必须等于从该节点流出的流量，即：

式中，

为速度矢量，

为散度符号；

b.动量方程

式中，ρ为密度，μ为粘度，g为重力加速度，p为压强，

为表面张力的体积力形式；

c.相函数方程

在贮液箱内，所有相体积分数综合为1，求解相函数方程是根据不同相之间容积比率的连续方程计算得到的，即：

式中，F为相函数，定义为液体燃料的体积与网格体积的比值；

由此，可以得到VOF模型中气液两相流动的密度和粘度为：

式中，ρ_g为气相密度，ρ_l为液相密度，μ_l为液相粘度，μ_g为气相粘度；

步骤S4详细步骤如下：

根据式(27)中的角速度变化曲线，

计算角加速度

并代入如下方程中：

采用Matlab中simulink方式进行数值仿真，从而求得液体晃动η。把η代入式刚液耦合方程中，从而获得采用等效力学方法建立的液体晃动模型对

产生的影响；

依据CFD仿真分析结果，将液体晃动对贮液箱产生的力矩代入(29)式，计算CFD方法获得的液体晃动模型对

产生的影响，进而比较两种方法对角加速度

产生的影响差异：

式中，M_CFD为CFD仿真获得的力矩；

因角加速度差值变化曲线变化较为随机，无法进行有效的函数拟合，因此对差值进行快速傅里叶变换处理，求取其主要影响频率，根据差值的最大振幅，拟合成一个正弦函数补偿，进而计算一个补充液体晃动力矩，代入刚液耦合方程中由此，最终获得刚液航天器动力学模型，如式(30)所示：

式中，M_l为液体晃动补充力矩；

S5、同时考虑柔性附件振动和液体晃动对刚体的耦合影响，建立大型柔性充液航天器动力学模型：

首先分析柔性航天器模型，假设p为柔性结构上的任意一个单元，其位置表示为：

r_p＝r_Op+r_d (31)

式中，柔性结构在静止状态下任意一点到航天器的质心的距离为r_Op＝r_Oo+r_op，假设柔性附件与航天器X轴重合，航天器质心到柔性附件与刚体连接处o点的距离为r_Oo＝[r_Oo0 0]^T，连接处o点到p点的距离r_op＝[x 0 0]^T，柔性附件的结构位移为r_d＝[0 0 w]^T，则p点的速度表示为：

柔性附件的动量矩H_p表示为：

因只对帆板进行横向积分，故只考虑r_Op的积分运算，因此得：

式中，

为r_p斜对称矩阵，定义

则上式写为：

因此，基于航天器刚体动量矩，刚柔航天器系统的总体动量矩表示为：

H_f＝H_p+H_m＝J_fω+h_p (36)

式中，

为刚柔航天器的转动惯量，由刚柔耦合航天器系统的动量矩守恒定理，可得：

式中，u(t)为输入力矩；基于式(32)及式(7)，柔性附件的动能表示为：

考虑航天器处于微重力环境及太阳光压小等因素，因此忽略航天器重力势能的影响与太阳光压引起的势能改变，只考虑弹性附件的应变能，柔性附件的势能表示为：

基于式(38)及(39)，建立拉格朗日函数L_p＝T_p-U_p，应用拉格朗日定理得：

式中，W_c为柔性结构的阻尼力，表示为：

式中，ε为柔性附件的阻尼系数；综上，柔性附件的振动方程表示为：

忽略上述方程的角速度二阶小量，可以写为：

把位移表达式(7)中w(x,t)代入运动学方程中，对方程两侧同时乘以

并对x积分可得：

式中，

为刚柔耦合矩阵，其中

n为第i阶模态，C_f＝diag{2ε_iΩ_i}为柔性附件的柔性矩阵，

为柔性附件的刚度矩阵，ε_i为第i阶模态的阻尼比，Ω_i为模态的固有频率；

然后分析充液航天器模型，根据等效弹簧质量的动量表达式(20)，等效弹簧质量的动量矩可以表示为：

H_s＝J_sω+h_s (46)

式中，

I为单位矩阵，

在旋转坐标系下，结合刚体航天器动力学方程，完成耦合系统模态影响分析根据动量矩守恒定理，充液航天器系统无扰动状态下的动力学方程为：

式中，

结合式(21)和式(47)，令δ_l＝-m_lir_li ^×，M_η＝m_li，C_l＝c_li，K_l＝k_li，考虑液体晃动为小幅晃动，省略上式二阶小量，并把前两阶晃动位移代入并合并为一个矩阵方程，得：

式中，

为刚液耦合矩阵，M_η＝[m_l1m_l1 m_l2 m_l2]^T为晃动液体质量矩阵，C_l＝[c_i1 c_i1 c_i2 c_i2]^T为晃动液体柔性矩阵，K_l＝[k_l1 k_l1 k_l2 k_l2]^T为晃动液体刚度矩阵，η为晃动液体模态值；

最后建立大型柔性充液航天器模型，根据刚柔航天器、刚液航天器模型分析，刚-柔-液耦合航天器动量矩表示为：

H＝H_m+H_p+H_s＝Jω+h_s+h_p (50)

根据系统总动量矩守恒得：

将式(44)，(45)，(30)，(49)代入式(51)得：

式中，

为刚柔耦合矩阵，C_f＝diag{2ε_iΩ_i}为柔性附件的柔性矩阵，

为柔性附件的刚度矩阵，式中，ε_i为第i阶模态的阻尼比，Ω_i为模态的振动频率，

为刚液耦合矩阵，M_η＝[m_l1 m_l1 m_l2 m_l2]^T为晃动液体质量矩阵，C_l＝[c_i1 c_i1 c_i2 c_i2]^T为晃动液体柔性矩阵，K_l＝[k_l1 k_l1 k_l2 k_l2]^T为晃动液体刚度矩阵，η为晃动液体模态值，M_l为液体晃动补充力矩。

本发明的特点及有益效果是：

发明对于大型柔性充液航天器动力学模型建模方法的研究和发展具有十分重要的意义。该项发明具有国际先进水平，它可以作为大型柔性充液航天器动力学模型建模方法的一种新的补充，进而有助于大型柔性充液航天器动力学模型建模技术的发展。该技术不仅有效提高了大型柔性充液航天器动力学模型准确性，同时为未来大型柔性充液航天器控制器设计奠定了良好的理论技术基础。

附图说明：

参见图1，柔性附件受力图。

参见图2，柔性附件微元受力图。

参见图3，带有椭圆形储液箱的刚液航天器模型图。

参见图4，晃动的液体部分等效模型图。

参见图5，贮箱网格划分图。

参见图6，液体随时间变化体积分数图。

参见图7，两种方法对的角加速度影响的差值图。

参见图8，刚柔液航天器简化模型。

具体实施方式

本发明的目的是提供一种大型复杂耦合航天器动力学模型建模方法，使建立的数学模型更加贴合实际卫星模型，并为航天器控制器的设计提供基础。本发明具体步骤如下：

S1、先根据太阳能帆板的结构尺寸，忽略了对航天器太阳能帆板的动力学影响微弱的小的杆件结构，建立航天器太阳帆板简化模型。考虑到太阳能帆板是一个典型的悬臂梁结构，为使分析更加直观，首先忽略刚体运动，只考虑刚体和柔性附件连接处力的作用，单独分析悬臂梁结构，采用假设模态法，将太阳能帆板简化为欧拉-伯努利悬臂梁。令P(x,t)是单位长度悬臂梁的横向外力分布，m(x)是悬臂梁的质量分布，EI(x)是悬臂梁的刚度分布，E是弹性模量，I(x)为悬臂梁在x处的惯性矩阵，w(x,t)是距离悬臂梁原点x处的截面在时间t时刻的纵向位移。M为每一段微元所受扭转力矩，F_s为微元所受的剪切力。对微元dx进行受力与力矩分析，可以得到如下力和力矩平衡方程：

式中，m(x)＝ρ_wingW_wingH_wing为单位长度质量。考虑到dx的二阶小量影响很小，故将其省略，可以得到：

梁的弯曲位移与扭转力矩M的关系可知：

将式(3)、(4)代入式(1)，可以整理为：

方程(5)即为悬臂梁的运动方程。基于该方程，后续将进行梁的固有振动特性分析，进而求得柔性帆板各阶模态的固有频率以及振型函数。

首先，考虑航天器运行在微重力环境下，受到重力很小，且几乎不受到其他外力的影响。因此，悬臂梁可以被视为自由运动模态，即悬臂梁的横向外力分布P(x,t)＝0，则基于式(5)可得：

对于式(6)中的w(x,t)，采用假设模态分析法可表示为：

式中，φ_n(x)为模态函数，χ_n(t)为广义坐标。

将式(7)代入式(6)，可以得到：

对于(8)式，左边对时间为t常数，右边对坐标x为常数，因此为保证(8)成立，其必须等于一个常数，可以记为Ω²，如式(9)所示，式中Ω为梁的固有振动频率。

利用分离变量方法，将(9)写成两个独立的常微分方程，为：

式中：

从而得到：

式(10)与式(11)是求解自由梁的标准方程。

为了获得自由梁的固有频率Ω以及结构振型，需求解方程(10)，(11)。(11)式的通解为：

χ(t)＝A₁sinωt+A₂cosωt (12)

式(10)的通解为：

φ(x)＝De^rx (13)

将式(13)代入式(10)可得：

r⁴-β⁴＝0 (14)

解得：

r_1,2＝±β r_3,4＝±iβ (15)

因此式(10)的通解可以表示为：

φ(x)＝D₁e^βx+D₂e^-βx+D₃e^iβx+D₄e^-iβx (16)

将式(16)转换为三角函数的形式：

φ(x)＝a_n[sin(βx)-sinh(βx)-α_n(cos(βx)-cosh(βx))] (17)

式中，

对式(17)进行归一化，从而求得系数a_n：

考虑到边界条件，由于悬臂梁一段自由，一段固定在航天器主体上，得到悬臂梁的边界条件为：

w(0,t)＝0，w′(0,t)＝0，w″(0,t)＝0，w″′(0,t)＝0

初始条件为：

将上述边界条件、初始条件代入式(7)，可得模态函数的边界条件及初始条件：

φ(x)|_x＝0＝0，φ′(x)|_x＝0＝0

φ″(x)|_{x＝L_wing}＝0，φ″′(x)|_{x＝L_wing}＝0

将模态函数的边界条件、初始条件代入(17)可得：

cos(βL_wing)·cosh(βL_wing)+1＝0 (19)

式(19)为超越方程，因此无法得到其精确的解，故在此用MATLAB程序来求解此方程，得到较为精确的数值解，从而获得柔性附件的固有频率及振型。

S2、应用动量矩守恒定理建立航天器系统的姿态动力学方程，参考NASA的贮液箱内液体晃动的研究成果，利用弹簧质量等效模型的前两阶模态对液体晃动部分进行等效，建立液体晃动动力学模型。

假设充液航天器贮液箱模型为椭圆形，令航天器质心坐标为O，储液箱液体分为两个部分，第一部分为不参与晃动的液体，其质量为m_l0，其质心与航天器质心距离为r_l0；第二部分为参与晃动的液体，将其等效为二阶弹簧质量模型并考虑其模态振动，各阶晃动模型参数为：晃动质量m_li，弹簧刚度k_li，阻尼c_li，平衡状态时各阶质量块距离航天器质心的距离为r_zi，晃动状态时各阶质量块距离航天器质心的距离为

η_i＝[η_i1 η_i2]^T，式中，η_i1为晃动质量沿着OX轴的晃动位移，η_i2为晃动质量沿着OY轴的晃动位移，i＝1,2为液体等效过程中取前两阶晃动位移。

等效弹簧质量的动量可以表示为：

式中，m_l为液体燃料的总质量，r_l为其到质心的距离。

在旋转坐标系下，等效液体晃动的弹簧质量动量守恒，可得液体晃动动力学方程为：

S3、流体力学分析(CFD)方法同样是研究液体晃动问题常用方法之一，可以通过模拟微重力环境下液体晃动，从而分析贮箱内流动参数的变化情况，较为深入地揭示携带的燃料对航天器姿态运动的影响。但CFD方法无法得到明确的航天器数学模型，难以开展后续控制工作，因此，本发明采用CFD与等效力学结合的方法，建立刚液航天器耦合模型。

采用CFD软件的流体体积函数(VOF)方法进行三维液体晃动仿真计算，得到椭球体贮箱内流动参数的变化规律，并进行深入分析，从而完善之前采用等效力学模型所得的液体晃动数学模型。

1、计算模型

采用基于有限体积法的VOF多相流模型模拟贮液箱内液体晃动现象。VOF模型中，对第i相流体的体积分数记为α_i，当α_i＝0时，贮液箱内不包含第i相流体；当α_i＝1时，贮液箱中充满第i相流体；当0<α_i<1则控制体处于相界面的位置，控制体中所有的体积分数之和等于1，对于航天器贮液箱内液体分布情况来说，一般为气液两项流：

α_l+α_g＝1 (22)

式中，α_l为液相体积分数，α_g为气相体积分数。

VOF方法的控制方程组包括：

a.连续性方程

所谓连续性方程，就是对贮液箱中任一点来说，流向该节点的流量必须等于从该节点流出的流量，即：

式中，

为速度矢量，

为散度符号。

b.动量方程

式中，ρ为密度，μ为粘度，g为重力加速度，p为压强，

为表面张力的体积力形式。

c.相函数方程

在贮液箱内，所有相体积分数综合为1，求解相函数方程是根据不同相之间容积比率的连续方程计算得到的，即：

式中，F为相函数，定义为液体燃料的体积与网格体积的比值。

由此，可以得到VOF模型中气液两相流动的密度和粘度为：

式中，ρ_g为气相密度，ρ_l为液相密度，μ_l为液相粘度，μ_g为气相粘度。

2.数值仿真

给定贮液箱x，y轴方向直径为1m,z轴方向直径为1.5m。贮液箱壁面处具有粘性阻尼，粘性阻尼是影响液体晃动主要因素之一，并且数值计算中采用了壁面函数。对贮箱进行网格划分。

采用有限体积法对结构化网格进行空间上离散，采用二阶全隐式格式进行时间上离散。数值计算中采用标准k-ε湍流模型，近壁区采用了标准壁面函数。充液液体为甲基肼(MMH)，充液比为50％，贮箱中气体为空气，压力为一个大气压，模型考虑表面张力的影响，贮箱所处重力环境为g＝0.2m/s²，液体在原点处绕着x轴负方向运动(右手定则)。其角速度变式(27)所示，假设航天器24秒进行60度的旋转。下表给出了流体的物理参数。

表1.流体的物理参数

S4、在航天器进行大机动姿态运动时，采用等效力学方法建立的液体晃动模型不能完全表达液体晃动的规律，而单纯通过CFD方法可以模拟微重力环境下液体晃动，但无法建立相关数学模型，难以进行控制器设计。因此，综合考虑两种液体建模方法，依据CFD分析软件获得的液体晃动内部参数结果对等效力学方法建立的液体晃动模型进行补充，以获得更贴近设计晃动规律的液体晃动数学模型。根据式(27)中的角速度变化曲线，计算角加速度

并代入如下方程中。

采用Matlab中simulink方式进行数值仿真，从而求得液体晃动η。把η代入式刚液耦合方程中，从而获得采用等效力学方法建立的液体晃动模型对

产生的影响。

依据CFD仿真分析结果，将液体晃动对贮液箱产生的力矩代入(29)式，计算CFD方法获得的液体晃动模型对

产生的影响，进而比较两种方法对角加速度

产生的影响差异。

式中，M_CFD为CFD仿真获得的力矩。

因角加速度差值变化曲线变化较为随机，无法进行有效的函数拟合，因此对差值进行快速傅里叶变换处理，求取其主要影响频率，根据差值的最大振幅，拟合成一个正弦函数补偿，进而计算一个补充液体晃动力矩，代入刚液耦合方程中由此，最终获得刚液航天器动力学模型，如式(30)所示：

式中，M_l为液体晃动补充力矩。

S5、在以上基础上，需要同时考虑柔性附件振动和液体晃动对刚体的耦合影响，建立大型柔性充液航天器动力学模型。

首先分析柔性航天器模型，假设p为柔性结构上的任意一个单元，其位置可表示为：

r_p＝r_Op+r_d (31)

式中，柔性结构在静止状态下任意一点到航天器的质心的距离为r_Op＝r_Oo+r_op，假设柔性附件与航天器X轴重合，航天器质心到柔性附件与刚体连接处o点的距离为r_Oo＝[r_Oo0 0]^T，连接处o点到p点的距离r_op＝[x 0 0]^T，柔性附件的结构位移为r_d＝[0 0 w]^T，则p点的速度可以表示为：

柔性附件的动量矩H_p可以表示为：

因只对帆板进行横向积分，故只考虑r_Op的积分运算，因此可得：

式中，

为r_p斜对称矩阵，定义

则上式可以写为：

因此，基于航天器刚体动量矩，刚柔航天器系统的总体动量矩可以表示为：

H_f＝H_p+H_m＝J_fω+h_p (36)

式中，

为刚柔航天器的转动惯量。

由刚柔耦合航天器系统的动量矩守恒定理，可得：

式中，u(t)为输入力矩。

基于式(32)及式(7)，柔性附件的动能可以表示为：

考虑航天器处于微重力环境及太阳光压小等因素，因此忽略航天器重力势能的影响与太阳光压引起的势能改变，只考虑弹性附件的应变能，柔性附件的势能可以表示为：

基于式(38)及(39)，建立拉格朗日函数L_p＝T_p-U_p，应用拉格朗日定理可得：

式中，W_c为柔性结构的阻尼力，表示为：

式中，ε为柔性附件的阻尼系数。

综上，柔性附件的振动方程可以表示为：

忽略上述方程的角速度二阶小量，可以写为：

把位移表达式(7)中w(x,t)代入运动学方程中，对方程两侧同时乘以

并对x积分可得：

式中，

为刚柔耦合矩阵，其中

n为第i阶模态。C_f＝diag{2ε_iΩ_i}为柔性附件的柔性矩阵，

为柔性附件的刚度矩阵，ε_i为第i阶模态的阻尼比，Ω_i为模态的固有频率。

然后分析充液航天器模型，根据等效弹簧质量的动量表达式(20)，等效弹簧质量的动量矩可以表示为：

H_s＝J_sω+h_s (46)

式中，

I为单位矩阵，

在旋转坐标系下，结合刚体航天器动力学方程，完成耦合系统模态影响分析根据动量矩守恒定理，充液航天器系统无扰动状态下的动力学方程为：

式中，

结合式(21)和式(47)，令δ_l＝-m_lir_li ^×，M_η＝m_li，C_l＝c_li，K_l＝k_li，考虑液体晃动为小幅晃动，省略上式二阶小量，并把前两阶晃动位移代入并合并为一个矩阵方程，可得：

式中，

为刚液耦合矩阵，M_η＝[m_l1 m_l1 m_l2 m_l2]^T为晃动液体质量矩阵，C_l＝[c_i1 c_i1 c_i2 c_i2]^T为晃动液体柔性矩阵，K_l＝[k_l1 k_l1 k_l2 k_l2]^T为晃动液体刚度矩阵，η为晃动液体模态值。

最后建立大型柔性充液航天器模型，根据刚柔航天器、刚液航天器模型分析，刚-柔-液耦合航天器动量矩可以表示为：

H＝H_m+H_p+H_s＝Jω+h_s+h_p (50)

根据系统总动量矩守恒可得：

将式(44)，(45)，(30)，(49)代入式(51)可得：

式中，

为刚柔耦合矩阵，，C_f＝diag{2ε_iΩ_i}为柔性附件的柔性矩阵，

为柔性附件的刚度矩阵，式中，ε_i为第i阶模态的阻尼比，Ω_i为模态的振动频率，

为刚液耦合矩阵，M_η＝[m_l1 m_l1 m_l2 m_l2]^T为晃动液体质量矩阵，C_l＝[c_i1 c_i1 c_i2 c_i2]^T为晃动液体柔性矩阵，K_l＝[k_l1 k_l1 k_l2 k_l2]^T为晃动液体刚度矩阵，η为晃动液体模态值，M_l为液体晃动补充力矩。

结合附图对本发明作进一步详述。

参见图1，柔性附件受力图。柔性附件简化为欧拉-伯努利悬臂梁，其中P(x,t)是单位长度悬臂梁的横向外力分布，m(x)是悬臂梁的质量分布，EI(x)是悬臂梁的刚度分布，E是弹性模量，I(x)为悬臂梁在x处的惯性矩阵，w(x,t)是距离悬臂梁原点x处的截面在时间t时刻的纵向位移。

参见图2，柔性附件微元受力图。图中M为每一段微元所受扭转力矩，F_s为微元所受的剪切力。对微元dx进行受力与力矩分析，可以得到力和力矩平衡方程。

参见图3，带有椭圆形储液箱的刚液航天器模型。假设充液航天器贮液箱模型为椭圆形，根据其形状，可以知道航天器质心坐标为O，储液箱液体分为两个部分，第一部分为不参与晃动的液体，其质量为m_l0，其质心与航天器质心距离为r_l0

参见图4，晃动的液体部分等效模型。该图为晃动液体等效图，各阶晃动模型参数为：晃动质量m_li，弹簧刚度k_li，阻尼c_li，平衡状态时各阶质量块距离航天器质心的距离为r_zi，晃动状态时各阶质量块距离航天器质心的距离为

η_i＝[η_i1 η_i2]^T，式中，η_i1为晃动质量沿着OX轴的晃动位移，η_i2为晃动质量沿着OY轴的晃动位移，i＝1,2为液体等效过程中取前两阶晃动位移。

参见图5，贮箱网格划分图。在CFD方法中对贮液箱网格划分是常用的手段，计算叠加贮液箱网格每个小微元的受力情况从而获得整个贮液箱的受力情况。贮液箱x，y轴方向直径为1m,z轴方向直径为1.5m。贮液箱壁面处具有粘性阻尼，粘性阻尼是影响液体晃动主要因素之一，并且数值计算中采用了壁面函数。

参见图6，液体随时间变化体积分数。根据仿真结果可知，液体在前10秒内对贮箱的作用力矩较大，在10秒之后，液体逐渐与贮箱共同运动并且绕贮箱壁均匀流动，从而对贮箱的力矩逐渐减小。不同时刻在x＝0平面剖面的液体体积分数如图6所示。从图中可以看出，当时间小于10秒时，液体晃动剧烈，对贮液箱壁产生较大的冲击作用，当时间大于10秒时，液体绕贮液箱壁近似均匀流动，从而对贮箱壁产生的冲击作用较小，当时间大于40秒时，液体慢慢下沉于贮箱底部。

参见图7，两种方法对的角加速度影响的差值。角加速度差值变化曲线变化较为随机，无法进行有效的函数拟合，因此对差值进行快速傅里叶变换处理，求取其主要影响频率，根据差值的最大振幅，拟合成一个正弦函数补偿，进而计算一个补充液体晃动力矩。

参见图8，刚柔液航天器简化模型。大型柔性充液航天器需要同时考虑柔性附件振动和液体晃动对刚体的耦合影响，对刚-柔-液耦合复杂航天器进行模型分析，模型简图8所示。

下面给出具体实例：

假设单块太阳能帆板长度为L_wing＝7.448m，太阳能帆板宽度为W_wing＝1.8m，太阳能帆板厚度为H_wing＝0.02m，太阳能帆板弹性模量E＝1.93GPa，航天器总共携带550kg推进剂液体。在姿态快速机动过程中仅剩余半箱燃料，不参与晃动的液体质量为m_l0＝0kg，其质心与航天器质心距离为r_l0＝1.137m，前二阶晃动液体质量为m_l1＝200kg，m_l2＝25kg，晃动液体质心与航天器的质心z方向距离为r_z1＝1.127m，r_z2＝0.994m，二阶晃动液体弹簧刚度为k_l1＝55.21N/m k_l2＝7.27N/m，二阶晃动液体阻尼为c_i1＝3.334N·s/m，c_i2＝0.237N·s/m。将上述参数代入式(2.58)和式(2.59)，获得耦合复杂航天器模型参数如下：

为刚液耦合矩阵，M_η＝diag([200 200 25 25])为晃动液体的质量矩阵，

C_l＝diag([3.334 3.334 0.237 0.237])为晃动液体的柔性矩阵，K_l＝diag([55.21 55.21 7.27 7.27])为晃动液体的刚度矩阵。

Claims (4)

1.一种大型复杂耦合航天器动力学模型建模方法，其特征是，步骤如下：

S1、先根据太阳能帆板的结构尺寸，忽略对航天器太阳能帆板的动力学影响微弱的小的杆件结构，建立航天器太阳帆板简化模型；

S2、应用动量矩守恒定理建立航天器系统的姿态动力学方程，利用弹簧质量等效模型的前两阶模态对液体晃动部分进行等效，建立液体晃动动力学模型；

S3、采用CFD与等效力学结合的方法，建立刚液航天器耦合模型，具体是采用CFD软件的流体体积函数方法进行三维液体晃动仿真计算，得到椭球体贮箱内流动参数的变化规律，并进行深入分析，从而完善S2采用等效力学模型所得的液体晃动数学模型；

S4、综合考虑S2、S3两种液体建模方法，依据CFD分析软件获得的液体晃动内部参数结果，对等效力学方法建立的液体晃动模型进行补充，以获得更贴近设计晃动规律的液体晃动数学模型；

S5、需要同时考虑S1中柔性附件振动和S2、S3中液体晃动对刚体的耦合影响，建立大型柔性充液航天器动力学模型；

S6、最后根据刚柔航天器、刚液航天器模型分析，建立大型柔性充液航天器模型。

2.如权利要求1所述的大型复杂耦合航天器动力学模型建模方法，其特征是，步骤S1详细步骤如下：

考虑到太阳能帆板是一个典型的悬臂梁结构，为使分析更加直观，首先忽略刚体运动，只考虑刚体和柔性附件连接处力的作用，单独分析悬臂梁结构，采用假设模态法，将太阳能帆板简化为欧拉-伯努利悬臂梁,令P(x,t)是单位长度悬臂梁的横向外力分布，m(x)是悬臂梁的质量分布，EI(x)是悬臂梁的刚度分布，E是弹性模量，I(x)为悬臂梁在x处的惯性矩阵，w(x,t)是距离悬臂梁原点x处的截面在时间t时刻的纵向位移,M为每一段微元所受扭转力矩，F_s为微元所受的剪切力,对微元dx进行受力与力矩分析，得到如下力和力矩平衡方程：

式中，m(x)＝ρ_wingW_wingH_wing为单位长度质量,考虑到dx的二阶小量影响很小，故将其省略，得到：

梁的弯曲位移与扭转力矩M的关系得到：

将式(3)、(4)代入式(1)，整理为：

方程(5)即为悬臂梁的运动方程，基于该方程，后续将进行梁的固有振动特性分析，进而求得柔性帆板各阶模态的固有频率以及振型函数：

首先，考虑航天器运行在微重力环境下，受到重力很小，且几乎不受到其他外力的影响，因此，悬臂梁被视为自由运动模态，即悬臂梁的横向外力分布P(x,t)＝0，则基于式(5)得：

对于式(6)中的w(x,t)，采用假设模态分析法表示为：

式中，φ_n(x)为模态函数，χ_n(t)为广义坐标；

将式(7)代入式(6)，得到：

对于(8)式，左边对时间为t常数，右边对坐标x为常数，因此为保证(8)成立，其必须等于一个常数，记为Ω²，如式(9)所示，式中Ω为梁的固有振动频率：

利用分离变量方法，将(9)写成两个独立的常微分方程，为：

式中：

从而得到：

式(10)与式(11)是求解自由梁的标准方程；为了获得自由梁的固有频率Ω以及结构振型，需求解方程(10)，(11)，(11)式的通解为：

χ(t)＝A₁sinωt+A₂cosωt (12)

式(10)的通解为：

φ(x)＝De^rx (13)

将式(13)代入式(10)得：

r⁴-β⁴＝0 (14)

解得：

r_1,2＝±βr_3,4＝±iβ (15)

因此式(10)的通解表示为：

将式(16)转换为三角函数的形式：

φ(x)＝a_n[sin(βx)-sinh(βx)-α_n(cos(βx)-cosh(βx))] (17)

式中，

对式(17)进行归一化，从而求得系数a_n：

考虑到边界条件，由于悬臂梁一段自由，一段固定在航天器主体上，得到悬臂梁的边界条件为：

w(0,t)＝0，w′(0,t)＝0，w″(0,t)＝0，w″′(0,t)＝0

初始条件为：

w(x,t)|_t＝0＝w(x,0)，

将上述边界条件、初始条件代入式(7)，得模态函数的边界条件及初始条件：

φ(x)|_x＝0＝0，φ′(x)|_x＝0＝0

φ″(x)|_{x＝L_wing}＝0，φ″′(x)|_{x＝L_wing}＝0

将模态函数的边界条件、初始条件代入(17)得：

cos(βL_wing)·cosh(βL_wing)+1＝0 (19)

式(19)为超越方程，因此无法得到其精确的解，故在此用MATLAB程序来求解此方程，得到较为精确的数值解，从而获得柔性附件的固有频率及振型。

3.如权利要求1所述的大型复杂耦合航天器动力学模型建模方法，其特征是，步骤S2详细步骤如下：

假设充液航天器贮液箱模型为椭圆形，令航天器质心坐标为O，储液箱液体分为两个部分，第一部分为不参与晃动的液体，其质量为m_l0，其质心与航天器质心距离为r_l0；第二部分为参与晃动的液体，将其等效为二阶弹簧质量模型并考虑其模态振动，各阶晃动模型参数为：晃动质量m_li，弹簧刚度k_li，阻尼c_li，平衡状态时各阶质量块距离航天器质心的距离为r_zi，晃动状态时各阶质量块距离航天器质心的距离为

η_i＝[η_i1 η_i2]^T，式中，η_i1为晃动质量沿着OX轴的晃动位移，η_i2为晃动质量沿着OY轴的晃动位移，i＝1,2为液体等效过程中取前两阶晃动位移；

等效弹簧质量的动量表示为：

式中，m_l为液体燃料的总质量，r_l为其到质心的距离；

在旋转坐标系下，等效液体晃动的弹簧质量动量守恒，得液体晃动动力学方程为：

4.如权利要求1所述的大型复杂耦合航天器动力学模型建模方法，其特征是，步骤S3详细步骤如下：

采用基于有限体积法的VOF多相流模型模拟贮液箱内液体晃动现象，VOF模型中，对第i相流体的体积分数记为α_i，当α_i＝0时，贮液箱内不包含第i相流体；当α_i＝1时，贮液箱中充满第i相流体；当0<α_i<1则控制体处于相界面的位置，控制体中所有的体积分数之和等于1，对于航天器贮液箱内液体分布情况来说，为气液两项流：

α_l+α_g＝1 (22)

式中，α_l为液相体积分数，α_g为气相体积分数；

VOF方法的控制方程组包括：

a.连续性方程

所谓连续性方程，就是对贮液箱中任一点来说，流向该节点的流量必须等于从该节点流出的流量，即：

式中，

为速度矢量，

为散度符号；

b.动量方程

式中，ρ为密度，μ为粘度，g为重力加速度，p为压强，

为表面张力的体积力形式；

c.相函数方程

在贮液箱内，所有相体积分数综合为1，求解相函数方程是根据不同相之间容积比率的连续方程计算得到的，即：

式中，F为相函数，定义为液体燃料的体积与网格体积的比值；

由此，可以得到VOF模型中气液两相流动的密度和粘度为：

式中，ρ_g为气相密度，ρ_l为液相密度，μ_l为液相粘度，μ_g为气相粘度；

步骤S4详细步骤如下：

根据式(27)中的角速度变化曲线，

计算角加速度

并代入如下方程中：

采用Matlab中simulink方式进行数值仿真，从而求得液体晃动η。把η代入式刚液耦合方程中，从而获得采用等效力学方法建立的液体晃动模型对

产生的影响；

依据CFD仿真分析结果，将液体晃动对贮液箱产生的力矩代入(29)式，计算CFD方法获得的液体晃动模型对

产生的影响，进而比较两种方法对角加速度

产生的影响差异：

式中，M_CFD为CFD仿真获得的力矩；

因角加速度差值变化曲线变化较为随机，无法进行有效的函数拟合，因此对差值进行快速傅里叶变换处理，求取其主要影响频率，根据差值的最大振幅，拟合成一个正弦函数补偿，进而计算一个补充液体晃动力矩，代入刚液耦合方程中由此，最终获得刚液航天器动力学模型，如式(30)所示：

式中，M_l为液体晃动补充力矩；

S5、同时考虑柔性附件振动和液体晃动对刚体的耦合影响，建立大型柔性充液航天器动力学模型：

首先分析柔性航天器模型，假设p为柔性结构上的任意一个单元，其位置表示为：

r_p＝r_Op+r_d (31)

式中，柔性结构在静止状态下任意一点到航天器的质心的距离为r_Op＝r_Oo+r_op，假设柔性附件与航天器X轴重合，航天器质心到柔性附件与刚体连接处o点的距离为r_Oo＝[r_Oo 00]^T，连接处o点到p点的距离r_op＝[x 0 0]^T，柔性附件的结构位移为r_d＝[0 0 w]^T，则p点的速度表示为：

柔性附件的动量矩H_p表示为：

因只对帆板进行横向积分，故只考虑r_Op的积分运算，因此得：

式中，

为r_p斜对称矩阵，定义

则上式写为：

因此，基于航天器刚体动量矩，刚柔航天器系统的总体动量矩表示为：

H_f＝H_p+H_m＝J_fω+h_p (36)

式中，

为刚柔航天器的转动惯量，由刚柔耦合航天器系统的动量矩守恒定理，可得：

式中，u(t)为输入力矩；基于式(32)及式(7)，柔性附件的动能表示为：

考虑航天器处于微重力环境及太阳光压小等因素，因此忽略航天器重力势能的影响与太阳光压引起的势能改变，只考虑弹性附件的应变能，柔性附件的势能表示为：

基于式(38)及(39)，建立拉格朗日函数L_p＝T_p-U_p，应用拉格朗日定理得：

式中，W_c为柔性结构的阻尼力，表示为：

式中，ε为柔性附件的阻尼系数；综上，柔性附件的振动方程表示为：

忽略上述方程的角速度二阶小量，可以写为：

把位移表达式(7)中w(x,t)代入运动学方程中，对方程两侧同时乘以

并对x积分可得：

式中，

为刚柔耦合矩阵，其中

n为第i阶模态，C_f＝diag{2ε_iΩ_i}为柔性附件的柔性矩阵，

为柔性附件的刚度矩阵，ε_i为第i阶模态的阻尼比，Ω_i为模态的固有频率；

然后分析充液航天器模型，根据等效弹簧质量的动量表达式(20)，等效弹簧质量的动量矩可以表示为：

H_s＝J_sω+h_s (46)

式中，

I为单位矩阵，

在旋转坐标系下，结合刚体航天器动力学方程，完成耦合系统模态影响分析根据动量矩守恒定理，充液航天器系统无扰动状态下的动力学方程为：

式中，

结合式(21)和式(47)，令

M_η＝m_li，C_l＝c_li，K_l＝k_li，考虑液体晃动为小幅晃动，省略上式二阶小量，并把前两阶晃动位移代入并合并为一个矩阵方程，得：

式中，

为刚液耦合矩阵，M_η＝[m_l1 m_l1 m_l2 m_l2]^T为晃动液体质量矩阵，C_l＝[c_i1 c_i1 c_i2 c_i2]^T为晃动液体柔性矩阵，K_l＝[k_l1 k_l1 k_l2 k_l2]^T为晃动液体刚度矩阵，η为晃动液体模态值；

最后建立大型柔性充液航天器模型，根据刚柔航天器、刚液航天器模型分析，刚-柔-液耦合航天器动量矩表示为：

H＝H_m+H_p+H_s＝Jω+h_s+h_p (50)

根据系统总动量矩守恒得：

将式(44)，(45)，(30)，(49)代入式(51)得：

式中，

为刚柔耦合矩阵，C_f＝diag{2ε_iΩ_i}为柔性附件的柔性矩阵，

为柔性附件的刚度矩阵，式中，ε_i为第i阶模态的阻尼比，Ω_i为模态的振动频率，

为刚液耦合矩阵，M_η＝[m_l1 m_l1 m_l2 m_l2]^T为晃动液体质量矩阵，C_l＝[c_i1 c_i1 c_i2 c_i2]^T为晃动液体柔性矩阵，K_l＝[k_l1 k_l1 k_l2 k_l2]^T为晃动液体刚度矩阵，η为晃动液体模态值，M_l为液体晃动补充力矩。

Images