不同流固耦合方法进行轮胎滑水性能仿真研究的对比分析

2020-07-20刘从臻

轮胎工业 2020年6期

苑阳，刘从臻*，钱浩

（1.山东理工大学交通与车辆工程学院，山东淄博 255049；2.江苏大学汽车与交通工程学院，江苏镇江 212013）

流固耦合方法涉及到流体与固体间的相互影响，其对流体域和固体域分别求解的特性可更加真实地模拟滑水现象的规律。近年来其在模拟轮胎滑水及深入解析滑水性能方面得到广泛应用。董斌等[1-3]使用混合花纹轮胎仿真模型，采用大型有限元仿真软件Fluent计算获得了不同工况下胎面所受的水流动压力及胎面接触区域的水流速度分布。王国林等[4-6]采用计算流体力学（CFD）的气-液二相流数值模型将轮胎滑水过程简化为流体动压力作用使轮胎产生附加变形，并将临界滑水速度仿真预测值与已有的经验公式计算值进行对比分析。

从计算原理上看，单向耦合算法仅建立在水流压力场对轮胎结构场影响的情况下，与实际滑水现象有差异，从本质上存在误差。

随着多场耦合技术和计算机软、硬件技术的发展，双向流固耦合算法的应用逐渐受到各行专业人士和学者的青睐。在轮胎滑水双向流固耦合计算方面，臧孟炎等[7]建立了基于Ls-Dyna软件的轮胎滑水模型，通过分析轮胎接地前缘形成的楔状水层逐渐被转动的车轮挤压进胎面排水沟槽的状态及其沿纵向排水沟槽喷出的状态等研究不同水膜厚度下车辆临界滑水速度的变化。H.Grogger等[8]以光面轮胎和纵向花纹轮胎为研究对象，在考虑结构、大变形和与路面相互作用等因素的基础上，提出一种基于弱耦合流体-结构相互作用的集成软件工具，建立了有效且准确传递基本变量的耦合模块和处理流体域网格以匹配轮胎外表面的变形模块，将流体压力分布与轮胎仿真模型相结合，重复迭代运算，获得准确的流体压力分布和上浮力。

目前部分仿真分析方法未考虑轮胎与水层流固耦合的双向作用，而且很少通过流体域分析轮胎的滑水特性。因此，通过单向、双向流固耦合方法进行仿真对比，并基于流体域分析轮胎的滑水特性，从而明晰单向耦合算法的不足是很有必要的。

本研究以某厂家提供的185/60R15轮胎为例，基于Workbench平台的轮胎滑水流固耦合仿真分析方法，建立精确的轮胎滑水水流模型，在相同工况下进行单向、双向流固耦合算法的仿真对比。

1 理论背景

1.1 流固耦合研究简介

随着计算科学以及数值分析方法的不断发展，流固耦合或流固交互作用研究从20世纪80年代以来受到了广泛关注，从数据传递角度出发，其可以分为单向流固耦合和双向流固耦合[9]。一般来说，对于大多数耦合作用现象，如果只考虑静态结构性能，单向流固耦合分析已经足够，但是如果要考虑结构的动力学特性，双向流固耦合分析就必不可少。双向流固耦合同时考虑流体与固体，计算精度高，但是设置复杂，计算周期长。相比之下，单向流固耦合计算所需时间短、资源较少，仅适用于固体结构形变量不足以显著影响流体域的情况。单向流固耦合方法将流体域计算所得结果作为负荷施加到结构域中，而没有将结构域分析结果传递给流体域，设置相对简单，计算周期较短。

1.2 水流控制方程

本研究轮胎滑水流固耦合模型中的水流流动满足质量守恒、能量守恒和动量守恒定律，这3个物理量守恒构成了水流的基本控制方程[10]。

质量守恒方程：

能量守恒方程：

动量守恒方程：

式中ρw——流体密度；

t——时间；

ui，uj——速度张量；

xi，xj——坐标张量；

E——动能；

σw，εw——流体柯西应力、应变；

p——流体静压；

τij——应力张量；

ρwgi——重力体积力；

Fi——其他体积力。

1.3 湍流模型

适合湍动能（k）-耗散率（ε）模型的流动类型比较广泛，包括旋均匀剪切流、自由流（射流和混合层）、腔道流动和边界层流动。该模型k和ε的方程[10]为

其中C1＝max[0.43，η/（η＋5）]，η＝Sk/ε

式中μ——动力粘性系数；

μt——湍流粘性系数；

Gk——由于平均速度梯度引起的湍动能；

Gb——由于浮力引起的湍动能；

YM——可压缩湍流脉动膨胀对总耗散率的影响；

η——用户修正系数；

S——长度；

C1ε，C2，C3ε——在Fluent软件中为默认值常数，C1ε＝1.44，C2＝1.90，C3ε＝0.09；

σk，σε——k和ε的湍流普朗特数，σk＝1，σε＝1.20。

1.4 滑水经验公式

一般而言，有限元模型建立及计算结果有效性最合适的验证方法是与试验结果进行对比，然而滑水试验对测试场地条件和设备要求很高，因此本仿真采用滑水经验公式对轮胎滑水模型的合理性进行验证。W.B.Horne等[11]通过分析大量轮胎滑水的试验结果，提出了著名的NASA滑水经验公式：

式中vp——临界滑水速度；

P——轮胎充气压力。

NASA滑水经验公式只适合于积水较薄且被水覆盖的路面状况，为此，D.Dunlap等[12]在公式（6）的基础上，提出了滑水速度与胎面花纹深度、胎面宽度和水膜厚度之间关系的滑水预测公式：

DT——胎面花纹深度，一般情况下取值范围为1.6～10.0 mm；

DW——水膜厚度；

T——胎面宽度。

本仿真分析中，P＝220 kPa，DT＝7 mm，DW＝10 mm，T＝185 mm，代入公式（7）中，计算pv′ 为102.59 km·h-1（28.50 m·s-1）。

2 滑水物理模型

以某厂家提供的185/60R15轮胎为例，通过参考文献[13]获得花纹沟排列方向及花纹类型等对滑水性能的影响，进而建立垂直花纹轮胎和路面三维壳体模型，并赋予轮胎和路面材料属性，定义轮胎-路面接触算法，进行网格划分及边界条件定义；在CFD模块中建立流体域模型并进行网格划分和边界条件定义[3]。

2.1 结构模型与边界条件设置

结构场计算的对象为轮胎和刚性路面。轮胎胎面与胎侧分开建模、划分网格，共产生64 783个单元和36 483个节点，如图1所示。

图1 轮胎网格模型

网格划分完成后对模型设置边界条件，边界条件包括结构负荷和约束，如图2所示。

图2 轮胎边界条件

2.2 流体域模型与边界条件设置

对分析水层建立六面体模型并与轮胎进行布尔减操作，剔除两者相交部分，而剩余部分即为流体域模型。然后进行网格划分和边界条件定义。在Fluent软件中对流体域采用补片协调方法的四面体网格划分，共产生196 350个四面体单元和40 333个节点，如图3所示。流体域模型的边界条件定义为速度入口、压力出口以及壁面3种形式，并且在速度入口施加Profile文件定义水的运动，如图4所示。

图3 流体域模型

图4 流体域边界条件

2.3 求解器设置

为了使单向、双向两种流固耦合算法的计算结果具有可比性，不仅两种算法的边界条件设置、网格划分完全一致，流体域、结构域对应的求解器设置也完全一致。其中，流体域求解器设置项目较多，主要设置项如下：基于压力传递求解，采用可实现的k-ε湍流模型、平顺可重塑的动网格和半隐式方法求解压力耦合方程组。

3 结果与分析

3.1 流体域

t＝1.5 s时的单向、双向流固耦合流体域的动水压力分布如图5所示。

图5 流体域动水压力云图

由图5可见，两种方法的压力峰值位置分布情况基本一致，但在单向流固耦合流体域中接地前端胎面最大动水压力作用区域比较集中，说明该处水流动能大部分转化为作用于胎面的压力能，容易对轮胎产生明显的径向抬升作用，而双向流固耦合模拟轮胎变形改变了周围水流的压力场，正、负压力均降低，使压力分布相对较为均匀，难以对胎面形成集中有效的径向抬升作用，从而降低了轮胎发生滑水现象的几率。

将流体域流固耦合面处的动水压力等于或大于轮胎所受负荷时定义为轮胎发生完全滑水的时刻。在相同工况下轮胎-干路面法向接触力及单向、双向流固耦合动水压力随时间的变化曲线如图6所示。

图6 在相同工况下轮胎-干路面法向接触力及单向、双向流固耦合动水压力随时间的变化曲线

由图6可见，单向流固耦合计算发生完全滑水时刻为t＝1.1 s，而双向流固耦合计算发生完全滑水时刻为t＝1.2 s。

从图6可直观地看出不同耦合方法模拟整个滑水过程流固耦合面处动水压力变化的异同，动水压力曲线变化趋势基本一致，动水压力变化大致分为3个阶段。

（1）充气加载阶段。该阶段流速设置为零。

（2）水流冲击阶段。待轮胎充气加载稳定后，水逐渐加速冲击轮胎，由于垂直花纹沟槽的排水能力有限，两种算法的胎面所受动水压力都呈现迅速增大的趋势，而双向流固耦合模拟轮胎变形改变了周围水流的压力场，正、负压力均降低，使相同时刻动水压力较小，因此该阶段单向流固耦合计算的动水压力先与轮胎所受法向负荷平衡，即单向流固耦合模拟轮胎滑水发生先于双向流固耦合模拟。

（3）完全滑水阶段。该阶段轮胎动水压力大于或等于所受法向负荷，轮胎处于完全滑水状态。

3.2 轮胎结构性能分析与对比

以接触面中心为参考点，t＝1.5 s时轮胎的变形云图如图7所示。

图7 轮胎变形云图

由图7可见，该时刻单向、双向流固耦合计算的轮胎最大变形分别为17.652和15.649 mm。轮胎最大变形发生在接地区域前后两侧，并且变形量分布具有对称性；由于轮胎接地前缘楔形区域的存在，轮胎与流体互相挤压，使该区域水流出现滞留，从而使胎面中间区域动水压力高于胎侧区域；相同时刻下，双向流固耦合计算的最大动水压力和最大变形量小于单向流固耦合计算结果。

为了方便对两种流固耦合方法的计算结果进行对比，生成相同工况下轮胎变形随时间的变化曲线，如图8所示。

由图8可见，轮胎发生变形也可大致分为3个阶段。