液体晃动对飞行器控制系统影响分析*

邓林啸 吕竹勇 王天舒 吕敬†

（1.北京航空航天大学航空科学与工程学院，北京 100191）（2.北京空天技术研究所，北京 100074）（3.清华大学航空航天学院 北京 100084）

引言

在飞行器中，除了部分导弹和少数运载火箭使用固体燃料以外，其余大部分的飞行器包括航天器和航空器都使用液体燃料.而随着飞行器航程和任务复杂度的增加，液体燃料在飞行器整体质量上的占比越来越高，则自由液面所导致的液体晃动对飞行器的影响也越来越不可忽视[1].从上世纪60年代起，相关的研究人员就开始研究储箱内的液体燃料对飞行器的姿态控制、飞行性能和稳定性的影响[2～4].

当前的研究中，对于飞行器整体耦合动力学建模过程中的晃动力计算通常采用CFD、SPH[5]、运动脉动球模型[6]、等效力学模型[7]等几种方式.但是受限于计算资源和内存，常用的CFD、SPH方法等在对充液储箱内的液体流场进行计算时用时太长，不利于整体耦合动力学方程的计算，所以通常可以利用简单的等效力学模型来代替复杂的晃动流场计算.岳宝增等[4]利用广义准坐标下的拉格朗日方法得到航天器刚体部分运动和液体燃料晃动的耦合动力学方程，采用复合控制方法对航天器姿态进行控制.张诗琪等[5]使用SPH方法计算晃动力，然后引入刚体约束方程得到固液耦合的动力学方程，分析了在周期激励情况下的固液耦合特性.

对于飞行器姿态机动过程中流固耦合问题，目前多数学者都聚焦于真空微重环境下对航天器的研究，鲜少有对大气中飞行的航空器的分析，故本文研究大气飞行环境下液体晃动对飞行器控制系统的影响.

本文针对稳态运动中的充液飞行器，此时液体晃动为小幅晃动，故而利用等效力学模型[8]描述液体晃动，将液体晃动问题转化为多刚体运动动力学问题；在利用虚功率原理推导出多贮箱充液系统动力学方程后，结合线性小扰动理论，建立起充液飞行器线性小扰动方程组，并推导出飞行器纵向模态传递函数.最后分析了横放圆柱贮箱尺寸和充液比对飞行器纵向典型传函系数的影响.

1 充液系统的增量方程

充液贮箱在小幅晃动情况下时，可以根据等效原则将液体晃动等效为等效摆力学模型，参考文献[8]中液体晃动等效力学模型推导方法，利用虚功率原理推导出携带nf个贮箱充液飞行器系统动力学方程和第i阶液体晃动方程如下：

上述中，Rc为飞行器在惯性系下位移，ωc为飞行器绕机体系转动角速度，q为单摆相对飞行器机体系摆角；m、S、J为系统质量、系统相对飞行器机体系原点静矩、转动惯量矩阵，Spi、Jpi为第i阶单摆相对单摆系原点静矩、转动惯量；g为重力加速度；ωpi为第i阶单摆相对飞行器机体系角速度；F、T为飞行器所受气动力和气动力矩；Qt、Qr分别为液体晃动导致的平动耦合力、转动耦合力；rci、为第i阶单摆相对单摆悬挂点矢量；Hpi为第i阶单摆摆动面法线方向；SpHi、JpHi、Jrsi为与晃动摆相关项具体表达式可参考文献[8].本文所用各坐标系示意图及其相互转换关系如图1所示，各坐标系定义如下：

图1 各坐标系关系示意图Fig.1 The diagram of the relationship of each coordinate system

地面坐标系Oxyz：坐标系原点取在飞行器发射点上，x轴沿弹道面与地面交线指向目标点处，z轴沿铅垂面向下，y轴与其它两轴垂直并构成右手系.

机体坐标系Ox1y1z1：原点在飞行器质心，Ox1轴与飞行器纵轴重合，指向头部为正，Oy1轴位于飞行器横向对称面内与Ox1轴垂直，Oz1轴垂直于另外两轴，方向按右手系确定.

弹道坐标系Ox2y2z2：原点在飞行器瞬时质心上，Ox2轴与导弹的速度矢量V重合，Oy2轴垂直于包含速度矢量V的铅垂面，Oz2轴与其他两轴垂直构成右手坐标系并位于包含速度矢量的铅垂面内.

速度坐标系Ox3y3z3：坐标系原点取在飞行器质心上，Ox3轴与飞行器质心的速度矢量V重合，Oy3轴位于飞行器横向对称面内与Ox3轴垂直，Oz3轴垂直于其它两轴并构成右手直角坐标系.

为了简化推导，引入下列假设.

假设一：飞行器携带轴对称贮箱，且贮箱几何中心放置于飞行器质心处，此时对于轴对称贮箱，其前两阶等效摆面法线分别与飞行器机体系x、y共线；

假设二：飞行器推力沿机体系x轴方向，由于飞行器本体姿态角变化很小，可以认为三个姿态通道是解耦的.

假设三：把速度看做时间的已知函数V(t)=V0(t)，且小扰动、未扰动运动的侧向参数及纵向角速度足够小.

将飞行器机体系原点放置于系统质心，将平动方程投影到弹道坐标系，转动方程投影到飞行器机体系、液体晃动方程投影到单摆本体系，并引入相关力与力系数，得到系统平动方程、转动方程、第i阶液体晃动方程增量方程如式（4）～（6）所示：

其中，θ、ψV为弹道倾角和弹道偏角；δy、δz为升降舵偏角和方向舵偏角；α和β为飞行迎角和侧滑角；γ、φ、ϑ为滚转角、偏航角、俯仰角；P、X、Y、Z为发动机推力、气动力在速度坐标系下分量；G、V为系统重力和飞行器速度；ωxωyωz为飞行器角速度在机体坐标系分量；JxJyJz为系统绕飞行器机体系轴转动惯量；HixHiyHiz为第i阶晃动摆相关项Jrsi在飞行器机体系下投影；PV为使飞行器速度V产生单位变化时所需推力增量，其余上标含义类似.

2 飞行器纵向传递函数推导

3 升降舵到俯仰角速度控制分析

3.1 开环性能分析

对于有翼式飞行器，若不计重力影响，即a33=0，同时略去下洗动力系数a'24，则不考虑液体晃动时，由式（24）略去液体晃动参数an可得飞行器升降舵到俯仰角开环传递函数，并表示为典型基本环节传递函数：

3.2 刚体飞行器算例

预设飞行器质量为1000kg，绕机体轴y轴转动惯量为1000kg·m2.未扰时平飞运动速度V=100m/s；飞行器相关动力系数取值为：a22=-0.071，a24=-3.84，a25=-8，a35=0.005，a33=0.带入上述系数，未受液体晃动扰动前飞行器升降舵到俯仰角速度开环传函如下所示

其传递系数为KM=-0.1614：对于阶跃升降舵输入，响应稳定值为-0.1614；时间常数TM=0.5099，其固有频率ωc=1/TM=1.961；相对阻尼系数ξM=0.0385；气动力时间常数T1=12.1359.根据式（33）做出其零极点、根轨迹、伯德图和阶跃响应曲线如图2所示：如图2所示，左上角为开环传函零极点分布图，可以看到开环极点实数小于零，所以开环系统稳定.右上角为比例系数的根轨迹图，可以看到负反馈闭环极点均在虚轴左侧，故而其闭环稳定.左下角为开环传函的伯德图，其相位裕量为79.4deg，增益裕量为15.8dB.左下角为其开环阶跃响应.可以看到开环系统稳定，但其阶跃响应超调量大，调节时间长，稳定误差大.

图2 开环响应曲线Fig.2 Open loop response curve

3.3 液体晃动对系统阶跃响应的影响

考虑图3所示横放圆柱贮箱，其长度为L，半径为R，选取L=1m，R=0.5m，充液比50%的贮箱.贮箱中心放置于飞行器质心处，轴向和飞行器机体轴重合；液体为航空煤油，此时液体重量为306.3kg；分析其对系统纵向传递函数的影响.式（27）和（32）分别为刚体飞行器和充液飞行器升降舵到俯仰角速度开环传递函数.对比两式，未充液系统（原始模型）在增加贮箱液体后，首先需要考虑飞行器质量和惯性张量变化对动力系数的影响，即贮箱液体的重量和转动惯量对传函的影响（质量惯量修正模型）；然后再考虑晃动参数ah对其的影响，即液体晃动对传函的影响（液体晃动模型），其示意图如图4所示.

图3 50%充液比贮箱内液体的有限元模型Fig.3 Finite element model of liquid in tank with 50% liquid filling ratio

图4 不同模型示意图Fig.4 Schematic diagram of different models

未充液飞行器模型如式（33）所示，考虑横放圆柱贮箱质量及惯量，对动力系数a22，a24，a25，a34，a35，a33进行修正，得到其修正模型开环传递函数为：

最后考虑液体晃动相关参数ah对其的影响，得到其液体晃动模型开环传递函数为：

做出上述三种模型开环传函零极点的图像如图5所示，可看到所有极点实部均小于零，即开环系统稳定；修正模型相比于原始模型而言，在附加液体质量和转动惯量的影响下，开环传函零极点均右移，其阻尼和自然频率均减小，阻尼频率变化不大，所以其阶跃响应波动增大，而周期基本不变；液体晃动模型相比于修正模型而言，其多出了两对偶极子，零点没有出现移动，主极点左移趋势明显，晃动频率基本不变，但阻尼却大幅降低，其阶跃响应振荡加剧，动态响应将进一步恶化.在后续分析中，忽略偶极子的影响，将液体晃动模型简化为二阶晃动模态进行分析，如式（35）所示.

图5 不同模型零极点分布图Fig.5 Distribution of poles and zeros of different models

不同模型的飞行器传递系数、时间常数、相对阻尼系数、气动力时间常数如图6所示，可以看到修正模型跟原始模型之间变化的主要是传递系数降低，即飞行器操纵性降低；气动力时间常数增大，即零点出现大幅偏移.而晃动模型相对于修正模型，相对阻尼系数出现很大的变化，即其阶跃响应超调量将会增加，飞行器操纵时过载将会增加.不同模型的阶跃响应曲线如图7所示，液体晃动对飞行器开环响应影响很大，其超调量和调节时间均大幅增加.同时引起飞行器操纵性降低，过载增加.

图6 不同模型系数变化示意图Fig.6 Variation diagram of different model coefficients

图7 开环阶跃响应Fig.7 Open loop step response

3.4 不同充液比对传递函数的影响

选取长度L=1，半径R=0.5m的横放圆柱贮箱，贮箱充液比变化时，分析其开环响应.充液比取0%～90%时，修正模型和晃动模型飞行器传递系数、时间常数、相对阻尼系数、气动力时间常数如图8所示，其充液0%时为无充液飞行器模型相关参数.随着充液比的增加，因为液体附加质量和惯性的影响，飞行器传递系数逐渐降低，操纵性降低；飞行器时间常数变化不大，即飞行器固有频率基本不变；相对阻尼系数逐渐降低，液体晃动对相对阻尼系数影响最为剧烈，即液体晃动会导致阶跃响应超调量和调节时间增加，动态稳定性降低；气动力时间常数逐渐增加，主要是由附加液体质量和惯量引起，即随着充液比的增加，零点绝对值逐渐减小.

3.5 不同长度L对传递函数的影响

保持圆柱半径R=0.5m，充液比为50%.L分别取1m，1.2m1.4m，1.6m，1.8m，分析其对系统控制的影响.加入横放圆柱贮箱后，随着贮箱长度的增加L的增加，相关参数的变化如图9所示，表中给出了不同长度下的修正模型和晃动模型飞行器传递系数、时间常数、相对阻尼系数、气动力时间常数的数值.可以看到随着贮箱长度L的增加，液体附加质量和惯量引起了传递系数降低，时间常数增加，相对阻尼比降低，气动力时间常数增加；分别代表飞行器操纵性降低，固有频率降低，相对稳定性变差，零点绝对值降低.而晃动模型相对修正模型，仍然是出现了相对阻尼系数大幅降低，阶跃响应超调量增加，调节时间变长，相对稳定性进一步变差.即对于阶跃响应而言，振荡峰值将会随着L的增加而加剧，阻尼振荡频率减小，振荡周期将增加.

图9 不同L下模型传函系数变化图Fig.9 Variation of model transfer coefficient with different L

3.6 不同半径R对开环响应的影响

保持横放贮箱长度L=1m，充液比为50%，充液液体为航空煤油，在常重g=10的情况下分析.半径R分别取0.1～0.5m，分析其对系统控制的影响.随着贮箱半径尺寸的变化，修正模型和液体晃动模型典型基本环节表示的传函相关参数如图10所示。如图10所示.随着R的增加，修正模型和晃动模型飞行器传递系数减小，操纵性降低，机动性能降低；飞行器时间常数轻微增加，飞行器固有频率稍稍降低；气动力时间常数增大，零点绝对值减小；相对阻尼系数主要还是由液体晃动引起，阶跃响应超调量和调节时间均大幅增加，相对稳定性降低.

图10 不同R下模型传函系数变化图Fig.10 Variation of model transfer coefficient with different R

4 小结

本文针对携带液体贮箱的飞行器，基于等效力学模型将稳态工作下的液体小幅晃动问题转化为多体系统动力学问题，通过虚功率原理推导了携带多贮箱的多阶耦合充液系统动力学方程，实现了液体、固体的实时耦合；并利用线性小扰动理论给出了充液飞行器流固耦合的线性小扰动方程，给出了纵向模态升降舵偏角到飞行器姿态传函.

最后针对携带横放圆柱贮箱的飞行器，利用主导极点概念，分别从贮箱内液体充液比和贮箱尺寸的角度分析了液体晃动对飞行器纵向典型传递函数系数的影响.液体晃动主要影响了典型传函中的相对阻尼系数，舵偏角阶跃响应超调量增加，调节时间变长，相对稳定性变差.