基于流固耦合的罐式半挂列车液体模型构建研究

周淑文，王春晖，张国胜，杨贵军，张思奇3

（1.东北大学机械工程与自动化学院，辽宁 沈阳，110819；2.交通运输部公路科学研究院，北京 100088；沈阳建筑大学交通与测绘工程学院，辽宁 沈阳，110168）

0 引言

目前，我国石油和化工产业的高速崛起与发展为我国工业生产制造、民生建设等各个方面提供了大量必需的能源和产业原材料以及药物制品，在基础产业建设、交通运输、医疗等相关领域做出了十分重要的贡献。数据显示，2020年我国危险品日运输量超300万吨，其中近80%通过道路运输由危货车完成。截至2020年底，全国共有危化品罐式运输车辆18.76万辆，占比51.51%。液体货物运输需求量大，这也带来了更多的相关交通事故，2012年至2018年发生车祸共2 860起，造成1 911人死亡，我国发生涉及危货品运输的事故以及事故人员伤亡统计呈逐年增加的趋势，且并未得到有效的控制[1]。

罐式半挂车重量大且质量分布不均，而且在行驶过程中液体会剧烈流动与车体发生碰撞，导致对该种液体运输车动力学建模十分困难，难以建立十分精确的动力学模型[2]。在罐式半挂车的研发过程中，建立一个准确的动力学仿真模型将极大提高相关数据计算和分析的便捷性和准确性。

针对罐内液体动力学，国内外已有众多学者和专家进行研究。国外学者对于罐内液体晃动的分析研究最早是以外形相对规则且体积有限的罐体为研究对象，且起步较早，根据相关文献的记载从20世纪50年代便开始有对罐内液体动力学特性的公式推导计算研究。采用的方法有流体动力学法、等效机械模型法、准静态方法以及试验和仿真法等。流体动力学法是将流体动力学理论和充液系统动力学方法相结合以解析式的方式求解液体晃动的问题。等效机械模型法是根据液体晃动的运动状态，搭建运动状态相似的机械模型来模拟其运动行为。2000年，Mohamed I.Salem[3]运用仿真动力学软件对罐车侧倾稳定性进行了分析，并且通过搭建机械钟摆模型模拟水平圆柱截面罐内液体晃动的运动状态。2010年，Casasanta[4]通过建立钟摆模型模拟椭圆截面罐内液体晃动，以此来研究罐车在罐体不同装载工况下行驶时车辆的侧倾稳定性。

国内对罐内液体晃动运动特性的研究起步较晚，最初起于20世纪80年代，所采用的理论方法与国外相仿。2015年，张海涛[5]建立非满载液体在晃动幅度较大的情况下的流固耦合振动模型，通过高精度算法拟合了耦合体系的液体晃动参数。2018年，吴相稷等[6]在前人的理论基础上，通过求解多质量椭圆摆模型推导出液体晃动时各参数之间的函数关系。

本文基于流固耦合理论，建立等效力学模型来近似模拟液体晃动，可以在保证模型准确描述罐内液体受力变化的前提下，简化计算过程。首先根据能量守恒定律和牛顿第二定律建立起流固耦合方程，并将其转化为矩阵，简化求解；然后通过对不同充液高度的非满载罐中液体的晃动和液体质心运动轨迹进行研究，建立椭圆摆的运动方程，得到罐内液体质心运动的等效机械运动模型的动态轨迹；最后采用MATLAB的ODE算法进行求解，得到振荡频率和振幅的关系。

1 流固耦合模型的建立

车辆运行过程中，车辆横向位置的改变导致罐内液体产生了一定幅度的晃动，由此对罐体产生的冲击力（非线性力）对行驶过程的稳定性会有较大的影响[7-9]。为研究液体晃动对罐式半挂车行驶中横向稳定性的影响，最有效的方式就是采用流固耦合理论[10]。

1.1 流体控制方程

由流体力学守恒理论可知，最重要的一个前提是满足质量守恒、动量守恒、能量守恒。恒温恒压下对于正常可压缩的流体，守恒方程为：

质量守恒方程：

动量守恒方程：

式中：t为时间；ff为体积力矢量；ρf为流体密度；υ为流体速度矢量；τf为剪切力张量。

式中：p为流体压力；μ为动力黏度；e为速度应力张量。

1.2 固体控制方程

由牛顿第二定律推导固体部分的守恒条件：

式中：ρs为固体密度；σs为柯西应力张量；fs为体积力矢量。

式中：λ为导热系数；SE为能量源。

受材料性质的限制，对固体分析时需考虑温度造成的影响：

式中：aT为与温度相关的热膨胀系数。

1.3 流固耦合方程

在耦合效应产生的过程中，满足各介质之间的守恒关系是流体部分和固体部分的耦合面形成的前提。所以根据守恒关系，建立方程如下：

式中：τ为应力值；q为热流量；d为位移；T为温度。

为了提高求解效率，可以采用矩阵求解方式，具体方程如下：

式中：k为迭代步长；Aff为流体部分系统矩阵；为流体部分所求结果；Bf为流场外部作用力；Ass为固体部分系统矩阵；为固体部分所求结果；Bs为固体的外部作用力；Asf为流体部分耦合矩阵；Afs为固体部分耦合矩阵。

2 液体晃动等效力学机械模型建立

2.1 液体晃动等效力学机械模型的数学形式

理论分析和试验研究表明，用自由液面振荡描述的一阶晃动模态是罐体内最重要的液体晃动模态。因此，我们从研究不同充液高度的非满载罐中液体的晃动和求解液体质心运动轨迹开始研究。

如图1所示，当a/b=1时，横截面为圆形；当a/b 〉 1时，横截面为椭圆形。图1中，a是罐体宽度的一半，b是罐体高度的一半，h0是液位与y轴的交点，Φ是液表面的倾斜角度。

图1 截面为圆形或椭圆形的非满载充液罐示意图Fig.1 Schematic diagram of non-full-loaded liquid-filled tank with circular or elliptical cross-section

定义液位高度与油箱高度的比值为油箱内的充液百分比，可表示为：

液体质心的轨迹可由下式求得：

罐内液体的横截面积可表示为：

当h0＜0时，

当h0＞0时，自由液面与y轴的交点（定义为h）随着自由液面倾角的不同而变化。因此，（x，y）坐标系下的自由液面可以表示为：

自由液面与罐壁的交点可以表示为：

液体的横截面积定义为S，其在x轴和y轴上的静力矩由式（12）和（15）可得。无论自由液面倾角如何，液体的横截面积保持不变。使h在给定的范围内以相当小的步长变化，并计算S及其在每个h值处的静态力矩。使用式（11）可得到液体的质心CG，确保满足|S-A|≤δ（δ是一个非常小的正值，取决于h步长）。

当自由液面倾角在适当范围内变化时，液体质心的轨迹如图2所示，这表明液体质心的轨迹保持与容器外围平行。

图2 液体质心的运动轨迹Fig.2 The motion of the GG of liquid

由此可见，使用液体质心的变化来表示罐内液体的晃动效果是完全合理且符合后续的公式计算要求的。

相比较其他简单机械而言，单摆和弹簧机构在描述横向液体晃动时准确性更高[11]。所以本文采用单摆和弹簧机构来描述横向液体晃动。摆线的长度和摆点的坐标位置是摆球运动轨迹的决定性因素，本文所研究的椭圆运动轨迹需要在此基础上对相关参数有所改变，计算摆线的长度和摆点位置的变化范围，并推导其变化规律，得到的模型便可以准确地反映出截面为椭圆的罐体内液体晃动时质心的运动路线。

2.2 罐内液体晃动椭圆摆运动方程

车辆在行驶过程中其运动形式是多变的，这也就导致了其所连接的罐内液体的运动状态并不是单一形式。所以通过建立椭圆摆的运动方程的方法来研究罐内液体晃动的速度变化、频率特性、简谐运动等相关运动特性[12]。

椭圆摆的振荡轨迹和基本参数如图3所示。假设摆的振荡轨迹与液体质心不重合，ap+bp为摆臂的长度，其中ap是摆振荡轨迹主轴的一半，bp轴是小轴的一半。acg是液体体积的重心椭圆轨迹长轴的一半，bcg是其小轴的一半。θ是摆角，α是连接原点和摆锤质量的直线与y轴的夹角。

图3 椭圆摆振荡轨迹示意图及基本参数Fig.3 Schematic diagram and basic parameters of oscillation track of elliptical pendulum

罐体周边、摆点摆动轨迹和液体重心均平行，可表示为：

摆运动分析如图4所示，其中xy为罐车坐标，xγ为接地固定坐标。x是XY原点和xγ原点之间的距离。

图4 摆运动分析示意图Fig.4 Schematic diagram of pendulum motion analysis

由图4可知，摆点质量的绝对位置可表示为：

因此，摆点质量的速度和加速度可表示为：

摆质量的动能定义为：

假设势能的零点位于摆平衡位置的表面。因此，摆的重力势能可以表示为：

根据式（20）和（21），利用拉格朗日函数可以得到摆系统的动力学方程，其表达式为：

椭圆摆系统的运动可以表示为：

将式（25）代入式（23）得到：

采用MATLAB的ODE算法求解式（26）。在求解过程中，我们使a/b在0.25步长下在1到2之间变化，摆幅在10步长下在10到180度之间变化。

图5显示了具有小振幅和大振幅的摆振的振荡频率和角速度。不同截面的罐体截面面积相同，当截面为圆形时（a/b= 1），a = b =0.360 2 m。

如图5（a）所示，摆振频率与摆臂长度和幅值有关，当a/b=1时，摆振频率随幅值的增加而减小。而其他摆的振荡频率随着振幅的增大而增大，当振幅达到一定值时频率达到最大值，之后频率减小。例如，当a/b= 2时，最大频率出现在90度的振幅处。对于所有的摆，当振幅低于某一值时，其振荡频率几乎保持不变。

如图5（b）～（d）所示，随着振幅的增加，摆的运动变得不规则，非线性程度增加，特别是当振幅大于170度时。所幸的是，为了避免车辆翻车，实际半挂车的侧向加速度小于0.45g，即自由液面倾角始终小于90度。因此，摆可以近似地假设为线性。

图5 椭圆摆的运动特性Fig.5 Motion characteristics of elliptical pendulum

3 结语

本文基于流固耦合理论推导出了罐内液体晃动等效机械模型，建立椭圆摆的运动方程，得到罐内液体质心运动的等效机械运动模型的动态轨迹。最后采用MATLAB的ODE算法进行求解，得到振荡频率和振幅的关系，发现在一定加速度范围内，摆可以近似地假设为线性。下一步工作将在此基础上对罐式半挂车侧翻特性及预防控制方法进行研究，以进一步提高罐式半挂车的抗侧翻性能。