极限运行工况下液罐双向流固耦合分析

蒋德飞，张竹林，邹彦冉

山东交通学院 汽车工程学院，济南 250357

0 引言

液罐车辆在道路运输过程中由于行驶状态的改变，罐体受到惯性力的作用引起液体晃动，表现出罐体与充装介质复杂的耦合效应[1]。液体的晃动导致罐车行驶不平稳，影响罐车的操纵稳定性及行车安全，同时晃动的液体产生瞬态载荷易导致罐体应力集中，造成罐体及防波板损坏。利用双向流固耦合方法对罐体内液体的晃动进行仿真分析，研究罐体的结构受力对保证罐车的安全稳定运行具有理论价值和实际意义。

近年来，国内外相关学者针对罐体内液体晃动对罐体的作用力进行了大量研究。何华等[2]采用Fluent软件仿真证明增加防波板数量可减小液体晃动，防波板数量增多时液体晃动的峰值载荷明显降低，同时对罐车不同制动初速度下罐体的瞬态应力进行对比分析，证明峰值应力存在时间与制动初速度有关；陈志伟[3]以充装液体介质的罐体为基础，建立液体晃动模型，分析刚性壁面应力分布和罐体内充装系数的大小及其对防波板和封头作用力的影响；刘忠亮[4]提出可通过对罐体V型支座的结构和位置的优化设计，使罐体封头和筒体中的应力最小；张韬等[5]对罐体制动工况液体的晃动进行流固耦合分析，得出罐车在达到极限制动工况0.52 s后液体晃动对罐体的压力达到最大；蔡昌全等[6]建立液化气体运输半挂车的整体模型，进行应力分析，对应力集中位置进行应力分类并进行强度校核，以确保满足运输安全要求；王若平等[7]按照文献[8]中的相关规定分别对纵向工况、横向工况、垂直向上工况、垂直向下工况进行应力、应变分析，指出罐体结构设计中的薄弱环节，为铝合金罐车的设计制造提供了可靠依据。

上述研究大多集中于采用单向流固耦合方式将仿真分析过程中的数据从流体传递给固体，该方法不易得到耦合过程中罐体的应力变化及分布规律，同时也忽略了罐体的结构变形对流体的影响，仿真结果相比双向流固耦合法求解精度较差。双向流固耦合方法流体与结构进行双向数据交换，固体以边界位移为载荷传递给流体，流体域将计算的流场数据返回到固体边界上，该方法能够得到结构域在流场提供的边界条件下的结构形变、应力分布等。

相对于其他行驶工况，罐车在制动和转向时罐体受力较复杂。根据文献[9]的要求，本文基于罐车在附着系数最大的沥青混凝土路面上制动时沿行驶方向不发生抱死拖滑可达到的最大制动减速度，以及转向时罐车沿行驶方向不发生拖滑的最大减速度和不发生侧翻的最大侧向减速度，采用双向流固耦合方法系统分析罐体在制动和转向2种极限工况下，不同充液比下的罐体所受应力随时间的变化趋势以及防波板的应力分布情况，对罐体的强度进行检算。

1 计算模型与计算方法

1.1 罐体模型和网格划分

某罐车罐体总长度为10 350 mm，宽度为2410 mm，罐体内部包含有6个防波板，将罐体内腔室分隔为7个舱室，相邻2个防波板上的过人孔交错布置，罐体结构如图1所示(图中单位mm)。

本文侧重于研究极限运行工况下罐内液体晃动产生的动载荷在液罐车罐体上的应力分布，采用三维软件SolidWorks构建罐体物理模型，将物理模型导入ANSYS-Workbench中，建立罐体有限元模型，采用ANSYS-Workbench中基于直接建模思想的重要几何前处理软件ANSYS-SpaceClaim对罐体在仿真过程中可能存在的尖角、小边及缝隙等进行检查和修补[10]，得到较精确的罐体有限元模型如图2所示。

a)罐体结构简图 b)防波板结构示意图 图1 罐体结构尺寸 图2 罐体有限元模型

根据图2的罐体有限元模型，通过ANSYS-SpaceClaim填充得到罐体的结构域与流体域模型[11]。通过ANSYS-Mechanical模块(该模块是ANSYS重要的前处理软件之一，主要功能包括网格智能生成和各种结果的数据处理)对罐体结构域和流体域模型进行网格划分，需要保证结构域和流体域网格尺寸与形状的一致性，网格过渡部分保证平滑且节点相互对应，可减少网格总数，提高耦合数据交互的便利性，避免在软件分析中产生模型的网格不匹配，影响仿真结果，提高计算效率。网格划分后，罐体结构域有限元模型包括277 264个单元，553 300个节点，罐体流体域模型包括352 429个单元，66 644个节点，结构域和流体域有限元模型如图3、4所示。

图3 罐体结构域网格模型 图4 罐体流体域网格模型

1.2 双向流固耦合

流体域模型可用于模拟罐体内部一定量的空气和水，根据相关约束和工况求解罐体结构域与流体域网格模型，流体域与结构域的数据相互交换，且在各个时间步都相互作用。采用双向流固耦合方法能够得到结构域在流场提供的边界条件下的时间历程，如结构形变、应力等变量的时间历程[12]。

双向流固耦合方法仿真精度相比单向流固耦合更符合实际运行工况，但对罐体的结构域和流体域模型进行仿真分析计算量较大，双向耦合分析具体流程如图5所示。进行瞬态计算时，罐体的结构域模型在ANSYS-Transient Structural中进行计算得到罐体的位移变形量；罐体的流体域模型在ANSYS-Fluent中仿真计算得到罐体在液体晃动下的压力、速度分布；结构域和流体域通过ANSYS-System Coupling完成耦合数据的交互传递。在耦合分析过程中，将总时间划分为n个时间步，在每个时间步内都要进行流体域与固体域的求解和数据传递，同时该时间步的结果作为下一时间步的初值进行计算，继而完成所有时间步求解[13-16]。

图5 双向流固耦合分析流程

1.3 两相流计算模型

本文对于流体域中流体运动采用ANSYS-Fluent中的流体体积模型(volume of fluid model,VOF)模拟，用来处理无相互穿插的多相流问题。

采用体积率函数表示流体自由面的位置和流体所占的体积，用每个单元内含有流体的体积和单元体积的比β表示流体的填充程度，当β=0时，该单元不含有流体；当β=1时，该单元充满流体；当0<β<1时，该单元含有自由液面。

VOF模型通过求解运动方程组和处理穿过区域的流体可以精确模拟气液分界面的变化情况和液体晃动过程[17-18]。多相流VOF模型的基本控制方程包括流体连续性方程与动量守恒方程。

1.3.1 流体连续性方程

式中：u、v、w为流体质点在x、y、z轴方向上的速度分量，ρ为流体质点处的密度，t为时间。

1.3.2 动量守恒方程

式中：Fx、Fy、Fz分别为流体质点在x、y、z方向上力的分量，p为流体质点处的压力。

1.4 求解器参数设置

罐体的流体域模型在ANSYS-Fluent中采用压力基(pressure-based)进行求解，鉴于罐体中气液两相的交界面随着运算时间发生变化，因此采用瞬态求解。在ANSYS-Fluent中湍流模型采用标准k-epsilon数值模型[19]。

流体域求解采用基于非迭代算法(pressure-implicit with splitting of operators，PISO)，该算法以压力基为理论基础，为瞬态不可压缩流动设计的一种求解方法，这种算法在求解瞬态问题时有明显优势，瞬态计算时，PISO算法能够在大时间步、动量及压力松弛因子(取值范围为0～1，通过改变压力松弛因子可控制收敛的速度和改善收敛情况)为1时仍能保持计算稳定，提高收敛速度[20-21]。

结构域与流体域有限元分析模型在仿真分析时需要定义材料属性。结构域材料采用铝合金，该材料屈服极限为325 MPa；为了安全起见，在对罐体进行应力校验时用屈服应力除以安全系数得到许用应力，根据文献[22-23]的规定，设计安全系数取1.5，计算得到罐体材料的许用应力为217 MPa。流体域中的气液两相模型通过ANSYS-Fluent中的Patch功能将流体域划分为气、液2部分，本文首先研究罐体充液比(罐体充装液体的体积与罐体容积的比值)为50%时的液体晃动，即50%的气体沿y轴正向分布，50%的液体沿y轴负向分布，气体和液体材料分别设置为空气和水；通过ANSYS-Fluent中Cell Registers建立不同充液比下的流体模型，从而实现多种充液比下罐体受力的仿真分析。

为防止流体域发生液面晃动后可能产生的网格变形、歪曲率过大或产生负体积网格[24]，本文动网格采取弹簧光顺网格控制方法，其计算控制方程为：

式中：Fi为网格节点i所受的力，Δxi、Δxj分别为网格节点i、j的位移，ni为与i相连的网格节点数，kij为i与j之间的弹簧刚度。

1.5 边界条件

双向流固耦合分析过程中罐体结构域内表面和流体域外表面需要进行耦合数据交互，在ANSYS-Transient Structural中可完成流体域耦合面的设置；结构域内表面命名为Solid_Move，流体域外表面命名为Fluid_Move，方便在ANSYS-System Coupling模块中创建数据交互通道，完成仿真过程中的双向数据交换。

设定罐体z方向和x方向对应的减速度用来模拟罐体制动工况和转向工况；流体域在Fluent中的边界条件设置为：给定流体域y轴负方向的重力加速度g=9.81 m/s2，不考虑传热影响，设置压力参考点，定义耦合面为流体区域的外表面，打开变形的相邻边界层，建立初始液体区域和气体区域，在液面处设置101 kPa，进行全局初始化。为了提高罐体双向流固耦合仿真精度，求解时间步长设置为0.005s，罐体结构域计算求解子步设置为3，总的模拟仿真时间为2 s，分析达到极限工况后2 s内流体域的运动参数和结构域的应力。

罐车行驶在不同路面时，在沥青或干混凝土路面上的附着系数最大，其滑动附着系数为0.75，制动减速度达到7.34 m/s2；峰值附着系数为0.8，制动减速度最高可达到7.8 m/s2，因此在模拟汽车制动工况时最大制动减速度为8.0 m/s2。

由于在良好路面上，车辆轮胎的附着系数最高可达0.8，即侧向减速度达到0.8g时，汽车开始侧滑；而罐车的侧翻阈值为0.4～0.6g，当侧向减速度达到0.4g，罐车尚未发生侧滑时，就已发生侧翻[25]，因此在模拟转向工况时，罐体z轴方向的最大减速度为8.0 m/s2，x轴方向的最大减速度为4.0 m/s2。

2 极限工况下罐体应力计算结果及分析

2.1 极限制动工况

在极限制动工况下对罐体进行瞬态分析，仿真分析完成后可在ANSYS-Fluent后处理软件CFD-Post中查看不同时刻的罐内气液2相分布图，观察仿真过程中的流体运动状况。气液2相图选取的平面为流体域沿y轴方向的纵切面；罐体内的液体晃动如图6所示，罐体后端的液体不断向罐体前端流动，同时空气通过防波板向罐体后端流动，最终罐体前端液体增多，空气与液体的交界面形成倾斜的近似斜面。

图6 极限制动工况下不同时刻罐体内的气液2相分布

罐体与防波板的应力分布分别如图7、8所示。由图7可知：罐体在晃动过程中最大应力为55.709 MPa，出现在罐体前封头与筒体接触的部位。由图8可以看出：防波板每处过人孔所受的应力沿汽车行驶方向逐渐增大，最大应力为41.596 MPa，出现在罐体沿行驶方向最前方的防波板直径最大的通液口处。

图7 极限制动工况下罐体应力云图 图8 极限制动工况下防波板应力云图

2.2 极限转向工况

在极限转向工况下对罐体进行瞬态分析，不同时刻罐体内气液2相分布情况如图9所示，以0.4 s为时间间隔分别列出罐车达极限转向工况后0、0.4、0.8、1.2、1.6、2.0 s时罐体的气液2相图，转向时罐体具有侧向加速度，气液2相图采用立体模型显示。从图9中可以看出：极限转向工况初始时刻液体域平面呈水平状，随着时间延长，罐体后端舱室的液体涌入前端舱室的同时，也有沿x方向的侧向运动。

图9 极限转向工况下不同时刻罐体内的气液2相分布

与极限制动工况相比，封头、防波板和筒体的峰值应力都明显增大，罐体最大应力为126.620 MPa，如图10所示，约为制动工况的2倍，防波板最大应力为89.073 MPa，如图11所示。与制动工况相比，罐体与防波板的承载应力区域位置发生了变化，罐体最大应力集中在前封头与筒体接触的上侧区域，防波板的最大应力发生在沿侧向加速度方向的防波板上部区域。

图10 极限转向工况下罐体应力云图 图11 极限转向工况下防波板应力云图

2.3 不同充液比下的仿真数值模拟

2.3.1 极限制动工况

为了分析同一运行工况下罐体充液比不同时罐车罐体的受力情况，选取充液比分别为30%、50%、80%，研究罐体及防波板的受力随时间的变化情况。极限制动工况下不同充液比时罐体及防波板的最大应力随时间的变化曲线如图12、13所示。

图12 极限制动工况下罐体最大应力时间历程 图13 极限制动工况下防波板最大应力时间历程

由图12、13可知：1)各防波板由于受到液体载荷冲击，应力均随时间呈波动性变化，其应力时间历程与罐体相似。2)随着充液比的增大，罐体和防波板所受应力逐渐增大，且应力增长速度加快，当充液比为80%时，罐体和防波板所受的应力最大，分别为94.433、70.846 MPa；随着充液比的增大，最大应力出现的时间延后，峰值应力出现在极限制动工况后0.07～0.12 s；随着时间继续增长，应力幅值逐渐减小。

2.3.2 极限转向工况

极限转向工况下，充液比为30%、50%、80%时罐体和防波板的最大应力变化趋势与制动工况相似，如图14、15所示。由图14、15可以看出：当充液比为80%时，罐体与防波板所受的应力最大，分别为223.12、148.76 MPa,出现在极限转向工况后0.05～0.10 s。

图14 极限转向工况下罐体最大应力时间历程 图15 极限转向工况下防波板最大应力时间历程

由极限工况下罐体应力计算结果分析可知：罐体与防波板的最大应力均发生在充液比为80%的极限转向工况，分别为223.12、148.76 MPa，分别位于罐体筒体与前封头的连接处与防波板沿侧向加速度方向的上部区域。

由于存在边缘应力的影响，罐体局部应力较大，属于一次应力中的局部薄膜应力，根据文献[22]计算局部薄膜应力强度的许用极限应力为325 MPa，因此该工况下罐体与防波板均满足强度要求。

3 结论

1)本文以某型号液体运输车辆罐体为研究对象，建立罐体结构域与流体域模型，采用双向流固耦合方法利用有限元软件ANSYS对罐车在制动和转向2种极限运行工况下罐体内液体的冲击晃动进行仿真分析，得到充液比为50%时2种极限工况下罐体和防波板的应力时间历程和气液2相分布图，从而得到2种工况下罐体与防波板的应力。

2)分析制动、转向2种极限运行工况充液比分别为30%、50%、80%时罐体与防波板的受力情况，得到峰值应力分别出现在极限制动工况后0.07～0.12 s、极限转向工况后0.05～0.10 s。罐体与防波板的最大应力均发生在充液比为80%的极限转向工况，分别为223.12、148.7 MPa。罐体的最大应力部位位于罐体筒体与前封头的连接处，防波板最大应力集中在防波板沿侧向加速度方向的上部区域，2处应力均小于许用极限应力325 MPa，故罐体与防波板均满足强度要求。