CN102306001A - 一种确定储箱内液体晃动作用力的仿真方法 - Google Patents

Abstract

一种确定储箱内液体晃动作用力的仿真方法，储箱内液体晃动所产生的沿储箱o_ty_t轴的作用力F_y采用公式F_y＝F_τ cos β+F_ne_n，j计算获得，储箱内液体晃动所产生的沿储箱o_tz_t轴的作用力F_z采用公式F_z＝F_τ sin β+F_ne_n，k计算获得，储箱内液体晃动所产生的沿绕储箱o_tx_t轴的以逆时针方向为正的晃动力矩M_x采用公式M_x＝(F_τ sin β+F_ne_n，k)Y-(F_τ cos β+F_ne_n，j)Z计算获得。本发明方法将储箱内的液体等效为质量集中在质心处的质点(质心点)，质心点只能在质心约束面上运动，通过求解质心点与质心面之间的相互作用来计算液体晃动时对储箱的作用力和力矩。该方法计算效率高，同时还克服了现有单摆或弹簧-振子等效力学模型法只能适用于小幅非线性晃动的缺点。

Description

一种确定储箱内液体晃动作用力的仿真方法

技术领域

本发明涉及一种确定液体晃动所产生作用力的仿真方法，可以应用于航天器化学推进剂贮箱内大幅液体晃动的动力学分析，从而确定储箱内化学推进剂晃动时产生的作用力和力矩。

背景技术

目前，对于航天器化学推进剂贮箱内的液体晃动问题，主要的分析方法有基于计算流体动力学(CFD)的数值仿真方法，以及工程上常采用的单摆或弹簧-振子等效力学模型方法两种。

CFD方法一般采用有限元或有限差分等数值分析方法，通过求解流体运动的Navier-Stokes方程来计算液体晃动作用力。该方法计算精度较高，但计算量很大，灵活性不足，特别是难以用于贮箱内液体晃动与航天器姿态的耦合分析，不便于进行航天器的控制系统仿真分析。

单摆或弹簧-振子等效力学模型是工程常用的简化分析模型，即用一单摆或弹簧-振子系统来代替储箱内的液体，常用于航天器控制系统的仿真分析。但由于该模型是基于小幅非线性晃动的假设建立的，因此不适用于大幅非线性晃动计算分析。

发明内容

本发明的技术解决问题是：克服现有技术的不足，提供了一种基于质心点及质心约束面的等效力学模型来确定储箱内液体晃动所产生作用力的方法，便于进行非线性大幅液体晃动时的仿真分析。

本发明的技术解决方案是：一种确定储箱内液体晃动作用力的仿真方法，储箱内液体晃动所产生的沿储箱o_ty_t轴的作用力F_y采用公式F_y＝F_τ cos β+F_ne_n，j计算获得，储箱内液体晃动所产生的沿储箱o_tz_t轴的作用力F_z采用公式F_z＝F_τ sin β+F_ne_n，k计算获得，储箱内液体晃动所产生的沿绕储箱o_tx_t轴的以逆时针方向为正的晃动力矩M_x采用公式M_x＝(F_τ sin β+F_ne_n，k)Y-(F_τ cos β+F_ne_n，j)Z计算获得；其中o_tx_ty_tz_t为储箱本体坐标系，原点o_t位于储箱的几何中心，z_t轴为储箱的对称轴，正向指向储箱的顶端点，y_t轴为储箱水平轴，x_t轴与y_t轴及z_t轴构成右手坐标系；β为质心点处的质心面切向单位矢量e_τ与o_ty_t轴所成的夹角，e_n，j和e_n，k分别为质心点处的质心面法向单位矢量e_n沿o_ty_t轴和o_tz_t轴的分量，Y和Z为质心点在储箱本体坐标系下的坐标值，

F_τ＝μ|V_t|，m_f为液体的质量，α为有效质量系数，V_t为质心点相对质心面运动速度的大小，ρ为曲面半径，A_a为质心面在惯性坐标系下的加速度矢量，μ为摩擦系数。

本发明与现有技术相比的优点在于：本发明方法在集中质量近似思想的基础上，结合航天器储箱内液体晃动运动的特点，通过质心点与质心约束面之间的相互作用来确定液体晃动的作用力。本发明对液体晃动的幅度不作限制，与现有的基于单摆或弹簧-振子的等效力学模型相比，适用范围更广。特别是对于航天器在星球上着陆，在轨交会对接，以及在轨快速机动这类工况激起的液体晃动问题非常适用。

附图说明

图1为本发明方法的计算分析原理图；

图2为液体沿储箱壁运动时的作用力分布示意图；

图3为本发明实施例中的仿真结果对比图。

具体实施方式

本发明方法的主要思想是：将储箱内的液体等效为质量集中在质心处的质点(质心点)，质心点只能在质心约束面上运动，通过求解质心点与质心面之间的相互作用来计算液体晃动时对储箱的作用力和力矩。

如图1所示，首先建立两个坐标系，其中o_ix_iy_iz_i为惯性坐标系(地心坐标系)，该坐标系的原点o_i位于地心，z_i轴的正方向指向北极，o_ix_iy_i与z_i轴垂直；o_tx_ty_tz_t为固连在液体储箱上的本体坐标系，原点o_t位于储箱的几何中心，z_t轴为储箱的对称轴，正向指向储箱的顶端点，y_t轴为水平轴，x_t轴与y_t轴及z_t轴构成右手坐标系。储箱内液体的质心约束面应为一三维封闭曲面，但考虑到一般航天器储箱都是轴对称结构，为便于计算，本方法将其简化为二维。即质心面如图1中所示的虚线组成的椭圆，当贮箱内充有一定量的液体时，在常重条件下缓慢转动贮箱，通过解析的方法计算每一时刻液体的质心位置，液体质心可能经过的所有位置点便构成了二维质心约束面。

另外，图1中φ为储箱本体坐标系原点到质心点的连线与o_ty_t轴的夹角，而e_τ为质心点处的质心面切向单位矢量，e_n为质心点处的质心面法向单位矢量。定义β角为矢量e_τ与o_ty_t轴所成的夹角，θ_x为储箱的安装角(即储箱对称轴z_t轴与竖直方向的夹角)。F_τ和F_n分别为质心点对质心面沿切向及法向的作用力，用于计算液体晃动的作用力。

1、质心点运动方程

在惯性坐标系下，针对液体质心点建立力平衡方程，由牛顿第二定律可得：

$m_f\frac{{\mathrm{dV}}_f}{\mathrm{dt}}=F_n+F_{\mathrm{\τ}}---\left(1\right)$

其中，m_f表示液体的质量，V_f为质心点在惯性坐标系下的速度矢量，F_n和F_τ表示质心面对质心点的作用力分别沿曲面法向和切向的分量(见图1)。

定义V_a为储箱(即质心面)在惯性坐标系下的速度矢量，V_t为质心点相对质心面运动速度的大小，同时根据质心点只能在质心面上运动这一约束条件可知，质心点在每一瞬时的速度方向必须与质心面的切向一致，因此有：

V_f＝V_a+V_te_τ    (2)

将上式对时间进行求导并代入式(1)，即可得与V_t相关的力平衡方程：

$A_a+{\stackrel{\·}{V}}_t{e}_{\mathrm{\τ}}+\frac{V_t^2}{\mathrm{\ρ}}((e_n\×e_{\mathrm{\τ}})\×e_{\mathrm{\τ}})=\frac{1}{m_f}(F_n+F_{\mathrm{\τ}})---\left(3\right)$

其中，ρ为曲面半径，由定义可知

A_a为储箱在惯性坐标系下的加速度矢量，对于航天工程，由于储箱与航天器是固连的，因此A_a即为航天器相对惯性坐标系的运动加速度。

为了求解V_t，将上式两边分别与矢量e_τ点乘可得：

${\stackrel{\·}{V}}_t=-\frac{F_{\mathrm{\τ}}}{m_f}-(e_{\mathrm{\τ}}\·A_a)---\left(4\right)$

根据β角的定义(见图1)可知β＝φ+90°，则：

$\stackrel{\·}{\mathrm{\β}}=\stackrel{\·}{\mathrm{\φ}}=\frac{V_t}{\mathrm{\ρ}}---\left(5\right)$

同时式(4)可改写为：

${\stackrel{\·}{V}}_t=-\frac{F_{\mathrm{\τ}}}{m_f}-A_j\cos \mathrm{\β}-A_k\sin \mathrm{\β}---\left(6\right)$

其中A_j和A_k分别为储箱加速度A_a沿o_ty_t轴和o_tz_t轴方向的分量。式(5)和式(6)为质心点的运动方程，可通过数值积分求得每一时刻V_t和β的大小。

2、作用力计算

将式(3)两边分别与矢量e_n点乘，根据各矢量的定义及矢量运算法则可得法向力F_n为：

$F_n=m_f[\frac{V_t^2}{\mathrm{\ρ}}-(e_n\·A_a)]---\left(7\right)$

上式中，定义F_n沿曲面法向向内为正(见图1)。

切向力F_τ采用简化模型来计算，认为F_τ跟质心点与质心面的相对速度V_t的大小成正比，其方向与质心点的速度方向相反：

F_τ＝μ|V_t|；  μ≥0    (8)

摩擦系数μ的单位为N·s/m，其取值范围通常为0-100。

法向力F_n和切向力F_τ的方向随着质心点在质心面上运动而不断变化，为了便于分析，需将它们转换成液体对储箱的作用力。根据作用力与反作用力的关系，可投影得到液体晃动对储箱的作用力和力矩的计算公式，即：

F_y＝F_τ cos β+F_ne_n，j                             (9)

F_z＝F_τ sin β+F_ne_n，k                             (10)

M_x＝(F_τ sin β+F_ne_n，k)Y-(F_τ cos β+F_ne_n，j)Z    (11)

其中，F_y，F_z和M_x分别表示液体对储箱的沿o_ty_t轴和o_tz_t轴的晃动力以及绕o_tx_t轴的晃动力矩(逆时针方向为正)；e_n，j和e_n，k分别为单位法向矢量e_n沿o_ty_t轴和o_tz_t轴的分量，Y和Z为质心点在储箱本体坐标系下的坐标值。

令二维椭圆形质心面的曲面方程为：

aY²+cZ²＝1                                         (12)

a，c为曲面常数，与储箱几何形状及充液比有关，在实际计算时，需要通过解析的方法来提前确定。根据矢量e_τ的定义可知：

aY cos β+cZ sin β＝0                             (13)

由式(12)和(13)可计算得到：

$Z={\left[\frac{1}{\frac{c^2}{a}{\tan }^2\mathrm{\β}+c}\right]}^{1/2}---\left(14\right)$

$Y={\left[\frac{1-{\mathrm{cZ}}^2}{a}\right]}^{1/2}---\left(15\right)$

综上，根据外部加速度以及储箱和液体的相关参数，可计算质心曲面常数a、c，以及曲面半径ρ。在每一时刻，由式(7)和(8)可分别计算法向力和切向力的大小，然后对式(5)和(6)进行积分计算V_t和β的值，并更新质心点位置。

需要注意的是，在前面分析时将液体等效为一质点，由此会造成一定的误差。其原因如图2所示，将真实分布的液体划分为若干部分，则当液体沿储箱壁运动时，每一部分液体都会对储箱产生垂直于各自运动速度方向的作用力，即图2中的F_i，而F_cm表示将液体等效为质点时对应的作用力。由于各部分液体的运动方向不同，因而各部分作用力的方向也不同，而真实液体对储箱的作用力为各部分力的矢量和，因此显然有：∑F_i＜F_cm。

为补偿上述误差，引入液体有效质量系数α，对法向力计算公式中与运动速度相关的计算项进行修正，即式(7)变为：

$F_n=m_f[\mathrm{\α}\frac{V_t^2}{\mathrm{\ρ}}-(e_n\·A_a)]---\left(16\right)$

显然，有效质量系数α的大小与储箱的充液比有关，同时也受储箱几何外形参数的影响。对于外形相似的储箱，只需针对某一个储箱，在不同的充液比下进行实验测出液体晃动作用力，然后将仿真结果与实验结果进行对比，便可确定有效质量系数α与储箱充液比之间的关系，在后续工程应用时可依据该结果来确定α的取值。

为了验证本发明的方法，将计算结果与国外落塔实验中的结果进行对比。此实验于上世纪70年代利用Martin Marietta公司的落塔实验装置进行，储箱如图1所示为柱冠形，半径为6.35cm，总高16.25cm，充液比为25％，液体密度为1.905g/cm³，储箱及支撑装置的质量为0.8kg，储箱的安装角θ_x＝45°(见图1)。此时取α＝0.23，μ＝0.875N·s/m。

国外实验中给出的实验测量结果包括了储箱及支撑机构的惯性力，为了便于跟实验结果比较，在式(9)和式(10)中加上相应的惯性力项，即：

F_y＝F_τ cos β+F_ne_n，j-A_jm_s    (17)

F_z＝F_τ sin β+F_ne_n，k-A_km_s    (18)

其中，m_s为液体储箱和支撑机构的质量。

仿真计算时，根据实验测量出的储箱加速度，利用线性插值的方法给出各时间步的加速度大小。计算得到的晃动力F_y和F_z与对应的实验测量结果的比较如图3所示，其中图3a为晃动力F_y的结果，而图3b为晃动力F_z的结果，图中的实线为本方法仿真计算得到的液体晃动作用力，而点为实验测量得到的作用力。从图中可以看出，考虑液体有效质量系数后的仿真计算结果与实验结果吻合的非常好，表明本发明方法可以比较准确的计算出液体晃动的作用力和力矩。

本发明说明书中未作详细描述的内容属本领域技术人员的公知技术。

Claims (1)

1.一种确定储箱内液体晃动作用力的仿真方法，其特征在于：储箱内液体晃动所产生的沿储箱o_ty_t轴的作用力F_y采用公式F_y＝F_τ cos β+F_ne_n，j计算获得，储箱内液体晃动所产生的沿储箱o_tz_t轴的作用力F_z采用公式F_z＝F_τ sin β+F_ne_n，k计算获得，储箱内液体晃动所产生的沿绕储箱o_tx_t轴的以逆时针方向为正的晃动力矩M_x采用公式M_x＝(F_τ sin β+F_ne_n，k)Y-(F_τ cos β+F_ne_n，j)Z计算获得；其中o_tx_ty_tz_t为储箱本体坐标系，原点o_t位于储箱的几何中心，z_t轴为储箱的对称轴，正向指向储箱的顶端点，y_t轴为储箱水平轴，x_t轴与y_t轴及z_t轴构成右手坐标系；β为质心点处的质心面切向单位矢量e_τ与o_ty_t轴所成的夹角，e_n，j和e_n，k分别为质心点处的质心面法向单位矢量e_n沿o_ty_t轴和o_tz_t轴的分量，Y和Z为质心点在储箱本体坐标系下的坐标值，

F_τ＝μ|V_t|，m_f为液体的质量，α为有效质量系数，V_t为质心点相对质心面运动速度的大小，ρ为曲面半径，A_a为质心面在惯性坐标系下的加速度矢量，μ为摩擦系数。

Images