航天器弹性充液贮箱双向流固耦合响应特性数值分析

杨鹏飞， 陈 刚， 薛 杰， 王 珺

(1.西安交通大学 航天航空学院，西安 710049； 2.气动噪声与控制重点实验室，四川 绵阳 621000；3.西安航天动力研究所，西安 710100)

液体晃动问题是指液体自由表面由于受到外加激励而产生的运动，这种问题广泛存在于工程实际应用的各个领域，如航空航天、石油化工等。就航天器贮箱来说，其受力及振动环境复杂，微小的外激励可能会造成内部液体剧烈晃动，引起附加惯性质量分布和流-固耦合刚度改变，严重时会影响结构振动模态。特别是当储液容器结构很薄时，内部液体晃动可能与贮箱弹性模态耦合，产生不可预料和控制的振动，这对航天器的安全性和稳定性产生了极大影响。

目前对液体晃动问题的研究也越来越受到研究者们重视。近几十年来，国内外学者对贮箱液体晃动问题已展开了深入的研究。Biswal等[1]在不同的圆环形阻尼板的尺寸、位置及数量情况下，通过有限元方法研究了液体晃动特性。Cho等[2]同样在有阻尼板作用下采用了有限元方法研究了二维矩形容器的液体晃动特性，探讨了不同阻尼板的尺寸及安装位置对晃动特性的影响。郝亚娟等[3]采用单一拉格朗日法,对理想流体贮箱中矩形隔层板的耦合振动问题进行了研究,探讨了上下层流体的高度、隔层板的几何尺寸、流体密度、流体晃动频率等因素对隔层板动力特性的影响。王琼瑶等[4]针对部分充液罐车转向或制动过程中罐体内液体的晃动问题，研究了防波板的几何参数(开孔的大小、形状及位置等)对液体瞬态晃动时的载荷转移及晃动力的影响。Chen等[5]利用改进 的SPH法模拟了二维水箱液面晃动问题，得到了液体自由液面形状和壁面处最大压力，并与实验结果进行了对比验证。程选生等[6]基于势流理论在微幅晃动范围内研究矩形贮液结构的动力响应，结果表明在微幅晃动范围内，液体最大晃动波高和地震加速度幅值成线性关系。张展博等[7]研究了低频大幅晃动条件下水平多孔挡板不同浸没深度对液面晃动的影响，发现了波峰和波谷变化的规律以及在特殊情况下挡板导致波峰高度增加的现象。

对于考虑流固耦合效应的液体晃动问题，Strand等[8]研究了具有柔性壁面的二维矩形水箱，分析了液体晃动与柔性壁面间的耦合关系。刘小民等[9]针对铁路罐车刹车过程中的液体晃动进行了研究，分析了不同材料下防波板受液体晃动冲击时的受力情况 ，同时揭示了充液比对罐体壁面及防波板受力的影响。得出了罐体端面受力随充液比的增加而增大。吴作伟等[10]以高速动车组车下水箱为研究对象，采用双向流固耦合方法，对车载水箱箱体内液体冲击晃动问题进行了建模与数值分析，结果表明:箱体所受最大应力基本与充液比成正比 ,且最大应力位置基本相同。姚灿等[11]针对薄壁贮箱的流固耦合问题，在理论分析基础上，采用有限元数值方法进行贮箱在充液不加压和充液加压两种工况下的模态分析，得到流固耦合下贮箱的含阻尼非线性固有频率。Yang等[12]采用双向流固耦合方法，研究了三维储罐的液面分布特征和容器的结构响应，结果表明了储罐结构的响应振幅随着地震的增加而增大，并主要受液体晃动的影响。

上述研究主要探讨了对挡板的几何参数与液体晃动的相互影响关系，大多没有采用双向流固耦合方法，仅有的少数采用双向流固耦合方法的研究也没有分析结构柔性对液体晃动和结构应力分布的影响。本文采用双向流固耦合方法，获取液贮箱液体晃动波高以及贮箱结构应力随时间的变化历程，研究了结构柔性对液体晃动和结构应力分布的瞬态响应的影响。分析了贮箱的材料、厚度以及箱内液体的充液比对弹性贮箱液体晃动及结构受力特性的影响规律和流固耦合效应，从而为利用结构柔性进行弹性贮箱液体晃动抑制提供设计依据和方法支撑。

1 数值方法

1.1 流动控制方程

假设箱内液体无旋、不可压缩，流体模型采用层流模型，忽略温度的影响，则贮箱内液体流动满足三个基本方程：

1.1.1 质量守恒方程

单位时间内流体微元中质量的增量，等于同一时间间隔内流入该微元的净质量。将质量守恒方程写成微分形式，表示如下：

(1)

式中:散度divu=∂u1/∂x+∂u2/∂y+∂u3/∂z,t、u、ρ分别为时间、速度矢量、流体密度，u1、u2和u3是速度矢量u在各方向上的分量。

1.1.2 动量守恒方程

在单位时间内流场域流体微元体中动量增加的量，必须等于同一单位时间内出入该微元体流体的净动量和净冲量的和，其中边界表面上的表面力和整个体积区域上的体积力会对微元体产生冲量。写成微分形式，并利用连续性方程化简得：

(2)

1.1.3 能量守恒方程

单位时间内流体微元中能量的增量，等于同一时间间隔内流入该微元的能量与外力做功的和。其微分形式为：

(3)

e表示单位质量流体的内能，体积力F=(Fx,Fy,Fz)。

1.2 结构动力学方程

由经典力学理论，动力学平衡方程表示如下：

(4)

1.3 流固耦合基本方程

综合固体与流体控制方程，流固耦合问题的有限元控制方程可以整理为下式:

(5)

式中，

式中：q表示来自流体的变量，fs表示流体结构的相互作用。

1.4 双向流固耦合策略

在很多耦合问题中，流体的作用力影响结构的变形，同时结构的位移又影响流场的形态，我们把这种类型的分析叫做“双向耦合”。双向流固耦合计算分为流体域和固体域两部分，流体和固体计算通过耦合界面传递数据，不考虑温度的影响，流固耦合界面满足运动学条件和动力学条件，即：

(6)

式中：df和ds分别表示流体和结构的位移，τf和τs分别为流体和结构的应力。

对于上述方程，当给定初始条件及边界条件后，可以通过建立耦合方程的通用形式进行统一求解，目前，主要有两种求解方法：直接耦合解法和分离法。本文采用直接耦合解法，将流体和固体的控制方程耦合至同一方程矩阵求解，方程表述如下：

(7)

求解过程概括如下：首先，将总计算时间划分为多个时间步，在每个时间步内，求解耦合系统的线性方程组并更新得到的解。计算应力和位移残量并和给定的迭代容差比较。如果解还不收敛，且没有到达流固耦合迭代的最大次数，则回到第1步继续下一个迭代。随后，该时间步的结果将作为下一时间步的初值进行计算，继而完成所有时间步的求解，如图 1 所示。

图1 双向流固耦合方法

2 计算模型及条件

2.1 计算模型

本文研究对象的结构尺寸如图2所示，贮箱的圆柱部分的长度为142 mm，箱体球形端半径为213 mm ，隔板位于圆柱端与球形端的交接处，圆柱端外部设有法兰结构，具体位置和尺寸如图所示。箱体厚度2 mm,法兰结构厚度5 mm。

图2 计算模型及结构尺寸示意图

对箱体和流体分开建模，建模时保持壳单元和流体单元接触面处节点一一对应。在几何空间位置上，流场和结构场的流固耦合边界应该重合，由于建模时需要考虑薄壁特点，选用壳单元来建立罐壁模型；设罐体外一点与罐体进行刚性连接，其方向为x轴正方向，罐体重力加载到x轴负方向；沿罐体轴向为y方向。流体域的网格数为45 080，固体域的网格数为5 958。贮箱整体结构(左)和流体(右)的有限元模型如图3。

图3 贮箱结构及流体网格

2.2 载荷加载

采用输入加速度时程曲线的方法对贮箱结构施加振动，该方法加载简单，且易于实现。选取的振动时程曲线如图4。定义时间步长为0.002 s,共计算500步，计算总时间1 s，加速度加载到y方向。

图4 加速度时程曲线

3 结果分析

3.1 不同贮箱材料对流固耦合响应的影响

本部分在充液比为0.5,箱内液体为酒精，箱体厚度为2 mm的情况下，采用双向流固耦合方法研究了不同材料的贮箱模型对箱内液体微幅晃动及结构应力变化的影响。箱体材料及其相应的材料属性，如表1所示。

表1 贮箱材料及其属性

图5为不同材料的贮箱模型液体晃动波高时程曲线。从图中可以看出，贮箱内液体晃动呈周期性变化，随着时间的增加，液体晃动愈来愈剧烈，在0.688 s时达到峰值。

图5 液体晃动波高时程曲线

图6为不同材料的贮箱在液体晃动波高峰值时刻云图，在0.688 s时刻，结构钢材料贮箱模型晃动波高最大，为18.23 mm；钛合金材料贮箱模型晃动波高次之，为14.28 mm;铝合金材料贮箱模型晃动波高最小，为10.85 mm。

图6 液体晃动波高峰值时刻云图

可以发现，随着材料杨氏弹性模量和密度的减小，贮箱内液体晃动波高峰值减小，在外加激励不变的情况下，晃动幅度逐渐趋于平缓。

图7为不同材料的贮箱模型等效应力随时间变化曲线，从图中可以看出，箱体等效应力呈不规则的周期振荡，随着材料杨氏模量和密度的增加，箱体所受的等效应力增大，且振荡幅度增加。

图7 贮箱等效应力时程曲线

图8为不同材料的贮箱模型应力峰值时刻云图，从图中可以看出，等效应力最大值均位于箱体靠近法兰位置，在0.004 s时刻，铝合金材料贮箱最大等效应力为523 599 Pa,在0.216 s时，钛合金材料贮箱和结构钢材料贮箱等效应力达到最大值，分别为584 793 Pa和698 608 Pa。可以发现，随着材料杨氏模量和密度的增加，贮箱模型等效应力增加。

图8 贮箱等效应力峰值时刻云图

3.2 贮箱厚度对流固耦合响应的影响

本部分采用双向流固耦合方法对2 mm、3 mm和4 mm厚度的贮箱液体晃动及其箱体结构所受应力的影响进行了对比研究。贮箱箱体材料为钛合金，箱内充酒精。

图9为不同厚度的贮箱模型液体晃动波高时程曲线。如图所示，贮箱内液体晃动呈周期性变化，随着时间的增加，液体晃动愈来愈剧烈，在达到峰值后晃动幅度逐渐减小。在0.688 s时刻，3种贮箱模型晃动波高均达到峰值，2 mm贮箱模型晃动波高最小，为14.28 mm; 3 mm贮箱模型晃动波高较大，为17.12 mm;在0.686 s时，4 mm贮箱模型晃动波高最大，为19.03 mm。可以发现，随着贮箱厚度的增加，液体晃动幅度增大。

图9 液体晃动波高时程曲线

图10为不同材料的贮箱模型应力峰值时刻云图。从图中可以看出，在0.216 s时刻，充2 mm厚度贮箱最大等效应力为584 793 Pa；在0.268 s时，3 mm厚度贮箱最大等效应力为418 996 Pa，4 mm厚度贮箱最大等效应力为342 455 Pa。可以发现，随着贮箱厚度的增加，贮箱模型等效应力减小。

图10 贮箱等效应力时程曲线

3.3 不同充液比对贮箱流固耦合响应的影响

本部分采用双向流固耦合方法研究不同充液比对箱内液体晃动和箱体的影响。仍然选取钛合金作为贮箱箱体材料，箱体厚度2 mm,贮箱内充入液体为酒精。

图11为不同充液比的贮箱模型液体晃动波高时程曲线。从图中可以看出，贮箱内液体晃动呈周期性变化，随着时间的增加，液体晃动愈来愈剧烈，在达到峰值后晃动幅度逐渐减小；在0.716 s时刻，充液比为0.2贮箱峰值波高最大，为18.91 mm；在0.688 s时，充液比为0.5贮箱峰值波高次之，为14.28 mm;在0.6 s时，充液比为0.8贮箱峰值波高最小，为9.743 mm。随着贮箱充液比的增加，液体晃动幅度减小，且晃动峰值时刻提前。

图11 液体晃动波高时程曲线

图12为不同材料的贮箱模型应力峰值时刻云图。从图中可以看出，在0.216 s时刻，充液比为0.2贮箱最大等效应力为352 640 Pa；在0.216 s时，充液比为0.5贮箱最大等效应力为584 793 Pa；在0.004 s时，充液比为0.8贮箱最大等效应力为806 876 Pa。可以发现，随着充液比的增加，贮箱模型等效应力增大。

4 结 论

本文针对薄壁贮箱的流固耦合问题，采用双向流固耦合方法对考虑流固耦合效应的充液贮箱进行了瞬态响应分析。通过改变贮箱的材料、充液比以及贮箱厚度，研究了贮箱结构柔性对液体晃动及结构受力特性的影响。得到以下结论：

(1)减小贮箱材料的杨氏模量和密度，箱内液体的晃动幅度减小，同时贮箱结构等效应力减小。可以发现，适当增加模型的柔性，可以有效抑制箱内液体晃动，同时减小液体晃动对箱体造成的强度破坏。

(2)增加贮箱结构的厚度，箱内液体晃动幅度增大，但贮箱结构受到的应力减小。

(3)通过改变贮箱内液体的充液比，可以发现，对本文所研究的贮箱模型，随着充液比的增加，箱内液体的晃动幅度减小，但箱体结构所受的应力增大。