CN107436981A - 车载金属油罐内部油料晃动数值模型的建立方法 - Google Patents

Abstract

本发明公开了一种车载金属油罐内部油料晃动数值模型的建立方法，包括油料晃动基本方程的建立，油料晃动的边界条件的确定，固体结构区域受力分析等步骤。运用VOF法和动坐标系建立了车载金属油罐内部油料晃动三维仿真模型，罐体、内部介质及防波板均采用实际尺寸建立三维模型，并可以实现到固体罐壁区域受力、变形的耦合分析。将罐车制动过程中制动减速度变化所对应的惯性力函数，加载到动量方程中求解，既节约了计算时间，又保证了较高的精度。本发明主要可以用于研究改善越野金属油罐车紧急转弯时的侧翻稳定性和紧急制动时罐体的结构强度。

Description

车载金属油罐内部油料晃动数值模型的建立方法

技术领域

本发明涉及金属油罐中油料晃动的分析领域，具体是一种车载金属油罐内部油料晃动数值模型的建立方法。

背景技术

车载金属油罐是运加油车的重要组成部分，需通过一定的方式固装于底盘车架上，分析其受力状态和计算结构强度非常困难，使用中常发生罐体开裂或者将罐体设计得很笨重，耗材增加，造成不必要的浪费等现象。其主要特点是：

(1)油罐安装在底盘上，车辆在行驶中罐体随底盘处于运动状态，油罐的受力情况不同于静态储油的卧式金属油罐，除承受气压、液压、罐体自重外，还要承受车辆紧急制动和道路颠簸引起的动载荷。

(2)受底盘载重量、牵引力及道路条件的制约，油罐额定容量较小，多数在20m³以下。油料在罐内存放时间不长，收发油料较频率。

(3)罐内油料在紧急制动或转弯时大幅晃动，对罐体结构强度及侧翻稳定性有较大影响。

油罐车在行驶当中，由于路面情况的变化，驾驶员的方向操纵，紧急情况的处理等，车辆行驶状态会发生改变，车辆行驶的不稳定与越野车载金属油罐内的液体油料产生相互作用，引起罐内液体晃动，对罐壁产生冲击，同时引起液体重心偏移，汽车质心变化，产生侧翻弯矩，造成轴荷分配不均，极易发生翻车事故。一旦晃动频率与底盘振动频率相交耦，还会产生共振，造成罐体破坏。因而越野金属油罐车在行驶中存在安全隐患，其侧翻时极易发生油品爆炸和泄露，产生巨大的危害和损失。国内外的大量报道都指出，充液油罐汽车翻车的概率比普通载货汽车大得多，例如：2009年10月28日，成绵高速一油罐车侧翻漏油，引起燃烧爆炸，造成数人遇难；2011年3日下午7点左右，一辆油罐车因交通事故发生侧翻造成200多人死亡。

当油罐车在行驶过程中紧急制动、突然加速、急转弯或者在凹凸不平的越野道路上颠簸时，车载金属油罐内部液体介质会发生晃动(Sloshing)，它产生的原因是外界对充液油罐施加的扰动，内部液体的晃动的形式也是多种多样的，如自上而下晃动、大幅晃动和混沌晃动等。车载金属油罐内部自由液面的位置在晃动过程中是时刻改变的，并且自由液面边界条件非常复杂，在晃动过程中可能产生碎波分离，这在实质上是复杂的非线性问题，对它的研究和求解非常困难。液体油料晃动的过程研究涉及动压力分布、作用力、力矩和自由液面固有频率等多个方面，这些参数直接影响到运加油车行驶的稳定性和车载金属油罐的性能，如可能导致油罐结构破坏，出现漏油现象。

对流体运动规律进行计算的学科即计算流体动力学(CFD)，它属于计算力学中的一个学科分支，CFD也是研究液体晃动的主要方法。CFD是根据液体流动的三大基本方程(质量守恒方程、动量守恒方程和能量守恒方程)，对液体流动情况进行数值模拟，其基本方法可以归纳为：用有限个离散点上的变量值的集合代替原来在时间域及空间域上连续的速度场和压力场等，根据这些量需要满足的基本物理关系，列出这些离散点上的代数方程组，然后对这些代数方程组进行求解以获得场变量的近似值。

经过四十多年的发展，CFD出现了多种数值解法，其分类方法有多种，按照对流体运动的坐标描述，可以分为欧拉法、拉格朗日法和任意拉格朗日-欧拉法；按照对控制方程的离散方式，可以分为有限差分法、有限元法、有限体积法和边界元法等；按照对流体自由面的计算方式，可以分为移动网格法、LEVEL-SET法、VOF(The Volume-Of-Fluid Method)法、和MAC(The Marker-And-Cell Method)法等。

①对流体运动进行描述时，根据我们采用的坐标网格是静止不动还是随流体运动，可以分为欧拉法、拉格朗日法和任意拉格朗日-欧拉法。

欧拉法又称为空间点法，其原理是设法在空间中的每一点上描述出流体运动随时间的变化情况，即流体的速度场、压力场、密度场等。采用欧拉法描述时，采用的是空间位置在流体运动过程中保持不变的网格。正是由于欧拉网格的位置、形状不会改变，当液体晃动时，对自由液面的追踪需要结合采用其他跟踪边界的方法，如VOF法和LEVEL-SET法等。

拉格朗日法又称为质点法，它是把流体域看做由许许多多的流体质点组成，只要知道每个流体质点在整个过程中的运动规律，那么，流动的规律也描述出来了。采用拉格朗日法描述时，可以无网格化，如果有空间网格，网格区域在流体区域，并且和流体一起运动。因此，初始时哪些网格节点所处位置是液体自由界面上，在运动过程中它们还是，可以比较精确地描述自由液面。采用拉格朗日法求解时，方程中不存在对流项，简化了数值求解过程，并且比较容易考虑表面张力效应。

采用拉格朗日法描述车载金属油罐油料晃动问题的缺点是当油料在晃荡过程中变形量较大时，网格容易产生畸变、扭曲，由此会增加数值计算的误差，甚至会出现求解失败，因此它一般只适用于液体的波动幅度较小，流动情况是低Reynold数流动以及无旋流动的问题，对于液体晃动幅度较大，流动范围很大的问题求解时，必须同时进行网格重构以改善网格的扭曲程度。

任意拉格朗日-欧拉法(ALE)是将拉格朗日法和欧拉法结合起来以充分利用二者的优势，可以处理大变形问题或移动的边界。在这种方法下，网格节点的运动不再与质点运动同步，计算网格在每一个时间步长上进行调整，既保持空间上的合理布置使计算域的形状不变，又要保证能追踪求解移动的自由界面边界。任意拉格朗日-欧拉法最先由Hirt和Donea进行研究；岳宝增，李笑天等人将其用于二维和三维自由液面大幅晃动问题的模拟；周宏用任意拉格朗日-欧拉法对充有液体的容器系统进行了流固耦合的运动学问题研究；王勇利用ALE方法分析了汽车燃油箱内的液体晃动问题；陈建平采用ALE方法分析了载液系统内的液体晃动问题；管延敏对水平激励下二维矩形箱体液体晃荡现象进行了数值计算，研究了不同布置形式的挡板对液体晃动情况的影响。胡振兴采用ALE有限元法对大型立柱式储液罐在地震载荷作用下的动力响应进行了分析。任意拉格朗日-欧拉法由于在计算过程中需要及时进行网格重构，计算效率不高。

②液体运动的控制方程是对时间和空间的微分方程，微分方程的数值求解需要对其进行时间和空间上的离散，按照在空间上离散方式的不同，主要方法有有限差分法、有限元法、有限体积法、边界元法和光滑粒子流体动力学(SPH)方法等：

有限差分法简单、使用简单，在网格划分较细密时，计算结果有较高的精度，主要适用于结构化的正交网格。

端木玉等采用有限差分法对不同液舱结构形式对晃荡所产生的影响进行了分析研究,计算分析了不同结构形式液舱的防晃效果；Jannette B.Frandsen用非线性有限差分法研究了罐体内液体晃动现象。

有限差分法可以应用到具有比较复杂边界的问题，但是，对于处理带有防波板的越野车载金属油罐内部液体晃动这类复杂边界的流场还比较困难。

有限元法是通过变分原理或拟变分原理将待分析的连续体分成一系列连续单元，在每个单元上通过所选择的形函数建立相关的物理量方程，然后组合所有单元上的方程形成方程组，最后对方程组求解，得到每个单元节点上的未知量。有限元法在结构力学中运用得十分广泛，由于流体力学中的微分方程大多是非线性的，不具有自伴性，有限元法在流体力学中的应用进展缓慢，直到上世纪七八十年代有限元法才逐渐运用于流体力学中。有限元法采用的是复杂的形函数，难于处理流体控制方程中的非线性项。有限元法的优点则在于可以使用任意划分的网格，其对网格的要求不如有限差分法严格，在分析流固耦合问题时可以进行直接双向耦合，但是对于大型问题的分析，计算量和计算时间都很长，并且，能更准确反映流体特性的单元模块还在不断研究更新中。

贾善坡等采用有限元法研究分析了圆柱形容器内三维液体晃动的固有频率和模态问题；S.Mitra和K.P.Sinhamahapatra用二维有限元法分析了充液容器的流固耦合问题，其中液体假设为无粘性、不可压流体；K.C.Biswal等采用有限元法分析了圆柱形罐体内柔性防波板对液体晃动的影响，其中液体假设为无粘性、不可压流体且无旋。

有限体积法的原理是在计算区域上划分网格，对网格节点分配相应控制体积，将待求微分方程在所划分网格上对每个控制体积积分得出离散方程，这些离散方程一般都能自动满足流体在控制体积内的质量、动量、和能量守衡的要求，因此对整个求解区域也满足守恒。有限体积法其实属于介于有限元法和有限差分法之间的一种方法。有限体积法是目前CFD方法中应用最广泛、计算效率和精度都较高的一种。

K.P.Thiagarajan采用有限体积法计算了二维空间内不同充液高度和不同频率外加晃动激励下的液体晃动现象及冲击力，充液高度变化范围：10％～90％，外加晃动激励周期变化范围：0.8～2.8s。

边界元法又称边界积分法，其基本思想是：将流体运动的控制方程转换为流体域边界上的积分方程，然后在边界上引入有限个积分单元，进行积分方程的离散，由于积分方程的离散仅在边界上进行，离散后的方程仅含有边界上的节点未知量，可以降低所需求解的问题的维数。边界元法应用时，可以直接对边界上流体速度等物理量求解，也可以在积分面和流域边界上布置某种类型的虚拟分布的奇点，然后求解相应未知奇点强度的封闭方程组，得到流体区域各点处的速度势和速度的反映。

边界元法的缺点是难以处理粘性流问题，由于粘性在准确分析流体运动情况时是不可忽略的因素，一些研究采取了引入人工粘性的方法来获得稳态解，但是在采用人工粘性时，粘性系数如何取值还没有较好的通用方法。边界元法在处理不需要考虑液体的粘性，容器内液位较深，且容器形状简单的液体晃动问题时能得到较好的计算结果。

王小贞采用了边界元法研究了容器中粘性、不可压缩液体的小幅晃动问题。在频域内考虑二维线性化Navier-Stokes方程，以问题的物理变量作为数值分析的未知函数，并推导了该问题分析的边界积分方程；HE Zhao采用多区域边界元法研究了带有孔钢梁的二维和三维罐体内的非线性液体晃动现象。

目前边界元方法还处在发展完善阶段，在边界元方法中所得方程的系数矩阵都是满秩的非对称矩阵，此外在三维晃动问题中还需考虑边界积分方程中的奇异积分处理问题，用边界元方法求解带自由液面的三维晃动问题还处在进一步的研究发展之中。

光滑粒子流体动力学(SPH)方法是采用粒子集合和插值核函数来估算空间函数及其导数，采用离散粒子的求和来代替方程的积分运算。

刘汉涛等应用SPH法求解了拉格朗日形式的N-S方程。对近岸水动力学中波浪传播、变形、与建筑物的相互作用问题进行了二维和三维数值计算；刘富等采用SPH方法对棱形液舱在外加激励作用下，舱内液体在不同充液系数情况下的晃动特性进行了三维数值模拟,发现液舱的充液系数及晃动周期对于液体的晃动特性及液体作用在舱壁上的压力有重要的影响。

SPH法主要集中于二维晃动的数值计算，目前关于液体三维晃动的研究较少，并且如何确定最优的插值核函数也还有大量研究要做。

③分析液体晃动自由液面十分重要，典型的方法有移动网格法、VOF法、LEVEL-SET法和MAC法等。

移动网格法属于界面跟踪方法，由于网格跟随界面移动，因此，它通常与拉格朗日型计算网格联用，边界上任意质点x_i的位置可以通过质点速度在时间域上的积分获得，即:

随着边界质点x_i的运动，网格节点的位置也会发生变化，这也是动网格法得名的由来。虽然用这类方法获得液体自由界面时计算较简单，但是，通常在数值计算时，需要采用很小的时间步长，并且如果液体运动变形较大，会导致网格随之发生大的变形，从而不能保证控制方程的离散精度。

后面三种方法则是属于界面捕捉方法，它是在固定不动的网格中进行，网格区域不限于流体域，自由表面的获得是通过对液体自由界面附近的网格中液体所占比例的计算得到。

VOF法和LEVEL-SET法在本质上是相似的，都是先在计算网格中定义一个具有以下性质的函数：

这里假设C₂<C₁，函数的解则通过求解偏微分方程(1.3)得到。

式中，是流体的速度。

Level-set方法中采用的构造函数C是连续的函数代入式(1.3)，可得：

对(1.4)式进行空间和时间上的离散，求解Level-set函数就可得到液体界面位置。标准的Level-set方法存在计算时质量不易保持守恒这一问题，Sussman将其与VOF方法联合运用，联合算法克服了Level-set方法单独使用时的质量不守恒问题，提高了计算精度。

VOF方法是最为著名的界面捕捉方法，它采用的构造函数是离散的体积比函数C，C函数定义为一个单元中液体所占体积的比值，单元网格全充满液体时其比值为1，单元网格中没有液体时其比值为0，因此边界上的单元其C函数的值在(0,1)区间内。把C函数代入式(1.3)并求解，可以得到C函数在流体边界单元中的情况，C的数值表示液体单元内的体积比值，它的微分值表示边界的法线，有了这两个数值，边界面在单元中的位置也就能就确定了。VOF方法在捕捉流体界面时十分灵活稳健，它能适应复杂的几何形状和液面翻转、碎波等情况。

Kang等人采用VOF模型研究了制动过程中罐车罐体内液体晃动问题，对不同充液比和不同防波板面积下的罐体的受力情况进行了比较分析，在论文中，

忽略了由防波板隔开的各腔室之间的相互作用，将罐车简化为单室模型进行研究；

Dongming Liu等采用VOF模型和大涡模型研究了带自由表面的三维液体晃动现象，通过实验验证了数值模型在激励强时比解析解更接近实验结果；Wemmenhove等采用改进的VOF模型计算了两维情况下液体晃动产生的动压力。

MAC方法是在求解的初始时刻把一些无质量的粒子(标记子)置于液体自由液面上，在整个流体运动的求解过程中，这些标记子所标记的即是流体的自由液面。它是以压强和速度为基本变量计算自由界面流动的界面捕捉方法，MAC法的优点在于基本概念十分简单，可以处理极其复杂的液面情况，例如碎波；但是，其计算量也相当大，尤其在三维计算的情况，因为它不仅要计算基本控制方程，还要对标记子进行追踪，计算存储量也很大。

油罐车在行驶过程中，伴随着制动、转弯等驾驶员的操纵行为，罐内液体产生晃动，对金属罐体施加额外的晃动载荷，影响整个车载系统的运动轨迹，常导致结构破坏或运行失稳。对于这类问题的分析主要有两大类方法，一类是采用等效力学模型，另一类是基于计算流体力学(CFD)理论来建立充液系统的液体晃动模型。

等效力学模型主要有质量弹簧模型和等效摆模型等，建立等效力学模型时应遵循的等效原则为：等效时，等效系统的质量与所需等效的液体质量相等，质心所在位置保持不变，等效系统的固有频率与液体晃动频率一致，所建立的等效系统形成的作用力和作用力矩与贮液容器内液体晃动产生的作用力与作用力矩分别相等。等效摆模型在航天领域应用较广泛，主要用来分析计算航天器所载液体燃料的小幅晃动问题。质量弹簧模型是用质量弹簧系统来模拟液体晃动时对容器壁的冲击力。

等效力学模型研究的关键在于如何建立由质量、弹簧、摆长以及阻尼等机械参数构成的等效力学模型。Dodge研究了矩形、圆柱形和椭圆形箱体中液体晃动等效力学模型的参数确定；Chu研究了能更准确地反映液体晃动特性的单摆等效模型；Unruh建立了一个比简单单摆模型更完善的等效的复合球摆模型，用以模拟旋转腔体内液体的晃动，该模型考虑了充液系统绕对称轴转动以及径向摆动两个自由度，并且能够计及系统旋转运动时的离心力；谢普根据势流理论建立了半挂式液罐车罐体内液体晃动的非线性模型和线性单摆阻尼模型，根据等效原则确定了单摆模型各项参数，并将此模型与液罐车刚体动力学模型结合，建立液罐车刚一液耦合系统动力学模型，分析了液体晃动对半挂液罐车稳定性的影响；卢军基于势流理论研究了具有任意充液比的水平放置圆柱罐体内的液体晃动特性，将流体作用在罐体上的力和力矩分为液体惯性项和液体晃动项，结合车辆动力学研究了液体晃动对铁路罐车蛇行运动稳定性的影响。

等效力学模型的优点在于将液体部分和固体结构部分作为一个整体进行综合考虑，降低了系统分析的复杂性，减少了数值计算的工作量，因此，可以采用相对比较简单的分析求解方法实现对原来复杂系统的处理，在实际工程中应用较多。但此种模型通常只能解释主导阶的流体动力学行为，对于较高阶频率的动力学行为可能较难把握。在运用等效力学模型时，弹簧刚度、摆长、质量、惯性系数等等效力学模型参数的求解对整个分析的准确性至关重要，这些参数的获得，既有运用传统势流理论和实验得出，又有根据液体晃动过程的数值模拟来确定的。

基于计算流体力学(CFD)理论来建立充液系统的液体晃动模型，把液体对罐壁的冲击力传递到固体区域求解，能更真实准确地反映液体的晃动对固体结构产生的冲击。在分析计算时，不需要再对液体建立等效力学模型和对等效力学模型参数求解，只需要利用数值方法直接求解描述液体晃动的控制方程即可。

Vaibhav Singal采用VOF多相流模型研究了车载煤油罐体内液体晃动对防波板及罐体的压力情况；M.Toumi等研究了载液罐车的稳定性问题；Rumold通过对不可压缩液体的流动方程和车辆动力学模型进行迭代求解，进行了车辆在转弯和紧急制动时的液体晃荡研究。

发明内容

本发明建立了越野车载金属油罐内部液体晃动的数值计算模型，能够准确反映油罐车在制动时的变加速度工况和转弯工况；建立了越野车载金属油罐流固耦合受力分析模型，用以计算得出越野金属油罐车转弯和制动时罐壁和内部防波板的应力和侧翻力矩随时间变化的规律。

为了实现上述目的本发明采用如下技术方案：车载金属油罐内部油料晃动数值模型的建立方法，包括以下步骤：

(1)建立油料晃动基本方程

建立固定于大地上的固定坐标系x₀,y₀,z₀，和以油罐中心为坐标原点与油罐一起运动的运动坐标系x,y,z；

建立随罐体运动的运动坐标系x,y,z下罐内液体的运动控制方程为：

连续性方程

动量方程

流体体积分数输运方程

式中ρ为流体密度，为固定坐标系x₀,y₀,z₀下流体的绝对速度；为单位质量上的流体质量力函数；μ为流体动力学粘性系数；表示在运动坐标系下的流体速度矢量；C为流体体积分数；为矢量算子，为单位质量流体表面力，为运动坐标系下的流体速度变化矢量，为单位质量的惯性力：

为运动坐标系的旋转角速度矢量；为运动坐标系原点到所研究的液体微团k的矢径，为运动坐标系相对固定坐标系的平动加速度；为运动坐标系相对固定坐标系的转动加速度，即

(2)确定油料晃动的边界条件

油料边界包括两个部分，一部分是油料与油罐交界的固壁Ω_g，另一部分是油料与空气交界的自由表面Ω_f；

其中在油罐壁上油料的流体速度满足：

式中，为罐体在运动坐标系下相对于运动坐标系坐标原点运动的速度，忽略油罐罐壁的变形量，则

油料与空气交界的自由表面Ω_f的表面应力条件为：

在Ω_f上

式中，p₀为液体表面大气压力；v_n为垂直于液体自由表面的法向速度，以流体域外法线方向为正；v_τ为切向速度；p为流体表面力；μ为流体动力学粘性系数；τ为剪应力；n为垂直于液面的单位矢量；

(3)固体结构区域受力

根据油料晃动基本方程和油料晃动的边界条件求解出罐内液体晃动冲击力后，把它加载到油罐罐壁上进行结构强度的计算，通过平衡微分方程、几何方程及本构方程在固体网格区域进行有限元法求解，得到油罐罐体的应力、应变、位移参数。

平衡微分方程、几何方程及本构方程分别为：

平衡微分方程：

式中，σ_x、σ_y、σ_z分别表示x、y、z方向的正应力分量，τ_xy、τ_yz、τ_zx分别表示x、y、z方向的剪应力分量，f_x、f_y、f_z分别表示x、y、z方向的体积力分量，τ_xz、τ_yx、τ_zy分别表示x、y、z反方向的剪应力分量。

几何方程：

式中，ε_x、ε_y、ε_z分别表示x、y、z方向的正应变，γ_xy、γ_yz、γ_zx分别表示x、y、z方向的剪应变，u、v、w分别表示x、y、z方向的位移分量；

本构方程：

式中，E表示弹性模量，υ表示泊松比，G表示剪切模量。

本发明运用VOF法和动坐标系建立了车载金属油罐内部油料晃动三维仿真模型，罐体、内部介质及防波板均采用实际尺寸建立三维模型，并可以实现到固体罐壁区域受力、变形的耦合分析。将罐车制动过程中制动减速度变化所对应的惯性力函数，加载到动量方程中求解，既节约了计算时间，又保证了较高的精度。本发明其主要应用是可以用于研究改善越野金属油罐车紧急转弯时的侧翻稳定性和紧急制动时罐体的结构强度。

本发明主要优点包括，通过建立的车载金属油罐内部油料晃动三维仿真模型，还可以采用c语言编写真实反映油罐车紧急转弯和紧急制动过程中所对应的惯性力函数程序，并加载到动量方程中求解罐体的应力、应变、位移和冲击力和侧翻力矩等参数。可以实现油罐罐壁区域受力、变形的耦合分析，研究越野金属油罐车在不同充液系数紧急转弯时的侧翻稳定性和紧急制动时罐体的结构强度，提高运加油车行驶稳定性和安全性。

附图说明

图1为车载金属油罐坐标图；

图2为50％充液系数流态图；

图3为50％充液系数油料对罐体冲击力；

图4为50％充液系数罐体最大压强；

图5为50％充液系数侧翻力矩；

图6为50％充液系数罐体最大应力；

图7为60％充液系数流态图；

图8为60％充液系数油料对罐体冲击力；

图9为60％充液系数罐体最大压强；

图10为60％充液系数侧翻力矩；

图11为60％充液系数罐体最大应力；

图12为70％充液系数流态图；

图13为70％充液系数油料对罐体冲击力；

图14为70％充液系数罐体最大压强；

图15为70％充液系数侧翻力矩；

图16为70％充液系数罐体最大应力；

图17为80％充液系数流态图；

图18为80％充液系数油料对罐体冲击力；

图19为80％充液系数罐体最大压强；

图20为80％充液系数侧翻力矩；

图21为80％充液系数罐体最大应力；

图22为94％充液系数流态图；

图23为94％充液系数油料对罐体冲击力；

图24为94％充液系数罐体最大压强；

图25为94％充液系数侧翻力矩；

图26为94％充液系数罐体最大应力；

图27为30km/h转弯各种充液系数侧翻力矩；

图28为30km/h转弯各种充液系数罐体最大应力。

具体实施方式

1.油料晃动基本方程

车载金属油罐三维模型及所建坐标系统如图1所示，图中坐标系x₀,y₀,z₀为固定于大地上的坐标系，坐标系x,y,z为建立在油罐中心并与油罐一起运动的随体坐标系，在油罐静止状态下，y₀和y的方向都为垂直于地面向上。

油罐内装载的油料可以认为是粘性不可压缩流体，密度ρ为常数，分析时不虑流体与外界热能的传递，在固定坐标系x₀,y₀,z₀下，流体运动控制方程为：

连续性方程

动量方程流体体积分数输运方程

式中，为固定坐标系x₀,y₀,z₀下流体的绝对速度；为单位质量上的流体质量力函数；μ为流体动力学粘性系数；ρ为流体密度；C为流体体积分数。为矢量算子，为单位质量流体表面力，为固定坐标系x₀,y₀,z₀下流体的绝对速度变化矢量。

车载金属油罐固定在罐车底盘上，在罐车匀速行驶时，其相对于固定坐标系是匀速运动，罐内液体与罐壁保持相对静止；在制动、转弯等工况下，油罐相对于固定坐标系的速度在改变，罐内液体相对罐体发生相对运动，产生冲击。如果采用定坐标系建立方程,流体区域划分网格离散后，其边界运动是跟随车体一起，也就是随着时间变化，流体边界位置会发生变化，流体区域的网格会发生变形，在计算过程中需要采用动网格。模拟罐车制动、转弯工况的过程中液体边界位置变化很大,为保证计算精度，应在网格畸变达到一定程度时采取重新划分网格。但是，在计算过程中，如果把控制网格变形的标准取得较严，那么在很小的时间步长内就需要重新划分网格，这对于分析小尺寸和结构简单的容器影响不大，对于所研究的全尺寸的内部防波板纵横交错、结构复杂的越野车载金属油罐的分析，计算量和耗时都巨大；如果变形标准取得较宽松，计算精度又会受影响。为解决这一问题，简化流体运动方程的求解，本发明将运动方程建立在随罐体运动的坐标系x,y,z下，令：

式中，表示在运动坐标系下的流体速度矢量，表示运动坐标系相对于静止坐标系的速度矢量，可表示为：

式中，为运动坐标系的平动速度矢量；为运动坐标系的旋转角速度矢量；为运动坐标系原点到所研究的液体微团k的矢径，由此得出在随罐体运动的坐标系下罐内液体的运动控制方程为：

连续性方程

动量方程

流体体积分数输运方程

其中，是随罐体运动的坐标系为非惯性坐标系统而引入的单位质量的惯性力：

式中，为运动坐标系相对静止坐标系的平动加速度；为运动坐标系相对静止坐标系的转动加速度，即

驾驶员紧急制动时，移动脚踩踏板需要一定的反应时间，随着驾驶员踩踏板的动作，踏板力迅速增大，但是，制动鼓与制动蹄片之间存在间隙，制动器制动力增长过程需经历一定的时间，司机松开制动踏板，制动力消除还需经历一定的时间，制动器的作用时间既受驾驶员踩踏板的速度的影响，又受制于制动器结构，油罐车制动过程中，为随时间变化的数值，为准确描述制动过程，本发明采用c语言编写了真实反映油罐车紧急转弯和紧急制动过程中3.47m/s²、6.37m/s²及8.82m/s²制动减速度变化所对应的惯性力函数程序，并加载到动量方程中求解，既节约了计算时间，又保证了较高的精度。其中，减速度为8.82m/s²时的程序如下：

2.油料晃动边界条件

车载金属油罐内部装载液体的区域为Ω，通常越野金属油罐车装载油料时，最大充液比为94％，机动过程中加油后，油料体积还会减少，因此，油罐内部始终有与空气相交的自由液面存在，在这种情况下，液体边界包括两个部分，一部分是液体与固体交界的固壁Ω_g，另一部分是液体与空气交界的自由表面Ω_f。

由于是粘性流体，在液体与固体交界的固壁Ω_g上本发明采用的是非滑移边界条件，即在固壁上流体速度满足：

式中，为罐体在运动坐标系下相对于坐标原点运动的速度，如果忽略油罐罐罐壁的变形量，则

在自由界面上的运动学条件要求界面上的流体法向速度与网格法向运动速度一致，也就是界面上没有质量通量，动力学条件必须满足表面应力条件，即自由表面两边的法向应力大小相等，方向相反，切向应力大小和方向相同。在本发明模型中，重力和大气压力的影响占主要地位，因此忽略表面张力，此时的表面应力条件为：

在Ω_f上 (11)

式中，p₀为液体表面大气压力；v_n为垂直于液体自由表面的法向速度，以流体域外法线方向为正；v_τ为切向速度；p为流体表面力、μ为流体动力学粘性系数、τ为剪应力、n为垂直于液面的单位矢量。

3.固体结构区域受力计算

根据油料晃动基本方程和油料晃动的边界条件求解出罐内液体晃动冲击力后，需要把它加载到固体罐壁上进行结构强度的计算，流体区域的力传递到固体壁面是通过在流固耦合面(FSI)上进行流体区域的网格节点力到固体区域网格节点力映射得到，加载了液体载荷的固体罐壁结构通过式13～15的平衡微分方程、几何方程及本构方程在固体网格区域进行有限元法求解，还采用c语言编写真实反映油罐车紧急转弯和紧急制动过程中所对应的惯性力函数程序，并加载到动量方程中，最终求解罐体的应力、应变、位移和冲击力和侧翻力矩等参数。

平衡微分方程：

式中，σ_x、σ_y、σ_z表示x、y、z方向的正应力分量，τ_xy、τ_yz、τ_zx表示x、y、z方向的剪应力分量，f_x、f_y、f_z表示x、y、z方向的体积力分量。τ_xz、τ_yx、τ_zy分别表示x、y、z反方向的剪应力分量。

几何方程：

式中，ε_x、ε_y、ε_z分别表示x、y、z方向的正应变，γ_xy、γ_yz、γ_zx分别表示x、y、z方向的剪应变，u、v、w分别表示x、y、z方向的位移分量。

本构方程：

式中，E表示弹性模量，υ表示泊松比，G表示剪切模量。

本发明对7000L车载金属油罐在半径为9m的弯道以30km/h速度行驶时产生的冲击力和侧翻力矩进行仿真研究，分别对充液系数为50％、60％、70％、80％及94％时紧急转弯油料晃动进行了模拟仿真。

(1)50％充液系数仿真结果及工况分析

从运加油车开始转弯至0.72s左右，油料由下向上向转弯外侧涌动，然后又从上往下涌动至1.2s左右，然后油料又由下向上向转弯外侧涌动，罐内油料流态如图2所示，油料对罐体的冲击力如图3所示(曲线从上到下为x、z、y方向)；油料对罐体在Z方向的冲击力较小，对X、Y方向的冲击力较大，在0.36s左右对罐体X方向充击力最大(约为34822N)；在0.24s左右对罐体Y反方向冲击力最大(约为36583N)。

罐体在0.3s左右在罐体转弯外侧与底部交界区域压强最大(约为14480Pa)，同时在油罐内侧前部产生真空(约为6497Pa)，如图4所示，罐体所受侧翻力矩在0.36s左右最大(约为37621N.m)，在0～0.3s左右，侧翻力矩为负(即无侧翻倾向)，侧翻力矩变化情况如图5所示；罐体所受最大应力约为290MPa，应力变化情况如图6所示。

(2)60％充液系数工况分析

从运加油车开始转弯至0.9s左右，油料自下而上向转弯外侧涌动，又自上而下涌动至1.5s左右，然后油料又由下向上向转弯外侧涌动，罐内油料流态如图7所示，油料对罐体的冲击力如图8所示；油料对罐体在Z方向的冲击力较小，对X、Y方向的冲击力较大，在0.36s左右对罐体X方向冲击力最大(约为40658N)，在0.24s左右对罐体Y反方向冲击力最大(约为41995N)。

罐体在0.3s左右在罐体转弯外侧与底部交界区域压强最大(约为18250Pa)，同时在油罐内侧前部产生真空(约为8722Pa)，如图9所示；罐体所受侧翻力矩在0.54s左右最大(约为36929N.m)，在0～0.15s左右，侧翻力矩为负(即无侧翻倾向)，侧翻力矩变化情况如图10所示。罐体所受最大应力约为244MPa，应力变化情况如图11所示。

(3)70％充液系数工况分析

从运加油车开始转弯至0.6s左右，油料由下向上向转弯外侧涌动，然后又从上往下涌动至1.5s左右,然后油料又由下向上向转弯外侧小幅涌动，罐内油料流态如图12所示，油料对罐体的冲击力如图13所示，油料对罐体在Z方向的冲击力较小，对X、Y方向的冲击力较大，在0.24s左右对罐体X方向冲击力最大(约为42600N)，在0.12s左右对罐体Y反方向冲击力最大(约为46855N)。

罐体在0.3s左右在罐体转弯外侧与底部交界区域压强最大(约为13930Pa)，同时在油罐内侧前部产生真空(约为2577Pa)，如图14所示；罐体所受侧翻力矩在0.36s左右最大(约为39713N.m)；在0～0.12s左右，侧翻力矩为负(即无侧翻倾向)，侧翻力矩变化情况如图15所示；罐体所受最大应力约为191MPa，应力变化情况如图16所示。

(4)80％充液系数工况分析

从运加油车开始转弯至0.9s左右，油料由下向上向转弯外侧涌动，然后又从上往下涌动至1.5s左右,然后油料又由下向上向转弯外侧涌动，罐内油料流态如图17所示。

油料对罐体的冲击力如图18所示，油料对罐体在Z方向的冲击力较小，对X、Y方向的冲击力较大，在0.24s左右对罐体X方向充击力最大(约为48925N)，在0.12s左右对罐体Y反方向冲击力最大(约为51045N)。罐体在0.3s左右在罐体转弯外侧与底部交界区域压强最大(约为14210Pa)，同时在油罐内侧前部产生真空(约为2409Pa)，如图19所示；罐体所受侧翻力矩在0.27s左右最大(约为45487N.m)，在0～0.1s左右,侧翻力矩为负(即无侧翻倾向)，侧翻力矩变化情况如图20所示。罐体所受最大应力约为297MPa，应力变化情况如图21所示。(5)94％充液系数工况分析

从运加油车开始转弯至1.26s左右，油料由下向上向转弯外侧涌动，然后就趋于稳定，小幅波动，罐内油料流态如图22所示，油料对罐体的冲击力如图23所示；油料对罐体在Z方向的冲击力较小，对X、Y方向的冲击力较大，在0.12s左右对罐体X方向冲击力最大(约为47706N)，在0.24s左右对罐体Y反方向冲击力最大(约为57364N)。

罐体在0.3s左右在罐体转弯外侧与底部交界区域压强最大(约为16420Pa)，同时在油罐内侧上部产生真空(约为1653Pa)，如图24所示。罐体所受侧翻力矩在0.18s左右最大(约为36270N.m)，侧翻力矩变化情况如图25所示；罐体所受最大应力约为244MPa，应力变化情况如图26所示。

(6)结果分析

充液系数在80﹪时侧翻危险性最大，其最大侧翻力矩约为45487N.m，其次是充液系数70﹪时最大侧翻力矩约为39713N.m，然后是50﹪充液系数时，最大侧翻力矩为37621N.m，60﹪充液系数时，最大侧翻力矩为36929N.m，如图27所示，94﹪充液系数时，最大侧翻力矩为36270N.m，充液系数94﹪(即额定容量)时，侧翻稳定性较好。每种充液系数的最大侧翻力矩都是在开始转弯初期，然后逐渐降低。

最大应力出现在充液系数80﹪时(约为297MPa)，属于安全状态，每种充液系数的最大应力均出现在转弯初期，如图28所示。

Claims (2)

1.车载金属油罐内部油料晃动数值模型的建立方法，包括以下步骤：

(1)建立油料晃动基本方程

建立固定于大地上的固定坐标系x₀,y₀,z₀，和以油罐中心为坐标原点与油罐一起运动的运动坐标系x,y,z；

建立随罐体运动的运动坐标系x,y,z下罐内液体的运动控制方程为：

连续性方程

动量方程

流体体积分数输运方程

式中ρ为流体密度，为固定坐标系x₀,y₀,z₀下流体的绝对速度；为单位质量上的流体质量力函数；μ为流体动力学粘性系数；表示在运动坐标系下的流体速度矢量；C为流体体积分数；为矢量算子，为单位质量流体表面力，为运动坐标系下的流体速度变化矢量，为单位质量的惯性力：

为运动坐标系的旋转角速度矢量；为运动坐标系原点到所研究的液体微团k的矢径，为运动坐标系相对固定坐标系的平动加速度；为运动坐标系相对固定坐标系的转动加速度，即

(2)确定油料晃动的边界条件

油料边界包括两个部分，一部分是油料与油罐交界的固壁Ω_g，另一部分是油料与空气交界的自由表面Ω_f；

其中在油罐壁上油料的流体速度满足：

式中，为罐体在运动坐标系下相对于运动坐标系坐标原点运动的速度，忽略油罐罐壁的变形量，则

油料与空气交界的自由表面Ω_f的表面应力条件为：

在Ω_f上

<mrow> <mi>p</mi> <mo>-</mo> <mn>2</mn> <mi>&amp;mu;</mi> <mfrac> <mrow> <mo>&amp;part;</mo> <msub> <mi>v</mi> <mi>n</mi> </msub> </mrow> <mrow> <mo>&amp;part;</mo> <mi>n</mi> </mrow> </mfrac> <mo>=</mo> <msub> <mi>p</mi> <mn>0</mn> </msub> </mrow>

式中，p₀为液体表面大气压力；v_n为垂直于液体自由表面的法向速度，以流体域外法线方向为正；v_τ为切向速度；p为流体表面力；μ为流体动力学粘性系数；τ为剪应力；n为垂直于液面的单位矢量；

(3)固体结构区域受力

根据油料晃动基本方程和油料晃动的边界条件求解出罐内液体晃动冲击力后，把它加载到油罐罐壁上进行结构强度的计算，通过平衡微分方程、几何方程及本构方程在固体网格区域进行有限元法求解，得到油罐罐体的应力、应变、位移参数。

2.根据权利要去1所述车载金属油罐内部油料晃动数值模型的建立方法，其特征在于：所述平衡微分方程、几何方程及本构方程分别为：

平衡微分方程：

<mfenced open = "" close = "}"> <mtable> <mtr> <mtd> <mrow> <mfrac> <mrow> <mo>&amp;part;</mo> <msub> <mi>&amp;sigma;</mi> <mi>x</mi> </msub> </mrow> <mrow> <mo>&amp;part;</mo> <mi>x</mi> </mrow> </mfrac> <mo>+</mo> <mfrac> <mrow> <mo>&amp;part;</mo> <msub> <mi>&amp;tau;</mi> <mrow> <mi>x</mi> <mi>y</mi> </mrow> </msub> </mrow> <mrow> <mo>&amp;part;</mo> <mi>y</mi> </mrow> </mfrac> <mo>+</mo> <mfrac> <mrow> <mo>&amp;part;</mo> <msub> <mi>&amp;tau;</mi> <mrow> <mi>x</mi> <mi>z</mi> </mrow> </msub> </mrow> <mrow> <mo>&amp;part;</mo> <mi>z</mi> </mrow> </mfrac> <mo>+</mo> <msub> <mi>f</mi> <mi>x</mi> </msub> <mo>=</mo> <mn>0</mn> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mfrac> <mrow> <mo>&amp;part;</mo> <msub> <mi>&amp;tau;</mi> <mrow> <mi>y</mi> <mi>x</mi> </mrow> </msub> </mrow> <mrow> <mo>&amp;part;</mo> <mi>x</mi> </mrow> </mfrac> <mo>+</mo> <mfrac> <mrow> <mo>&amp;part;</mo> <msub> <mi>&amp;sigma;</mi> <mi>y</mi> </msub> </mrow> <mrow> <mo>&amp;part;</mo> <mi>y</mi> </mrow> </mfrac> <mo>+</mo> <mfrac> <mrow> <mo>&amp;part;</mo> <msub> <mi>&amp;tau;</mi> <mrow> <mi>y</mi> <mi>z</mi> </mrow> </msub> </mrow> <mrow> <mo>&amp;part;</mo> <mi>z</mi> </mrow> </mfrac> <mo>+</mo> <msub> <mi>f</mi> <mi>y</mi> </msub> <mo>=</mo> <mn>0</mn> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mfrac> <mrow> <mo>&amp;part;</mo> <msub> <mi>&amp;tau;</mi> <mrow> <mi>z</mi> <mi>x</mi> </mrow> </msub> </mrow> <mrow> <mo>&amp;part;</mo> <mi>x</mi> </mrow> </mfrac> <mo>+</mo> <mfrac> <mrow> <mo>&amp;part;</mo> <msub> <mi>&amp;tau;</mi> <mrow> <mi>z</mi> <mi>y</mi> </mrow> </msub> </mrow> <mrow> <mo>&amp;part;</mo> <mi>y</mi> </mrow> </mfrac> <mo>+</mo> <mfrac> <mrow> <mo>&amp;part;</mo> <msub> <mi>&amp;sigma;</mi> <mi>z</mi> </msub> </mrow> <mrow> <mo>&amp;part;</mo> <mi>z</mi> </mrow> </mfrac> <mo>+</mo> <msub> <mi>f</mi> <mi>z</mi> </msub> <mo>=</mo> <mn>0</mn> </mrow> </mtd> </mtr> </mtable> </mfenced>

式中，σ_x、σ_y、σ_z分别表示x、y、z方向的正应力分量，τ_xy、τ_yz、τ_zx分别表示x、y、z方向的剪应力分量，f_x、f_y、f_z分别表示x、y、z方向的体积力分量，τ_xz、τ_yx、τ_zy分别表示x、y、z反方向的剪应力分量；

几何方程：

<mfenced open = "" close = "}"> <mtable> <mtr> <mtd> <mrow> <msub> <mi>&amp;epsiv;</mi> <mi>x</mi> </msub> <mo>=</mo> <mfrac> <mrow> <mo>&amp;part;</mo> <mi>u</mi> </mrow> <mrow> <mo>&amp;part;</mo> <mi>x</mi> </mrow> </mfrac> <mo>,</mo> <msub> <mi>&amp;gamma;</mi> <mrow> <mi>x</mi> <mi>y</mi> </mrow> </msub> <mo>=</mo> <mfrac> <mrow> <mo>&amp;part;</mo> <mi>v</mi> </mrow> <mrow> <mo>&amp;part;</mo> <mi>x</mi> </mrow> </mfrac> <mo>+</mo> <mfrac> <mrow> <mo>&amp;part;</mo> <mi>u</mi> </mrow> <mrow> <mo>&amp;part;</mo> <mi>y</mi> </mrow> </mfrac> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>&amp;epsiv;</mi> <mi>y</mi> </msub> <mo>=</mo> <mfrac> <mrow> <mo>&amp;part;</mo> <mi>v</mi> </mrow> <mrow> <mo>&amp;part;</mo> <mi>y</mi> </mrow> </mfrac> <mo>,</mo> <msub> <mi>&amp;gamma;</mi> <mrow> <mi>y</mi> <mi>z</mi> </mrow> </msub> <mo>=</mo> <mfrac> <mrow> <mo>&amp;part;</mo> <mi>w</mi> </mrow> <mrow> <mo>&amp;part;</mo> <mi>y</mi> </mrow> </mfrac> <mo>+</mo> <mfrac> <mrow> <mo>&amp;part;</mo> <mi>v</mi> </mrow> <mrow> <mo>&amp;part;</mo> <mi>z</mi> </mrow> </mfrac> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>&amp;epsiv;</mi> <mi>z</mi> </msub> <mo>=</mo> <mfrac> <mrow> <mo>&amp;part;</mo> <mi>w</mi> </mrow> <mrow> <mo>&amp;part;</mo> <mi>z</mi> </mrow> </mfrac> <mo>,</mo> <msub> <mi>&amp;gamma;</mi> <mrow> <mi>z</mi> <mi>x</mi> </mrow> </msub> <mo>=</mo> <mfrac> <mrow> <mo>&amp;part;</mo> <mi>u</mi> </mrow> <mrow> <mo>&amp;part;</mo> <mi>z</mi> </mrow> </mfrac> <mo>+</mo> <mfrac> <mrow> <mo>&amp;part;</mo> <mi>w</mi> </mrow> <mrow> <mo>&amp;part;</mo> <mi>x</mi> </mrow> </mfrac> </mrow> </mtd> </mtr> </mtable> </mfenced>

式中，ε_x、ε_y、ε_z分别表示x、y、z方向的正应变，γ_xy、γ_yz、γ_zx分别表示x、y、z方向的剪应变，u、v、w分别表示x、y、z方向的位移分量；

本构方程：

<mfenced open = "" close = "}"> <mtable> <mtr> <mtd> <mrow> <msub> <mi>&amp;epsiv;</mi> <mi>x</mi> </msub> <mo>=</mo> <mfrac> <mn>1</mn> <mi>E</mi> </mfrac> <mo>&amp;lsqb;</mo> <msub> <mi>&amp;sigma;</mi> <mi>x</mi> </msub> <mo>-</mo> <mi>&amp;upsi;</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;sigma;</mi> <mi>y</mi> </msub> <mo>+</mo> <msub> <mi>&amp;sigma;</mi> <mi>z</mi> </msub> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> <mo>,</mo> <msub> <mi>&amp;gamma;</mi> <mrow> <mi>x</mi> <mi>y</mi> </mrow> </msub> <mo>=</mo> <mfrac> <mn>1</mn> <mi>G</mi> </mfrac> <msub> <mi>&amp;tau;</mi> <mrow> <mi>x</mi> <mi>y</mi> </mrow> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>&amp;epsiv;</mi> <mi>y</mi> </msub> <mo>=</mo> <mfrac> <mn>1</mn> <mi>E</mi> </mfrac> <mo>&amp;lsqb;</mo> <msub> <mi>&amp;sigma;</mi> <mi>y</mi> </msub> <mo>-</mo> <mi>&amp;upsi;</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;sigma;</mi> <mi>x</mi> </msub> <mo>+</mo> <msub> <mi>&amp;sigma;</mi> <mi>z</mi> </msub> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> <mo>,</mo> <msub> <mi>&amp;gamma;</mi> <mrow> <mi>y</mi> <mi>z</mi> </mrow> </msub> <mo>=</mo> <mfrac> <mn>1</mn> <mi>G</mi> </mfrac> <msub> <mi>&amp;tau;</mi> <mrow> <mi>y</mi> <mi>z</mi> </mrow> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>&amp;epsiv;</mi> <mi>z</mi> </msub> <mo>=</mo> <mfrac> <mn>1</mn> <mi>E</mi> </mfrac> <mo>&amp;lsqb;</mo> <msub> <mi>&amp;sigma;</mi> <mi>z</mi> </msub> <mo>-</mo> <mi>&amp;upsi;</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;sigma;</mi> <mi>x</mi> </msub> <mo>+</mo> <msub> <mi>&amp;sigma;</mi> <mi>y</mi> </msub> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> <mo>,</mo> <msub> <mi>&amp;gamma;</mi> <mrow> <mi>y</mi> <mi>x</mi> </mrow> </msub> <mo>=</mo> <mfrac> <mn>1</mn> <mi>G</mi> </mfrac> <msub> <mi>&amp;tau;</mi> <mrow> <mi>z</mi> <mi>x</mi> </mrow> </msub> </mrow> </mtd> </mtr> </mtable> </mfenced>

式中，E表示弹性模量，υ表示泊松比，G表示剪切模量。