基于VOF法的球罐内液体晃荡数值模拟及载荷计算

刘岑凡 杨智荣 于哲敏 王国军

(1.中国特种设备检测研究院 北京 100029)(2.大连理工大学 大连 116024)

特种设备中许多压力容器是液体的储存装置，如各种储罐、球罐和槽车等。当发生地震或其他载荷激励时[1，2，3]，这些设备会对时间相关的加速度载荷产生响应，并引起其盛装的液体介质晃荡，液体随即又以动压力的形式反作用于这些设备，从而改变设备的振动和变形状态。液体晃荡的作用不仅在设备自身振动时产生，且振动消失后仍将持续一段时间。设备振动时导致的液体晃荡，不仅受容器本身的振动控制，还会受到惯性和重力的作用[2]；液体晃荡会对限制其运动的固体表面通过动压力的方式产生晃荡力，从而反过来影响容器的运动。显然，与不含液体介质的特种设备相比，此类设备在加速度载荷下的结构运动特性会发生一定的改变，与此同时受到结构振动和液体晃荡的影响，含液体结构在地震或其他载荷激励作用下的动力响应会非常复杂。

受加速度载荷激励引起的液体晃荡附加力对于设备安全有很大的影响。当地震发生时，如储罐、球罐、水塔等，地震作用导致的水体晃动而产生的晃荡力会是导致其结构破坏的主要原因之一[1]，其中破坏形式主要为支座破坏，从而导致设备倾倒；但同时液体晃荡可以吸收振动能量，达到减震的目的，如高层建筑的减震水箱[2]，为了满足使用功能要求，其容量高达上百吨，可能会增加结构的惯性力导致结构震害增加[3]，但若设计合理水箱可作为调频液体阻尼器（Tuned Liquid Damper, TLD）的一部分，利用液体晃荡过程中产生的动侧压力提供减振力[4]。因此，从保障特种设备安全及震后设备安全评价角度看，为了控制或利用液体晃荡，优化特种设备抗震结构，充分利用水体晃荡力来削减地震产生的惯性力，抑制地震对设备产生破坏，十分有必要对振动作用下液体的动力响应进行研究。

液体在振动作用下的动力响应由液体晃荡和液体晃荡力与支撑结构动力耦合作用两部分组成，而如何精确获取液体的晃荡及其产生的晃荡力是其中的重点和难点。液体晃荡的研究方法大致可以归为理论研究、实验研究和数值计算三大类[5-9]。如Chen等[5]通过直接数值求解欧拉方程的理论方法求解储罐大幅晃荡过程；祁江涛[6]和刘文夫等[7]通过VOF方法结合动网格理论对LNG船液舱晃荡进行了研究；大连理工大学岳前进课题组[8]通过搭建大型振动台对LNG储液舱的晃荡载荷分布规律进行了研究。目前理论研究对复杂构体或粘性较大流体的晃荡难以获得精确的解析解；通过振动台实验可以模拟地震载荷下储罐水体晃荡，但投资巨大，且存在尺度效应；随着计算机技术的高速发展，数值计算方法越来越受到重视，已经逐渐成为现代研究液体晃荡问题的重要手段[10-13]。

VOF（Volume of Fluid）方法[10，11，13]是一种成熟的自由液面追踪的欧拉数值方法，它以追踪流体体积所占网格单元体积分率的途径描述自由表面变化。本文基于VOF方法对中国特种设备检测研究院载荷响应实验室充液球罐晃荡实验进行数值仿真，仿真所得球罐晃荡附加力与实验吻合良好，验证了VOF方法在特种设备晃荡研究中的可行性，为后续研究打下基础。

1 VOF方法

本文利用ANSYS Fluent®17.0软件，采用VOF方法模拟水体在球罐中的晃荡过程，并结合分段线性化界面计算方法（Piecewise-liner Interface Calculation,PLIC）确定流体间的界面并捕捉气泡的存在。湍流模型使用了k-ε模型封闭湍流粘度项。在振动过程中水和空气的总体积基本保持不变，密度几乎不变，流体的压缩性对晃荡作用的影响可忽略，故采用不可压缩模型[2]；认为晃荡过程为等温过程不涉及能量输运方程的计算；本模拟的两相流动中的介质为空气和水，不存在反应，因此不存在质量源相和相间交换项。则根据以上模型简化，建立主要的控制方程如下：

体积分率输运方程：

体积分率之和为1：

式中：

f ——体积分率；

ρ ——密度；

v ——速度；

q ——第q相。

求解过程中仅求解n-1相的输运方程，而主相的体积分率由体积分率的限制条件即体积分率之和为1确定。

VOF方法中仅求解一个动量输运方程，认为所有相的速度都是相同的，这种假设的引入在相界面处速度差异较大的情况下会造成较大的误差。

动量输运方程：

式中ρ，μ为控制体内的流体密度和粘度，其表达式为：

式中：

μt——湍流粘度，本研究中选用k-ε模型对此项进行封闭；

g ——重力加速度；

F ——外部体积力。

2 球罐振动实验与模拟

杨智荣等[14]利用中国特种设备检测研究院载荷响应实验室的多通道动态载荷实验台对不同充装液位球罐在不同频率、振幅外载激励下的响应规律进行了研究。

实验系统如图1、图2所示，实验用球罐直径为1m，液体为水并添加染色剂便于数字摄影分析。本次实验的振动输出为正弦曲线输出，振幅15mm，频率1.5Hz；通过加速度传感器和三维力传感器测量振动过程中的加速度和球罐水平受力，传感器分布见图2。其中传感器1和3测量主振动方向加速度（X方向），传感器2测量垂直主振动方向加速度（Y方向）。传感器4、传感器5、传感器6、传感器7分别安装在四个支腿底部，可测量四个支腿在X和Y方向的受力，通过加和支腿上的受力再去除装置的惯性力可以得到晃荡附加力。

图1 球罐振动实验装置与球罐数值模拟网格划分

图2 球罐水体晃动实验测点分布

本文选取三个不同液位的实验过程进行了模拟，盛装量分别为30%、50%和70%，采用VOF方法耦合k-ε湍流模型，通过用户自定义函数（User De fine Function，UDF）添加动量源相的方式实现振动的模拟，输入的加速度为传感器1和2的实测加速度，其中测得的加速度（50%液位）如图3所示。CFD网格采用ANSYS ICEM软件划分，全部为六面体网格，总共约245000网格单元，经网格无关性测试表明此网格分辨率可以满足模拟需求。模拟采用非稳态模拟，时间步长为0.001s，计算时长12s，振动持续10.5s。罐体表面设置为无滑移壁面边界条件。

图3 装载量为50%液位时测得的X和Y方向的加速度

表1 球罐晃荡数值模拟的物性参数和模型设置总结

ANSYS Fluent中详细的物性和模型等的设置见表1。求解在中国特种设备检测研究院数值仿真实验室的高性能图形工作站上进行，采用20线程并行计算，单次模拟耗时6h左右。

3 模拟结果与讨论

图4展现了50%液位时球罐晃荡模拟不同时刻的液面演变，可以看到在晃荡初期液体晃荡速度较小，液面整体倾斜保持基本平直；随着振动时间增加，速度增大，液面会呈现出波浪状推进；且在波浪的推进过程中会出现液体与液面的分离拍打，形成气泡等现象。

图4 50%液位球罐晃荡模拟不同时刻的液面演变

不同液位沿振动方向（X方向）的液体晃荡所引起的附加力与实验对比，如图5所示，Force_X_exp和Force_X_cal分别代表实验和模拟得到的晃荡附加力。从图5中可以看到随着液位的升高，液体晃荡所引起的附加力逐渐增大，但并非随质量线性增加；当液面从50%涨到70%时相比液面从30%涨到50%，虽然球罐内增加的液体质量一致，但晃荡附加力的增加幅度明显加大；结果显示附加力大小和频率与实验值吻合较好，证明此模型可用于振动过程中的球罐结构分析和内部晃荡频率预测；在实验中发现当液位较低时，水体晃荡过程中湍动强度较大，产生了较多的气泡和水花[14]，但在模拟中由于湍流模型的缺陷和VOF方法均匀化的限制，并未很好地捕捉到这一现象，这可能是造成水位较低（30%液面）时实验与模拟结果虽然仍在可接受范围内(＜15%)，但明显大于水位较高（70%液面）时的误差。在未来的研究中将更多地去考察湍流带来的影响，以及如何在模拟中准确的进行描述。振动输入在10.5s处停止，在图中此时刻之后仅显示模拟结果，可以看到在X方向上晃荡附加力在振动源相消失之后，短时间内仍会保持相似的幅值与频率。

图5 不同液位沿振动方向（X方向）液体晃荡模拟得到的附加力与实验结果对比

图6 不同液位沿垂直振动方向（Y方向）液体晃荡模拟得到的附加力与实验结果对比

不同液位沿垂直振动方向（Y方向）的液体晃荡所引起的附加力与实验对比，如图6所示，Force_Y_exp和Force_Y_cal分别代表实验和模拟得到的晃荡附加力。与X方向受力相比，由于振动源的高频特性，其模拟结果也呈现高频状态，实验和模拟结果基本吻合，实验测得的液体晃荡力略大于模拟值；Y方向晃荡附加力的幅值与X方向相比显然是不容忽视的，由此可见在实际工程问题中不可忽略垂直主振动方向的液体晃荡附加力。与X方向类似，随着液位升高，液体晃荡所引起的附加力逐渐增大，但并非随质量线性增加，从50%到70%液面的晃荡附加力增幅明显更大。在振动源相消失后，Y方向的晃荡附加力表现出与X方向附加力不同的响应规律，其幅值会立即下降，这可能是由于Y方向的液体晃荡呈现高频的特点造成的。图5中也能体现这一规律，图3中X方向的振动输入加速度，可以认为是由两个波叠加而生成的，一个幅值大频率低的长波和一个类似于Y方向加速度的高频短波，图5中停止振动后液体晃荡力的曲线立即变得平滑，可以认为是高频短波的影响在振动消失之后很快衰减。

为考察液体晃荡附加力的频率是否与振动频率一致，对模拟得到的晃荡附加力进行了正弦函数拟合，例如50%液位的拟合结果见图7。振动源的频率为1.5Hz，而30%、50%和70%液面通过拟合分析得到的频率分别为1.497Hz、1.485Hz和1.503Hz，可见液体晃荡产生的附加力的频率与振动频率基本一致。

图7 50%液面时模拟得到的晃荡附加力正弦函数数值拟合结果

4 结论与展望

本文通过VOF方法耦合界面重构模型和k-ε湍流模型对不同液位下的球罐晃荡行了模拟，对球罐壁面压力进行积分得到了晃荡引起的附加力。模拟过程实现了实际测量加速度与计算流体软件的耦合，有利于本模型的推广。根据主振动方向和垂直其方向的晃荡附加力的实验和模拟结果，发现在考察球罐抗震安全性时垂直振动方向的高频振动的影响不可忽略；但液面较低时（30%液面）湍流强度较大，模拟与实验结果的偏差较大，湍流模型在球罐晃荡模拟中的作用需要进一步研究；从对比结果看出VOF方法可以较好预测振动状态下的液体晃荡附加力，证明该模型可用于大型球罐振动晃荡过程的研究。未来可通过模拟对大型储液设备的晃荡附加力的影响因素进行更加深入的研究，并引入流固耦合仿真对设备进行含液体晃荡力的结构分析，为储液设备达到本质安全发挥重要的作用。

参考文献

[1]Sezen H， Livaoglu R， Dogangun A．Dynamic Analysis and Seismic Performance Evaluation of Above-ground Liquid-containing Tanks [J]．Engineering Structure， 2008，30(03)： 794-803．

[2]尚春雨．带自由水面的结构地震反应分析方法及应用[D]．上海：上海交通大学，2008．

[3]刘云贺，张俊发，田洁，等．顶部设置水箱的建筑结构控震研究[J]．西安理工大学学报，1999，15(04)：71-75．

[4]荣峤．高层结构利用调谐液体阻尼器（TLD）控制水平地震反应的研究[D]．北京：北京交通大学，2009．

[5]Chen W， Haroun M A， Liu F．Large amplitude liquid sloshing in seismically excited tanks [J]．Earthquake Engineering and Structural Dynamics， 1996， 25： 653-669．

[6]祁江涛，顾民．LNG船液舱晃荡的数值模拟[J]．中国造船，2007（b11)： 541-549．

[7]刘文夫，薛鸿祥，唐文勇，等．独立B型LNG棱形液舱晃荡载荷数值分析[J]．船舶工程，2015，37(07)：22-25+72．

[8]卫志军，阮诗伦，岳前进，等．超大型储液舱内晃荡载荷实验研究[J]．海洋工程装备与技术，2014，1(01)：55-61．

[9]卫志军，岳前进，张文首，等．大尺度储舱液体晃荡砰击压力测量方法研究[J]．中国科学：物理学力学 天文学，2014，44（07）：746-758．

[10]ANSYS Inc．， ANSYS FLUENT 17．0 Documentation [Z]．2016．

[11]Youngs D L．Time-Dependent Multi-Material Flow with Large Fluid Distortion， Numerical Methods for Fluid Dynamics [C]．Morton K W and Baines M J， editors．1982， Academic Press．

[12]Hasheminejad M S， Mohammadi M M， Jarrahi M．Liquid Sloshing in Partly-filled Laterally-exited Circular Tanks Equipped with Baffles [J]．Journal of Fluids and Structures， 2014， 44： 97-114．

[13]龚国毅，赵成璧，唐友宏，等．基于VOF和浸入边界法的黏性二相流模型对LNG液舱晃荡的数值模拟[J]．科学技术与工程，2013，13(10)：2630-2634．

[14]杨智荣．大型球罐结构抗震安全保障中的关键技术研究[R]．北京：2016．