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

邓林啸; 吕竹勇; 王天舒; 吕敬; Deng Linxiao; Lv Zhuyong; Wang Tianshu; Lv Jing

网刊加载中。。。

使用Chrome浏览器效果最佳，继续浏览，你可能不会看到最佳的展示效果，

确定继续浏览么?

复制成功，请在其他浏览器进行阅读

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

- ORCID：
邓林啸 ¹
- ORCID：
吕竹勇 ²
- ORCID：
王天舒 ³
- ORCID：
吕敬 ¹
✉

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

最近更新：2021-11-08

DOI：10.6052/1672-6553-2021-017

摘要

对于小幅晃动的充液飞行器流固耦合问题，可以采用等效力学模型来描述液体晃动.然而在研究液体对刚体的影响上，暂未有人从频域的角度进行建模.本文在已知液体各阶模态等效力学模型的基础上，利用虚功率原理和线性小扰动原理给出充液飞行器纵向模态舵偏角到姿态传函，并根据主导极点概念将充液飞行器高阶系统近似为二阶系统典型传递函数，以携带横放圆柱贮箱飞行器为例，分析了不同充液比和贮箱尺寸对典型传函系数的影响.结果表明，液体质量和转动惯量主要影响飞行器传递系数、时间常数，使得飞行器机动性和纵向模态自然频率降低；而液体晃动主要影响典型传函中相对阻尼系数，使得飞行器纵向相对稳定性降低.

关键词

小幅晃动; 流固耦合; 等效力学模型; 典型传递函数

2020-11-17收到第1稿，2021-03-10收到修改稿.

引言

在飞行器中，除了部分导弹和少数运载火箭使用固体燃料以外，其余大部分的飞行器包括航天器和航空器都使用液体燃料.而随着飞行器航程和任务复杂度的增加，液体燃料在飞行器整体质量上的占比越来越高，则自由液面所导致的液体晃动对飞行器的影响也越来越不可忽视^［

1］.从上世纪60年代起，相关的研究人员就开始研究储箱内的液体燃料对飞行器的姿态控制、飞行性能和稳定性的影响^{［参考文献 2~4}2~4］.

当前的研究中，对于飞行器整体耦合动力学建模过程中的晃动力计算通常采用CFD、SPH^［

5］、运动脉动球模型^{［参考文献 6

百度学术}6］、等效力学模型^{［参考文献 7

百度学术}7］等几种方式.但是受限于计算资源和内存，常用的CFD、SPH方法等在对充液储箱内的液体流场进行计算时用时太长，不利于整体耦合动力学方程的计算，所以通常可以利用简单的等效力学模型来代替复杂的晃动流场计算.岳宝增等^{［参考文献 4

百度学术}4］利用广义准坐标下的拉格朗日方法得到航天器刚体部分运动和液体燃料晃动的耦合动力学方程，采用复合控制方法对航天器姿态进行控制.张诗琪等^{［参考文献 5

百度学术}5］使用SPH方法计算晃动力，然后引入刚体约束方程得到固液耦合的动力学方程，分析了在周期激励情况下的固液耦合特性.

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

本文针对稳态运动中的充液飞行器，此时液体晃动为小幅晃动，故而利用等效力学模型^［

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

1 充液系统的增量方程

充液贮箱在小幅晃动情况下时，可以根据等效原则将液体晃动等效为等效摆力学模型，参考文献［

8］中液体晃动等效力学模型推导方法，利用虚功率原理推导出携带n_f个贮箱充液飞行器系统动力学方程和第i阶液体晃动方程如下：

m {\ddot{R}}_{c} - S \times {\dot{ω}}_{c} - \sum_{j = 1}^{n_{f}} \sum_{i = 1}^{n} (S_{p i} \times {\dot{ω}}_{p i}) = F + m g - Q_{t}

（1）

\begin{array}{l} S \times {\ddot{R}}_{c} + J \cdot {\dot{ω}}_{c} + \\ \sum_{j = 1}^{n_{f}} \sum_{i = 1}^{n} [J_{p i} \cdot {\dot{ω}}_{p i} - r_{c i} \times (S_{p i} \times {\dot{ω}}_{p i})] = S \times g + T - Q_{r} \end{array}

（2）

S_{p H i}^{T} \cdot {\ddot{R}}_{c} + J_{r s i}^{T} \cdot {\dot{ω}}_{c} + J_{p H i} δ {\ddot{q}}_{i} + C_{p i} \cdot ω_{p i} + Q_{i i} = - (S_{p H i}^{T} \cdot g)

（3）

上述中， $R_{c}$ 为飞行器在惯性系下位移， $ω_{c}$ 为飞行器绕机体系转动角速度，q为单摆相对飞行器机体系摆角；m、S、J为系统质量、系统相对飞行器机体系原点静矩、转动惯量矩阵， $S_{p i}$ 、 $J_{p i}$ 为第i阶单摆相对单摆系原点静矩、转动惯量；g为重力加速度； $ω_{p i}$ 为第i阶单摆相对飞行器机体系角速度；F、T为飞行器所受气动力和气动力矩； $Q_{t}$ 、 $Q_{r}$ 分别为液体晃动导致的平动耦合力、转动耦合力； $r_{c i}$ 、为第i阶单摆相对单摆悬挂点矢量； $H_{p i}$ 为第i阶单摆摆动面法线方向； $S_{p H i}$ 、 $J_{p H i}$ 、 $J_{r s i}$ 为与晃动摆相关项具体表达式可参考文献^［

8］.本文所用各坐标系示意图及其相互转换关系如图1所示，各坐标系定义如下：

图1 各坐标系关系示意图

Fig.1 The diagram of the relationship of each coordinate system

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

机体坐标系 $O x_{1} y_{1} z_{1}$ ：原点在飞行器质心， $O x_{1}$ 轴与飞行器纵轴重合，指向头部为正， $O y_{1}$ 轴位于飞行器横向对称面内与 $O x_{1}$ 轴垂直， $O z_{1}$ 轴垂直于另外两轴，方向按右手系确定.

弹道坐标系 $O x_{2} y_{2} z_{2}$ ：原点在飞行器瞬时质心上， $O x_{2}$ 轴与导弹的速度矢量 $V$ 重合， $O y_{2}$ 轴垂直于包含速度矢量 $V$ 的铅垂面， $O z_{2}$ 轴与其他两轴垂直构成右手坐标系并位于包含速度矢量的铅垂面内.

速度坐标系 $O x_{3} y_{3} z_{3}$ ：坐标系原点取在飞行器质心上， $O x_{3}$ 轴与飞行器质心的速度矢量 $V$ 重合， $O y_{3}$ 轴位于飞行器横向对称面内与 $O x_{3}$ 轴垂直， $O z_{3}$ 轴垂直于其它两轴并构成右手直角坐标系.

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

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

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

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

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

\{\begin{array}{l} m \frac{d Δ V}{d t} - \sum {\underset{̲}{S_{p h i}}}_{x} Δ {\ddot{q}}_{i} = (P^{V} - X^{V}) Δ V + \\ (- P α - X^{α}) Δ α - G c o s θ Δ θ \\ m V \frac{d Δ θ}{d t} - \sum {\underset{̲}{S_{p h i}}}_{z} Δ {\ddot{q}}_{i} = (P^{V} α + Y^{V}) Δ V + \\ (P + Y^{α}) Δ α + G s i n θ Δ θ + Y^{δ_{y}} Δ δ_{y} \\ - m V c o s θ \frac{d Δ ψ_{V}}{d t} - \sum {\underset{̲}{S_{p h i}}}_{y} Δ {\ddot{q}}_{i} = (- P + Z^{β}) Δ β + \\ (P α + Y) Δ γ_{V} + Z^{δ_{z}} Δ δ_{z} \end{array}

（4）

\{\begin{matrix} J {}_{x}\cdot Δ {\dot{ω}}_{x} + \sum_{i = 1}^{n} H_{i x} \cdot Δ {\ddot{q}}_{i} = Δ T_{x} \\ J {}_{y}\cdot Δ {\dot{ω}}_{y} + \sum_{i = 1}^{n} H_{i y} \cdot Δ {\ddot{q}}_{i} = Δ T_{y} \\ J {}_{z}\cdot Δ {\dot{ω}}_{z} + \sum_{i = 1}^{n} H_{i z} \cdot Δ {\ddot{q}}_{i} = Δ T_{z} \end{matrix}

（5）

H_{i x} \cdot Δ {\dot{ω}}_{x} + H_{i y} \cdot Δ {\dot{ω}}_{y} + H_{i z} \cdot Δ {\dot{ω}}_{z} + J_{p H i} \cdot Δ {\ddot{q}}_{i} + C_{p i} \cdot Δ {\dot{q}}_{i} = - Δ \underset{̲}{(S_{p H i}^{T} \cdot g)}

（6）

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

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

对于轴对称贮箱，等效力学模型前两阶晃动质量已经能够描述液体运动，根据假设一，影响飞行器纵向传函的主要是第二阶单摆，即在飞行器纵向对称面内摆动的单摆.将式（4）第一式，式（5）第二式和式（6）略去横侧向参数 $Δ ψ_{v}$ ， $γ$ ， $φ$ ， $q_{1}$ 的相关项，并引入几何关系，写出纵向线性小扰动方程组如式（7）所示：

\{\begin{array}{l} \begin{array}{l} \frac{d Δ V}{d t} = \frac{(P^{V} - X^{V})}{m} Δ V + \frac{(- P α - X^{α})}{m} Δ α - \frac{G c o s θ}{m} Δ θ \\ \frac{d Δ θ}{d t} = \frac{(P^{V} α + Y^{V})}{m V} Δ V + \frac{(P + Y^{α})}{m V} Δ α + \\ \frac{G s i n θ}{m V} Δ θ + \frac{Y^{δ_{y}}}{m V} Δ δ_{y} \\ Δ {\dot{ω}}_{y} + \frac{H_{2 y}}{J {}_{y}} \cdot Δ {\ddot{q}}_{2} = \frac{M_{y}^{V}}{J {}_{y}} Δ V + \frac{M_{y}^{α}}{J {}_{y}} Δ α + \frac{M_{y}^{δ_{y}}}{J {}_{y}} Δ δ_{y} + \\ \frac{M_{y}^{ω_{y}}}{J {}_{y}} Δ ω_{y} + \frac{M_{y}^{\dot{α}}}{J {}_{y}} Δ \dot{α} \end{array} \\ \begin{array}{l} H_{2 y} \cdot Δ {\dot{ω}}_{y} + J_{p H 2} Δ {\ddot{q}}_{2} + C_{p 2} Δ {\dot{q}}_{2} = - m_{s 2} g l_{2} Δ q \\ Δ α = Δ ϑ - Δ θ \\ \frac{d ω_{y}}{d t} = \ddot{ϑ} \end{array} \end{array}

（7）

引入与《导弹飞行力学》^［

9］一书相同纵向动力系数，将式（7）写为如下形式：

\{\begin{array}{l} \begin{array}{l} \frac{d Δ V}{d t} - a_{11} Δ V - a_{13} Δ θ - a_{14} Δ α = 0 \\ \frac{d^{2} Δ ϑ}{d t^{2}} - a_{21} Δ V - a_{22} \frac{d Δ ϑ}{d t} - a_{24} Δ α - a_{24}^{'} \frac{d Δ α}{d t} \\ + a_{26} \cdot Δ {\ddot{q}}_{2} = a_{25} Δ δ_{y} \\ \frac{d Δ θ}{d t} - a_{31} Δ V - a_{33} Δ θ - a_{34} Δ α = a_{35} Δ δ_{y} \end{array} \\ \begin{array}{l} H_{2 y} \cdot Δ {\dot{ω}}_{y} + J_{p H 2} Δ {\ddot{q}}_{2} + C_{p 2} Δ {\dot{q}}_{2} = - m_{s 2} g l_{2} Δ q \\ Δ α = Δ ϑ - Δ θ \end{array} \end{array}

（8）

其中 ${a_{i}}_{j}$ 定义与文献［

9］一致，

a_{26} = H_{2 y} / J_{y}

，

m_{s 2}

l_{2}

为沿x方向摆动单摆质量和摆长；

H_{2 y}

为沿x方向摆动单摆相关量

J_{r s i}

在飞行器机体系y轴上分量.将式（8）做拉氏变换后写为象函数，经整理后写为下列矩阵

\begin{array}{l} [\begin{matrix} s - a_{11} & 0 & - a_{13} & - a_{14} & 0 \\ - a_{21} & s^{2} - a_{22} s & 0 & - a_{24} - a_{24}^{'} s & a_{26} \\ - a_{31} & 0 & s - a_{33} & - a_{34} & 0 \\ 0 & - 1 & 1 & 1 & 0 \\ 0 & H_{2 y} s^{2} & 0 & 0 & J_{p H 2} s^{2} + C_{p 2} s + m_{s 2} g l_{2} \end{matrix}] \\ [\begin{matrix} Δ V (s) \\ Δ ϑ (s) \\ Δ θ (s) \\ Δ α (s) \\ Δ q_{2} (s) \end{matrix}] = [\begin{matrix} 0 \\ a_{25} \\ a_{35} \\ 0 \\ 0 \end{matrix}] Δ δ_{y} (s) \end{array}

（9）

根据克莱姆定理算出每一项为：

$Δ V = Δ_{V} (s) / Δ (s), Δ ϑ = Δ_{ϑ} (s) / Δ (s), Δ θ = Δ_{θ} (s) / Δ (s), Δ α = Δ_{α} (s) / Δ (s), Δ q_{2} = Δ_{q_{2}} (s) / Δ (s)$ 其中 $Δ (s)$ 为主行列式，Δ_V（s），Δ_ϑ（s），Δ_θ（s），Δ_α（s），Δq₂（s）为伴随行列式，伴随行列式可通过方程（9）右端列带入主行列式的各列得到.主行列式由齐次方程的系数组成

Δ (s) = |\begin{matrix} s - a_{11} & 0 & - a_{13} & - a_{14} & 0 \\ - a_{21} & s^{2} - a_{22} s & 0 & - a_{24} - a_{24}^{'} s & a_{26} \\ - a_{31} & 0 & s - a_{33} & - a_{34} & 0 \\ 0 & - 1 & 1 & 1 & 0 \\ 0 & H_{2 y} s^{2} & 0 & 0 & J_{p H 2} s^{2} + C_{p 2} s + m_{s 2} g l_{2} \end{matrix}|

（10）

整理后得

Δ (s) = P_{1} s^{6} + P_{2} s^{5} + P_{3} s^{4} + P_{4} s^{3} + P_{5} s^{2} + P_{6} s

（11）

其中系数P₁-P₆由主行列式 $Δ (s)$ 决定.根据传递函数定义，算出升降舵偏角 $Δ δ_{y} (s)$ 为输入，以 $Δ V (s)$ ， $Δ ϑ (s)$ ， $Δ θ (s)$ ， $Δ α (s)$ ， $Δ q_{2} (s)$ 为输出的传递函数.

\begin{array}{l} G_{δ_{y}}^{V} (s) = \frac{Δ V (s)}{Δ δ_{y} (s)} = \frac{Δ_{V} (s)}{Δ (s) Δ δ_{y} (s)} \\ = \frac{A_{1} s^{4} + A_{2} s^{3} + A_{3} s^{2} + A_{4} s}{P_{1} s^{6} + P_{2} s^{5} + P_{3} s^{4} + P_{4} s^{3} + P_{5} s^{2} + P_{6} s} \end{array}

（12）

\begin{array}{l} G_{δ_{y}}^{ϑ} (s) = \frac{Δ ϑ (s)}{Δ δ_{y} (s)} = \frac{Δ_{ϑ} (s)}{Δ (s) Δ δ_{y} (s)} \\ = \frac{B_{1} s^{4} + B_{2} s^{3} + B_{3} s^{2} + B_{4} s}{P_{1} s^{6} + P_{2} s^{5} + P_{3} s^{4} + P_{4} s^{3} + P_{5} s^{2} + P_{6} s} \end{array}

（13）

\begin{array}{l} G_{δ_{y}}^{θ} (s) = \frac{Δ θ (s)}{Δ δ_{y} (s)} = \frac{Δ_{θ} (s)}{Δ (s) Δ δ_{y} (s)} \\ = \frac{C_{1} s^{4} + C_{2} s^{3} + C_{3} s^{2} + C_{4} s}{P_{1} s^{6} + P_{2} s^{5} + P_{3} s^{4} + P_{4} s^{3} + P_{5} s^{2} + P_{6} s} \end{array}

（14）

\begin{array}{l} G_{δ_{y}}^{α} (s) = \frac{Δ θ_{α} (s)}{Δ δ_{y} (s)} = \frac{Δ_{α} (s)}{Δ (s) Δ δ_{y} (s)} \\ = \frac{D_{1} s^{4} + D_{2} s^{3} + D_{3} s^{2} + D_{4} s}{P_{1} s^{6} + P_{2} s^{5} + P_{3} s^{4} + P_{4} s^{3} + P_{5} s^{2} + P_{6} s} \end{array}

（15）

\begin{array}{l} G_{δ_{y}}^{q_{2}} (s) = \frac{Δ q_{2} (s)}{Δ δ_{y} (s)} = \frac{Δ_{q_{2}} (s)}{Δ (s) Δ δ_{y} (s)} \\ = \frac{E_{1} s^{4} + E_{2} s^{3} + E_{3} s^{2} + E_{4} s}{P_{1} s^{6} + P_{2} s^{5} + P_{3} s^{4} + P_{4} s^{3} + P_{5} s^{2} + P_{6} s} \end{array}

（16）

其中 $A_{i}, B_{i}, C_{i}, D_{i}, E_{i}$ 系数由其伴随行列式决定.当设计飞行器时，只研究扰动运动短周期阶段.为了简化制导系统的分析，需对导弹的纵向扰动运动现象加以简化，即只讨论速度偏量 $Δ V$ 可以忽略的短周期扰动运动，即在假设三的前提下，

建立短周期扰动运动近似方程组.由扰动运动方程组（8）中去掉描述 $Δ V$ 变化的第一个方程，在其余的方程式中，令 $Δ V = 0$ ，于是，得到简化后的扰动运动方程组，再将其拉普拉斯化，写成矩阵形式：

\begin{array}{l} [\begin{matrix} s^{2} - a_{22} s & 0 & - a_{24} - a_{24}^{'} s & a_{26} \\ 0 & s - a_{33} & - a_{34} & 0 \\ - 1 & 1 & 1 & 0 \\ H_{2 y} s^{2} & 0 & 0 & J_{p H 2} s^{2} + C_{p 2} s + m_{s 2} g l_{2} \end{matrix}] \\ [\begin{matrix} Δ ϑ (s) \\ Δ θ (s) \\ Δ α (s) \\ Δ q_{2} (s) \end{matrix}] = [\begin{matrix} a_{25} \\ a_{35} \\ 0 \\ 0 \end{matrix}] Δ δ_{y} (s) \end{array}

（17）

方程组（17）主要是描述导弹的角运动，只适用于不超过几秒钟的短暂时间.将式（17）第四式单独拿出，有

Δ q_{2} (s) = \frac{- H_{2 y} s^{2} Δ ϑ (s)}{J_{p H 2} s^{2} + C_{p 2} s + m_{s 2} g l_{2}}

（18）

将式（18）带入（17）第一式得到：

[(1 - \frac{H_{2 y} a_{26}}{J_{p H 2} s^{2} + C_{p 2} s + m_{s 2} g l_{2}}) s^{2} - a_{22} s] Δ ϑ (s) - (a_{24} + a_{24}^{'}) Δ α (s) = a_{25} Δ δ_{y} (s)

（19）

定义晃动相关系数

a_{h} = 1 - \frac{H_{2 y} a_{26}}{J_{p H 2} s^{2} + C_{p 2} s + m_{s 2} g l_{2}}

（20）

则有

[\begin{matrix} a_{h} s^{2} - a_{22} s & 0 & - a_{24} - a_{24}^{'} s \\ 0 & s - a_{33} & - a_{34} \\ - 1 & 1 & 1 \end{matrix}] [\begin{matrix} Δ ϑ (s) \\ Δ θ (s) \\ Δ α (s) \end{matrix}] = [\begin{matrix} a_{25} \\ a_{35} \\ 0 \end{matrix}] Δ δ_{y} (s)

（21）

主行列式为：

\begin{matrix} Δ (s) = |\begin{matrix} a_{h} s^{2} - a_{22} s & 0 & - a_{24} - a_{24}^{'} s \\ 0 & s - a_{33} & - a_{34} \\ - 1 & 1 & 1 \end{matrix}| \\ = a_{h} s^{3} + P_{1} s^{2} + P_{2} s + P_{3} \end{matrix}

（22）

其中

\{\begin{array}{l} P_{1} = a_{h} (- a_{33} + a_{34}) - a_{22} - a_{24}^{'} \\ P_{2} = a_{22} (a_{33} - a_{34}) - a_{24} + a_{33} a_{24}^{'} \\ P_{3} = a_{24} a_{33} \end{array}

（23）

利用克莱姆定理，得到传递函数为

\begin{array}{l} G_{δ_{y}}^{ϑ} (s) = \frac{Δ J (s)}{Δ δ_{y} (s)} = \frac{Δ_{ϑ} (s)}{Δ (s) Δ δ_{y} (s)} = \\ = \frac{(- a_{35} a_{24}^{'} s + a_{25}) + a_{25} (a_{34} - a_{33}) - a_{35} a_{24}}{a_{h} s^{3} + (a_{h} (- a_{33} + a_{34}) - a_{22} - a_{24}^{'}) s^{2} + (a_{22} (a_{33} - a_{34}) - a_{24} + a_{33} a_{24}^{'}) s + a_{24} a_{33}} \end{array}

（24）

\begin{array}{l} G_{δ_{y}}^{θ} (s) = \frac{Δ θ (s)}{Δ δ_{y} (s)} = \frac{Δ_{θ} (s)}{Δ (s) Δ δ_{y} (s)} \\ = \frac{a_{h} a_{35} s^{2} - a_{35} (a_{22} + a_{24}^{'}) s + a_{34} a_{25} - a_{24} a_{35}}{a_{h} s^{3} + (a_{h} (- a_{33} + a_{34}) - a_{22} - a_{24}^{'}) s^{2} + (a_{22} (a_{33} - a_{34}) - a_{24} + a_{33} a_{24}^{'}) s + a_{24} a_{33}} \end{array}

（25）

\begin{array}{l} G_{δ_{y}}^{α} (s) = \frac{Δ θ_{α} (s)}{Δ δ_{y} (s)} = \frac{Δ_{α} (s)}{Δ (s) Δ δ_{y} (s)} \\ = \frac{- a_{h} a_{35} s^{2} + (a_{35} a_{22} + a_{25}) s - a_{25} a_{33}}{a_{h} s^{3} + (a_{h} (- a_{33} + a_{34}) - a_{22} - a_{24}^{'}) s^{2} + (a_{22} (a_{33} - a_{34}) - a_{24} + a_{33} a_{24}^{'}) s + a_{24} a_{33}} \end{array}

（26）

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

3.1　开环性能分析

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

\begin{array}{l} G_{δ_{y}}^{ϑ} (s) = \frac{a_{25} s + a_{25} (a_{34} - a_{33}) - a_{35} a_{24}}{s^{2} + (a_{34} - a_{22}) s - a_{22} a_{34} - a_{24}} \\ = \frac{K_{M} (T_{1} s + 1)}{s (T_{M}^{2} s^{2} + 2 T_{M} ξ_{M} s + 1)} \end{array}

（27）

其中 $K_{M}$ 为飞行器传递系数，代表稳态时输出变量与输入变量之比： $K_{M} = \frac{- a_{25} a_{34} + a_{34} a_{24}}{a_{22} a_{34} + a_{24}} (s^{- 1})$ ； $T_{M}$ 为飞行器时间常数， $1 / T_{M}$ 为纵向模态自然频率. $T_{M} = \frac{1}{\sqrt[]{- a_{24} - a_{22} a_{34}}} (s)$ ； $ξ_{M}$ 为飞行器相对阻尼系数， $ξ_{M} = \frac{- a_{22} + a_{34}}{2 \sqrt[]{- a_{24} - a_{22} a_{34}}}$ ； $T_{1}$ 为飞行器气动力时间常数， $T_{1} = \frac{a_{25}}{a_{25} a_{34} + a_{35} a_{24}} (s)$ .同理可以得到其他传函典型基本环节传递函数：

G_{δ_{y}}^{θ} (s) = \frac{K_{M} [1 + T_{1} \frac{a_{35}}{a_{25}} s (s - a_{22})]}{s (T_{M}^{2} s^{2} + 2 T_{M} ξ_{M} s + 1)}

（28）

G_{δ_{y}}^{α} (s) = \frac{K_{M} T_{1} [1 - \frac{a_{35}}{a_{25}} (s - a_{22})]}{(T_{M}^{2} s^{2} + 2 T_{M} ξ_{M} s + 1)}

（29）

根据式（27）~（29）可知， $G_{δ_{y}}^{{\dot{θ}}_{y}} (s)$ 、 $G_{δ_{y}}^{\dot{θ}} (s)$ 和 $G_{δ_{y}}^{α} (s)$ 特征方程式为

T_{M}^{2} λ^{2} + 2 T_{M} ξ_{M} λ + 1 = 0

（30）

式（30）根为：

\begin{array}{l} λ_{1,2} = - \frac{ξ_{M}}{T_{M}} \pm \frac{\sqrt[]{ξ_{M}^{2} - 1}}{T_{M}} \\ = - \frac{1}{2} (a_{34} - a_{22} - a_{24}^{'}) \pm \\ \sqrt[]{\frac{1}{4} (a_{34} - a_{22} - a_{24}^{'})^{2} + (a_{22} a_{34} + a_{24})} \end{array}

（31）

飞行器具有纵向静稳定性的条件为 $a_{22} a_{34} + a_{24} < 0$ ；下面仅限于讨论满足静稳定性的飞行器.此外，从式（27）~（29）可知，只要分析出升降舵到俯仰角典型基本环节参数，就能分析其余所有传函，故而以下对 $G_{δ_{y}}^{{\dot{θ}}_{y}} (s)$ 典型基本环节参数进行分析.充液飞行器 $G_{δ_{y}}^{ω_{y}} (s)$ 略去下洗动力系数 $a_{24}^{'}$ 且不计重力影响，再将式（20）带入式（24）后写可写为式（32）形式，因 $H_{2 y} a_{26}$ 为小量，故而略去公因式（即偶极子）后表达成典型环节传函.

\begin{array}{l} G_{δ_{y}}^{ω_{y}} (s) = \frac{a_{25} s + a_{25} a_{34} - a_{35} a_{24}}{s^{2} + (a_{34} - a_{22} / a_{h}) s + (- a_{22} a_{34} - a_{24}) / a_{h}} \cdot (\frac{J_{p H 2} s^{2} + C_{p 2} s + m_{s 2} g l_{2}}{J_{p H 2} s^{2} + C_{p 2} s + m_{s 2} g l_{2} - H_{2 y} a_{26}}) \\ \approx \frac{a_{25} s + a_{25} a_{34} - a_{35} a_{24}}{s^{2} + (a_{34} - a_{22} / a_{h}) s + (- a_{22} a_{34} - a_{24}) / a_{h}} = \frac{K_{M} (T_{1} s + 1)}{(T_{M}^{2} s^{2} + 2 T_{M} ξ_{M} s + 1)} \end{array}

（32）

3.2　刚体飞行器算例

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

G_{δ_{y}}^{{\dot{θ}}_{y}} (s) = \frac{- 8 s - 0.6208}{s^{2} + 0.151 s + 3.846}

（33）

其传递系数为 $K_{M} = - 0.1614$ ：对于阶跃升降舵输入，响应稳定值为-0.1614；时间常数 $T_{M} = 0.5099$ ，其固有频率 $ω_{c} = 1 / T_{M} = 1.961$ ；相对阻尼系数 $ξ_{M} = 0.0385$ ；气动力时间常数 $T_{1} = 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）分别为刚体飞行器和充液飞行器升降舵到俯仰角速度开环传递函数.对比两式，未充液系统（原始模型）在增加贮箱液体后，首先需要考虑飞行器质量和惯性张量变化对动力系数的影响，即贮箱液体的重量和转动惯量对传函的影响（质量惯量修正模型）；然后再考虑晃动参数 $a_{h}$ 对其的影响，即液体晃动对传函的影响（液体晃动模型），其示意图如图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）所示，考虑横放圆柱贮箱质量及惯量，对动力系数 $a_{22}$ ， $a_{24}$ ， $a_{25}$ ， $a_{34}$ ， $a_{35}$ ， $a_{33}$ 进行修正，得到其修正模型开环传递函数为：

G_{δ_{y}}^{ω_{y}} (s) = \frac{- 7.768 s - 0.4614}{s^{2} + 0.1302 s + 3.733}

（34）

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

\begin{array}{l} G_{δ_{y}}^{ω_{y}} (s) = \frac{- 174.2 s^{3} - 12.5 s^{2} - 4971 s - 295.3}{22.43 s^{4} + 1.72 s^{3} + 723.2 s^{2} + 40.45 s + 2389} \\ = \frac{- 7.7664 (s + 0.05941) (s^{2} + 0.01235 s + 28.54)}{(s^{2} + 0.06124 s + 3.736) (s^{2} + 0.01544 s + 28.51)} \end{array}

（35）

做出上述三种模型开环传函零极点的图像如图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%时为无充液飞行器模型相关参数.随着充液比的增加，因为液体附加质量和惯性的影响，飞行器传递系数逐渐降低，操纵性降低；飞行器时间常数变化不大，即飞行器固有频率基本不变；相对阻尼系数逐渐降低，液体晃动对相对阻尼系数影响最为剧烈，即液体晃动会导致阶跃响应超调量和调节时间增加，动态稳定性降低；气动力时间常数逐渐增加，主要是由附加液体质量和惯量引起，即随着充液比的增加，零点绝对值逐渐减小.

图8 不同充液比下模型传函系数变化图

Fig.8 Variation of model transfer coefficient with different filling ratios

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 小结

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

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

参考文献

苗楠，王天舒，李俊峰.微重环境下液体晃动研究进展.力学与实践，2016，38（3）：229~236 [百度学术]

Miao N， Wang T S， Li J F. Research progress of liquid sloshing in micro gravity environment. Mechanics and Practice， 2016，38 （3）： 229~236（in Chinese） [百度学术]

Baiter， H.F. Stability Boundaries of Liquid-Propelled Space Vehicles with Sloshing. AIAA Journal， 1963， 1（7）：1583～1589 [百度学术]

祝乐梅，岳宝增.充液航天器姿态的自适应非线性动态逆控制.动力学与控制学报，2011，9（4）：321~325 [百度学术]

Zhu L M， Yue B Z. Adaptive nonlinear dynamic inverse control of attitude of liquid filled spacecraft. Journal of Dynamics and Control， 2011，9（4）：321~325（in Chinese） [百度学术]

岳宝增，于嘉瑞，吴文军.多储液腔航天器刚液耦合动力学与复合控制.力学学报，2017，49（2）：390~396 [百度学术]

Yue B Z， Yu J R， Wu W J. Rigid liquid coupling dynamics and composite control of multi cavity spacecraft. Acta Mechanica Sinica，2017，49（2）：390~396（in Chinese） [百度学术]

张诗琪，刘锦阳.基于SPH的充液刚体刚-液耦合动力学建模.动力学与控制学报，2020，18（4）：86~94 [百度学术]

Zhang S Q， Liu J Y. Rigid liquid coupling dynamic modeling of rigid body filled with liquid based on SPH. Journal of Dynamics and Control， 2020，18（4）：86~94（in Chinese） [百度学术]

邓明乐.液体大幅晃动等效力学模型及航天器刚—液—柔—控耦合动力学研究［博士学位论文］.北京理工大学，2017（Deng M L. Equivalent mechanical model of liquid sloshing and rigid liquid flexible control coupling dynamics of spacecraft ［Ph. D Thesis］. Beijing University of technology， 2017（in Chinese）） [百度学术]

李青，马兴瑞，王天舒.非轴对称贮箱液体晃动的等效力学模型.宇航学报，2011，32（2）：242~249 [百度学术]

Li Q， Ma X R， Wang T S. Equivalent mechanical model of liquid sloshing in non axisymmetric tank. Acta Astronautics Sinica， 2011，32（2）：242-249（in Chinese） [百度学术]

李青. 充液挠性系统动力学分析及在航天工程中的应用研究［博士学位论文］.北京：清华大学，2010.（Li Q. Dynamic analysis of fluid filled flexible system and its application in aerospace engineering ［Ph. D Thesis］. Beijing： Tsinghua University， 2010 （in Chinese）） [百度学术]

钱杏芳，林瑞雄，赵亚男.导弹飞行力学.北京：北京理工大学，2006：175~229（Qian X F， Lin R X， Zhao Y N. Missile flight mechanics. Beijing： Beijing University of Technology， 2006：175 ~ 229（in Chinese）） [百度学术]

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

摘要

关键词

引言

1 充液系统的增量方程

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