基于FEM/DEM的轮胎-沙地相互作用的仿真*

赵春来　臧孟炎

(华南理工大学 机械与汽车工程学院， 广东 广州 510640)

基于FEM/DEM的轮胎-沙地相互作用的仿真*

赵春来臧孟炎†

(华南理工大学 机械与汽车工程学院， 广东 广州 510640)

摘要：介绍了离散元与有限元耦合方法(FEM/DEM)的基本原理，提出了对应轮胎-沙地接触计算的Slave区域自适应定义方法.基于室内土槽实验，建立了轮胎-沙地相互作用的三维仿真模型(不考虑轮胎变形)，其中，沙地使用球形离散单元(DE)模拟，土槽、轮胎使用六面体有限单元(FE)描述.然后，采用基于FEM/DEM的自主开发软件CDFP对光面和人字形沟槽两种胎面结构轮胎在沙地的行驶行为进行仿真，通过光面轮胎仿真结果与相应实验结果的对比，验证了该方法的可行性和有效性.对比两种胎面结构轮胎的仿真结果(法向反作用力、挂钩牵引力等)，解释了越野轮胎开凿沟槽的必要性.

关键词：固体力学；有限元；离散元；行驶行为

越野车辆通常在沙地类的松软路面上作业.轮胎与沙地间的相互作用带来的沙粒流动、轮胎空转和沉陷等是影响汽车沙地行驶行为的主要因素，也关系到车辆的作业能力、作业效率和能源消耗等.因此，深入研究轮胎-沙地间的相互作用，对越野汽车的设计、性能评价等具有重要指导意义.近年来，随着计算机技术的飞速发展，数值分析方法被广泛地应用于轮胎-沙地相互作用的仿真分析.离散元有限元耦合方法(FEM/DEM)是该领域研究中比较有前景的方法.目前，该方法在研究机械与土壤相互作用方面已经有很多应用，例如工程机械与土壤(David[1])相互作用、岩石切削(Oate[2])、土壤耕作(Takashi[3])、鱼雷与海床(Rosenil[4])相互作用等.Nakashima等[5-7]采用二维FEM/DEM研究了轮胎-沙地相互作用，采用离散元模拟沙地，有限元模拟轮胎，并开发了相应的计算程序.但二维方法在研究轮胎行驶行为时存在很多局限性，例如无法描述复杂的轮胎结构特性、沙粒的侧向流动及轮胎侧向力等.文中介绍了一种作者自主开发的、适用于越野车轮胎-沙地相互作用仿真分析的三维FEM/DEM方法，提出了对应轮胎-沙地间接触计算的Slave区域自适应定义方法，仿真分析了三维轮胎在沙地的行驶行为和不同花纹轮胎的特点.

1FEM/DEM算法介绍

1.1　单元运动方程

在FEM/DEM动力学算法中，离散单元和有限单元的运动方程为

mi(d2ui/dt2)=Fi

(1)

Ii(d2θi/dt2)=Mi

(2)

式(1)用于描述平动(适用于有限单元和离散单元)，式(2)用于描述转动(仅适用于离散单元);i为单元编号，mi为单元质量，Ii为转动惯量，ui为位移，θi为转角，Fi、Mi分别表示合外力和转矩.运动方程采用显式中心差分法求解.

1.2　单元间相互作用

接触算法是计算离散单元间、离散单元与有限单元间相互作用力的前提.离散单元间的接触判断采用基于格子搜索的“C-grid”方法[8].离散单元与六面体有限单元间的接触判断采用笔者所在课题组开发的“点-面”接触算法[9-10]，将离散单元视为接触从节点(Slave Node)，有限单元的接触面视为接触主面(Master Surface).其主要步骤为：

步骤1将接触对所在的三维空间分成合适大小的盒子区域并编号;

步骤2在事先定义的接触区域找出某个盒子包含的离散单元，并将盒子的编号赋给这些单元;

步骤3根据接触主面上的节点所在的盒子区域，给每一个节点赋一个盒子编号;

步骤4在离散单元所在盒子及与其相邻的26个盒子中，找出距离该单元最近的节点的编号，该最近节点周围的面即为潜在的接触面;

步骤5通过离散单元质心位置与潜在接触面间的相对位置，确定接触是否发生，并计算穿透量.

需要指出的是,显式算法的时间步长通常比较小，可以认为当前时步的“潜在接触对”相对于上一个时步没有变化[11].因此，上述搜索过程中的步骤2、3、4(全局搜索过程)不需要每个时步都执行，其搜索时间间隔ΔT与单元运动速度、搜索盒子尺寸和临界尺寸等有关，具体可参考文献[10].

单元间的接触作用模型如图1所示.其中，hij为单元间的重叠量，vi、vj、ωi、ωj分别为单元i和j的平动速度和角速度.Fn、Fs分别为法向力和切向力，μ为摩擦系数.

Fn=Fn,e+Fn,v

(3)

(4)

式中，Fn,e、Fn,v、Fs,e、Fs,v分别为单元间的法向弹簧力、法向阻尼力、切向弹簧力和切向阻尼力，通过Hertz-Mindlin理论求解[12],

(5)

Fn,v=-γnmijvn,ij

(6)

(7)

Fs,v=-γsmijvs,ij

(8)

式中，nij、sij分别为法向和切向单位向量,vn,ij、vs,ij分别为单元间法向和切向相对速度,Eij、Gij、mij、Rij分别为

(9)

(10)

mij=(mimj)/(mi+mj)

(11)

Rij=(RiRj)/(Ri+Rj)

(12)

图1　单元接触模型Fig.1　Contact models of elements

1.3　轮胎-沙地相互作用的FEM/DEM描述

用FEM/DEM分析轮胎-沙地相互作用的二维示意图如图2所示.假设轮胎沿x正方向行驶.则轮胎的挂钩牵引力N、路面竖直方向反作用力P及滑转率s分别为

N=G-|R|

(13)

P=∑fy

(14)

s=(1-V/(rω))

(15)

式中,G=∑fx+，R=∑fx-, f为轮胎有限单元与沙粒离散单元间的接触力,V为轮胎质心水平速度,s为轮胎的滑转率,r为轮胎的滚动半径,ω为轮胎的角速度.

图2　轮胎与沙地间的相互作用示意图Fig.2　Schematic diagram of tire-sand interactions

2Slave区域自适应定义方法

常用的接触定义方法对Slave和Master区域的定义在建模时一次完成，计算过程中不能发生变化.越野车辆沙地行驶性能评价需要车辆在不同条件沙地行驶长达80 m的距离，如果仿真分析采用现行的接触定义方法，所有离散单元沙地区域都需要定义为与轮胎接触的Slave区域(如图3所示)，轮胎-沙地的接触搜索计算将需要耗费大量时间.考虑到计算过程中轮胎在较大范围移动的特点，文中在离散元与有限元“点-面”接触算法的基础上，提出了Slave区域自适应定义法，即初始建模时定义一个较小的Slave区域，而且这个区域随着轮胎的移动而移动，以减少参与接触判断的离散单元数目，提高计算效率.具体描述如下.

(1)初始Slave区域确定如图3所示，以轮胎质心O为参考点，定义轮胎-沙地接触搜索的初始Slave区域.

图3　Slave区域Fig.3　Slave search region

图中r为轮胎半径，B为轮胎宽度，l1为Slave区域沿轮胎行驶方向的扩充量，l2为Slave区域沿轮胎行驶方向相反方向和沿轮胎两侧的扩充量，l3为Slave区域沿轮胎下陷量方向的扩充量.Slave区域大小为x：2r+l1+l2,y：2l2+B,z：r+l3.初始Slave区域确定后，开始如图4(a)所示的轮胎-沙地相互作用仿真计算.其中，t为计算时间，t0为初始时刻.

(2)Slave区域自动更新 假设ΔT为有限元轮胎和离散元沙地“点-面”接触算法中全局搜索时间间隔(详见1.2节).经过时间ΔT后，轮胎行驶至图4(b)所示的初始Slave区域左侧(轮胎质心移动距离为ΔT×V).在新的全局搜索之前，更新Slave区域：以新的轮胎质心坐标为参考点，建立新的“点-面”搜索Salve区域，如图4(c)所示.

图4　Slave区域自适应定义Fig.4　Concept of the adaptive Slave region

为确保与轮胎相接触的沙地在轮胎-沙地相互作用计算中能够有效地参与接触计算，Slave区域的扩充量应满足：

(16)

式中，Dmax为最大离散单元直径，Vz为轮胎z方向速度.

将以上Slave区域自适应定义方法在自主开发软件CDFP[14]中实现，即可保证轮胎-沙地接触计算中，将表征沙地的Slave区域限定在以轮胎中心为参考点的局部区域，该区域随轮胎中心的移动(包括水平方向的行驶和竖直方向的下陷)而自适应更新，可以显著减少轮胎-沙地接触搜索时间.

3数值算例

3.1　仿真模型

为确认三维FEM/DEM方法研究轮胎-沙地行驶行为的可行性，基于室内土槽实验[15]，建立了光面轮胎在沙地土槽中的行驶仿真模型，如图5所示.

图5　土槽实验三维仿真模型Fig.5　3D simulation model of the soil-bin experiment

测试实验中，轮胎通常由硬路驶入松软路面.因此，模型中路面由硬路和沙地两部分组成.硬路用六面体有限单元离散，选用弹性材料，不考虑其破坏；沙地采用土槽中随机分布的球形离散单元模拟，单元半径范围为5~7 mm；土槽和轮胎采用六面体有限单元描述，不考虑其变形.模型仿真参数如表1、2所示.

表1　材料参数Table 1　Material parameters

表2　仿真计算参数Table 2　Simulation calculation parameters

3.2　仿真流程

仿真过程如下：(1)散体离散单元在自重下达到平衡状态；(2)轮胎在自重与外载荷(1 295 N)作用下达到平衡状态；(3)轮胎质心加载角速度5 rad/s及相应水平速度，分析不同滑转率下轮胎的行驶行为.程序执行流程图如图6所示，其中t为计算时间，Tter为计算终止时间，Δt为显式算法的时间步长.

4仿真结果及讨论

4.1　光面轮胎仿真结果及分析

图7为光面轮胎在30%滑转率下的行驶轨迹云图，云图显示沙粒在竖直方向(z方向)的位移(为便于观察，隐去了硬路和土槽).由图可知：仿真分析可获得明显的车辙，轮胎下的离散单元被压实下沉、位移为负值，轮胎两侧的单元受轮胎挤压而“隆起”、位移为正值.

图6　程序流程图Fig.6　Program flowchart

图7　光面轮胎在30%滑转率工况下的车辙Fig.7　Rut of smooth tire under 30% slip ratio condition

改变轮胎行驶速度，分析光胎在不同滑转率下的行驶行为，并将挂钩牵引力随滑转率的变化趋势与文献[15]中的实验结果对比，如图8所示.

图8　挂钩牵引力随滑转率变化趋势Fig.8　Variation trend of drawbar pull with slip ratio

图8中挂钩牵引力值取轮胎在沙地进入稳定状态后的均值.由图可知，当滑转率小于25%时，尽管仿真结果与实验值有较大差别，但整体趋势一致.即在滑转率小于25%前，挂钩牵引力随着滑转率的增加而明显增大；大于25%后，挂钩牵引力几乎不随滑转率的增加而变化.

图9(a)为轮胎在30%滑转率工况下的沙粒流动速度矢量图，由图可知沙粒流动可分为两个区域，前区顺时针方向和后区逆时针方向.该现象与图9(b)所示的实验结果[16]一致.

图9　沙粒流动趋势Fig.9　Flow trend of the sand particles

4.2　不同胎面花纹仿真结果比较

为进一步确认离散元与有限元耦合方法在越野车辆沙地行驶性能仿真评价方面的有效性，文中在原有光面轮胎的基础上，分析了一款同规格、人字形沟槽花纹轮胎的沙地行驶性能，花纹沟槽数目为12，与轮胎纵向的夹角为60°，沟槽宽度为23.4 mm.图10为人字形沟槽花纹轮胎在30%滑转率下的行驶轨迹仿真结果.图11为两种胎面结构轮胎在相同滑转率下所受的路面法向反作用力时间历程.由图10、11可知：由硬路面进入沙地时，两种轮胎所受的路面法向反作用力均出现较大波动，随后在1295 N(轮胎自重及外载荷之和)附近波动.

图10　花纹轮胎在30%滑转率工况下的车辙Fig.10　Rut of groove tire under 30% slip ratio condition

图11　滑转率为30%时路面法向反作用力时间历程Fig.11　Time history of road vertical reaction force under 30% slip ratio

图12为两种轮胎的挂钩牵引力时间历程.由图可知，当轮胎由硬路进入沙地后，挂钩牵引力先出现较大波动，随后迅速减小，并进入相对稳定状态.沙地上花纹轮胎的挂钩牵引力均值约为155 N，大于光面轮胎的仿真结果110 N.这是因为，轮胎花纹沟槽增大了轮胎与沙地的接触面积，导致驱动力增加，也说明了越野车轮胎开凿花纹的必要性.

图12　滑转率为30%时挂钩牵引力时间历程Fig.12　Time history of drawbar pull under 30% slip ratio

图13为两种轮胎的下陷量时间历程.由图可知，当轮胎由硬路面进入沙地后，下陷量迅速增加，并很快趋于稳定状态.沙地上花纹轮胎的下陷量值约为28 mm，大于光面轮胎的24 mm.这一结果是由于花纹轮胎的切向作用力将沙粒“扒开”，导致其下陷量值较光面轮胎大.

图13　滑转率为30%时的轮胎下陷量时间历程Fig.13　Time history of tire sinkage under 30% slip ratio

以上仿真分析结果充分说明了离散元与有限元耦合分析方法评价越野车辆沙地行驶性能的有效性.

5结语

文中在介绍离散元与有限元耦合方法的基础上，针对轮胎有限单元与沙地离散单元接触算法，提出了Slave区域自适应定义方法以大幅提升轮胎与沙地间全局搜索效率.然后，使用基于FEM/DEM方法的自主开发软件CDFP仿真分析了轮胎-沙地间的相互作用.光面轮胎的仿真结果与对应土槽实验的挂钩牵引力和轮胎周边沙粒运动趋势、人字形花纹轮胎和光面轮胎在挂钩牵引力和轮胎下陷量仿真结果的差异，充分说明了FEM/DEM方法研究此类问题的有效性.

尽管FEM/DEM方法描述轮胎-沙地间的相互作用的有效性已经被初步验证，但仍存在如下问题：首先，沙粒采用的是球形离散单元，这与实际沙粒形状存在差别，即如何通过调整离散单元间的细观参数实现对实际沙粒宏观性能的等效；其次，考虑到模型规模和计算时间，文中沙地模型中的离散单元尺寸大于实际沙粒尺寸，这种假设对沙地的宏观力学性能如何影响，以及如何考虑这种影响以服务于工程分析；第三，离散单元方法的计算效率较低，这也是该方法工程应用的一大瓶颈，基于GPU并行计算的开发有望为该方法的工程应用提供有力支撑.

参考文献：

[1]Horner D,Peters J,Carrillo A.Large scale discrete element modeling of vehicle-soil interaction [J].Journal of Engineering Mechanics,2001,127(10):1027-1032.

[3]Okayasu T,Fukuda T,Tsuchiya K,et al.Development of soil tillage simulator using fem-dem quasi-coupling method [C]∥Proceedings of the 7th International Symposium on Machinery and Mechatronics for Agriculture and Biosystems Engineering(ISMAB).Taiwan：the Chinese Institute of Agricultural Machinery,2014:1-6.

[4]Mendes R B,Alves J L D,Silva C E.Simulation of torpedo pile launching by coupled discrete and finite element analisys [C]∥Procedings of the XXVII Iberian Latin American Congress on Computacional Methods in Enginee-ring (XXVII CILAMCE).Para-BRAZIL：the Brazilian Association of Computational Methods in Engineering,2006:1-19.

[5]Nakashima H,Takatsu Y,Shinone H.Analysis of tire tractive performance on deformable terrain by finite element-discrete element method [J].Journal of Computational Science and Technology,2008,4(2):423-434.

[6]Nakashima H,Takatsu Y,Shinone H.FE-DEM analysis of the effect of tread pattern on the tractive performance of tires operating on sand [J].Journal of Mechanical Systems for Transportation and Logistics,2009,2(1):55-65.

[7]Nakashima H,Oida A.Algorithm and implementation of soil-tire contact analysis code based on dynamic FE-DE method [J].Journal of Terramechanics,2004,41:127-137.

[8]Williams J,Perkins E,Cook B.A contact algorithm for partitioning N arbitrary sized objects [J].Engineering Computations,2004,21(2/3/4):235-248.

[9]高伟，臧孟炎.基于内聚力模型的夹层玻璃梁冲击破坏过程仿真 [J].华南理工大学学报：自然科学版，2013,41(5):139-144.

Gao Wei,Zang Meng-yan.Simulation of impact fracture process of laminated glass beam based on cohesive model [J].Journal of South China University of Technology：Natural Science Edition,2013,41(5):139-144.

[10]Zang M Y,Gao W,Lei Z.A contact algorithm for 3D discrete and finite element contact problems based on pe-nalty function method [J].Computational Mechanics,2011,48(5):541-550.

[11]Benson D J,Hallquist J O.A single surface contact algorithm for the post-buckling analysis of shell structures [J].Computer Methods in Applied Mechanics and Engineering,1990,78(2):141-163.

[13]Han K,Peric D,Crook A,et al.A combined finite/discrete element simulation of shot peening processes-part I：studies on 2D interaction laws [J].Engineering Computations,2000,17(5):593-620.

[14]臧孟炎，雷周，高伟.显式离散元与有限元耦合软件CDFP V1.0[CP/CD].著作权登记号：2011SR057163,2011-08-12.

[15]Shinone H,Nakashima H,Takatsu Y.Experimental analysis of tread pattern effects on tire tractive performance on sand using an indoor traction measurement system with forced-slip mechanism [J].Engineering in Agriculture,Environment and Food,2010,3(2):61-66.

[16]庄继德.计算汽车地面力学 [M].北京：机械工业出版社,2002:88-90.

Simulation of Tire-Sand Interactions Based on FEM/DEM

ZhaoChun-laiZangMeng-yan

(School of Mechanical and Automotive Engineering,South China University of Technology, Guangzhou 510640,Guangdong,China)

Abstract:In this paper, the basis of the combined finite and discrete element method(FEM/DEM) is introduced, and an adaptive Slave region method is proposed according to the FEM/DEM simulation of tire-sand interactions. Moreover, the three-dimension simulation models of tires running on sand are constructed based on a soil bin experiment with the tire deformation being neglected, where the sand is modeled by using spherical discrete elements, the soil bin and the tires are discretized into hexahedron finite elements. Then, the running behaviors of the smooth and groove tires on sand at different slip ratios are simulated by using the self-developed software CDFP, which is an explicit combined discrete-finite element program. The simulation results of the smooth tire are compared with the corresponding experimental ones, thus verifying the feasibility and effectiveness of the proposed method. In addition, the comparison of the simulation results(such as road reaction force and tire drawbar pull) of the smooth and groove tires explains the necessity of the grooves of off-road vehicle tires.

Key words:solid mechanics; finite element; discrete element;running behavior

中图分类号：U 461.5+1

doi:10.3969/j.issn.1000-565X.2015.08.012

文章编号：1000-565X(2015)08-0075-07

作者简介：赵春来(1988-)，男，博士生，主要从事汽车CAE技术研究.E-mail： chunlai0401@163.com†通信作者： 臧孟炎(1961-)，男，教授，博士生导师，主要从事汽车系统动力学、汽车CAE技术研究.E-mail： myzang@scut.edu.cn

*基金项目：国家自然科学基金资助项目(11172104)；广东省科技厅对外科技合作资助项目(2013B051000012)

收稿日期：2015-02-12

Foundation items: Supported by the National Natural Science Foundation of China(11172104) and the International S&T Cooperation Program of Guangdong Province(2013B051000012)