液罐车内弹性挡板对液体晃动的响应分析*

2023-12-13王琼瑶林国勉

机电工程技术 2023年11期

徐凯，王琼瑶，林国勉

（五邑大学智能制造学部，广东江门 529000）

0 引言

液体晃动现象常出现在部分充液容器中，并且涉及到工程领域的各个方面，包括航空航天领域、陆运领域、海洋工程领域以及储罐等。液体的剧烈晃动可能导致充液容器的结构破坏、液罐车侧翻事故的频发以及航天器发射的失败等。因此，研究如何抑制充液容器内液体的晃动对工程应用具有十分重要的意义。

自20 世纪60年代开始，Abramson 等[1]比较全面地开展了常重力环境下各种形状的挡板的液体晃动实验研究。Hasheminejad 和Mohammadi[2]使用线性势理论的二维动力学方法研究了液体在无黏性和不可压缩的条件下，水平圆柱挡板容器的3 种挡板安置方式对液体晃动的影响。结果表明：容器在高充液比情况时，底部安装垂直挡板不是抑制部分填充容器内液体晃动的最有效的手段。Yu和Chu[3]研究了部分充满液罐车的液体晃动的非线性运动，研究结果表明：在充液比在90%时，瞬时侧向力最小。Yan 和Rakheja[4]研究了部分充液的液罐车在有、无挡板的情况下其直线制动性能特征。研究结果显示当罐体内有挡板时，液罐车的制动性能明显提高。Wang 等[5]基于VOF （Volume of Fluid）模型研究了无挡板，传统刚性挡板和组合挡板在横、纵方向上对液罐车运行稳定性的影响。研究结果表明：组合挡板能很好地降低液体对罐体结构的冲击强度。虽然刚性挡板能有效地抑制容器内液体的晃动，但是刚性挡板会增大罐车结构的质量，还会使罐内的清理工作变得麻烦。随着复合材料的不断发展，人们发现弹性材料也可以制作成挡板。不仅可以起到抑制罐内液体晃动的作用，而且还可以减轻罐车结构质量。Duan 等[6]采用流固耦合方法分析了在加速过程中的液体晃动现象、挡板的变形及受力情况，研究了在充液比为0.5 和0.8 时水平圆柱三维模型的液相分布。研究表明：充液比为0.8 的储罐液体分布更稳定。Zhang等[7]运用SPH （Smoothed Particle Hydrodynamics）法研究了二维矩形罐多弹性挡板及其组合挡板对液体晃动的影响。研究表明：垂直或T 形挡板与水平挡板的组合在抑制晃动方面要优于单个挡板。Saghi 等[8]提出了一种水弹性模型，研究由具柔性壁面的梯形和矩形储罐上的摇摆引起的晃动载荷。结果表明：梯形储罐比矩形储罐有更好的抑制液体晃动的性能。Yu 等[9]采用有限元分析和任意拉格朗日-欧拉（ALE）方法研究了挡板渗透率对液体晃动的影响。研究发现一个最有效的抑制晃动的临界渗透系数，并且带有孔洞挡板的储罐内会出现旋涡。Cho 等[10]基于速度势的非线性有限元法研究在水平激励下二维罐内液体的大幅晃动，研究表明：在罐体内液体的晃动量与防波板的设计参数有很强的相关性。Meng等[11]提出了一种半解析数学模型来研究矩形刚性容器与弹性挡板相互作用下液体的晃动特性。研究表明薄挡板能提高液体晃动的固有频率。包文红等[12]采用Fluent 软件研究液罐车在不同加速度激励下，液体晃动对罐壁的冲击力随时间变化的规律。研究结果表明液罐车在减速阶段挡板受到的冲击力最大。徐晓美等[13]对液罐车在不同充液比和不同激励时的稳定性问题进行了研究。研究结果显示液罐车充液比为50%时，液体晃动对车辆侧倾稳定性影响最大；而且随着充液比的增大，罐体所受的力和力矩也随之增大。Pozzetti 等[14]在VOF 算法的基础上，提出了一种全新的三相流多尺度方法。该方法在考虑不同长度尺度的影响的同时还能大幅度减少计算负担，研究结果显示：三相流多尺度方法比流体体积法提供了更高的精度。

基于上述研究，本文在传统刚性挡板的基础上探索研究弹性挡板的阻尼特性，研究液罐车以横向斜坡阶跃制动时罐内液体晃动对挡板产生的冲击，探讨液罐车在这种冲击下弹性挡板的相对长度（H/R）和厚度对液体晃动的影响。此项研究将对罐车内液体晃动的抑制研究提供一种新的思路。

1 计算模型

1.1 流体域的控制方程

本文使用开源软件OpenFOAM（Open Field Operation And Manipulation）对罐体内的液体晃动现象进行模拟。采用两相流模型对液体晃动过程中出现的气相（气体）和液相（液体）进行模拟。由于涉及液体晃动的液体流速很小，相应的雷诺数很小，晃动模型采用层流模型即可。涉及液体晃动的控制方程为不可压缩流体的连续性方程和Navier-Stokes动量方程为：

式中：ρ为罐体内液体密度；U为速度矢量；t为时间；p为压力；τ为剪切应力。

在两相流模型中，涉及到捕捉气液交界面的问题，本研究采用VOF （Volume of Fluid）方法来捕捉自由液面。VOF 法引入F（x，y，t）函数，主要是跟踪自由面的流体体积分数，表示流体占所在网格体积的百分数。当F=0 时，网格内全为气体域；当F=1 时，网格全处于液体域；当0＜F＜1 时，表示网格一定含有自由液面。因此，函数F的求解，就可以跟踪自由液面的位置。函数F随时间变化的控制方程为：

式中：θ为部分单元体参数，其值在0 到1 之间；u、v分别为x、y方向上的流体速度分量。

由于圆柱壁面边界处流体和固体没有相对运动，所以采用无滑移边界条件；在运动中通量修正保证弹性挡板面通量为0，所以采用运动的壁面边界条件。

1.2 弹性挡板的控制方程

由结构力学，弹性挡板的控制方程可以表达为：

式中：ρ为弹性挡板密度；d为位移；f为体积力矢量；σ为柯西应力张量。

使用定常剪切模量参数输入超弹模型Neo-Hookean，其弹性应变能函数表达式为：

式中：μ为剪切模量；I1为第一偏应变常量；J为弹性变形梯度的行列式；D1为材料的不可压缩参数，详细定义为：

式中：γ为泊松比。

1.3 数值模型建立

选取液罐横向截面，建立一个半径R=1.015 m、厚度为0.05 m的二维圆形刚性罐体（罐体的厚度为一个网格尺寸大小）。如图1 所示，在罐体中心处放置高为H、宽为L的弹性挡板，采用了无量纲的方法设置了8 种弹性挡板的相对长度，其具体参数设置如表1 所示，弹性挡板的具体材料参数如表2 所示，力矩的参考点在罐体的底部点O′，罐体充液比为50%。二维圆形模型的构建、网格的划分和流体域数值模拟计算均在OpenFOAM 软件中进行，其中网格是在blockMeshDict 里构建的，固体域数值模拟计算在Deal.II软件中进行的。

表1 计算参数设置

表2 弹性挡板的材料设置

图1 罐体和挡板结构图

1.4 晃动响应分析

在充液比为50%时，如果液罐车发生变道制动、转向以及换道运动，其自由液面晃动最为剧烈，更易发生侧翻现象[15]。故本文选取罐体的充液比为50%。液罐车在转向时的离心加速度用斜坡阶跃的加速度进行模拟，加速度随时间变化的关系如式（7）所示。

式中：a(t)为加速度激励；k为加速度斜坡的上升率；a是稳态加速度的幅值，由于模拟的是液罐车的转向运动，幅值取0.25g；β为圆弧过度参数，一般取0.2。加速度曲线如图2所示。

图2 圆滑斜坡阶跃加速度激励

液体载荷转移量是通过动态液体质心（cg）与静止状态下液体质心之间的变化来评估的，其相应的计算公式为：

式中：xcg和ycg分别为x轴和y轴的液体载荷转移分量；x0和y0为静止状态下液体质心坐标；V为液体总体积；Va是液体单元a的体积，xa和ya是单元a的质心坐标；ψ为液体域，如图1所示。

对于任意单元b，用单元中心处压力与单元边界面面积的乘积来计算任意单元b对壁面的作用力，将该作用力在整个湿周（液体与壁面接触的部分）进行积分来计算作用在壁面的总晃动力，其计算公式为：

式中：F为液体晃动对边界面的总作用力；p为任意单元b的压力；s为单元与边界接触的面积。

根据力矩的定义，任意单元b相对于参考点O′的矢径与单元b产生的作用力的叉乘即为该单元对参考点O′的力矩。对该力矩在整个湿周进行积分即为晃动液体作用在壁面的总晃动力矩。相应公式为：

式中：M为作用在罐体壁面的总力矩；l为单元到参考点O′的位置矢量。

2 网格无关性验证

一般来说网格划分越细密，其仿真的结果越精确。但是实际上当网格加密到一定的程度时，仿真的结果变化就很小，而且会导致仿真时间成本大幅增加。故为了综合效率和精度一般认为，在满足库朗数限制条件下，当仿真数据随网格尺寸的变小而只产生微小变化时，就达到网格无关性了。网格尺寸对液体晃动仿真结果影响较大，故在仿真中选取了3种不同数量的网格。

针对充液比为50%的流体域，在OpenFOAM 中使用blockMeshDict 对模型进行分块划分及O 型块处理，每种网格的模拟是在ax= 0.25g的圆滑过渡的斜坡阶跃加速度下进行的，以模拟罐车转向时的转向加速度。设置了3 种不同数量的网格（31 528、44 038 和61 946），分别代表不同精细程度的网格。如图3 中3 种网格之间的计算结果相差微小，为了较小的仿真时间选取网格数为31 528为计算网格。

图3 网格无关性验证

3 二维横向模型晃动分析的结果和讨论

3.1 弹性挡板形变对液体晃动的影响

图4 为弹性挡板中心处水平位移的最大形变量随H/R的值和挡板厚度变化的趋势图。由图可知，弹性挡板中心处的水平位移峰值随着H/R的比值的增大而增大。从图4 中还可以看出，随着弹性挡板厚度的增加，其中心处对应的水平位移峰值减小；并且在同一H/R值时，弹性挡板的厚度增加5 mm 时，其水平位移峰值下降约10 mm。原因是当杨氏模量相同的时候，刚度随弹性板厚度的增大而逐渐增大，弹性挡板抵抗变形的能力相应增强。

图4 弹性挡板水平位移峰值

为研究弹性挡板中心处的水平位移量和压力变化与液体晃动的关系。图5 为10 mm 厚度的弹性挡板在板的相对长度（H/R）值分别为1.1、1.3、1.5、1.7 时，挡板中心处的水平位移量和压力的时间历程曲线。由图可知，弹性挡板中心处位移随时间历程曲线与正弦曲线相似，这表明在液体冲击作用下，弹性挡板在其平衡位置来回振动，并且其振动的幅值随着H/R值的增大而增大。

从图5 中还可以看出，弹性挡板中心处的压力峰值主要集中在0.7 s 附近，即液体晃动的第一个周期内，这表明在斜坡阶跃加速度作用下，液体晃动产生的冲击载荷在第一个周期内最大，此后冲击载荷逐渐减小。对比不同长度的弹性挡板下液体冲击载荷的峰值可以发现弹性挡板的长度对其压力峰值的影响较小。另外，观察压力的时间历程曲线可以看出，“双峰”现象较为明显，表明液体晃动运动出现了非线性现象。由于液体与弹性挡板的相互作用，晃动液体在较短的时间内多次与弹性挡板发生碰撞，弹性板的形变吸收了一部分液体晃动动能，使得压强峰值分化，避免出现过大的液体冲击力。

3.2 弹性挡板对液体载荷转移量的影响

液体载荷转移量直接影响作用在罐车上整体的力和力矩大小，是评价液罐车运行稳定性的重要指标。仿真中罐体的激励为ax= 0.25g，充液比为50%，设置了3 种厚度弹性挡板，即10 mm，15 mm，20 mm。弹性挡板的相对长度（H/R）分别设置为1.0、1.1、1.2、1.3、1.4、1.5、1.6、1.7 共8 个值。为研究弹性挡板对液体载荷转移量的影响，表3 列出了罐体内无挡板和设置不同弹性挡板的情况下，其液体载荷转移量峰值及其下降率。

表3 不同弹性挡板液体载荷转移量的峰值及与无挡板横纵载荷转移量峰值相对应的下降率

无挡板的横纵载荷转移量峰值分别为x0=0.193 m、y0=0.046 m。根据表3显示，相比较于无挡板的罐体，带弹性挡板的罐体的横向载荷转移量峰值降低19.9%～25.9%，并且其相应的纵向载荷转移量峰值降低32.4%～45.2%。从表3还可以看出：H/R值从1.0增大到1.6时，随着弹性挡板厚度的增大，其对应纵向和横向载荷转移量峰值的下降率呈减小趋势；当H/R值增大到1.7时，情况正好相反，但载荷转移量的变化幅度不大，在0.6%左右。该结果表明：在选择弹性挡板作为防晃装置时，减小挡板的厚度可以提高其对应载荷转移量峰值的下降率，即提高挡板抑制液体晃动的效果。

3.3 力和力矩

图6 为横向加速度激励下，作用在罐体上力的峰值随H/R值的变化趋势图。由图可知，和无挡板模型相比，弹性挡板可以很好地减少圆形罐体横向晃动力的峰值。原因是当水波撞击弹性挡板时，弹性挡板可以很好地把液体的动能转化为弹性势能，当弹性势能积蓄到一定的时候，会与后续的水波发生对冲消耗。并且在图5 的“双峰”现象也显示了弹性挡板在降低液体对罐体冲击力上有很好的效果。

图6 横向晃动力峰值

图7 为罐体内无挡板和放置弹性挡板时，液体晃动产生的侧倾力矩峰值随板的相对长度的变化趋势图。对于放置弹性挡板的罐体，罐体内液体晃动产生的侧倾力矩的峰值明显的小于罐体内无挡板时对应的侧倾力矩峰值，带弹性挡板的罐车具有更好的运行稳定性。

图7 侧倾力矩峰值

从图6～7 还可以看出，弹性挡板的厚度变化对横向晃动力和侧倾力矩的峰值影响不大。故在后续使用弹性挡板作为阻尼装置时，建议选取较小的厚度，这样既可以在性能上满足需求还可以实现挡板的轻量化。

3.4 罐体内的速度矢量和压强分析

图8 是H/R值为1.0、挡板厚度为10 mm 时，不同时刻流体域的速度矢量分布和压强分布图。从图中可以看出，在自由液面附近出现了破碎波的现象，由于液体与气体之间存在粘性摩擦力，液体表面速度越快，与气体之间的耗散就越大，其动能的损失转换为内能，并且不增加压力。即液体与气体边界交换的耗散会影响流体的势能。

图8 液体表面破碎的速度矢量分布和压强分布（H/R=1，板厚10 mm）

图9 为带弹性挡板的圆形罐体内的流体域的速度矢量分布和压强分布图，其中时间为1.36 s，板厚10 mm，板的相对长度分别为1.0、1.4和1.7。在图5中1～2 s内压强的峰值随着弹性挡板的相对长度的增加，压强却渐渐减小并趋近于零。此数据可以由图9 中的现象解释，从图中可以发现：随着弹性挡板相对长度的增加，其两侧液面的高度差逐渐减小，此时挡板中心处的左右压强差也在减小，此外，在挡板右下角处出现了涡流，并且随着弹性挡板相对长度的增加，涡流明显增大，底部的压强也随之增加。如图8～9 中液体产生的涡流会使内摩擦力做功，造成进一步的耗散。