储液容器内液体晃荡的非线性动力学分析＊

李文盛，赵友清，贾善坡，2，3，王 凯，谭继可

（1.长江大学城市建设学院，湖北 荆州434023；2.山东大学岩土与结构工程研究中心，山东 济南250061；3.中国科学院武汉岩土力学研究所岩土力学与工程国家重点实验室，湖北 武汉430071）

液体晃荡问题广泛存在石油、核工程、船舶、化工等领域，具有非常复杂的流体运动现象，呈现很强的非线性和随机性。当外界激励频率接近储液系统的固有频率时，内部储存液体晃动将变得更剧烈并导致容器破坏，不仅造成重大经济损失及人员伤亡，还造成更严重的次生灾害。因此，对液体晃荡和控制晃荡的研究引起了广泛关注［1－2］。

目前，对液体晃荡问题的数值研究方法主要有MAC法、VOF法、有限单元法（混合插值FEM、分布FEM、ALE有限元方法）及边界元方法等，并取得了一些重要成果。A.A.Amsden等［3］采用MAC方法研究了带有自由液面的大幅晃动问题，但这种方法计算存储量大、稳定性差，V.Armenio［4］提出了一种改进的 MAC法；包光伟等［5］、J.H.Jung等［6］采用VOF方法对液体晃动问题进行了数值计算；刘永涛等［7］基于VOF法得到液体晃荡运动的数值计算模型；周宏等［8］利用任意拉格朗日欧拉法（ALE）实现了自由液面的跟踪，并对晃动问题进行了分析讨论；R.Sygulski［9］假设流体为无黏、不可压缩、微幅晃动的理想流体，采用边界元方法研究了三维液体的晃动问题；C.Z.Wang等［10］基于流体速度势理论对流体非线性晃荡问题进行了有限元分析。

1 非线性波动理论模型

使用流体动力学非线性波动理论求解储液系统流体自由液面大幅波动问题，将三维模型简化为二维几何非线性平面应变问题，本文中运用有限元软件ABAQUS［11］，采用状态方程进行求解。在非线性波动理论模型有限元网格中，流体部分采用4节点平面应变单元CPE4R，每个节点具有2个位移自由度；刚体部分采用R2D2单元；刚体与流体相互作用使用法向接触，其中刚体表面为主面，流体表面为从面；储液容器壁高足够大以保证晃动时液体没有溢出。

为了得到稳定的数值解，假设流体材料为几乎不可压缩的，且黏性很小。流体体积响应通过us－up线性状态方程和牛顿黏性剪切模型建模。在储液系统模型中，为确定材料体积强度和压力提供了一个水动力模型作为液体密度和能量（质量内能）的函数。Mie－Grüneisen状态方程中，利用激波速度us和粒子速度up的线性关系，有：

式中：ζ＝1－ρ0／ρ，ρ0是扰动前流体密度，ρ是扰动引起变化的流体密度。

模型中，非线性通过体积应变计算，由于水的体积弹性模量很大（2.07GPa），几乎不可压缩，根据文献［11］，取体积弹性模量比真实值小2～3个数量级，即2.07MPa，也满足了介质的不可压缩性。水的剪切黏度ν应足够小，因为水无黏性，一个大的黏度将导致过硬反应。由牛顿流体Navier－Poisson定律，有：

式中：s为偏应力，ep为偏应变率。

在液体晃动问题中，重力作为恢复力是必不可少的，因此在初始状态将重力赋于所有流体单元。晃动波高满足：

液体晃动动态平衡方程：

罗四强说着拉起阿里的手，带着他朝花坛走去。刚走几步，另一侧的悼念厅又有哀乐响起。阿里又是一个寒噤，大声道：“姆妈醒了!”说话间，甩开罗四强的手，又朝着侧厅奔跑。罗四强急追好几步，才抓到他。

式中：üM为节点相对加速度矢量；MNM为对角化的一致质量矩阵，MNM＝∫VρNNNMdV，NN为单元形函数矢量；IN＝∫V（βN：σ）dV，σ为单元应力；PN＝∫SNMfds＋∫VNNFdV，f为面力，F 为体力。

当储液系统受到外加激励¨Xg时，上式可写为：

利用体积黏性压力抑制最高频率振荡单元，储液系统体积黏性线性表达式为［12］：

式中：b1为阻尼系数；cd为疏密波的速度；Le为单元特征长度；在方程（5）时，体积黏力也应包含在内。

2 算 例

模型为刚性矩形容器，底部固定于刚性基础。模型相关参数如下：底部截面宽度2a＝3.0m，刚壁高度H＝4.5m，容器内液体深h＝4.0m，密度ρ＝1 000kg／m3，声速c0＝1 438.75m／s（体积模量2.07GPa），黏滞系数0.001Pa·s，重力加速度g＝9.81m／s2，如图1所示。

图1 非线性波动理论模型Fig.1 The 2Dmodel considering nonlinear wave theory

二维储液容器自由晃动问题中，湿面在x、y方向是相互独立的，对Navier－Strokes方程分离变量即可获得流体自由晃动的固有频率。由于容器内流体材料为水，黏性足够小，可视为不可压缩、无黏、无旋的理想流体，则流体晃动的固有频率［13］：

自由液面波动形状：

式（8）～（9）表明，液体自由液面波动形状为正弦或余弦曲线。

2.1 数值验证

为了证明本文中所述有限单元法的准确性，将水平正弦波激励作用下二维矩形容器内液体晃动问题作为验证实例。在底部平行x方向施加水平正弦波激励¨ug＝Agsinωit，其中A为激励幅值（实例中取为0.01）。对每次施加的加速度频率ωi，获得其前4个周期内对应的液体自由表面晃动波高的最大值ηmax，如图2所示。由于所施加的外激励频率ωi与系统固有频率相等时系统即发生共振现象，而流体的共振主要表现为自由液面晃动时产生最大波高，通过波高峰值对应的频率可确定系统的固有频率。所以分别找出前三阶固有频率：第1阶，固有频率ω1＝3.20rad／s，η1，max＝0.201m；第2阶，固有频率ω2＝4.51rad／s，η2，max＝0.025m；第3阶，固有频率ω3＝5.60rad／s，η3，max＝0.029m。

为验证上述所得结果，通过式（7）计算前3阶液体自由晃动时固有频率解析解（见表1）。由表1可知，本文中采用的有限单元法具有有效性和可靠性。

图2 储液系统晃动波高Fig.2 Surface wave displacement for the tank－liquid system

表1 液体自由晃动频率Table 1 Frequencies of liquid sloshing

2.2 水平正弦激励下的液体晃荡

由图2可知，在不同频率的正弦激励作用下引起的液体晃动最大波高，其中第1阶固有频率引起的液体晃动最大波高是第2、3阶固有频率的6～8倍。可见，储液容器液体晃动响应主要由第1阶固有频率决定。图3为第1阶固有频率和第3阶固有频率激励作用时自由液面点B的波高曲线，可见，激励共振频率决定了液体晃动波形形状。

图3 不同谐频下自由液面点B的波面响应Fig.3 Free surface elevation of poit B for harmonic frequencies

2.2.1 第1阶激振频率时动力响应分析

算例中，A＝0.01，第1阶的固有频率ω1＝3.20rad／s，总时间为10s，分析步长Δt＝0.01s，水的黏滞系数足够小，取ν＝0.001。图4（a）为点A、O和B在第1阶激振频率液体晃荡时自由液面波高曲线，A、B处液面波动随着周期数增多而逐渐加强，意味着系统在该频率作用下发生了明显的低频共振现象，施加激励的频率即为系统固有频率值。液面中心点O处液体随时间变化几乎没有明显波动；左壁面点A液体晃动波高在t＝8.82s时出现最大值ηmax＝0.183m；右壁面点B液体晃动波高在t＝9.80s时出现最大值ηmax＝0.201m。观察点A、B处的波高曲线，两点的波动性态具有相反性，且液体晃动具有明显的非线性，主要表现在液体晃动的幅值向上的要大于向下的。

为了研究激励作用下不同液深处流体的晃动响应，在流体右界面（流体与刚壁交界面）由自由液面向下每隔0.5m依次取节点。图4（b）给出计算获得的点B、C、D、E、F、G液体晃动波高曲线，由图可以看出，波高幅值随液体深度的增加逐渐递减，这说明了水平正弦激励下底部固定的储液容器的振动是基本的梁型振动。

图4 第1阶频率作用下液体晃动波高曲线Fig.4 Variation in time of the surface wave in the first sloshing mode

图5给出了非线性分析时流体第1阶晃动模态图及对应液体晃动位移矢量图，由图可以看出，自适应网格技术使非线性流体动力响应分析保持很好的结构网格，整个计算过程中单元网格没有出现过大的畸形，保证了使用该方法计算结果的稳定与准确。

图5 第1阶液体晃动模态图和液体晃动位移矢量图Fig.5 Liquid shapes corresponding to the first sloshing mode and displacement vector plot

2.2.2 不同幅值激励时动力响应分析

取水平正弦激励¨ug＝0.01gsinω1t和¨ug＝0.005gsinω1t，图6分别给出了两种不同幅值激励作用下自由液面点A、B处液体晃动响应曲线。由图可知，液体晃动波高曲线与水平正弦激励频率基本相同，说明在不同幅值激励作用下液体晃动响应频率不随幅值变化；同时液体晃动上下波动的幅值出现明显不对称性，向上的波动高于向下的波动，并且相差程度随着幅值的增加不断增大。这些现象恰恰证明了液体晃动具有强的非线性。

图6 不同幅值激励下液体晃动的波高曲线Fig.6 Free surface elevation of liquid for different amplitudes

2.3 EI地震波作用下的液体晃荡

采用地震反应分析最常用的EI Centro地震波，最大峰值加速度为35cm／s2，加速度时间间隔为0.01s，总时长20s，图7（a）为该地震波加速度曲线，图7（b）为EI Centro地震波作用下液体自由液面相对两点A、B的晃动波高曲线。从图可以看出，晃动时点A处的波高在t＝17.80s时出现最大值ηmax＝9.3cm；晃动时点B处的波高在t＝15.10s时出现最大值ηmax＝6.9cm；在开始晃动的前12s，液面晃动波高幅值较小，随后液面晃动波高幅值逐渐变大并在17.80s后减小，可见，地震波对储液容器液体晃动反应需要一个激发过程，液体晃动波高幅值并不同随地震加速度幅值的变大而变大、减小而减小，最大液面波高幅值往往发生在激励加速度最大幅值以后，这是因为输入的地震波激励加速度虽减小了，但储液系统中输入的地震波能量还在增加。同时还可以看出，点A、B的波动性态具有相反性，且液体晃动具有明显的非线性，主要表现在液体晃动的幅值向上的大于向下的。

图7 EI Centro地震波和液面点A、B液体晃动的波高曲线Fig.7 EI Centro seismic wave and free surface displacement curve at points Aand B

3 结 论

使用有限元法，基于非线性波动理论建立储液容器液体晃动数学模型，通过对其施加水平简谐激励得到液体晃动的固有频率和模态，研究了二维刚性储液系统非线性晃动响应特性。该方法亦可适用于任意复杂形状三维弹性容器内液体的晃荡问题。由数值实例获得的结果表明，该方法求得容器内液体晃动的固有频率具有较高的精确度，外加激励频率及幅值对液体晃荡影响较大，表现了明显的非线性。

［1］Hasheminejad S M，Aghabeigi M.Sloshing characteristics in half－full horizontal elliptical tanks with vertical baffles［J］.Applied Mathematical Modelling，2012，36（1）：57－71.

［2］Xue Mi－an，Lin Peng－zhi.Numerical study of ring baffle effects on reducing violent liquid sloshing［J］.Computers＆ Fluids，2011，52：116－129.

［3］Amsden A A，Harlow F H.A simplified MAC technique for incompressible fluid flow calculations［J］.Journal of Computational Physics，1970，6（2）：322－325.

［4］Armenio V.An improved MAC method（SIMAC）for unsteady high－reynolds free surface flows［J］.International Journal for Numerical Methods in Fluids，1997，24（2）：185－214.

［5］包光伟，王振强，张天翔.火箭推进剂液体晃动关机响应的数值仿真［J］.宇航学报，2002，23（2）：84－88.Bao Guang－wei，Wang Zhen－qiang，Zhang Tian－xiang.Numerical simulation of slosh of the propellant fuel during the period the shut－down of the rocket［J］.Jourual of Astronautics，2002，23（2）：84－88.

［6］Jung J H，Yoon H S，Lee C Y，et al.Effect of the vertical baffle height on the liquid sloshing in a three－dimensional rectangular tank［J］.Ocean Engineering，2012，44：79－89.

［7］刘永涛，马宁，顾解忡.各种激励作用下液舱内液体晃荡的计算与分析［J］.船海工程，2009，38（5）：7－12.Liu Yong－tao，Ma Ning，Gu Xie－chong.Calculation and analysis of liquid sloshing loads in tanks under different kinds of stimulations［J］.Ship and Ocean Engineering，2009，38（5）：7－12.

［8］周宏，李俊峰，王天舒.用 ALE有限元模拟的网格更新方法［J］.力学学报，2008，40（2）：267－272.Zhou Hong，Li Jun－feng，Wang Tian－shu.Mesh update algorithm in ale finite method within free surface flow［J］.Chinese Journal of Theoretical and Applied Mechanics，2008，40（2）：267－272.

［9］Sygulski R.Boundary element analysis of liquid sloshing in baffled tanks［J］.Engineering Analysis with Boundary Elements，2011，35（8）：978－983.

［10］Wang C Z，Khoo B C.Finite element analysis of two－dimensional nonlinear sloshing problems in random excitations［J］.Ocean Engineering，2005，32（2）：107－133.

［11］ABAQUS／explicit user’s manual：Version 6.4［M］.Rhode Island：Hibbit，Karlsson and Sorensen Inc，2002.

［12］ABAQUS／theory user’s manual：Version 6.4［M］.Rhode Island：Hibbit，Karlsson and Sorensen Inc，2002.

［13］梅强中.水波动力学［M］.北京：科学出版社，1984.