基于拉格朗日描述的罐车内液体晃动模拟*

石怀龙 王勇 邬平波

(西南交通大学 牵引动力国家重点实验室,成都 610031)

引言

液体晃动会影响货物运输安全,如公路、铁路和航空运输中的液罐车以及海底运输管道等,主要是由于液体自由表面发生非线性运动,对罐体施加动态载荷.液体晃动行为研究起源于航空航天领域,主要难点为容器和液体接触面之间的运动边界条件模拟、液体自由表面的非线性运动模拟等.

Abramson等人基于线性弹簧振子模型研究了航天器内储油箱的液体振动,分析具有不同罐体形状内的液体晃动频率和振型,并给出接触作用力解析表达式,指出液体第一阶晃动频率对系统影响最大[1].Ranganathan等人建立了单摆-固定质量块模型,用单摆模拟液体一阶晃动效应,用固定质量块模拟未参与晃动的部分液体,并应用于公路槽罐车制动过程中的液体晃动行为[2,3].文献[4,5]详细阐述了等效模型研究进展,分析了不同建模方法和充比液时重载槽罐车的倾覆性能.然而,这些典型的等效物理模型过于简化,不能描述液体惯量和自由表面的连续性变化.Celebi基于Navier-Stokes方程建立了连续非线性模型,研究带有隔板的矩形贮箱内的液体晃动[6].Ibrahim[7,8]总结了二十一世纪初之前的液体晃动研究进展,明确了液体晃动的被动控制方法以及其在结构振动抑制中的应用,讨论了液体的非线性属性和结构非线性振动等编程实现的困难.Aliabadi[9]对比了连续流体模型与弹簧振子和单摆系统构成的等效模型的结果差异,分析了液体自由液面的变化和动态惯量特性,这些都是等效模型所不能实现的.文献[10]总结了近期液固耦合的建模手段和数值方法,并与试验数据对比分析各自优缺点.Bogomaz[11]和Vera等人[12]采用单摆等效模型模拟液体,建立考虑液体晃动的三维车辆系统动力学模型,分析车辆在启动、制动和连挂以及过曲线时的车钩振动、车体加速度变化情况,指出液体晃动对车辆承载和振动响应产生显著影响. Wang等[13]和Wei等[14]较早地采用连续介质力学方法模拟矩形贮箱内液体的晃动效应,可模拟液体表面的连续性变化.

虚拟质量法将液体考虑为附加质量代入系统质量阵,借助有限元软件计算当前液体质量,然后通过插值得到任意时刻的液体质量,即做了线性叠加假设[15].刘富等[16]通过光滑粒子流体动力方法研究矩形容器内隔板设置问题,合理设置隔板可有效抑制液体晃动.等效力学模型需借助流体动力学和有限元仿真等手段推算或拟合模型参数,并基于线性假设组建车辆系统三维晃动模型,可提高计算效率,但仅适用于基本规律性研究,不能描述液体自由表面和惯量的非线性变化[17,18].

综上,有必要提出可适用于描述液体惯量连续变化的非线性模型,并易于与车辆多体系统动力学模型集成.Shabana提出绝对节点坐标方法(Absolute Nodal Coordinate Formulation,ANCF),在全局坐标系下用节点位失的梯度代替传统有限元中的转角坐标,形成的柔性体质量矩阵为常数阵,在大变形问题研究中具有高精度等特点[19,20].

本文基于采用绝对节点坐标方法的柔性多体系统动力学理论,提出一种液体晃动模拟新方法,并应用于液罐车内的液体晃动问题.首先,基于流体力学基础理论推导液体粘性方程和体积不可压条件方程,然后,采用绝对节点坐标方法描述的实体单元对液体进行网格划分,并通过罚函数描述液体与罐体的接触作用,最后组建罐体-液体耦合多体系统动力学模型,仿真研究液罐车内液体的横向和纵向晃动行为.

1 基于绝对节点坐标方法的柔性多体系统

1.1 运动学描述

1.2 惯性力表达

(1)

1.3 弹性力表达

(2)

其中,单元刚度阵K为节点坐标的强非线性阵；格林-拉格朗日应变张量ε定义为：

(3)

图1 单元各构型定义及相互关系Fig.1 ANCF deformable body configurations

1.4 广义力表达

假设单元任意一点处受外力矢量作用f,则外力对应的虚功为fTδr,其中r为单元上的力作用点处的位置矢量,δr为r的虚变量.为获得与系统广义坐标相关的广义外力,将广义力写成节点坐标虚变量形式为：

(4)

则Qf=STf对应外力f的广义力.

1.5 运动微分代数方程组

由于绝对节点坐标方法的位置矢量定义在全局坐标系下,因此其运动约束方程描述和雅克比矩阵形式较简单.因此,这里仅给出含非线性约束方程的增广式多体系统运动微分代数方程组：

(5)

式中,Qe为系统广义外力,Qd为与约束方程相关的广义约束力,λ为拉格朗日乘子,Cq为与约束方程相关的雅克比矩阵.

2 流体力学基础理论

2.1 Navier-Stokes方程

流体力学研究方法通常可分为两类,欧拉法关注空间网格点处的流体流动情况,而拉格朗日法则追踪流体物质点运动轨迹.当流体粘性属性只与温度场相关,且内部应力与应变线性相关,称该类流体满足牛顿型特征,适用牛顿第二定律描述其运动.假设流体为不可压的牛顿类,连续性方程写为：

(6)

式中,ρ为流体密度(kg/m3),t为时间(s),r和v分别为位置(m)和速度矢量(m/s),为散度算子.假设流体不可压缩,即流体密度保持不变,则连续性方程式(6)简化为·v=0,即不可压条件.流体运动偏微分方程为：

(7)

σ={-p+λtr(D)}I+2μD

(8)

将公式(8)代入(7)得流体单元状态平衡方程：

(9)

(10)

式中,J=|J|；Cr=JTJ为右柯西-格林应变张量.

2.2 流体不可压条件

质量守恒条件为dm=ρodV=ρdv,其中ρo和dV分别为流体在初始构型时密度和体积,而ρ和dv分别为当前构型即发生形变后的密度和体积.流体变形前后的体积关系为[22]：

dv=rxdx·(rydy×rzdz)

=(rx·(ry×rz))dxdydz=JdV

(11)

上式两边同时乘以ρoρ,再利用质量守恒关系可得出流体发生形变前后的密度关系为ρo=ρJ.若流体不可压条件为J=1,表明流体发生形变前后的体积和密度均保持不变.

3 基于绝对节点坐标方法的液体运动描述

3.1 六面体实体单元

(12)

式中,a、b和c分别为单元沿轴x、y和z方向的长度；ξ=x/a,η=y/b,ξ=z/c,其中,ξ、η和ζ∈[0,1],ξk、ηk和ζk为节点k的自然坐标.单元j上任意一点的位置矢量表示为：

(13)

式中,I为3×3单位阵,Sj和ej分别为单元的形函数矩阵和节点坐标矢量.六面体实体单元含8个节点,共计96个自由度,其位置矢量及梯度矢量可以表达液体形状及其惯量变化.

3.2 体积不可压条件

(14)

(15)

3.3 单元内部粘性力

其中,vj为单元当前构型的体积,Vj为单元初始构型的体积.结合公式(10),液体单元的内部粘性力虚功表达式可简化为：

(16)

(17)

3.4 液体运动方程

根据液体单元的质量矩阵、外力矢量(式(4))、保证体积不可压条件的约束力(式(14)和(15))和粘性力(式(17)),则其运动微分方程组为：

(18)

若液体被划分为多个单元,通过组装单元运动方程来组建系统运动方程,利用数值仿真求得系统速度和位置矢量.

4 基于罚函数描述的液体与罐体接触关系

(19)

式中,K和C为罚函数刚度和阻尼系数.

由于正压力fn应垂直于罐体上接触点处切平面,因此需要计算接触点处的法向量v.已知罐体结构两端为半椭球体,中间部分为圆柱体,如图2所示,可根据液体空间位置计算出罐体上对应的接触点位置.摩擦力fc=μfn沿接触点处切平面方向,μ为摩擦系数.正压力和摩擦力同时作用在液体和罐体上,幅值相等但方向相反.

图2 液体与罐体接触关系的坐标系定义Fig.2 Interaction definition between the liquid and tank

(20)

(21)

4.1 中部圆柱体区域接触关系

(22)

式中：r为圆柱体半径,假设罐体为刚体.接触点处穿透速度定义为：

(23)

作用在液体上的力为：

(24)

(25)

4.2 端部椭球体区域接触关系

(26)

5 液罐车内液体晃动行为仿真

以液罐车装载类似水的液体为例,密度ρ=1000kg/m3,重力系数g=-9.81m/s2,粘性系数取μ=9.3×10-4Pa[14].体积不可压罚函数系数kIC=1.0×106N/m和cTD=1.0×104N·s/m.罐体长度12m,直径2.8m.模拟液罐车半载情况,液面高度1.4m,则液体质量约37.5t.

5.1 横向晃动模拟

在罐体上施加横向简谐激励f=0.5Hz、A=0.1m,仿真时间4s,模拟液罐车内液体的横向往复晃动过程,如图3和图4所示.t=0.0s时,半圆柱形液体在重力作用下开始与罐体内壁发生接触,最终构型与罐体内壁贴合；t=1.5s、2.5s时,液体发生纵向、横向和垂向变形以贴合罐体内壁,但不同断面处的液体自由表面形状不完全相同,即液体自由表面的运动呈现非线性；t=3.0s、4.0s时可以观察到罐体内液体的显著晃动行为,液体上表面位置呈非线性变化.仿真计算时仅采用了1个单元,参考第1.3节和第4节定义单元初始外形,以尽可能保证液体初始形状与罐体内壁贴合,从而节省计算工时.通过增加单元数,可以获得更加连续的液体自由表面变化形状,但基本上与1个单元情况相似,因此不再列写多个单元时的仿真结果.

5.2 纵向晃动模拟

对罐体施加纵向简谐激励f=0.1Hz、A=1m,模拟液罐车内液体的纵向晃动过程,仿真时间20s,结果如图5所示.可见,液体的纵向晃动必然引起液体质心的纵向移动.仿真计算表明其质心最大位移可达0.5m以上,按此值估算液罐车前后转向架轴重转移量可达15%,从而影响液罐车的动力学性能和车轮磨耗.液罐车采用长大编组方式运行,制动过程中的液体纵向晃动必然引起车钩/缓冲器受力的突变或制动效能损失.

图3 液罐车内液体的横向晃动(主视图)Fig.3 Liquid sloshing along lateral direction (axial view)

6 结论

1)本文提出采用绝对节点坐标方法的液体晃动模拟新方法,用于模拟液罐车内液体自由表面的连续性变化.该方法描述液体自由表面形状的连续性变化,并适用于研究具有复杂结构或形状容器的内部液体晃动问题.

2)基于流体力学牛顿类体基础理论,推导粘性方程和满足体积不可压缩的条件方程；采用绝对节点坐标方法描述的实体单元进行液体网格划分.

3)采用罚函数方法描述液体与罐体的接触作用,给出液体与罐体中部圆柱体或端部椭球体的接触关系,组建罐体-液体耦合系统动力学模型.

图4 液罐车内液体的横向晃动(俯视图)Fig.4 Liquid sloshing along lateral direction (side view)

图5 液罐车内液体的纵向晃动Fig.5 Liquid sloshing along longitudinal direction in tank

4)仿真计算液罐车内液体的横向和纵向晃动过程,发现液体自由表面形状呈非线性变化,不同断面处的高度和形状不同；液体的横向、纵向晃动必然引起液罐车的质心发生改变,从而影响轮重减载和车钩力变化.

5)本文限于研究罐体和液体之间的耦合运动,暂未考虑整车模型,也未考虑罐体内含隔板情况.

1Abramson H N. The dynamic behavior of liquids in moving containers-with applications to space vehicle technology. Technical Report NASA SP-106, 1966. San Antonio,TX: Southwest Research Institute

2Ranganathan R,Ying Y, Miles J B. Analysis of fluid slosh in partially filled tanks and their impact on the directional response of tank vehicles. SAE Technical Paper,(No. 932942), 1993

3Ranganathan R,Ying Y, MilesJ B. Development of a mechanical analogy model to predict the dynamic behavior of liquids in partially filled tank vehicles,SAE Technical Paper,(No. 942307), 1994

4Salem M I. Rollover stability of partially filled heavy-duty elliptical tankers using trammel pendulums to simulate fluid sloshing[Ph.D Thesis]. West Virginia: West Virginia University, 2000

5Zheng X, Li X, Ren Y. Equivalent mechanical model for lateral liquid sloshing in partially filled tank vehicles.MathematicalProblemsinEngineering, 2012:123～132

6Celebi M S, Akyildiz H. Nonlinear modeling of liquid sloshing in a moving rectangular tank.OceanEngineering, 2001,29(12):1527～1553

7Ibrahim R A, Pilipchuk V N, Ikeda T. Recent advances in liquid sloshing dynamics.AppliedMechanicsReviews, 2001,54(2):133～199

8Ibrahim R A. Liquid Sloshing Dynamics: Theory and Applications. Cambridge University Press, 2005

9Aliabadi S, Johnson A, Abedi J. Comparison of finite element and pendulum models for simulation of sloshing.Computers&Fluids, 2003,32(4):535～545

10 Rebouillat S, Liksonov D. Fluid-structure interaction in partially filled liquid containers: a comparative review of numerical approaches.Computers&Fluids, 2010,39(5):739～746

11 Bogomaz G I, Markova O M, Chernomashentseva Y G. Mathematical modelling of vibrations and loading of railway tanks taking into account the liquid cargo mobility.VehicleSystemDynamics, 1998,30(3-4):285～294

12 Vera C, Paulin J, Suarez B, et al. Simulation of freight trains equipped with partially filled tank containers and related resonance phenomenon.ProceedingsoftheInstitutionofMechanicalEngineers,PartF:JournalofRailandRapidTransit, 2005,219(4):245～259

13 Wang L, Octavio J R J, Wei C, et al. Low order continuum-based liquid sloshing formulation for vehicle system dynamics.JournalofComputationalandNonlinearDynamics, 2015,10(2):021022

14 Wei C, Wang L, Shabana A A. A total Lagrangian ANCF liquid sloshing approach for multibody system applications.JournalofComputationalandNonlinearDynamics, 2016,10(5):051014

15 马驰骋,张希农,柳征勇等. 变质量贮箱类流固耦合系统的振动响应及时频特性分析. 振动与冲击, 2014,33(21):166～171 (Ma C C, Zhang X N, Liu Z Y, et al. Dynamic responses and time-frequency feature analysis for a fluid-structure coupling system with a variable mass tank.JournalofVibrationandShock, 2014,33(21):166～171 (in Chinese))

16 刘富,童明波,陈建平. 储箱内液体晃动动力学分析及防晃结构优化. 计算机应用与软件, 2011,28(12):202～205 (Liu F, Tong M B, Chen J P. Inside-container liquid sloshing dynamics analysis and sloshing suppression structure optimization.ComputerApplicationsandSoftware, 2011,28(12):202～205 (in Chinese))

17 王勇. 考虑液体晃动的三大件转向架罐车耦合系统动力学性能研究[博士学位论文]. 成都:西南交通大学, 2004 (Wang Y. Study on the coupling system dynamics for three-piece bogie tank cars taking into account the liquid sloshing[Ph.D Thesis]. Chengdu, Southwest Jiaotong University, 2004 (in Chinese))

18 卢军. 任意充液比油罐车液体晃动及整车横向稳定性研究[硕士学位论文]. 成都:西南交通大学, 2009 (Lu J. Liquid sloshing and transverse stability analysis of tank truck with arbitrary proportional liquid[Master Thesis]. Chengdu:Southwest Jiaotong University, 2009 (in Chinese))

19 Shabana A A, Hussien H A, Escalona J L. Application of the absolute nodal coordinate formulation to large rotation and large deformation problems.JournalofMechanicalDesign, 1998,120(2):188～195

20 田强,刘铖,李培等. 多柔体系统动力学研究进展与挑战. 动力学与控制学报, 2017,15(5):385～405 (Tian Q, Liu C, Li P, et al. Advances and challenges in dynamics of flexible multibody systems.JournalofDynamicsandControl, 2017,15(5):385～405 (in Chinese))

21 Olshevskiy A, Dmitrochenko O, Kim C W. Three-dimensional solid brick element using slopes in the absolute nodal coordinate formulation.JournalofComputationalandNonlinearDynamics, 2014,9(2):021001

22 Shabana A A. Computational continuum mechanics,3rd ed. Cambridge, UK: Cambridge University Press, 2017