正弦激励下矩形容器内液体晃动特性研究

杨 琦，李德友，常 洪，王洪杰

（哈尔滨工业大学能源科学与工程学院，黑龙江 哈尔滨 150001）

液体晃荡是指在有限的空间内，带有自由表面的流体因受到外界因素激励或外加扰动而产生的液体运动。当外界激励的频率接近液体的固有频率时会引起共振，导致液体剧烈晃动，并对容器壁产生强烈的冲击，这是造成储液容器结构破坏的重要因素之一[1]，因此液体晃荡问题的研究对飞行器油箱、海洋核动力平台稳压器、航天液体推进器等储液系统的安全运行具有重要意义。载液容器在外界载荷激励下晃动时，容器内液体自由表面会产生不同的晃动波形，较常见的波形有驻波、行进波、水跃和组合波[2]，如图1 所示，其中图片的横轴与纵轴分别对应液面的宽度与高度。在外界激励频率远低于液体固有频率的晃荡中，波面只做周期性振动而不向前传播，此时的自由液面呈驻波形式；激励频率增大后，波面发生改变，会形成一系列波长较小的行进波，但其对容器壁面的冲击力要比驻波大很多；当激励频率增大到接近于液体固有频率时，对容器施加任一微小扰动都会使得自由液面出现跳跃，即产生了水跃现象；当载荷的振荡频率已经处于液体固有频率的微小邻域内时，共振现象发生，此时液体自由表面会形成幅度较大的驻波且伴随着行进波及水跃现象，这些波相互叠加形成了组合波。

图1 不同形式波形示意图[2]

有关液体晃荡问题的研究首先出现于航空航天领域，随着计算机技术和研究手段的快速发展，复杂的液体晃动问题吸引了大批国内外研究者投身其中。对于液体晃荡的研究方法大致可以归为理论分析、试验研究和数值模拟3 大类。理论分析是最早研究液体晃动问题的主要方法，以线性势流理论为基础，Abramson[3]对圆柱形、球形、矩形等规则几何形状容器内的液体晃动问题进行了研究，并且探讨了液体晃动产生的压力对贮箱结构所造成的影响。Chen 等[4]应用势流理论计算了车辆突然加速时考虑阻尼影响的液罐内自由液面的水动力特性。卢军等[5]基于势流理论对平放圆柱形贮箱内任意充液比的晃动问题进行研究，并用Galekin 方法对该边界问题进行了求解。理论分析方法通过建立描述流体运动的方程研究液体晃动问题，对简单模型的小幅晃动问题可准确地给出解析解，但当贮箱形状复杂或液体晃动幅度增大时，无法通过理论研究得到有效可靠的解析解。

实验方法一直是研究液体晃荡问题的重要手段。Wu 等[6]运用相似理论建立了模型试验台，研究了低液载率耦合激励下弹性罐内的晃动冲击载荷特性。Das 等[7]通过自由振动和强制振动试验对U 形容器作为高层建筑系统中的液体阻尼器的有效性进行了研究。Jin 等[8]以自来水和甘油为实验液体，通过对照实验探讨了流体黏度对中深度比水平激振引起的晃动的影响。Zhang[9]对倾斜底部调谐液体阻尼器的晃动模态特性进行了数值和实验研究。卫志军等[10]采用大尺度模型试验方法，对矩形舱内液体晃荡产生的冲击压力进行了研究，分析在不同外激振幅、频率和不同载液率下，液舱内晃荡冲击载荷特性。实验分析最能反映液体晃动的实际情况，但需投入大量的人力、物力和财力，实验周期较长，而且并不是所有物理量都能直接获得，此时对液体晃动问题进行数值模拟计算显得尤为重要。

数值分析方法有着理论分析和实验研究无法比拟的优势，不仅计算快、精度高、能获得大量的计算成果，而且计算成本低，能处理液体大幅晃动的非线性问题。Pandit 等[11]对部分填充斜底容器的晃动特性进行了数值研究，并采用Galer-kin 有限元方法推导了二维槽内液体晃动的运动模态方程。Guan 等[12]基于边界元法研究了三维折流板储罐的非线性晃动问题，系统地分析了挡板对三维罐内液体晃动的影响。Zhang 等[13−14]给出了影响挡板形状的因素及各种不同形状的挡板，通过数值计算方法探究最优的挡板形式。龙飞等[15]利用CFD方法仿真计算，研究了脉动风载荷作用下产生的液体晃动及外部激励共同作用于容器时的动力响应。翁羽等[16]使用CFD 和FEM 耦合的数值分析方法获得了部分充液容器在典型地震历程下的自由液面变化和结构响应特性。李骏等[17]利用FLUENT 有限元仿真软件对制动工况下货车燃油箱进行数值模拟，得到流固耦合作用下燃油箱的受力情况。李奇等[18]基于OpenFOAM 自行开发的水动力学求解器naoe-FOAM-SJTU 研究了船舶运动与液舱内液体晃荡间的相互耦合作用。

由此可见，数值分析方法已成为分析液体晃动问题的重要方法，而在数值模拟过程中，VOF 法是对自由液面的捕捉或追踪最常用的方法。尚春雨等[19]提出采用VOF 方法耦合k−ε湍流模型，通过UDF（用户自定义函数）加载振动加速度载荷的方法模拟液箱内液体晃动问题。刘岑凡等[20]采用同样的方法模拟了不同充液量下球罐内液体的晃动过程。此方法相当于将加速度值作为质量力反向加给流体，以实现容器运动条件下非惯性坐标系的转化，但多用于载液容器平动的情况，对于容器转动的情况，需要考虑转动牵连加速度、转动牵连向心加速度及科氏加速度，不容易理解且计算比较麻烦。

在液体晃动工程问题中，多数情况下液箱的位置及速度变化是可以确定的。本文采用VOF（Volume of Fluid）方法耦合湍流模型，通过UDF（User Define Function）直接加载液箱运动速度随时间的变化规律，给定液箱的刚体运动。该方法理解方便且计算简单，通过与实验结果的对比，验证了该方法的有效性。本文还使用该方法模拟了外界载荷激励下矩形载液容器的晃荡过程，得到容器内自由液面变化情况及液体对容器底部的冲击压力大小。通过改变外界激励载荷的频率与幅值，观察载液容器内液体的晃荡情况，分析各因素对液体晃动特性的影响，以期为储液系统的安全运行提供参考。

1 数值方法

1.1 VOF 法

本文采用VOF 法（流体体积法）来模拟带自由液面流体的运动。VOF 法通过研究网格单元中流体和网格的体积比函数F来确定自由液面，追踪流体变化[21]。其中流体体积函数F的控制微分方程为

利用连续方程∇ ·u=0，可得到其守恒形式：

该式决定了流体体积函数的输运过程。

定义函数F=F(x,y,z,t)表示区域内流体体积对计算区域体积的占有比率，即对于某个计算单元而言，F=1 表示单元被流体所充满；F=0 表示它是一个空单元；F介于0～1 之间，表示单元被流体部分充满，即自由边界存在于该单元中。通常，F=0.5 的网格会被提取出来作为气液两相的交界面[18]。F的梯度表示自由边界的法线方向，在得到各单元的计算值及其梯度后，根据F值和边界法向作一直线或平面切割单元，便可近似表示界面的位置，并以此边界设置边界条件[22]。VOF 法是一种非常流行的液面追踪方法，具有存储量小、适用于三维运算等优点，应用领域相当广泛，是目前处理液体晃动问题一种实用、稳定的方法。

1.2 运动规律设定

用户自定义函数（UDF）是一个在C 语言基础上来扩展FLUENT 特定功能的编程接口[19]。即当FLUENT 标准界面不能满足模拟的需求时，用户可以自己用C 语言编写扩展程序代码并与FLUENT动态链接，加载到FLUENT 环境中供其使用，以满足用户的各种需求。在研究载液容器内液体的晃荡问题时，我们经常会很容易地确定容器的位置及其速度变化规律。本文研究采用UDF 直接给定液箱刚体运动的方法（即通过UDF 提供每个时间步长的速度来指定Fluent 中某一特定区域的运动，Fluent 利用这些速度值来更新模型的节点位置）来模拟载液容器内液体晃动过程，观察其晃动情况。

1.3 矩形水箱固有频率

矩形液箱内液体晃动的固有频率可以由下式计算得到[7]：

式中：n为内部晃动模态数；L为液箱在运动方向上的长度；g为重力加速度；h为液箱的填充深度。

2 液体晃动仿真

2.1 研究对象

本文研究对象采用文献[23]中的调谐液体阻尼器（TLD）水箱模型，该容器的简化模型如图2 所示。水箱为矩形，水平放置，在外界非线性激励作用下周期性振动，水箱激励方向的长度为966 mm、宽为200 mm、高为360 mm，充液深度为119 mm。图3 为载液容器的计算模型及其网格划分。

图2 载液容器简化模型

图3 载液容器计算模型及其网格划分

为验证网格无关性，本文针对水箱模型划分了5 套不同数目的网格（41 万、55 万、77 万、93 万、108 万）进行计算，选取一个晃动周期内底面中心点处压力峰值为监测量。图4 为网格无关性验证结果。可以看出当节点数高于93 万时，压力幅值基本不再随网格数的增加而大幅变动，说明其受节点数的影响不大。综合考虑计算精度及计算时间等因素，决定选取节点数为108 万的网格进行后续计算。

图4 网格无关性验证

2.2 数值计算设定

本文利用ANSYS Fluent 软件，采用VOF 方法模拟液体在矩形容器中的晃荡过程。湍流模型使用k−ω模型封闭湍流黏度项。设置流动介质为空气和水，气液两相均匀连续介质，由于流体的压缩性对晃荡的影响可忽略不计，故采用不可压缩模型。求解器采用速度压强耦合求解器Coupled，压力修正方程的离散格式采用Body Force Weight 格式。容器壁面为不可滑移边界，不考虑传热的影响。非定常计算时间步长设置为0.004 s。在计算之前，分别设置非定常计算时间步长为0.002、0.004、0.008 s 进行无关性验证，选取了底面中心点处压力峰值作为监测量，验证结果如图5 所示。可以看出时间步长选为0.004 s 满足数值计算对计算速度与精度的要求。

图5 时间步长无关性验证结果

2.3 数值模拟算例与结果分析

为了验证本文所采用数值模型的准确性，选取了文献[23]中Tait 等的调谐液体阻尼器水箱振动实验，并在相同工况下对其进行数值计算，将本文所用数值模拟方法计算得到的结果与文献[23]中的实验结果进行比较。在0～60 s 内，载液容器以vx=A·fe·2πcos(2πfet)的正弦速度规律沿图2 所示x方向运动。其中，液箱振幅A=2.5 mm，振动频率fe=0.545 Hz。t=60 s 时刻容器突然静止，随后液体在惯性力的作用下自由晃动。

图6 为容器内液体晃动稳定后，距左侧壁48、193、388 mm 处液面随时间的变化曲线。从图6可以看出，数值计算结果与实验测量结果吻合良好，液面呈周期性变化，且晃动频率基本接近于外界激励频率。

图7 展示了一个晃动周期内不同时刻的液面演变过程。可以看出在该激励条件下，容器内液体的晃动幅度较小，液面以行进波的形式进行演变。随着容器的振动，液面峰值由最右侧壁面位置处移动到最左侧，接着又回到右侧。因此在一个晃动周期内，容器任何位置都会有液面波峰两次经过的情况，对应图6 中的液面变化曲线在一个周期内应该存在两个峰值，其中图6(c)最为明显，图6(b)隐约也有两波峰出现，只是193 mm 处更靠近容器左侧，液面峰值来回两次经过该处所需时间较短，且在该段时间内即使液面波峰已经移动到左侧壁处，由于距离很短（远小于液面波形的波长），该处液面并不会达到波谷状态，液面变化幅度很小，对应于图6(b)中液面曲线在该段时间内几乎是平的。同理，图6(a)中，距左侧壁48 mm 处双峰值现象就更加不明显，由于计算精度及图片展示的关系，使得我们观察到该处液面变化曲线似乎只有一个峰值。

图6 晃动稳定后不同位置处液面变化曲线[20]

图7 一个周期内不同时刻的液面演变

图8 为t=60 s 时刻容器突然静止后，液体在容器内晃动的自由衰减过程。从图8 可以看出：容器静止后，液体在惯性的作用下继续晃动；随着时间的增加，晃动幅值逐渐减小，最后趋于稳定。原因是受液体自身黏度及其与容器壁面之间的摩擦碰撞影响，液体晃动的能量被逐渐耗散，晃动幅度变小直至停止，而衰减频率同样接近于外界激励频率。

图8 容器静止后距左侧壁48 mm 处液面变化

本文在计算过程中监测液体晃动对容器壁的冲击压力大小，如图9 所示，设置D1、D2 两处压力监测点，分别为容器底面中心点及最右侧边中点。图10 为晃动稳定后，两监测点处的压力变化曲线。从图10 可以看出液体对容器底的冲击压力呈周期性变化，对于D1 点，在一个周期内压力值存在两个峰值，原因是液面峰值两次经过容器中心，该处液位高，对容器底部的压力就大；同理则D2 点压力值在一个周期内呈单峰值变化。

图9 压力监测点设置

图10 晃动稳定后监测点处压力变化曲线

图11 为D1、D2 两点压力的时域变化经FFT 变换后得到的频谱图。从图11 可以看出，当频率为0.55 Hz 时，压力幅值最大，即第一主频为0.55 Hz，说明在外界激励频率fe=0.545 Hz 时，容器内液体晃动频率为0.55 Hz，二者非常接近，相对误差仅为 ε=0.9%。图11(a)中1.09 Hz 为谐波频率，与第一主频呈两倍关系。

图11 监测点处压力频谱图

3 液体晃动仿真

3.1 不同激励频率下容器内液体晃动

对于本文所研究载液容器的晃动，填充深度h=119 mm，运动方向上液箱的长度L=966 mm，则载液容器一阶固有频率为0.566 Hz，显然上面计算中所选振动频率fe=0.545 Hz 已接近固有频率。为研究外界激励频率对液体晃动特性的影响规律，本文尝试改变激励频率大小（分别选取fe=0.272、1.09 Hz），对液箱模型进行数值模拟计算，观察液体的晃动情况。

图12、13 分别为不同频率下液体晃动过程中不同位置处液面变化对比及对容器底面冲击压力变化对比。可以看出，当外界激励频率（fe=0.272、1.09 Hz）远离载液容器的一阶固有频率（f=0.566 Hz）时，容器内液体的晃动幅度很小，自由液面仅在几毫米范围内变化，液体对容器底部的冲击压力也几乎为零，相比于外界激励频率fe=0.545 Hz 的情况，液体的晃动相当不明显。由此可见，外界激励频率对液体晃动特性的影响主要在于其接近载液容器固有频率的程度。当激励频率处于固有频率的微小邻域内时，将会发生共振现象，此时很小的激励幅值便可使容器内液体产生较大的晃动。

图12 不同频率下不同位置处液面变化对比

图13 不同频率下监测点处压力变化对比

图14 为不同外界激励频率下各监测点处的压力频谱图。从图14(a)(b)可以看出，当外界激励频率为0.272 Hz 时，液体晃动频率为0.27 Hz，接近于激励频率。而观察图14(c)(d)可以发现：当外界激励频率为1.09 Hz 时，第一主频为0.546 Hz，仅为激励频率的一半。分析其原因应该是，激励频率增大后，容器晃动过快，箱内液体的能量不足以支撑其跟随容器以较快的频率晃动，导致容器振动一个周期结束回到原处时，箱内液体的晃动刚进行到一半，因此液体晃动频率仅为激励频率的一半。为验证该分析，本文提取出激励频率为1.09 Hz 时D1 点处的压力变化曲线，如图15 所示。可以看出，此时液体晃动周期T=4.66–2.83=1.83 s，即晃动频率f=0.546 Hz，与FFT 分析结果一致。这说明当外界激励频率超过液体的固有频率时，载液容器内液体晃动频率不再接近于激励频率，而是会逐渐小于激励频率。

图14 不同频率下监测点处压力频谱图

图15 fe=1.09 Hz 时D1 点压力曲线

3.2 不同激励幅值下容器内液体晃动

从以上结果可以看出，虽然外界激励频率在固有频率附近，但液体晃动幅度仍较小。本文尝试增大外界激励幅值（取液箱振幅A=25 mm），观察载液容器内液体的晃动情况。图16 展示了增大振幅后载液容器晃动不同时刻的液面演变过程。可以看出液体晃动幅度明显增大，液面变化呈现组合波的形式，出现了翻转、水跃等现象，此时共振现象更明显。

图16 增大振幅后不同时刻的液面演变

图17 为不同振幅下容器内不同位置处液面变化对比。从图17 可以看出：振幅增大后，液面变化幅值增大，但晃动频率基本没有改变；相比于A=2.5 mm 时的情况，曲线变化趋势大致相同，但在某些时刻出现了小幅波动，如图17(a)中b–c 段所示。为分析其产生原因，在后处理过程中截取了容器Z方向中心截面处一个周期内不同时刻的液面演变过程，如图18 所示。从图18 可以看出，左侧壁附近液面在0.2T时刻达最高，之后液面逐渐降低（对应图17(a)中a–d 段），到0.9T时刻该处液面达最低。而在0.4T～0.6T时间段内，液面出现上下浮动的现象而非一直下降，此即为图17 中的b–c 段曲线出现小幅波动的原因。而该现象的产生应该是因为外界激励振幅增大后，液面以组合波的形式变化，即驻波、行进波和水跃等相互叠加的形式，因此即使在0.4T～0.6T时间段内液面最高峰已经移动到容器右侧位置，在左侧壁附近仍然会出现较小的波峰，使得该处液面上下浮动。图17(b)中的e–f 段与图17(c)中的g–h 段也是同样的道理。

图17 振幅增大后不同位置处液面变化对比

图18 z 方向中间截面上液面演变过程

图19 为不同振幅下D1、D2 两监测点处压力值变化对比。可以看出外界激励振幅增大后，液体晃动对容器底部的冲击压力也增大。这是因为液体晃动幅度增大，液面幅值增高，而静压与液位高度线性相关。从图17 也可以看出，振幅增大后，波峰处的液面高度增加，即液体堆积增多，对容器底部的压力也就增大；波谷处液面高度降低，液体量减少，对容器底部的冲击压力减小。此外，从图19 可以看出，压力变化曲线也出现了小幅波动的现象，经分析可知该现象同样是由该段时间内液面的上下浮动所引起。

图19 振幅增大后监测点处压力变化对比

图20 为振幅A=25 mm 时D1、D2 两点处压力频谱图。可以看出，外界激励幅值增大后，液体晃动变得剧烈，出现多个幅值较小的频率，该频率均为基频的整数倍，即谐波频率。此外，从图20(b)可以看出，液体晃动的频率为0.542 Hz，与外界激励频率0.545 Hz 相差很小，即增大激励幅值并未对晃动频率产生影响，说明液体晃动频率与外部激励幅值无关。

图20 振幅增大后监测点处压力频谱图

4 结论

本文以矩形载液容器为研究对象，通过对周期性边界激励下容器晃动过程进行数值模拟，验证了数值计算方法的准确性，分析了容器壁面受力及自由液面晃动与外部激励频率、激励幅值之间的变化关系，主要结论如下。

1）对静止的载液容器施加周期性边界激励，液体会出现较为强烈的晃动波，且晃动同样呈现周期性；撤掉外界激励后，液体晃动受阻尼作用影响逐渐衰减，而衰减频率接近于外界激励频率。

2）液体晃动频率取决于外界激励频率，与激励幅值无关。当激励频率小于或处于液体固有频率附近时，液体晃动频率接近于激励频率；当激励频率超过液体固有频率时，容器晃动过快而液体能量不足，导致晃动频率减慢，小于外界激励频率。

3）同一激励频率下，液体晃动幅度随激励幅值的增大而增大，液面也对应呈现不同的波形。当激励幅值较小时液面呈现为行进波，激励幅值增大时，液面转变为组合波的形式，伴随着翻转、水跃等现象。

4）当激励频率接近于液体固有频率时，共振发生，此时很小的激励便可使容器内液体产生较大的晃动；而当激励频率远离（远小于或大于）固有频率时，即使外界激励幅度很大，容器内液体晃动并不明显，液面几乎没有变化。这说明液体晃动可以吸收振动能量，达到减震的目的。该结论可用于液体阻尼器或减震器的研究中，对液体晃动的抑制同样具有参考意义。

本文阐明了外界激励各因素对矩形载液容器晃动的影响，可以为工程实际中一些近似周期性载荷激励情况下，液体晃动动态响应的预测提供有益参考，例如：海浪作用下船舶的周期性摆动、地震波作用下地面储液罐的晃动等。基于该数值方法可对三维复杂运动规律下液体晃动问题进行数值模拟研究。