唐浩等

摘 要：随着航空航天技术的发展，研究液体晃动显得越来越重要。液体晃动对飞行器燃油供给、飞行状态以及机体结构变形产生较大影响。如何准确模拟液面变化、计算液体对壁面的冲击力是目前存在的难点。针对这些问题，简要介绍了SPH方法的基本原理，对导弹油箱内燃油的晃动进行了仿真计算，模拟出液面变化。计算了液体重心位置的变化曲线，同时得到油箱与燃油之间的相互作用力，并对结果做出了分析。对导弹油箱内燃油晃动研究提供了一定的理论和工程参考价值。

关键词：光滑粒子动力学；燃油晃动；数值仿真；接触力

贮箱内液体的晃动是指液体自由表面由于受到外加扰动或激励而产生的运动。液体的晃动问题在航空航天、船舶、石油化工以及核动力等领域较为常见，从上世纪50年代起，各国工程师对这一问题已经给予了广泛关注[1]。储箱内带自由液面的液体晃动不仅会对贮箱壁产生冲击载荷，引起结构的疲劳破坏。对于载液比重大的载液系统，如航天器、大型民机以及油罐车等，在运动过程中，短时间内液体的大幅晃动引发的液体重心变化还会对这些载液系统的运动轨迹或飞行姿态产生不利的影响。依据相应的要求，油箱的几何尺寸应避免在飞行中使燃油产生剧烈晃动，使燃油晃动造成的重心变化在可接受范围内；由于飞机姿态变化所带来的燃油重心移动也处于可以接受的范畴[2]。在航空领域，飞机在起飞、着陆及飞行过程中，因气动紊流、发动机转子的不平衡，武器的发射等原因，会引起油箱内燃油的晃动。研究晃动产生的冲击力、晃动力矩所引发的油箱结构强度问题、液体晃动重心的变化对全机稳定性的影响是十分必要的[3]。

液体晃动是一种特殊的流动现象，拥有自由表面流动的相关特点，同时还与约束液体运动的壁面有着较为复杂的相互作用，此类问题的难点就在于对液面大幅晃动的准确模拟以及对液体与壁面相互作用的处理如何准确模拟液面大幅晃动以及处理[4]。因此，对液体晃动问题的求解极具挑战性。

光滑粒子力学（Smoothed Particle Hydrodynamic， SPH）方法广泛用于求解流体动力学问题，可以很好的适应自由边界的变化，并且对于自由表面大变形流动问题的处理效果也比较好。对于液体的边界条件处理，本文考虑了自由液面以及液体与固体壁面两种情况。自由边界问题，文章采用Dirichlet方法处理，根据表面粒子密度来判断特定粒子是否位于自由边界上。而对于液体与固体壁面的相互作用的区域，则在固体壁面上布置一组虚粒子，用来对内部粒子产生排斥力，从而对粒子与固体壁面之间的相互作用进行模拟。

文章基于SPH方法数值模拟导弹油箱受到激励时液体的大幅晃动，模拟了燃油液面变化，计算了燃油晃动过程中重心坐标的变化以及油箱和燃油相互作用力的变化，并对结果作了分析。

1 理论简介

1.1 SPH方法基本原理

SPH的理论基础是粒子方法[5-7]，使用一系列离散点或者粒子来描述连续的系统，系统中的离散的粒子将承载所有的物理信息。

对于粒子i的任意一个函数值都可以应用光滑函数对其支持域内所有粒子相应的函数值进行加权平均后获得，支持域是一个半径为kh的圆，k是确定粒子影响半径的系数，如图1所示，其中S是影响区域的边界。

通过粒子近似式，可以将连续积分形式下的函数转换成在支持域的粒子的离散化求和，通过这种方式实现无网格化。

1.2 核函数

本文研究选择的核函数为Monghan和Lattanzio提出的三次样条函数[8]，也称B样条函数。

为了满足一致性要求，上式中的系数αd在一维、二维和三维空间中分别取为 和 。

1.3 控制方程

拉格朗日描述下的粘性流体控制方程就是Navier-Stokes方程。

用希腊字母上标α、β表示坐标系的方向，可以通过张量理论中的指标法用来表示方程的叠加，并且可以在变化的拉格朗日空间内引入对总时间的导数，将N-S方程作粒子近似，并进行空间离散，得到粒子近似的N-S方程。

连续方程：

动量守恒方程：

式中，σ为总应力张量，由各向同性压力p和粘性应力 τ两部分组成，v为速度，ρ为密度，x为距离，σ为总应力张量。

1.4 边界条件

在液体和固体壁面的相互作用的区域内，液体不可以穿透贮箱壁，同时也不能脱离贮箱壁存在。因此在固壁边界上布置一组虚粒子，对邻近粒子施加排斥力，这样可以防止内部粒子穿透固壁边界，如图2所示。排斥力的方向和两个粒子的中心连线相同，排斥力的大小可以根据公式（5）计算得到[4]。

式中的n1和n2分别取值为12和4；r0是排斥力的作用范围，其值为离散后粒子的初始间距；D为常数，大小与最大速度的平方相近。

处理自由表面边界，则采用Dirichlet边界条件[11]，根据粒子i的粒子密度来判断是否位于自由边界，若粒子 密度满足式（6），那么认定其在自由表面上，并使该粒子密度等于流体表面的实际密度。

式中，ρ0是自由表面流体的实际密度，参数β的大小在0.85-0.98之间，本文将其取值为0.95。

2 模型的建立

2.1 模型说明

贮箱结构取自某型导弹的油箱[9-10]，几何模型如图3所示，左图为贮箱结构图，油箱的长约为830mm，宽约为650mm，高度约为560mm。油箱内部液体为煤油，液体体积约为0.148m3。

为方便将油箱内液体离散成均匀的SPH粒子，采用体积等效的方式，将油箱中的燃油用体积与其相等的形状简单的几何体代替。计算得到油箱内燃油的体积约为0.148m3，为保证体积等效的几何体完全包含于油箱内部，设计如图4所示长方体形状的水柱，燃油划分为5460个SPH单元，油箱划分为壳单元，不考虑油箱的变形，系统采用右手坐标系。

2.2 计算工况说明

油箱晃动的俯仰角度随时间变化的曲线，如图5所示。考虑到燃油从图4所示的液柱状态运动到平稳状态需要一段时间才能稳定，因此采用图5的角度激励曲线，初始的0-3.5s内俯仰角度为0，液柱在重力作用下运动。3.5s后燃油达到基本稳定，俯仰角度从负方向的-45度线性变化到正方向的45度，变化时间为3s，6.5s后俯仰角度重新变为0。在计算过程中，竖直方向始终作用地面的重力加速度G，值取9.8m/s2。

2.3 仿真流程说明

文章的仿真计算基于PAM Crash软件平台。完成了模型相关参数的设置后，将其导入Crash Pam Solvers求解器中进行计算，时间步长取0.0001s，仿真计算10s后停止。

3 仿真计算结果

3.1 油箱内液面的变化

在不同时刻油箱内燃油晃动的液面变化情况如图6所示，其中左侧为文章采用的数值方法计算得到的液面位置。

前3.5s燃油在重力作用下趋于稳定，液面平稳。3.5s-7.5s时，箱体做俯仰运动，液面发生波动，SPH方法能够比较真实地模拟出液面随油箱运动的变化情况，晃动过程中产生的波浪也能比较准确的模拟出来。

3.2 燃油的重心位置变化

油箱中燃油在晃动过程中，粒子的运动是相当无序的。由于液体已经离散成一系列有质量的SPH单元，因此重心位置可以根据划分得到的SPH单元加权求和得到。

燃油重心在X、Z三个方向的变化曲线，如图7所示。由于模型是关于XZ平面对称，同时贮箱系统也在XZ平面内运动，因此，重心在Y方向的变化不大，主要表现为X方向和Z方向的运动。

从以上重心位置的时间历程曲线可知，0-3.5s内燃油从液柱状态下落运动到稳定状态的过程中，X方向的坐标基本保持不变，Z方向坐标也基本保持不变。在3.5-6.5s内，油箱在外界激励下运动，使得液体X方向的坐标有较大的变化，近似于一个与所给激励类似的正弦变化，Z方向的坐标有一个比较小波动变化。6.5s之后，油箱停止运动，由于是绕Y轴的俯仰运动，箱体也对称于XZ平面，因此停止运动后，重心在X和Z方向坐标基本无变化。

3.3 油箱与燃油相互作用力变化

油箱内燃油晃动过程中，燃油对油箱有不同程度的冲击作用力，而且两者之间的相互作用力是不断变化的。由于对称面XZ两侧燃油对油箱的冲击作用力可以相互抵消，因此Y方向的作用力相对很小。X方向和Z方向相互作用力较大，变化曲线如图8所示。

图中可以看出，前3.5s内，燃油液柱运动到稳定状态的过程中，X方向的冲击力除了一开始液体下落时一直在很小的量级波动，而Z方向的冲击力逐渐增大。晃动过程中，X方向作用力波动明显，Y方向变化始终较小，Z方向也有较小波动。6.5s后X方向趋于0，Z方向逐渐减小到1130N左右，接近燃油重力。

4 结语

文章从贮箱内液体晃动问题研究的必要性出发，讨论了研究的意义，简要阐述了SPH方法的基本理论。随后针对某型导弹油箱中的燃油晃动问题展开了仿真计算和分析。

通过文章的相关研究工作可以得出如下结论：

⑴油箱内燃油晃动的结果较好，并且晃动过程中的波浪也能比较准确地描述出来，对油箱内防晃板的布置优化和燃油供给方式的选择提供一定参考。

⑵计算得到了液体晃动过程中重心坐标的变化，有利于对飞行器飞行姿态进行准确的控制调整。

⑶由于油箱燃油体积为0.148m3，通过计算可知重力为1138N，与仿真得到的稳定值1130N差别很小，验证了仿真的可靠性。

文章的研究可以应用在航天器液体推进剂储箱晃荡或者大型LNG船储箱液体晃荡问题，对飞机燃油晃动等研究也可以提供一定的参考价值。

[参考文献]

[1]刘富.贮箱内液体晃动动力学分析及结构防晃技术研究[D].南京：南京航空航天大学，2010.

[2]冷飞.基于SPH方法的飞机油箱燃油晃荡研究[D].南京：南京航空航天大学，2009.

[3]夏恒新，王为，宝音贺西，等.带球形底的圆柱容器中液体小幅晃动的实验研究和有限元分析[J].动力学与控制学报，2009，7（1）：66-70.

[4]孙云舫.SPH方法的改进及其在水波数值模拟中的应用[D].天津：天津大学机械工程学院，2009.

[5]Jiannong Fang，Aurèle Parriaux，Martin Rentschler，Christophe Ancey.Improved SPH methods for simulating free surface flows[J].Applied Numerical Mathematics，2009（59）：251-271..

[6]Sudarshan Tiwari，J?rg Kuhnert.Modeling of two-phase flows with surface tension by finite pointset[J].Journal of Computational and Applied Mathematics，2007（203）：376–386.

[7]Youn Kyung Song，Kuang-An Chang，Yonguk Ryu，Sun Hong Kwon.Experimental study on flow kinematics and impact pressure in liquid sloshing[J].Exp Fluids，2013（54）：1592-1611.

[8]Monaghan J J，Lattanzio J C.A refined particle method for astrophysical problems[J].Astronomy and Astrophysics， 1985，149（1）：135-143.

[9]Takehiro Himeno，Toshinori Watanabe，Satoshi Nonaka，Yoshihiro Naruo，Yoshifumi Inatani，Hiroshi Aoki.Numerical and experimental investigation on sloshing in rocket tanks with damping devices[J].AIAA，2007-5557：1-17.

[10]Sudarshan Tiwari，Axel Klar，Steffen Hardt，Alexander Donkov.Coupled solution of the Boltzmann and Navier–Stokes equations in gas–liquid[J].Computers & Fluids， 2013，10（71）：283-296.

[11]Koshizuka S，Nobe A，Oka Y.Numerical analysis of breaking waves using the moving particle semi-implicit method[J]. International Journal for Numerical Methods in Fluid，1998，26（7）：751-769.