闫玉龙，申云峰，岳宝增

(1. 太原理工大学数学学院， 太原 030024； 2. 郑州工业应用技术学院信息工程学院， 郑州 451150； 3. 北京理工大学宇航学院， 北京 100081)

0 引 言

随着现代航天事业的快速发展，航天器需要装配柔性附件和一定数量的液体推进剂以完成复杂的航天任务。例如，为实现航天器的长期在轨运行，航天器通常安装大型的太阳能帆板以提供必要的能量来源。因此，板类柔性附件是航天器系统的重要组成部分。航天器在轨过程中柔性附件会不可避免地产生弹性振动，对航天器的姿态和轨道造成严重的影响。此外，在轨航天器往往处于微重力状态，在该种情况下表面张力对液体推进剂的动力学行为起着不可忽略的影响，毛细效应使得储腔液体静液面发生大幅弯曲呈现新月形[1-3]。同时，微重力环境使得液体晃动的固有频率显著降低，与柔性附件的振动频率相接近，容易产生共振现象，对航天器的动力学行为造成显著影响。因此，对含有板类柔性附件的充液航天器动力学研究具有十分重要的理论意义和工程价值。

在储腔容积相同的情况下，球形储腔的表面积最小，能够有效地降低航天器的质量。众多学者对含有球形充液储腔航天器系统的液体晃动以及刚-液耦合动力学和控制问题进行深入研究[4-8]。由于该类曲壁轴对称储腔的壁面不平行于储液腔的对称轴，难以给出势函数的解析表达式，因此对于该类储腔液体晃动的解析研究是较为困难的。通过变分法和最小势能原理，能够获得凸轴对称储液腔静液面的控制微分方程以及稳态构形[9]。基于曲壁轴对称储腔的几何特点，文献[10]采用特殊坐标系得到凸轴对称储腔液体常重力下的晃动频率，且该方法适用于任意充液深度。Utsumi[11-14]通过引入特殊的球坐标系对于凸轴对称储腔的液体晃动动力学进行研究，采用Gauss超几何级数获得晃动势函数和液面波高的解析形式，从而得到球形储腔中的液体在横向和轴向激励下的晃动频率、波高、晃动力和力矩等晃动参数，该种方法可用于研究常重力和不同微重力环境的工况下液体晃动行为。对含有球形储腔的充液航天器，可采用Hamilton-Ostrogradskiy变分原理推导动力学方程，得到晃动液体受到任意方向外激励的等效力学模型[15]。对含有多个轴对称充液储腔的航天器，可通过复合控制方法对航天器的姿态和轨道机动控制进行研究。研究结果表明，若在设计姿态和轨道控制器时没有充分考虑液体燃料晃动的影响，则系统将会出现刚-液-控耦合问题并导致航天器姿态不稳定[16-18]。

Kane等[19]讨论了含有柔性附件刚体的平动和转动相互耦合的情况，首次提出由于柔性附件的几何非线性变形与刚体运动相互耦合导致的动力刚化这一概念。对由中心体和柔性附件组成的柔性航天器，可通过混合坐标法，采用准坐标形式的Lagrange方程，建立柔性航天器的混合坐标动力学方程[20]。文献[21-23]研究由刚性平台和若干相对于刚体发生转动的柔性附件组成的航天器系统，采用基于准坐标下的Lagrange方程构建系统的动力学模型。通过空间离散化和截断，将系统的偏微分方程转换为便于进行数值仿真和计算的非线性离散状态方程。对于中心刚体-旋转梁的系统，梁的振动频率和模态会随着刚体的转动惯量和梁的转速而发生改变[24]。文献[25]对附着在空间运动体上柔性悬臂梁的动力刚化问题进行系统研究，考虑梁横向二维振动和纵向一维振动的情况，采用假设模态法对柔性梁的运动进行离散，研究横向变形和纵向变形的耦合效应。此外，航天器柔性附件在瞬态热效应下产生热变形，对航天器的姿态产生扰动[26-27]；在航天器编队飞行问题中，也存在柔性航天器和刚性航天器的相对姿态动力学和主动控制问题，需要设计姿态反馈控制器实现航天器编队飞行[28]。由此可见，航天器动力学和控制的研究中广泛存在着刚-柔耦合动力学问题。

柔性航天器通常含有液体燃料，液体燃料的晃动与刚体运动、柔性附件振动相互耦合，会产生复杂的动力学现象和控制问题[29-30]。对于该类复杂充液柔性航天器系统，通常采用多体动力学的研究方法进行分析。文献[31-32]将晃动液体和柔性附件分别等效为两个质量-弹簧模型和悬臂板模型，并将重力梯度力矩作为航天器的主要扰动力矩，设计PID控制器进行航天器姿态镇定，并对耦合情况进行分析。吕敬等[33]对带弹性附件以及矩形贮箱的刚-液-弹耦合系统进行研究，得到外力矩作用下俯仰运动的耦合动力学模型，分析液体和弹性体对航天器主刚体运动的影响。

板类柔性附件具有较大的长宽尺寸，但是厚度较小。基于航天器板类柔性附件的几何尺寸和材料特点，可采用薄板模型对板类柔性附件进行建模和动力学行为研究。本文采用Kirchhoff-Love薄板模型对太阳能帆板的变形进行研究，通过变分原理推导微重力环境下凸轴对称充液储腔的液体晃动控制方程，获得刚-液-柔耦合航天器系统的动力学模型，研究了耦合航天器系统各部分的耦合效应，并讨论柔性附件的装配位置对耦合航天器系统动力学行为的影响，以期为航天器总体设计提供参考依据。

1 耦合航天器系统的动力学模型

考虑如图1所示的刚-液-柔耦合航天器系统，由航天器主刚性平台、板类柔性附件和轴对称充液储腔组成。在研究中做出如下基本假设：1)储液腔的液体运动是无黏、无旋、无阻尼、不可压缩；2)储液腔为刚性；3)液体表面位移相对于静平衡液面足够小，在研究过程中可用线性化理论进行处理；4)采用Kirchhoff-Love板模型分析板类柔性附件的变形，板弯曲时中面各点只有垂直中面的位移w，没有平行中面的位移。

图1 刚-液-柔耦合航天器示意图Fig.1 Diagram of rigid-liquid-flex coupled spacecraft system

1.1 柔性航天器动力学

考虑如图1所示的耦合航天器系统。航天器刚体平台质量为m0，为研究航天器位置和姿态运动，引入如下坐标系统：OeXeYeZe为惯性坐标系，以航天器刚性平台的质心O为原点建立航天器本体坐标系Oxyz，其中Ox,Oy,Oz轴为惯性主轴，航天器主刚体的转动惯量为J0=diag(J11,J22,J33)。令θ=[θx,θy,θz]T为本体坐标系相对于惯性系的Euler角。为方便起见，定义坐标变换矩阵按照Ox→Oy→Oz的顺序，从惯性坐标系OeXeYeZe到本体坐标系Oxyz的转换矩阵为

(1)

其中，sk=sinθk,ck=cosθk,k=x,y,z(下同)，rOeO为O相对于原点Oe的矢径，本体系在惯性系下的角速度为ω=[ω1,ω2,ω3]T，原点O在惯性系下的平动速度为v=[v1,v2,v3]T，则

(2)

式中:

(3)

板类柔性附件的一边与航天器主刚体相固连，其余三边为自由端，该类柔性附件也称为悬臂板。板的中性面位于航天器本体系的Oxy平面。如图2所示，点B位于板中性面的边缘，该点在航天器本体系的矢量为rOB=[xb,yb,0]T。板沿Ox,Oy,Oz轴的长度分别为a,b,h，即板的长度、宽度和厚度，板的密度为ρ，弹性模量为E，泊松比为μ，质量为ma=ρabh。

图2 板类柔性附件示意图Fig.2 Diagram of plate-type appendage

考虑板类柔性附件中性面上的任意点P，则有rBP=[x,y,0]T，其中x∈[0,a],y∈[0,b]，从而点P在本体系Oxyz的矢径为

rOP=rOB+rBP=[xb+x,yb+y,0]T

(4)

薄板单位元的惯性力为-ρhap3，根据D’Alembert原理，可得

(5)

其中,cz为板的结构阻尼系数。下面采用假设模态法求解以上偏微分方程。由于悬臂板振动的准确模态函数是未知的，因此需要寻找满足几何边界条件的特征函数，参考文献[34]，可将悬臂板模态函数表示为:

W(x,y)q(t)

(6)

式中:φm(x)为悬臂梁的第m阶模态函数，ψn(y)为自由-自由梁的第n阶模态函数，其中ψ1,ψ2分别对应于自由-自由梁的平动模态、转动模态。模态函数的表达式为

(7)

其中，

(8)

式中:当n为奇数时，对应于板的对称模态；当n为偶数时，对应于板的反对称模态。将式(6)代入式(5)，两边乘以WT，并对全板积分，可得

(9)

式中:

分别表示板的质量和刚度矩阵，方程(9)称为板的振动方程。根据板上任意点加速度的表达式(4)，并对板进行积分，可得到板的惯性力F，以及对本体系原点惯性力矩M=[Mx,My,Mz]T的解析形式

(10)

刚体航天器的Lagrange函数为

(11)

对主刚体采用准坐标下的Lagrange方程

(12)

其中,F0,M0分别为航天器刚体受到的非保守力和非保守力矩。结合方程(9)和(12)写成矩阵形式，可得到含有柔性附件的航天器的动力学方程

(13)

式(13)中质量矩阵的非对角部分M2表明航天器刚性平台和柔性附件的耦合效应。由于篇幅限制，质量矩阵和广义外力的具体形式不再给出。耦合航天器系统的动力学方程(13)和刚体运动学方程(2)构成含有太阳能帆板的航天器系统控制方程。

1.2 曲壁轴对称储液腔的液体晃动动力学

在轨航天器往往处于微重力环境，液体表面张力在储腔的液体晃动中将起着不可忽略的作用。当Bond数较小时储液腔的静液面会产生严重弯曲，此外，曲壁轴对称储腔具有壁面不平行于储腔对称轴的特点，储腔充液深度的变化导致静液面的构形发生改变，对储腔液体晃动动力学行为造成影响。下面考虑如图3所示的曲壁轴对称充液储腔。

图3 轴对称储腔示意图Fig.3 Diagram of axisymmetrical tanks

如图3所示，储腔中的虚线和实线分别表示静液面和扰动液面，分别用M,F表示；静液面与储液腔壁面相交的部分称为接触线。通过引入球坐标系统O1R1θ1φ1对储腔液体晃动动力学进行研究，其中原点O1为圆锥的顶点，圆锥的侧面与储腔壁面在接触线处相切。液体表面位移ζ的运动方向为R1方向。通过球坐标系，静液面、扰动液面、储腔壁面分别表示为：

以球腔的底部顶点o1为原点，建立直角坐标系o1x1y1z1，其中o1z1轴为储腔的对称轴，方向向上，平面o1x1z1为球腔垂直方向剖面，即球坐标系中的O1R1θ1平面。接触线上的任意点在直角坐标系o1x1y1z1的坐标记为(xc,yc,zc)。由于储腔是轴对称的，接触线任意点到o1x1y1平面的距离都相同，即接触线上点的zc值不发生变化。令R0为储腔垂直方向的半径，若zc>R0，球坐标O1的原点位于储腔上方；若zc

对于耦合航天器系统，考虑储腔为刚性，且液体为小幅线性晃动。储腔液体的速度势函数φ可通过叠加原理分解为相对晃动势函数φr和牵连晃动势函数φe。储腔中任意点p的速度为vp=v+ω×·(rOo1+ro1p)，则牵连晃动势函数的表达式为

φe=(v1+ω2z0-ω3y0)R1sinθ1cosφ1+(v2+

ω3x0-ω1z0)R1sinθ1sinφ1+(v3+ω1y0-

ω2x0)(h1-εR1cosθ1)

(14)

对于储腔液体晃动动力学，可通过变分原理[35]得到相应的控制方程

(15)

其中,pl,pg分别为液体和气体压强，σ,σ1,σ2分别为汽-液分界面F、固-液分界面W1、汽-固分界面W2的单位面积表面能，可推导得到液体相对晃动势函数φr以及自由液面晃动波高函数ζ[11]

(16)

式中特征函数Θmk的解析形式为

Θmk(θ1)=sinmθ1·G(m-αmk,αmk+m+1,

m+1,(1-cosθ1)/2)

其中,G(α,β,γ,x)为Gauss超几何级数，当|x|<1时，级数为收敛的。qm,pm为液体晃动的模态系数，amk,bmk,cmk是通过晃动频率决定的系数，la,lb是用于提升收敛速率的系数，αmk是由储腔边界条件确定的常数。对式(16)关于模态系数qm,pm进行变分，代入液体晃动的控制方程(15)，根据变分的独立性，可得到模态系数的控制方程

(17)

1.3 耦合航天器系统动力学

储液腔的液体晃动产生的晃动力和晃动力矩会对耦合航天器系统的动力学行为造成影响，同时耦合航天器系统的位置和姿态的变化也会影响液体的晃动行为。两者相互作用，构成耦合航天器动力学系统。曲壁轴对称储腔的液体晃动力Fslosh和晃动力矩Mslosh是由于液体的动压强以及表面张力动接触线的非平衡拖动产生的，由于篇幅限制，晃动力和晃动力矩的具体形式不再给出。将晃动力和晃动力矩代入方程(13)，可得

(18)

结合耦合航天器系统的运动学方程(2)和动力学方程(18)，以及液体晃动模态系数的控制方程(17)，构成含有板类柔性附件和轴对称储腔的耦合航天器系统的控制方程。

2 仿真校验

下面通过数值仿真，研究耦合航天器模型的动力学行为和各部分的耦合效应。首先验证薄板模型的正确性，与相关文献[36]的结果进行对比：考虑薄板绕其固定端旋转时，薄板的变形响应。其中物理参数为：薄板的长、宽和厚度分别为a=1.8288 m,b=1.2192 m,h=0.00254 m，薄板的弹性模量为E=70 GPa，密度为ρ=2000 kg/m3，泊松比为ν=0.3。薄板定轴转动的运动规律为

式中:T=30 s,Ω=0.2π rad/s，薄板末端角点弹性变形的响应如图4所示，其中Case 1为本文计算的结果，Case 2为文献的相关数据，两者具有较好的一致性。

图4 板的末端角点变形Fig.4 The response of corner deformation of plate

为验证晃动模型的有效性，考虑球腔液体在横向简谐激励下的晃动力响应。本文晃动模型的数值仿真结果与文献[37]的CFD数值仿真和实验结果进行对比。参数为：球腔半径为R0=0.148 m，充液比为0.5，简谐激励的幅值为0.93 mm，频率为0.98倍液体晃动的基阶频率。液体晃动力的响应如图5所示，其中Case 1为本文给出晃动模型的数值计算结果，Case 2和Case 3分别为文献[37]的CFD数值仿真结果和实验数据。根据本文晃动模型的数值仿真结果与文献[37]的实验和CFD计算结果相比，在频率和相位上有较好的吻合效果，幅值上略小于文献[37]的结果。

图5 球腔液体的晃动力响应Fig.5 The liquid sloshing force response of the spherical tank with the transverse excitation

下面通过数值仿真研究含轴对称充液储腔和板类柔性附件的航天器系统动力学行为和耦合机理。板类柔性附件是由两块铝合金板面板和蜂窝状铝合金内芯通过树脂粘合而成，其材料参数为[27]：长、宽、厚度分别为a=9 m,b=3 m,h=0.0262 m，刚度E=4.45 Gpa，材料密度ρ=94.5 kg/m3，泊松比ν=0.3，阻尼系数cz=0.001。航天主刚体的参数为：质量m0=5000 kg，关于航天器本体系的惯量矩阵J0=diag(1250,1250,1500) kg·m2，板的连接点位置xb=1.0 m,yb=-1.5 m。液体推进剂的参数为：储液腔为球形储腔，半径R0=0.4 m,充液比0.4，液体密度ρf=1000 kg/m3,表面张力σ=0.0725 N/m,重力加速度g=0.02 m/s2，储腔在本体系的位置(-0.1,-0.15,-0.4) m。在该参数下，液体晃动的前两阶固有频率为0.2864 rad/s, 0.4673 rad/s，柔性附件的最低阶固有频率为2.3616 rad/s，考虑耦合航天器系统受到如下形式的零角动量的外激励

图6(a)、图6(b)分别为航天器相对于惯性系的速度、角速度响应。如图6所示，当航天器系统受到外激励时，航天器的速度和角速度发生变化，航天器的速度分量v1,v3和角速度分量ω2不为零，而速度和角速度的其他分量为零；在外激励结束后，航天器的速度和角速度均不为零，即航天器仍然存在平动和转动。这主要是由于航天器刚性平台与推进剂晃动、柔性附件振动的相互耦合效应引起的，若忽略推进剂晃动和帆板振动的影响，则航天器在零角动量的外激励下为先加速后减速的过程，最后状态为静止。图7为储腔液面波高的响应ζx，液体晃动以一阶反对称模态为主，高阶模态也会对液体晃动造成影响。液体仅存在沿着Ox方向的晃动，不存在沿Oy方向的晃动，这与外激励作用形式和柔性附件装配位置有关。图8为帆板末端弹性变形的响应。当外力矩作用于航天器系统时，帆板发生振动，弹性变形的响应与外激励的形式、作用时间有关；在外激励作用结束后，帆板仍然存在弹性振动，由于帆板结构阻尼较小，帆板的响应为衰减振动。振动以一阶对称模态(即弯曲变形)为主，高阶对称模态也会对帆板的振动造成影响。

图6 航天器刚性平台的响应Fig.6 The response of rigid platform of spacecraft

图7 液体液面波高的响应Fig.7 The response of liquid surface displacement

图8 柔性附件末端角点的弹性变形Fig.8 The response of corner elastic deformation of flexible appendage

下面考虑帆板装配位置的改变对于耦合航天器系统动力学行为的影响。令连接点B的位置发生变化，参数yb变为-1.25 m，其余参数均不发生变化，航天器系统受到的外激励的形式不变，图9(a)、图9(b)分别为航天器相对于惯性系的速度、角速度的响应。由于帆板装配位置的改变，导致航天器主刚体存在更为复杂的动力学响应，航天器的速度和角速度的各个分量均不为零。液体晃动参数的响应如图10所示，图10(a)、图10(b)分别为储腔液面波高在Ox,Oy方向的响应。与前一种工况相比，液体存在两个方向的晃动，柔性附件的装配位置的改变使得液体晃动的动力学行为更为复杂，即柔性附件与储腔液体晃动之间的相互耦合效应。图11为柔性附件末端弹性位移的响应，帆板的响应类似于前一种工况，但相位和幅值略有不同。

综上所述，帆板装配位置的改变深刻影响航天器动力学和液体晃动动力学，反应了航天器主刚体运动、柔性附件振动和储腔液体晃动之间的复杂耦合效应。当帆板的对称轴位于本体系惯性主轴时，帆板振动对于航天器主刚体和液体晃动动力学的影响较小，即航天器的速度和角速度的某些分量为零，且液体仅在某个方向发生晃动。

图9 航天器刚体平台的响应Fig.9 The response of rigid platform of spacecraft

图10 液体液面波高的响应Fig.10 The response of liquid surface displacement

图11 柔性附件末端角点的弹性位移Fig.11 The response of corner elastic deformation of flexible appendage

3 结 论

本文对含有板类柔性附件和轴对称充液储腔的耦合航天器系统进行动力学建模和刚-液-柔耦合效应研究。通过数值仿真，对含有帆板和球形充液储腔的航天器系统的动力学行为进行研究。研究结果表明，耦合航天器在外激励作用下存在复杂的动力学行为，主刚体、柔性附件、储腔液体各部分之间存在明显的耦合效应。在零角动量阶跃驱动下，航天器主刚体的运动会引起储腔推进剂晃动和帆板的振动，同时液体晃动和柔性附件的振动对于航天器姿态和位置有着不可忽略的影响；当跃驱动结束后，刚-液-柔耦合航天器系统仍然存在平动和转动，若只考虑航天器主刚体运动时，刚体航天器在驱动完成后为静止状态；在驱动施加和撤除的瞬间，液体晃动力和晃动力矩存在不连续的变化，液体晃动和帆板振动的参数幅值存在较大的改变。柔性附件的振动以一阶对称模态为主，高阶对称模态也会对帆板的振动造成影响；帆板的装配位置改变影响系统动力学行为，当帆板的对称轴位于本体系惯性主轴时，帆板振动对于航天器主刚体和液体晃动动力学的影响较小，航天器的速度和角速度某些分量为零。这可为耦合航天器系统的总体设计和姿态-轨道控制器设计提供一定参考。