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

1.1　单元运动方程

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

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

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

式中,Ei、Ej、Gi、Gj、i、j分别为单元i、j的弹性模量、剪切模量和泊松比,mi、mj、Ri、Rj分别为单元i、j的质量和半径.求解离散单元与有限单元的接触力时，将有限单元视为半径无穷大的离散单元[13].

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

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

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

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

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

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