含隔板球形贮箱液体晃动激励频率影响分析

马 亮， 刘 昊， 魏 承， 汤 亮， 赵 阳

(1.哈尔滨工业大学 航天工程系，哈尔滨 150001； 2.北京控制工程研究所，北京 100190)

光滑粒子流体动力学方法(Smoothed Particle Hydrodynamics,SPH)自1977年由Lucy等[1-2]首次提出，在经历了40年的发展和改进过程后，现已在航空航天[3-5]、船舶海洋[6-8]、公路运输[9-11]等多个领域取得了广泛应用。

在航天领域，随着现代航天器对高精度载荷及长寿命飞行的深入需求，航天器贮箱内液体晃动所产生的晃动力对系统动力学行为的影响已经成为总体设计重点考虑的问题。

基于SPH方法的流体动力学问题研究经历了逐步完善的过程。Monaghan[12]于1994年首次将SPH方法引入流体动力学问题的求解，随后多位学者基于此方法对钝体绕流[13]、剪切流动[14]、溃坝[15]等问题进行了仿真分析。同时，SPH方法在计算精度[16]、计算稳定性[17]及边界处理[18]等方面经过不断地改进和修正，使其成为一种解决流体动力学问题的成熟建模方法。

现代航天器贮箱多为球形结构，而针对球形贮箱的液体晃动问题，现有研究多采用等效模型法或有限元法：黄华等[19]通过建立三维质心面等效模型，利用质心动力学方程及接触碰撞模型对椭球形贮箱的晃动力进行求解；赵阳等[20]采用任意拉格朗日-欧拉法(Arbitrary Lagrange-Euler,ALE)，对常重力及微重力环境下球形贮箱内的液体晃动现象进行了仿真分析；王天舒等[21]提出适用于微重力环境的球形贮箱液体晃动弹簧-质量等效模型，同时利用虚功率原理对晃动力进行了求解；苗楠等[22]建立了静平衡面垂直于等效重力的等效模型，从而将液体大幅晃动分解为等效重力的整体运动及在此基础上的小幅晃动；岳宝增等[23]提出网格移动算法，在ALE方法基础上对球形贮箱大幅晃动进行了数值模拟。

由于航天器在执行任务过程中，伴随着频繁的轨道及姿态机动，使得航天器贮箱存在大范围运动或大载荷受力，从而导致贮箱内液体出现大幅度晃动。现有研究所采用的等效模型方法[24]或有限元方法均无法适用于大幅晃动下液体出现破碎的现象，而针对航天器球形贮箱的液体大幅晃动问题又亟待解决。

通过以上分析，本文开展了基于SPH方法的球形贮箱液体晃动问题研究：首先根据SPH方法的基本方程，给出N-S方程的求解形式，同时确定人工黏度项及压力项方程，并采用动态边界以提高边界计算精度；然后基于SPH方法建立球形贮箱液体晃动模型，采用两种激励模式及三种激励频率进行仿真分析，分别记录舱壁受力和观测点压强情况；最终设计一种带隔板的球形贮箱结构，并分析隔板对液体晃动的抑制作用。

1 基于SPH方法的N-S方程求解

光滑粒子流体动力学方法的基本思想是通过两步近似，将空间函数及空间函数导数表示为核函数支持域内各离散点处函数值与核函数及核函数导数乘积求和的形式。有利于简化空间函数及空间函数导数求解的难度，使得该方法适用于求解N-S方程。

由Gingold等的研究可知，在经历核近似、粒子近似两步近似过程后，空间函数f(x)在粒子i处的函数值可以表示为

(1)

式中:W为核函数;h为光滑长度;N为核函数支持域内粒子总数。

空间函数f(x)的导数在粒子i处的函数值可以表示为

(2)

核函数的支持域由光滑长度h及光滑因子κ决定，支持域半径为κh，如图1所示。本文采用三次样条核函数[25]，其具体表达式为

式中:R=|x-x′|/h;光滑因子取为κ=2。

图1 核函数支持域及粒子近似Fig.1 Support domain of kernel function and the particle approximation

采用SPH方法的基本方程，可将N-S方程组中的连续性方程、动量方程表示为

(4)

(5)

式中:vij=vi-vj;Wij=W(xi-xj,h);fe为单位质量外力;α和β为坐标轴方向。

(6)

在弱可压缩模型下[27]，压强项可表示为

(7)

式中:ρ0为流体参考密度；压强-密度刚性系数γ=7;声速cs取决于流体的最大速度，取为cs=10Vmax。

SPH方法存在边界问题：在距边界一个κh长度的范围内(包括边界处)会出现核函数支持域与问题域交错的现象。如图2所示，方形粒子代表边界粒子，三角形粒子代表发生截断的流体粒子，圆形粒子代表未发生截断的流体粒子。为使这一范围内的粒子计算保持精度，需对边界问题进行特殊处理。

为解决SPH方法的边界问题，本文采用动态边界法[28]，该方法的思想是直接构造与流体实粒子性质相同的边界虚粒子，在计算过程中将边界虚粒子代入计算，解决边界附近(一个κh长度范围内)流体粒子核函数支持域截断问题。但边界虚粒子自身参数不随计算结果发生变化，而是保持固定或随规划运动轨迹进行运动，该方法下虚粒子厚度与核函数支持域半径相关。

图2 支持域与问题域截断导致边界问题Fig.2 Boundary problem caused by the overlap of support domain and computational domain

2 仿真工况设计

采用SPH方法对N-S方程进行离散化，并给出人工黏度项、压力项的经验方程，针对边界问题提出动态边界法，从而实现对N-S方程的完整求解。在此基础上完成球形贮箱液体晃动仿真模型的建立，主要包括贮箱及液体建模、激励形式及工况设计。

2.1 贮箱及液体

贮箱设计分为无隔板模型及有隔板模型，两种结构的唯一区别在于隔板结构，其他参数均相同。球形贮箱半径r=0.125 m，液体深度h=0.1 m，隔板高度hb=0.05 m，隔板位于xoz平面内。为测量液体晃动对舱壁产生的压力，设置了4个测量点，测量点均位于yoz平面的舱壁面上，不同测量点的高度不同。其中，p2点高度与隔板高度相同，p3点高度与液面高度相同，具体位置如图3所示，贮箱与隔板均设计为刚体。

图3 贮箱及液体建模(m)Fig.3 Model of spherical tank and fluid(m)

2.2 激 励

由于航天器贮箱在实际飞行过程中存在大范围运动或大载荷受力，仿真设计两种类型激励模式：力模式、运动模式。力模式的加载对象为液体，控制量为液体所受加速度；运动模式的加载对象为贮箱，控制量为贮箱位移。两种模式下均采用简谐激励形式：x=Asin(2πf·t+φ)，力模式下x表示液体y方向所受加速度，运动模式下x表示贮箱y方向位移，仿真的重力环境选为常重力。设计6组仿真工况，激励模式及简谐运动参数如表1所示。

表1 仿真工况设计Tab.1 Design of simulation and parameters

3 激励频率影响分析

图4所示为两种贮箱模型下，4个时刻点处晃动液体的侧视图。每个时刻下有隔板贮箱的液体晃动幅度均小于无隔板情况，为了定量分析隔板影响，分别从舱壁受力、测点压强两方面进行对比分析。

3.1 舱壁受力

力模式激励下液体的加载方向为y方向，液体在yoz平面的晃动行为明显，因此主要分析舱壁y,z两方向的受力情况。

舱壁y方向初始受力为零，随着液体受到简谐激励，开始出现周期性晃动力。不同的激励频率造成舱壁受力不同的曲线形式，当激励频率f=1.5 Hz时晃动力幅度达到最大，而三种激励频率下隔板的存在均使得晃动幅度减小，如图5(a)和图5(b)所示。

舱壁z方向初始受力为负值，这是由液体重力造成的，随着液体受到简谐激励，同样出现周期性晃动力。不同的激励频率对舱壁z方向晃动力影响巨大，且同样在激励频率f=1.5 Hz时晃动力幅度达到最大，隔板的存在同样使得晃动幅度减小，如图5(c)和图5(d)所示。可以看出，由于存在液体晃动现象，舱壁在y,z两方向的受力均增大，且为周期性受力。

表2给出了各工况下，舱壁受力稳定后的平均晃动幅值，体现出的规律与曲线图一致，同时可以看出隔板对激励频率大的工况晃动抑制效果更明显。

图4 晃动过程中各时刻下液面形态Fig.4 Shape of the fluid at different moment in sloshing

图5 舱壁受力随振动频率变化情况Fig.5 Force applied on bulkhead varies with vibrational frequency

表2 舱壁受力平均晃动幅值Tab.2 Average force amplitude applied on bulkhead

3.2 测点压强

运动模式激励下贮箱的加载方向同样为y方向，液体在yoz平面的晃动行为明显，因此在yoz平面设置了4个压力测量点：p1点高度低于隔板高度，p2点高度与隔板高度相同，p3点高度与液面高度相同，p4点高度高于液面。

对于p1点，在无隔板工况下，不同的激励频率使得压力曲线呈现出截然不同的形式：当f=1 Hz时，晃动行为与激励形式一致，表现为简谐波动；当f=1.5 Hz时，出现“双波峰”现象，前压力峰值略高于后压力峰值，这是由液体撞击舱壁后沿球形内轮廓上升后回落所形成；当f=2 Hz时，也会出现“双波峰”现象，但两次压力峰值相差较大，这是由于液体撞击舱壁后的上升幅度较小，没有形成翻卷、破碎，如图6(a)所示。在有隔板工况下，三种激励频率下的测点压力值均略有下降。当f=1 Hz时，压力曲线形式同样表现为简谐波动；当f=1.5 Hz时，“双波峰”的前后压力峰值更加接近；当f=2 Hz时，后压力峰值反高于前压力峰值，如图6(b)所示。这是由于p1点位于隔板下方，隔板的存在会对液体首次冲击舱壁存在阻碍作用，而当液体回落时隔板同样阻碍了液体的回流，正是这种作用使得前压力峰值有所减小而后压力峰值明显增大。但隔板的存在使得压力峰值达到“平均化”效果，这样不会使舱壁局部在某时刻的压力过大，而导致结构变形或损坏。

对于p2点，其高度高于p1点，但仍位于液面之下，因此在有无隔板及三种激励频率下的压力曲线形式与p1点都相似，但局部现象存在不同之处。在无隔板及激励频率为1.5 Hz,2 Hz的工况下，没有出现明显的“双波峰”，而变为前波峰加“后压力平面”的组合,如图6(c)所示。这是由于p2点高度增加，液体不存在回落或回落幅度极小；在有隔板及激励频率为1.5 Hz,2 Hz的工况下，前压力峰值高于后压力峰值，如图6(d)所示。这是由于p2点高度与隔板高度相同，隔板对液体首次冲击舱壁及液体回落回流的阻碍作用都减小，但隔板的存在同样使得测点压力值有所下降。

对于p3点，其高度位于液面处，压力曲线形式与前两个位于液面下的测点存在较大差异。在三种激励频率下压力曲线均呈现“单波峰”加零压力组合。当不存在隔板时，激励频率越大压力峰值越大,如图6(e)所示；当存在隔板时，随着激励频率的增加隔板对测点处压力的抑制效果越明显，如图6(f)所示。

对于p4点，其高度位于液面之上，压力曲线形式与p3点相似。在激励频率为1 Hz工况下已检测不到压力值，而在激励频率为1.5 Hz,2 Hz的工况下，压力曲线也变为“脉冲波”加零压力组合。隔板的存在同样对激励频率越大的工况测点处压力抑制越明显,如图6(g)和图6(h)所示。

图6 测点压强随振动频率变化情况Fig.6 Pressure of measuring point varies with vibrational frequency

4 结 论

通过建立有无隔板两种球形贮箱模型，针对两种激励模式及三种激励频率进行仿真分析，最终得到以下结论：

(1) 力模式激励下，无论贮箱有无隔板均在激励频率f=1.5 Hz时，得到舱壁受力的最大峰值，而隔板的存在会减弱舱壁所受晃动力。

(2) 力模式激励下，由于液体晃动现象的产生，舱壁在y,z两方向的受力均增大，且为周期性受力。

(3) 运动模式激励下，位于液面下的舱壁测点p1,p2表现出一致性，液体冲击舱壁及回落作用导致“双波峰”现象的出现，测点高度的不同又会使得该现象在两测点处存在差异。而隔板的存在对舱壁所受压力具有“平均化”效果。

(4) 运动模式激励下，位于液面处及液面以上的舱壁测点p3,p4表现出一致性，出现“单波峰”或“脉冲波”加零压力组合形式，且隔板的存在对激励频率越大的工况抑制效果越明显。