液体晃荡现象在工程的各个领域广泛存在，从广义的范畴讲它是指容器内介质自由表面波动引起的介质运动。除空载和满载外，可以说任何一个载液系统在一定条件下都可能发生晃荡现象，较为常见的载液系统有油罐车的储油罐、货轮的液舱以及飞机的油箱等。

对于飞机结构起飞、降落以及做各种机动时，油箱里的燃油冲击油箱壁，对箱体结构和飞行安全造成威胁。显然飞机的运动状态变化是导致机翼整体油箱中燃油发生晃荡的直接原因，在滚转^[1]、俯仰^[2]以及半滚倒转机动^[3]下，燃油在油箱内运动非常剧烈，强烈冲击箱体结构，在飞机设计^[4]中有相关规范要求。飞机降落是事故高发段，不同于飞机机动下的油箱晃荡，着陆时油箱的运动是惯性力引起的，主要运动形式为垂直方向的垂荡运动，因运动幅度较小未得到关注，但仍可能对箱体结构的安全性构成威胁。

理论解析方法和试验方法^[5-7]早期获得了许多研究成果。但随着载荷激励越来越复杂、容器几何越来越不规则，理论和试验方法逐渐被飞速发展的计算机模拟技术所替代。在油箱晃荡问题中，不同的数值计算方法被应用，贾善坡等^[8]应用泛函极值原理建立了任意充液容器内液体晃荡的有限元分析模型，确定了充液容器内液体固有频率、模态及地震波对液体晃荡的影响，对比理论解，证明了该方法的有效性。杨尚霖等^[3]以半滚倒转机动行为下的某型飞机整体油箱为研究对象，基于VOF建立ABAQUS与Star-CCM+联合仿真流固耦合方法进行分析，得到了晃荡燃油的形态变化，流体域压力以及油箱部件的应力时程曲线。杨瑞^[9]通过ALE有限元法准确地模拟出飞机整体油箱中的液体晃荡，得到了油箱运动、充液率和挡板对液体晃荡的影响。黄愉太^[10]采用ABAQUS有限元软件对比了SPH和CEL 2种方法在弹性体贮箱充液晃荡模型中的结果后，选择了更适合模拟大变形的SPH方法对复杂飞机整体油箱模型下的液体晃荡进行了分析。这些研究主要针对液体晃荡，未关注运动幅较小的垂荡运动，研究液体的垂荡运动对计算方法和计算精度提出了更高的要求。

本文以机翼油箱垂荡为背景研究，先利用基于粒子的流体动力学仿真软件nanoFluidX计算油箱模型中的流体运动，然后提取运动流体作用于油箱的受力，再利用有限元软件ABAQUS对油箱壁进行应力分析，获得了一套精度高、计算效率高的分析方法，可为飞机油箱晃荡问题的研究提供参考。

1 研究方法 1.1 模型与边界条件

通常机翼整体油箱由多个翼盒组成，在飞机降落并与地面接触的瞬间，油箱的垂向晃荡是其主要运动形式，燃油受迫垂向运动，冲击各个翼盒的箱体结构。在这一过程中，几何尺寸最大的翼盒下蒙皮所受弯矩最大，最易发生变形，所以选择几何如图 1所示的矩形油箱进行研究。

飞机着陆时通常垂直方向的下降率不超过5 m/s，以保证人体和飞机结构的安全。在载荷和结构惯性力的作用下，整体油箱的运动导致箱内燃油受迫晃荡，燃油晃荡带来的冲击效应和燃油的充液率有关。着陆瞬间整体油箱的主要运动形式为垂荡，在某型号飞机的一次着陆过程，测得起落架开始与地面接触的0.56 s内，机翼整体油箱处加速度变化如图 2所示。假设油箱初始静止，按图 2的加速度运动，其速度及位移变化如图 3所示。

从位移变化图中可以看出机翼整体油箱在着陆瞬间的运动过程，由于惯性作用在前段时间内开始下沉，油箱位移在0.3 s达到最大，待下沉达到最低点后又开始向上回弹，0.5 s回到初始位置并继续向上运动。

1.2 SPH-FEM联合仿真介绍

SPH方法通过流体动力学仿真软件nanoFluidX实现。在nanoFluidX中，油箱视为刚体，计算得到晃荡下箱体的压力分布。然后将压力分布作为载荷条件，输入到有限元软件ABAQUS中，利用FEM方法分析箱体受到的应力。

nanoFluidX是Altair研发的基于SPH方法的流体动力学仿真工具，相较于其他求解器有以下优势。在包含不规则几何的计算中，前处理不需要花费过多精力设计网格，导入几何模型，选择单元生成粒子即可。有不同类型的运动选项，来规定实体的运动，能够在实体单元与周围流体相互作用时，较准确获取结构所受的力和力矩。采用了GPU并行技术，可节省了计算时间。

以Delorme等^[11-12]进行的油箱晃荡试验为参照，验证nanoFluidX在油箱晃荡数值模拟中对油箱受力计算的准确性。试验中矩形液箱尺寸为90 cm×5 cm×58 cm(长×宽×高)，充液高度为9.3 cm。液箱绕其底边中轴进行一种正弦形式的周期旋转运动，角位移，振幅θ_max=4°，周期T=1.92 s。在nanoFluidX计算中，液箱视为刚体，加入水和空气，进行3个周期的计算。

模拟得到液箱内水的晃荡情况与Delorme等^[11-12]的试验照片进行对比(如图 4所示)，可见波行进时两者的形状相同。随着液箱晃荡，水夹杂着空气与箱壁碰撞后，波向相反的方向行进，在波的行进过程中先发生水跃现象，从而导致了射流以及波的破碎。由于试验照片拍摄角度以及数值模拟中使用的水和空气相关参数与真实试验的差异，波形存在一定差异。

在Delorme等的试验中，通过装在液箱内壁初始液面高度处的压力传感器，测量3个周期内水的晃荡对箱壁产生的压力。图 5为该位置动态压力试验结果和数值模拟结果的对比，可见模拟结果与试验结果吻合良好，冲击压力出现和结束的时间完全一致。由于波浪形态具有一定的差异及压力测量位置不完全一致，压力值存在较小偏差，但仍可证明nanoFluidX在模拟中压力计算的准确性，特别是对压力多次峰值的模拟与试验吻合良好。

1.3 空气相关性验证

如果机翼整体油箱内没有气体进入，随着泵将油吸入引擎，将会在油箱内形成真空，以致发生不输油等故障，严重危害飞行安全。因此，为了保证飞机燃油供应系统的正常运行，通常机翼油箱有通气孔，或通过压气机等增压装置产生增压空气。即机翼整体油箱在着陆发生垂荡时，箱内介质除燃油外还有空气的存在。

使用简化的飞机整体油箱刚体模型，粒子间距取6 mm，壁厚为5层粒子，总粒子数为287 663，油的材料属性见表 1。分别在充液率为25%, 50%, 75%情况下，将箱内剩余空间中充满空气，在1.1节中的实测油箱垂向加速度激励下，通过nanoFluidX进行计算。图 6显示了加入空气后对油箱下壁板受力产生的影响。可以看到加入空气前后，2条曲线大部分吻合较好，只在峰值处存在一定差异，峰值处的最大误差为18%，所有峰值的平均误差约为8%。所以可认为，此加速度激励下空气的影响可以忽略，本文后续计算忽略空气影响，以提高计算效率。如果发生剧烈冲击现象，空气的影响就不容忽视。

2 着陆瞬间油箱垂荡数值模拟及分析 2.1 模拟结果分析

取充液率为50%，油粒子数为388 275，燃油的相关参数见表 1。图 7为模拟过程中，油箱内燃油的形态，以及不同时刻燃油的速度。在整个过程中，燃油整体相较于箱体没有较大的位移，只是在某些时刻，油箱侧壁的燃油会沿壁向上运动一小段距离，并且由图 8所示燃油在整个运动过程中平均速度基本和油箱保持一致。燃油重心与油箱下壁板距离的时间变化历程见图 9，在0.027 s有最大距离为96.39 mm，0.419 s有最小距离95.63 mm，可见燃油重心相较于下壁板的距离最大变化量为0.76 mm，燃油在油箱中的相对位置基本上没有发生变化。

燃油晃荡同时又对油箱壁面产生冲击，油箱下壁板受燃油作用力，作用力峰值时程曲线见图 10，0~0.25 s油箱加速度变化的峰峰值较小，下壁板受力峰值都小于1 500 N，0.25~0.45 s因油箱加速度变化峰峰值增大，受力峰值上升到2 000 N附近，0.286 s时达到最大值2 315 N，0.45~0.576 s油箱加速度峰峰值再次减小，所以下壁板受力峰值基本上又都降到1 500 N以下。

为了探究飞机着陆瞬间充液率对油箱垂荡的影响，进一步分析充液率为25%, 75%, 95%情况下油箱晃荡情况。不同充液率下，燃油SPH粒子间距是相同的，燃油单元的属性和参数不变。图 11为不同充液率下的油箱下壁板受力时间历程图，由图可见高充液率下受力峰值出现的次数要小于低充液率，也就是说在相同垂荡加速度下，高充液率的燃油相较于低充液率的燃油更为稳定。不同充液率相应时刻载荷峰值的大小，随着充液率的增加而增大，并且整体上充液率越大，峰值出现的时间越晚。在4种充液率下，受力峰值较大的时段都处在加速度变化较为剧烈的阶段(0.25~0.45 s)。

2.2 模拟结果分析

本文所研究的是油箱垂荡问题，油箱只有垂直方向的位移，假设箱体侧壁与燃油之间没有相互作用力，燃油只受到来自箱体的竖直方向的载荷，且作用在燃油质心。可做以下等效简化：①将燃油等效为集中在初始质心位置的质量体，其质量M=, ρ是燃油密度，Σ为燃油域, V是燃油的体积。②该质量体与燃油具有相同的动量，，是质量体的速度，v表示燃油速度(时间和位置的函数)。③则该质量体的运动方程为ρ= ，g为重力加速度，f为表面力。根据假设燃油与侧壁没有相互作用力，表面力f只是箱体下壁板对其竖直方向的力，所以可将油箱下壁板对燃油的约束简化为固定在质心正上方与质量体连接的细绳，下壁板对燃油作用下的表面力f等效于细绳的拉力T，T最小为0，质量体只有竖直方向的位移。

通过简化的等效模型可以定性分析2.1节中实测加速度下油箱内燃油的运动状态。在等效模型中，质量体只受到重力和拉力，运动方程为T-Mg= ，因拉力T≥0，所以a≥-g。当a=-g时，质量体不受拉力T作用。油箱加速度突然发生变化时，燃油相对油箱的运动取决于其是否有能力和油箱运动状态保持一致。在油箱垂荡的过程中，根据油箱运动状态可分以下3种情况来讨论(油箱和油的初始速度相同并处于稳定状态):

1) 当油箱加速度A≥0时，质量体加速度a能取到和A相同的值，此时T≠0;

2) 当油箱加速度-g≤A＜0时，质量体同样可以和油箱有相同的运动状态;

3) 但当油箱加速度A＜g时，燃油加速度最大为-g，无法和油箱运动状态一致。此时T=0，质量体和油箱发生相对位移。

综合等效模型的3种运动情况，得出结论为油箱向下的加速度值大于g时，拉力T为0，且和下壁板的距离会逐渐增加。

结合2.1节实测得到的加速度变化图，向下加速度的最大值为9 m/s²，小于g，所以在0.56 s内燃油始终没有和油箱下壁板发生分离，燃油相对油箱发生线性垂荡，这与数值模拟结果相符。飞机着陆时机翼油箱处的加速度很大程度上取决于飞机垂直方向的速度，即下降率越大，机翼油箱处加速度也越大，有可能会超过9.8 m/s²。这种加速度下的燃油晃荡将不同于本文所采用的实测加速度，燃油会发生伴随砰击现象的非线性垂荡。

综合数值计算结果总结下壁板所受合力随实测加速度和充液率的变化关系。垂晃时燃油的运动方程为，在问题域内积分得到燃油与下壁板接触处的表面力，油箱下壁板受到的压力和燃油的表面力互为相反力，且燃油密度保持不变，则T为，其中S表示油箱下壁板，V₀表示油箱总体积，k为燃油充液率，燃油域Σ=V₀k。

当充液率为25%和50%时，燃油总量相对较少，燃油和油箱的运动形式保持一致，所以取为油箱的实测加速度，即计算所用的加速度边界条件，将加速度值代入并与计算结果进行对比(图 13)，2条曲线峰值大小及曲线走势极为吻合。但充液率为75%和95%时，燃油总量较多，黏滞性和可压缩性相对明显，不能用油箱加速度代替燃油加速度，但仍可用上式来计算加速度不超过9.8 m/s²时线性垂荡情况下下壁板的压力。

2.3 表面应力分布

将nanoFluidX油箱垂荡数值模拟结果中蒙皮的压力分布，作为ABAQUS应力计算的加载条件，通过结构响应对结构的安全性进一步评估。取充液率50%情况下0.286 s时刚体模型中下壁板的最大压力分布，如图 14a)所示。飞机机翼整体油箱下壁板为铝制机翼蒙皮，厚度为1.12 mm，长宽同样取为1 000 mm×450 mm，材料属性见表 2。机翼蒙皮在该压力作用下应力的分布如图 14b)所示，在长边中点处有最大值224.1 MPa，小于铝板的屈服强度。说明在充油率为50%时，实测油箱加速度下，单次箱内燃油晃荡不能对箱体结构造成显著破坏。

在充液率75%的情况下，0.403 s油箱下壁板受力取得最大值3 767 N，相应的蒙皮应力分布最大值为346.3 MPa，超过了铝的屈服强度325 MPa。即充液率为75%时，着陆瞬间造成的油箱晃荡会对机翼下蒙皮造成破坏。由此得出结论，在本文所用的实测加速度下，充液率较低时，单次着陆造成的油箱晃荡不会对整体油箱的结构安全产生威胁；而当充液率不低于75%时，着陆时油箱垂荡，导致机翼下蒙皮的应力超过了安全值，对油箱结构构成威胁。

本文选择的模型虽然比较简单，但建立的nanoFluidX与ABAQUS联合仿真方法对于复杂情况仍然适用，包括线性非线性加载晃荡、箱体结构复杂曲面等。针对这些复杂情况，本文所建立的仿真方法同样可得到较为准确的数值计算结果。

3 结论

本文介绍了一种高精度的油箱垂晃数值模拟方法，能够准确计算幅值较小的加速度激励下，油箱结构所受的冲击压力。结合理论和数值结果分析，得出以下结论：

1) 通过对油箱垂荡进行数值仿真，验证了本文建立的nanoFluidX与ABAQUS联合仿真方法的精确性，指出此方法可适用于复杂几何及不同加载等复杂情况。

2) 针对某实际着陆载荷，发现充液率越大燃油的黏滞性越明显，下壁板受燃油作用力的峰值也越大，但峰值出现时间会推迟。当充液率为0.75及以上时，燃油晃荡造成了较薄的机翼下蒙皮应力超过材料屈服强度、发生明显变形的现象，很好地解释了工程中出现的问题。

3) 针对油箱垂荡提出的单摆模型，可以分析油箱加速度对箱内燃油运动的影响，定性判断了箱内燃油的晃荡形态，确定了加速度是否超过9.8 m/s²为非线性垂荡和线性垂荡的判断标准。

4) 对于加速度不超过9.8 m/s²的实测机翼油箱加速度垂荡载荷下，总结了下壁板受燃油作用力T与充液率k以及燃油速度v之间的变化关系T=，可用于估算此类情况下壁板受力的情况。