微重力环境下部分充液贮箱内液体晃动特性分析林柯成;宋晓娟;吕书锋【摘要】针对微重力环境下充液航天器液体燃料晃动的动态特性为研究对象,采用液体晃动等效力学模型对充液航天器贮箱内的晃动液体进行等效,推导出微重力环境下贮箱内液体的边界条件,得出贮箱所受力和力矩的计算公式,并计算出了力和力矩的大小.应用有限元数值仿真法计算出贮箱内液体晃动的力和力矩,通过与解析法得到的结果进行比对,验证了数值模拟算法的可靠性.【期刊名称】《内蒙古工业大学学报（自然科学版）》【年(卷),期】2018(037)006【总页数】7页(P424-430)【关键词】液体晃动;等效模型;有限单元法;晃动力;晃动力矩【作者】林柯成;宋晓娟;吕书锋【作者单位】内蒙古工业大学机械工程学院,内蒙古呼和浩特 010051;内蒙古工业大学机械工程学院,内蒙古呼和浩特 010051;内蒙古工业大学理学院,内蒙古呼和浩特 010051【正文语种】中文【中图分类】U297.91充液航天器在执行航天任务的过程中需携带大量液体燃料推进剂.航天器在由发射轨道进入预定轨道的过程中,由于受到火箭动力、天体引力和大气阻力的影响,液体燃料会不可避免地发生晃动.微重力作用下液体晃动特性比常重力下更为复杂.航天器在执行航天任务的过程中,贮箱内的液体推进剂会不断被消耗为航天器提供持续动力.贮箱内的充液比会逐渐降低,液体的晃动特性也会跟着改变.近年来,大量学者对液体晃动的研究取得了很多新进展.文献[1]通过拉格朗日方程得到航天器刚体部分运动和液体燃料晃动的耦合动力学方程,提出了复合控制器,并通过数值模拟来验证控制器的精度和效率.控制器分别用于控制单燃料腔和四燃料腔航天器的轨道机动和姿态机动.文献[2]提出了一种新型网格移动算法,使用任意拉格朗日-欧拉有限元方法对球形贮箱中的三维液体大幅晃动进行数值模拟.模拟结果与实验数据吻合度很高.文献[3]针对四储箱充液航天器进行姿态机动控制研究.结果表明四储箱的三种空间布局对航天器姿态机动过程中的角速度、液体晃动力矩和控制力矩有不同影响.文献[4]建立了一种新的等效模型,首次提出液体的静平衡表面垂直于等效重力方向的假设,将液体的大幅运动分解为跟随等效重力的整体运动和在此基础上的小幅晃动.文献[5]借助有限元软件求解容器内三维液体晃动的固有频率和模态,对圆柱形容器内液体晃动的相关结果进行了总结,得出了该类型容器内液体晃动的一般规律.文献[6]创建了飞机机身油箱的有限元模型,壁板采用四节点四边形壳单元,箱内的液体域采用八节点六面体单元,单元算法采用多物质ALE方法,分析了液体最大晃动幅度和挡板长高比之间的关系.文献[7]用MSC-PATRAN创建了圆柱形贮箱的有限元模型,研究了贮箱半径和充液比对液体固有频率的影响.文献[8]改进了ALE有限元网格更新方法,使流-固耦合的分析结果更接近理论值.文献[9]将圆柱形贮箱内的液体晃动等效为弹簧-阻尼模型,通过有限元分析得到弹簧刚度系数和阻尼系数随充液部分高宽比的变化规律.圆柱形贮箱是形状最简单的贮箱,也是航天器上采用最多的贮箱.本文首先将圆柱形贮箱内的液体晃动等效为弹簧-质量模型,用弹簧和质量块的各参数表示液体对贮箱的作用力和力矩,之后选用ADINA和Fluent模拟该算例,验证数值模拟仿真的正确性和可靠性.为复杂贮箱及实际工程应用提供一定的理论基础.1 液体的边界条件假设圆柱形贮箱内的液体是理想液体,不可压缩,无黏度,无旋,则满足如下条件ρ=C(1)μ=0(2)▽·v=0(3)其中ρ为液体的密度,C表示常数,μ为液体的粘度,v为液体的速度.在自由晃动的条件下,液体的晃动速度可表示v=▽φ(4)其中φ为速度势,▽为Nabla算子,满足▽(5)液体的边界条件可表示为▽2φ=0(6)(7)定义波高函数z=η(x,y,t),则速度势和波高在自由液面上可有如下关系在z=η(x,y,t)(8)速度势函数和波高函数也可以写成下面的形式(9)(10)其中i、j分别为径向和周向的半波数,qij(t)是时间函数,φij(x,y,z)和Hij(x,y,z)是空间函数,其应满足如下方程▽2φij=0(11)(12)2 液体的晃动方程航天器在轨运行时,由于液体质量大,惯性大,因此在推进力的作用下,其晃动幅度通常不超过贮箱半径的15%,可视为小幅晃动,速度势满足方程(11)和(12).液体做小幅晃动时,根据文献[10],可将其等效为弹簧-质量模型,如图1所示.第n阶晃动质量距离贮箱底的高度可由式(13)计算.(13)图1 液体晃动等效弹簧-质量模型Fig.1 Spring-mass equivalent model of liquid-sloshing其中,hn为第n阶液体晃动质心距离贮箱底面的高度,h为充液高度,R为贮箱底面半径,λn是一阶Bessel函数导数的第n个根(见表1).表1 前6阶λn值Tab.1 Values of first 6λni123456λn1.8413.0543.8324.2015.3185.331将液体小幅晃动等效为弹簧-质量模型,等效原则为:(1)等效系统的总质量等于液体晃动的质量.(2)等效系统的各阶固有频率等于液体的各阶固有频率.(3)等效系统的力和力矩等于液体晃动受到的力和力矩.等效原则第(1)条用公式表示为(14)其中,m0为不晃动的液体质量,mn为第n阶液体晃动的质量,n为晃动阶数.等效原则第(2)条用公式表示为(15)其中,kn为第n阶液体晃动的弹簧刚度,ωn为液体晃动的第n阶固有频率,g为重力加速度.对于刚性圆柱形贮箱,当重力沿波高的方向时,各阶液体晃动的固有频率可由式(16)求得(16)各阶液体晃动质量可由式(17)求得(17)其中,ml为液体总质量.单个弹簧-质量系统的动力学方程为(18)其中,等号左面前二项为平动部分,后二项为转动部分.ηn为第n阶弹簧伸缩位移. 单个弹簧-质量系统对贮箱的反作用力为knxn,因此液体对贮箱的合作用力F可以表示如下:(19)其中,rn为第n阶晃动质量的回转半径,x0为不晃动液体的位移,xn为第n阶晃动液体的位移.将坐标原点定为贮箱底面圆心,在考虑转动的情况下,每个弹簧质子作用在原点的力矩Mn为(20)等效系统作用于贮箱上的合力矩为(21)若贮箱只做平动,则和转动相关的项均为0,动力学方程(18)可进一步简化为(22)合力方程(19)可简化为(23)力矩方程(21)可简化为(24)当液体做小幅晃动时,液体的质心位置可近似认为不变,因此有(25)当贮箱的速度为正弦函数时,其位移为余弦形函数,由式(22)可解得(26)(27)将式(26)、(27)带入式(23)、(24)中,求得等效系统作用于贮箱上的合力为(28)合力矩表达式为(29)3 仿真模拟用ADINA创建刚性圆柱形贮箱的有限元模型.贮箱内部底面半径为1m,高为4m.贮箱壁材料为铝合金,壁厚0.01m.液体推进剂为水,充液比为50%.贮箱所受重力加速度为0.1m/s2,方向沿z轴负方向.提取前6阶模态振型和固有频率.如图2所示为晃动液体燃料的前6阶振型.图2 液体晃动的前6阶振型Fig.2 First 6 vibration shapes of liquid-sloshing 将计算出来的各阶固有频率及晃动质量与式(16)、(17)求得的解析解相比较,比较结果如表2、表3所示.可见ADINA的仿真解与解析解吻合得较好,证明仿真法可靠.表2 各阶固有频率(单位:rad/s)Tab.2 Values of each level natural frequency (Unit: rad/s)第1阶第2阶第3阶第4阶第5阶第6阶仿真解0.4300.5540.6230.6520.7360.739解析解0.4290.5530.6190.6480.7260.730 表3 各阶晃动质量(单位:kg)Tab.3 Values of each level sloshing mass (Unit: kg)第1阶第2阶第3阶第4阶第5阶第6阶仿真解1423247119894342解析解1426247120904343提取液体晃动的其他参数,如波高、速度、应力等.用Fluent模拟液体的晃动,提取波高、速度、贮箱壁所受的力、力矩等变量.根据万有引力公式和黄金替换公式,求得g=0.1m/s时速度的平均值为2513m/s.假设贮箱受到的推进力F0为余弦力,有效值为106N,频率f为50Hz,则速度函数为v=0.63sin(100πt)+2513.本例利用有限单元法模拟5个完整的受力周期液面的晃动情况,将各变量绘制成曲线图,如图3～8所示.分析结果如下:(1)由图3、图4可知,左壁面液面高度先上升,右壁面液面高度先下降.这是因为贮箱运动方向向右,刚开始时液体由于惯性向左运动,液面向左倾斜.两侧壁液面高度变化周期和激振力周期相等.随着时间的推移,两侧壁液面高度均呈增大趋势,总体上是左侧壁大于右侧壁.而中轴线液面高度几乎不变.这是由于液体的浸润性,使得液面呈凹槽形.图3 贮箱两侧壁波高Fig.3 Wave heights on two side walls图4 贮箱中轴线波高Fig.4 Wave height on center line(2)液体晃动速度由图5所示,晃动速度总体上随着时间的推移而减小,最后趋于稳定.速度大致呈正弦形分布,周期和激振周期相当.正向速度大于负向速度,液体质心呈升高趋势.(3)图6、图7为液体对贮箱的作用力,在水平方向上比垂直方向上大2个数量级.这是因为贮箱所受的激振力沿水平方向.贮箱所受合力最大值为843727N.(4)图8为液体对贮箱的力矩,贮箱所受力矩的最大值为854650N·m.取前6阶液体晃动,将表2、表3中的解析解带入式(25)和式(26),可求得合力的最大值为855895N,力矩的最大值为816086N·m.与解析解相比,合力的误差为1.4%,力矩的误差为4.1%.仿真法的计算结果与解析解相比十分接近,误差小于5%,证明数值模拟方法的可靠性与正确性.图5 贮箱液体质心的垂直速度Fig.5 Velocity of liquid’s mass cente r图6 液体对贮箱的水平方向作用力Fig.6 Force on horizontal direction图7 液体对贮箱垂直方向作用力Fig.7 Force on vertical direction图8 液体作用在贮箱上相对于原点的力矩Fig.8 Torque on tank relative to origin point4 结论本文以在轨运行中携带部分充液圆柱形贮箱的充液航天器为研究对象,将晃动液体等效为弹簧-质量模型,推导了圆柱形贮箱内液体晃动的边界条件,以及液体对贮箱的力和力矩计算式.针对平动算例,用软件对解析公式进行了验证.数值模拟的结果与解析公式计算的结果吻合较好,验证数值仿真模拟的可靠性,为后续分析形状复杂的贮箱内的液体晃动特性提供了一定的参考价值.参考文献:【相关文献】[1] 岳宝增,于嘉瑞,吴文军.多液腔航天器刚液耦合动力学与复合控制[J].力学学报,2017,49(2):390～396.[2] 岳宝增,唐勇.球形贮箱中三维液体大幅晃动数值模拟[J].宇航学报,2016,37(12):1279～1284.[3] 邓明乐,岳宝增,黄华.液体大幅晃动类等效力学模型研究[J].宇航学报,2016,37(6):631～638.[4] 苗楠,李俊峰,王天舒.横向激励下液体大幅晃动建模分析[J].宇航学报,2016,37(3):268～274.[5] 贾善坡,许成祥.储液容器内液体自由晃动的有限元分析[J].船舶力学,2012,16(1～2):21～26.[6] 杨瑞.基于ALE有限元法的飞机整体油箱燃油晃动特性研究[D].哈尔滨:哈尔滨工业大学,2016.[7] 赵剑波.有加速度时轴对称贮箱内液体晃动的数值模拟[D].大连:大连理工大学,2008.[8] 周宏.液体晃动数值模拟及刚一液耦合动力学研究[D].北京:清华大学,2008.[9] M.Farid,N.Levy,O.V.Gendelman.Vibration Mitigration in Partically Liquid-filled Vessel Using Passive Energy Absorbers[J].Journal of Sound and Vibration,2017,406:51～73. [10] H.Abramson.The Dynamic Behavior of Liquids in Moving Containers[M].Washington DC,USA: NASA,1966.106.。
