部分充液罐车内液体晃动的瞬态响应分析

王琼瑶， 蒋开洪， Subhash Rakheja, 上官文斌

(1.华南理工大学 机械与汽车工程学院，广州 510640；2.宁波拓普集团股份有限公司，宁波 315800)

随着我国物流业和交通运输业的迅猛发展，液罐车已成为石油、液化天然气、危险化学品运输的主要交通工具。由于液体货物密度的不同以及最大载荷限制等原因，液罐车通常处于部分充液状态。当罐车在制动或者转向时，液体货物未受限制的自由面将会发生显著的晃动现象。液体的大幅晃动将会改变车轴的载荷分布，产生附加的力和力矩，从而对罐车的行驶性能造成严重的不利影响[1]。事故统计表明，相比较于其他的载货车辆，液罐车更容易发生交通事故，且造成严重的人员伤亡和财产损失[2]。因此，研究液罐车内液体的晃动具有重要的现实意义。

对液体晃动现象的研究可追溯到1960年，然而研究主要集中在对航天器、船舶以及地面固定容器内液体晃动进行分析。对公路罐体内液体晃动的研究方法主要分为三类：准静态方法，等效力学模型法以及液体晃动动力学模型法。准静态方法主要用来预测液体晃动达到稳态时的动力学响应，且该方法只适用于不带防波板的罐体[3]。然而，研究表明，液体晃动瞬态时的晃动力明显大于稳态时的作用力[4]。等效力学模型可以用来分析液体的瞬态晃动，而且模型可以很方便地与车辆模型进行耦合[5]。包光伟[6]以油罐车为背景，采用Galerkin方法建立了贮箱平动激励下离散的液体受迫晃动方程,建立了三维受迫晃动响应的等效力学模型。结果表明,纵向晃动固有频率随贮箱长度而下降,晃动质量随贮箱长度而增加。等效力学模型的不足之处是模型参数需要从线性晃动理论或者试验获得，而且只限于对不带防波板罐体内的液体晃动进行分析。因此，尽管防波板是一种有效的抑制液体晃动的装置，但准静态方法和等效力学模型法均不能分析防波板在抑制液体晃动方面的影响。液体晃动动力学模型法又可以细分为线性动力学方法以及CFD方法。这两种方法均可对带/不带防波板的罐体内的液体晃动行为进行分析。线性动力学方法只限于分析液体的小幅度晃动，且需假设液体为不可压缩，无黏性，无旋转的流体。计算流体动力学(CFD)方法一般用于分析液体的大幅度晃动，由于该方法需要大量的计算量，关于其与车辆动力学模型耦合分析的研究较少[7-8]。

横向防波板一般用来限制液罐车内液体货物的晃动，尤其是在罐车制动的过程中。这类防波板通常中心开有一个圆孔，底部开有一个近似半圆形的通油孔。然而，用线性晃动理论和CFD方法对液体在带防波板的罐体内的晃动进行分析的研究较少。Kolaei等[9]运用边界元方法分析了水平圆柱形罐体内液体的晃动现象，研究并考虑了不同类型的纵向防波板。研究结果表明，固定于罐体顶部的部分防波板在常见的充液比的情况下能最有效地抑制液体的晃动，固定于罐体底部的防波板只在充液比较少的情况下才有效。刘小民等[10]以带有单个防波板的运动罐体模型为基础，对运动罐车在刹车过程中罐体内液体的晃动现象进行了数值分析，研究结果表明：罐体端面受力随着充液比的增加而增大，当充液比(液体体积与罐体总体积的比值)为0.8时，罐体的总体受力最大。然而，以上研究中的防波板板均为平板，未考虑防波板的曲率对液体晃动的影响，也未考虑防波板的其他几何参数的影响，比如防波板的倾斜角，开孔的形状、位置和大小等。

本文建立了液罐车罐体内液体晃动的三维CFD模型，并应用此模型研究了防波板的几何参数(曲率，开孔的大小、形状及位置，防波板的倾斜角等)对液体瞬态晃动时的载荷转移及晃动力的影响。模型的有效性通过与准静态模型以及试验结果的对比分析进行了验证。用横向和纵向载荷转移量，作用在罐体端面、防波板上的力的稳态值和峰值以及晃动产生的俯仰力矩和侧倾力矩评估了防波板的设计参数的影响。此外，对防波板与罐体端面以及相邻防波板间的空气压力形成机理进行了分析，并研究了空气压力对载荷转移以及晃动力产生的影响。

1 液罐车液体晃动的动力学模型

本文通过FLUENT求解器对液体晃动的动力学问题进行分析，如图1所示。数值计算采用两相流模型，即液相和气相。气相以空气为研究对象，由于液体晃动过程中会出现气体压缩的现象，空气采用可压缩的理想气体。气液两相密度差异较大，忽略了液体黏性力的影响。液体与罐体壁面、防波板间的相互作用可能会引起液体的大幅晃动甚至是液体飞溅现象，故流体域采用湍流计算模型。在气液两相流的数值计算中，两相之间存在随时间变化的运动界面，界面呈现多种多样的形状，本文采用的流体体积法(Volume of Fluid, VOF)能够很好地分析液体自由面的变形，有效追踪液体自由面的瞬态位置

▽·(uλ)=0

(1)

式中：u是流体速度矢量；▽是梯度算子；λ是流体体积分数。当λ等于0或者1时，表示网格单元分别为气相和液相占据，当λ介于0和1之间时，则表示气相和液相的交界面。

图1 带横向防波板的液体晃动模型Fig.1 The liquid slosh model in a tank with baffles

给定流场Z方向重力加速度以模拟重力场。给定流场X方向和Y方向一定的加速度以分别模拟油罐车在制动和转向时的工况，加速度随时间的变化关系如下

(2)

式中：g(t)表示加速度；δ表示上升时间；G表示加速度的稳态值；β是使加速度值从上升期到平稳期平滑过渡的系数，取0.2；Ta是过渡时间。加速度曲线如图2所示。

图2 圆滑过渡的斜坡阶跃加速度激励Fig.2 Rounded ramp-step acceleration excitation

液体晃动时会使得液体的重心实时发生变化，液体的重心坐标(Xcg,Ycg,Zcg) 相对于罐体几何中心的变化可由下式表示

(3)

记液体晃动时液体与罐体壁面接触的区域为A1，空气与罐体壁面接触的区域为A2，如图3所示。对于单元e，其对罐体壁面产生的作用力可由单元中心的压强Pe与单元的面积矢量Ae求出。因而，液体晃动产生的合力为

图3 液体的载荷转移量，力以及力矩的计算原理图Fig.3 Schematic diagram of computation of load shift, force and moment

(4)

式中：Fi为晃动合力沿i(i=x,y,z)轴的分量；Aei为单元e的面积矢量沿i轴的投影面积。俯仰力矩，侧偏力矩以及侧倾力矩的参考点均为罐体几何中心点在罐体底部的投影点O′,通过单个单元的力矩在整个湿周的积分求得

(5)

式中：Fe为单元e产生的晃动力矢量；le为单元相对于点O′的位置矢量；M为力矩矢量。

液体晃动会对罐体壁面、端面和防波板产生作用力。因此，晃动产生的合力可以看成作用在罐体壁面的力Ft和作用在防波板上的力Fb的和，即

F=Ft+Fb

(6)

式(4)～(6)计算得到的力和力矩只考虑了液相部分，忽略了气相(空气)对罐体产生的作用力。当充液比较大时，防波板的孔将完全浸没在液体中，如图1(b)所示，在加速度激励较大的情况下，防波板与罐体端面间或相邻两个防波板间空气压力将迅速增大。若考虑空气压力影响，液体晃动产生的总合力为

FT=F+Fa

(7)

式中：FT为总合力矢量；Fa为空气产生的作用力矢量，其求解方式与液体晃动产生的作用力相似，为气相单元e产生的作用力在区域A2的积分

(8)

在加速度激励作用下，晃动力的稳态值与加速度的大小成线性关系。然而，晃动力的瞬态值远大于其稳态值。所以在分析防波板在抑制液体晃动的作用时，有必要同时考虑力、力矩的瞬态值和稳态值。防波板的不同设计因素的作用可以通过瞬态放大因子来进行评估，其定义为作用在罐体和防波板上的力或力矩的最大值和相应的稳定值的比值，即

(9)

2 模型有效性检验

2.1 模型验证方法

提出的CFD模型能够有效计算液体在横向或纵向加速度激励作用下晃动产生的力和力矩。然而模型的有效性需要通过相关的理论或者试验来验证。常见的的验证方法主要有准静态模型法、频率验证法以及试验验证法。

2.1.1 准静态模型法

在准静态模型中，液体的自由面假设为平面，见图4。

图4 准静态模型Fig.4 Quasi-static model

在给定的加速度激励下，作用在罐体上的力和力矩为恒定值

(10)

2.1.2 频率验证法

液体自由晃动的频率由罐体的形状、大小以及充液比有关。Kobayashi等[11]1989年提出了圆柱形罐体在不同充液比时液体晃动频率f的计算表达式

(11)

式中：g为重力加速度；L为液体自由面长度；h为静态时液面的高度。在试验中，给罐体一个横向或者纵向简谐加速度激励，液罐内液体将会在加速度激励的作用下晃动。当加速度激励停止后，液体开始自由晃动；在CFD模型中，在流场的纵向和横向分别给定一个持续时间为0.06 s的三角波脉冲加速度激励，液体分别在俯仰平面和侧倾平面自由晃动。通过对比试验中和模型中液体自由晃动的频率，即可进一步验证模型的有效性。

2.1.3 试验验证法

为了进一步验证模型的有效性，用0.264∶1的试验模型(不带防波板)进行验证。试验用的罐体由有机玻璃材料制成，罐体壁厚为10 mm，目的是尽量减小液体晃动过程中罐体壁面的振动。试验是在振动试验台上进行。试验的结构见图5，罐体通过支撑架固定在上支撑板上，上支撑板通过四个传感器与下支撑板连接，下支撑板固定在振动台上。支撑板具体的结构尺寸见表1。

图5 试验模型结构图Fig.5 The set-up diagram of the experimental model表1 支撑板的尺寸Tab.1 The dimension of the support plate

参数d/mw/mh/m数值1.580.70.03

试验用的液体为水，充液比为50%，对罐体模型在纵向和横向分别施加ax=0.5sin(2π×0.37t)以及ay=1.0sin(2π×0.95t)的简谐激励，罐体内的水在简谐激励的作用下来回晃动，见图6。

(a) 瞬态:纵向晃动(b) 瞬态:横向晃动

图6 罐体模型试验图

Fig.6 The experiment diagram

将试验测得的晃动力与力矩分别与仿真得到的晃动力与力矩进行对比，验证模型有效性。

2.2 验证结果及讨论

图7(a)为不带防波板的罐体在gx=0.6g的斜坡阶跃加速度激励下以及充液比分别为50%、70%和90%时，由本文的CFD模型获得的纵向力稳态值与准静态模型获得的纵向力稳态值的对比。

(a) CFD模型与准静态模型获得的纵向晃动力稳态值得对比

(b) CFD模型与式(12)获得的自由晃动频率的对比图7 CFD模型验证Fig.7 The validation of CFD model

由图7(a)可知，两者差别非常小，误差主要由网格单元计算液体的总体积时出现的偏差造成的。图7(b)为本文的CFD模型获得的液体晃动自由频率与式(12)获得的结果的对比，由图可知，当充液比为30%～90%时，纵向晃动的频率范围为0.16～0.26 Hz，横向晃动的频率范围为0.56～0.81 Hz。纵向晃动频率和横向晃动频率均随着充液比的增大而增大。且CFD模型获得的结果与理论值十分接近。

图8(a)和(b)分别为纵向激励和横向激励下，试验结果与仿真结果的对比。由图可知，在纵向激励和横向激励的分别作用下，仿真值和试验值都基本吻合。

在纵向晃动中，由仿真获得的纵向力峰值比试验值大，主要原因是在试验中，液体自由面在晃动过程中出现分离(液体飞溅)现象，而在仿真中这一现象并未出现或液面分离程度较低。分离的液体下落到罐体内，见图6(a)，使得试验获得的垂向力峰值出现双峰现象。由于俯仰力矩不仅与纵向力和垂向力有关，还与液体的纵向载荷转移量有关，双峰现象在俯仰力矩中出现较少(仅在大约t=5 s时刻)；

在横向晃动中，由试验获得的横向力和垂向力峰值均比仿真获得的峰值大，主要原因是在试验中，罐体的壁面为有机玻璃面，在液体晃动过程中，罐体壁面也会发生相应的小幅振动，而在仿真中，罐体壁面假设为刚体，故不会出现罐体壁面振动的情况。由于侧倾力矩不仅与横向力和垂向力有关，还与液体的横向载荷转移量有关，试验结果与仿真结果差异较小。

试验验证结果说明模型是有效的，可以有效预测三维液罐内液体晃动产生的力和力矩。

(a) 纵向激励下的晃动响应对比

(b) 横向激励下的晃动响应对比图8 试验数据与仿真数据对比Fig.8 The comparison between the experimental data and simulation data

3 带防波板的液罐模型的仿真分析

以实际的液罐车的罐体尺寸为依据，建立了1∶1相应的简化几何模型。简化的几何模型尺寸如图9所示。罐体的横截面为圆形，半径为1 015 mm。罐体的圆柱部分沿X轴方向的长度为7 550 mm。罐体端面的曲率半径为2 030 mm。三个防波板沿水平方向固定于罐体的中间位置。防波板的曲率半径为2 030 mm, 孔的直径为D, 底端通油孔的半径r为144 mm，其面积占罐体横截面积的1%。坐标原点O取在罐体的几何中心。罐体内的液体为燃油，其密度为850 kg/m3。本文中除特殊说明外，防波板的中孔的形状默认为圆形，中孔位于防波板的几何中心，孔径为884mm。仿真的初始条件为：罐体的充液比为50%或70%；给定流体域横向(Y方向)或者纵向(X方向)斜坡-阶跃加速度激励(如图2所示)，以及垂向(Z方向)的重力加速度激励。边界条件为：罐体的壁面以及防波板均为刚性壁面，且壁面为不可滑移边界。

图9 带三防波板的液罐几何图Fig.9 Geometry of the tank with three baffles

3.1 网格和时间步长无关性检验

液体晃动的非线性相应受网格大小以及时间步大小的影响很大。选用三种不同大小的网格数来检验网格无关性，其对应的网格数分别为36 190,53 040以及207 148。选用了三种时间步长，分别为0.02 s，0.01 s和0.002 s。图10为在加速度激励为gx=0.3g，gy=0.25g，充液比为50%的情况下，纵向力和横向力随时间的变化曲线。由图可知，当网格数由36 190增大到53 040时，纵向力和横向力的响应曲线无论是峰值还是稳态值均会略微增大，而当网格数从53 040增大到207 148时，纵向力和横向力的响应曲线无论是峰值还是稳态值均变化很小，对应的曲线几乎重合，当网格数为53 040左右时，得到的结果可以接受。同理可得出时间步长为0.01 s较为合适，此处不再赘述。

3.2 罐体内空气压力的变化

在不带防波板的罐体内，液体晃动过程中，空气压力产生的作用力通常非常小，可以忽略[12]。然而在带防波板的罐体内，当充液比较大防波板的孔完全浸没在液体中时，防波板与罐体端面间或相邻两个防波板之间的空气压力将会非常大，如图1(b)所示。例如，在防波板1的左右两侧的临近区域分别取一个空气单元，如图11(a)所示，当加速度激励为gx=0.6g，gy=0.25g，充液比为70%时，监测两个空气单元在液体晃动过程中的气压随时间的变化，见图11(b)。由图11可知，晃动过程中，防波板1左右两侧的气压均发生了大幅震荡，且左侧气压值始终大于右侧气压值，表明左侧空间空气的压缩程度相对更高。左侧气压的峰值约为27 kPa，稳态值约为16 kPa，气压变化频率为1.85 Hz，远大于液体晃动频率0.58 Hz。

(a) 横向力

(b) 纵向力图10 网格数对瞬态晃动力的影响Fig.10 Effect of grid size on transient slosh force

(a) 防波板1左右两侧的空气单元

(b) 防波板1两侧空气单元气压变化的时间历程图11 防波板1两侧的空气单元及其气压变化的时间历程Fig.11 The air element on both sides of baffle 1 and the time history of air pressure of the two elements

罐体端面与防波板之间以及相邻防波板之间气压的大幅变化将会影响作用在单个防波板上的合力。图12为当加速度激励为gx=0.6g，gy=0.25g，充液比为70%时，作用在防波板1上的纵向力合力以及液体边界单元、气体边界单元分别作用在防波板1上的纵向力。由图12可知，由气压产生的作用在防波板1上的力甚至比液体晃动作用在防波板1上的力还要大，且两者方向相反，从而减小作用在防波板上的合力，表明气压的影响可以限制液体的晃动。

图12 分别考虑气相和液相时防波板1上承受的纵向力Fig.12 The longitudinal force exerted on baffle 1 considering air phase and liquid phase, respectively

3.3 防波板的设计参数及其对液体晃动的影响分析

3.3.1 防波板的设计参数

研究了防波板的不同设计参数对液体晃动的影响，这些设计参数主要有：防波板中孔的直径大小D(图13(a))、防波板中孔的位置高度H(图13(b))以及防波板中孔的形状S(图13(c))等，如表2所示。

表2 防波板的设计参数Tab.2 The design parameters of baffles

表2中的粗体部分表示防波板设计参数中的参考选项，为了研究这些设计参数对液体晃动的影响，固定其中两项，变化第三项。与表2中设计参数对应的防波板的几何形状见图13，其中ξ为防波板的曲率。

(a) 中孔的直径大小D(b) 孔的位置高度H

(c) 孔的形状

图13 防波板的设计参数

Fig.13 The design parameter of baffles

3.3.2 防波板的设计参数对液体晃动的影响

研究表明，液体在纵向加速度激励下，其载荷转移量以及晃动力与防波板的开孔面积关系密切，但这些研究涉及的防波板均是平板。本文研究了弯曲型防波板的开孔面积的大小对液体晃动响应的影响。由于开孔面积对横向载荷转移以及横向晃动力影响较小，讨论仅限于纵向载荷转移和纵向力。

图14(a)为在加速度激励为gx=0.6g，gy=0.25g，充液比为50%和70%时，纵向力放大因子随防波板圆孔直径变化的关系。防波板的圆孔直径分别为D=642 mm，780 mm以及884 mm，对应的防波板的开孔面积(包括通油孔面积)分别占罐体横截面积的10.5%，15.8%和20%。由图可知，随着防波板的开孔面积增大，纵向力放大因子逐渐增大，与充液比的大小无关。主要原因是：随着防波板的开孔面积的增大，纵向载荷转移量增大，如图14(b)所示，表明液体晃动的平均速度增大，使得晃动的剧烈程度增大，从而增大了纵向力的大小。

(a) 纵向力放大因子

(b) 纵向力载荷转移量图14 防波板的开孔面积对纵向晃动响应的影响Fig.14 The effect of opening area of baffles on liquid slosh response

除了防波板的开孔面积，防波板孔的形状也会对液体晃动产生影响。图15为在加速度激励为gx=0.6g，gy=0.25g，充液比为50%和70%时，防波板的开孔形状对纵向力放大因子以纵向力载荷转移量的影响。孔的三种形状分别为圆形，椭圆形以及十字架形，它们的开孔面积相等，与底部的通油孔一起均占罐体横截面积的20%。由图15(a)可知，与圆形孔和椭圆形孔的防波板相比，防波板孔为十字架形时，纵向力的峰值最小，与充液比无关。主要原因是：在十字架形孔周围会形成涡流，涡流的形成会减小液体晃动的动能[13]，从而减小晃动的剧烈程度。然而，当充液比为70%时，十字架形孔的防波板将使液体纵向载荷转移量更大，见图15(b)。主要原因是：十字架形孔的上沿位置相比于圆形孔和椭圆形孔的上沿位置更高，如图15(b)中的小图所示。在中等或高充液比的情况下，比较不容易形成封闭空间，空气压力的变化相对较小，其对液体晃动的抑制作用相对较小。椭圆形孔由于其孔的上沿位置最低，容易形成封闭空间，液体晃动过程中，空气压力的变化较大，其对晃动的抑制作用也因此变大。因此对椭圆形孔的防波板而言，其引起的纵向载荷转移量最小，如图15(b)所示。

(a) 纵向力放大因子

(b) 纵向力载荷转移量图15 防波板的开孔形状对纵向晃动响应的影响Fig.15 The effect of opening shape of baffles on liquid slosh response

受防波板孔的形状对液体晃动影响的启发，研究了防波板孔的位置高度对液体晃动的影响。图16为在加速度激励为gx=0.6g，gy=0.25g，充液比为50%和70%时，防波板孔位置高度对纵向力放大因子以及纵向载荷转移的影响。由图可知，当充液比为50%时，随着孔的位置高度的增大，纵向力峰值减小，而纵向载荷转移量却随之增大。当充液比为70%时，随着孔的位置高度的增大，纵向力峰值以及纵向载荷转移量均随之增大。主要是因为孔的位置高度将影响防波板的有效面积以及罐体内空气压力的形成。当充液比为50%时，孔的位置高度增大，防波板的有效面积增大，因此纵向力峰值减小，同时，孔的位置较高时，罐体内防波板与罐体端面或相邻防波板间的封闭空间不易形成，空气压力的变化很小，纵向载荷转移此时较大；当充液为70%时，若孔的位置高度为H=715 mm，则孔完全浸没于液体中，罐体端面与相邻的防波板之间的封闭空间在晃动的一开始便形成了，液体在加速度激励的作用下向罐体端面聚集并压缩封闭空间，封闭空间的气压会突然增大(如图11所示)。增大的气压会阻碍液体进一步向端面聚集，从而减小液体的纵向载荷转移量以及液体晃动的幅值。随着孔的位置高度逐渐增大，孔逐渐露出液面，封闭空间在晃动初始阶段并未形成，气压几乎不变(一个标准大气压)且其对晃动无影响，直至液体向罐体端面聚集过程中封闭空间形成。因此，当充液比为70%时，随着孔的高度的增大，纵向力峰值以及纵向载荷转移量均随之增大。

(a) 纵向力放大因子

(b) 纵向力载荷转移量图16 防波板孔的位置高度对纵向晃动响应的影响Fig.16 The effect of opening height of baffles on liquid slosh response

4 结 论

(1) 建立了部分充液罐车流体动力学模型，并通过准静态模型和试验方法等验证了模型的有效性。

(2) 研究了防波板的几何参数对液体晃动的影响。研究表明，防波板的孔的大小、形状以及位置高度对液体晃动的影响很大。通过减小防波板孔的面积以及降低孔的位置高度，均可提高防波板抑制液体晃动的能力。

(3) 提出了当防波板的孔的位置高度较低或充液比较大时，防波板与罐体端面或者相邻防波板间空气压力的形成机制，研究表明空气压力可以较小液体的载荷转移从而降低液体晃动的剧烈程度。因此空气压力的形成能有效提高防波板抑制液体晃动的能力。