越野车辆通常在沙地类的松软路面上作业.轮胎与沙地间的相互作用带来的沙粒流动、轮胎空转和沉陷等是影响汽车沙地行驶行为的主要因素，也关系到车辆的作业能力、作业效率和能源消耗等.因此，深入研究轮胎-沙地间的相互作用，对越野汽车的设计、性能评价等具有重要指导意义.近年来，随着计算机技术的飞速发展，数值分析方法被广泛地应用于轮胎-沙地相互作用的仿真分析.离散元有限元耦合方法(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区域，而且这个区域随着轮胎的移动而移动，以减少参与接触判断的离散单元数目，提高计算效率.具体描述如下. (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区域的扩充量应满足： 式中，Dmax为最大离散单元直径，Vz为轮胎z方向速度. 将以上Slave区域自适应定义方法在自主开发软件CDFP[14]中实现，即可保证轮胎-沙地接触计算中，将表征沙地的Slave区域限定在以轮胎中心为参考点的局部区域，该区域随轮胎中心的移动(包括水平方向的行驶和竖直方向的下陷)而自适应更新，可以显著减少轮胎-沙地接触搜索时间. 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材料单元个数杨氏模量/MPa泊松比密度/(mg/mm3)轮胎3522.0 0.491.8土槽5507.5×1040.202.4沙地500007.5×1040.202.4硬路面1687.5×1040.202.4 表2　仿真计算参数Table 2　Simulation calculation parameters对象法向阻尼/(N·s/m)切向阻尼/(N·s/m)摩擦系数轮胎-沙地50450.45沙地-土槽50450.40 3.2　仿真流程 仿真过程如下：(1)散体离散单元在自重下达到平衡状态；(2)轮胎在自重与外载荷(1 295 N)作用下达到平衡状态；(3)轮胎质心加载角速度5 rad/s及相应水平速度，分析不同滑转率下轮胎的行驶行为.程序执行流程图如图6所示，其中t为计算时间，Tter为计算终止时间，Δt为显式算法的时间步长. 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 以上仿真分析结果充分说明了离散元与有限元耦合分析方法评价越野车辆沙地行驶性能的有效性. 文中在介绍离散元与有限元耦合方法的基础上，针对轮胎有限单元与沙地离散单元接触算法，提出了Slave区域自适应定义方法以大幅提升轮胎与沙地间全局搜索效率.然后，使用基于FEM/DEM方法的自主开发软件CDFP仿真分析了轮胎-沙地间的相互作用.光面轮胎的仿真结果与对应土槽实验的挂钩牵引力和轮胎周边沙粒运动趋势、人字形花纹轮胎和光面轮胎在挂钩牵引力和轮胎下陷量仿真结果的差异，充分说明了FEM/DEM方法研究此类问题的有效性. 尽管FEM/DEM方法描述轮胎-沙地间的相互作用的有效性已经被初步验证，但仍存在如下问题：首先，沙粒采用的是球形离散单元，这与实际沙粒形状存在差别，即如何通过调整离散单元间的细观参数实现对实际沙粒宏观性能的等效；其次，考虑到模型规模和计算时间，文中沙地模型中的离散单元尺寸大于实际沙粒尺寸，这种假设对沙地的宏观力学性能如何影响，以及如何考虑这种影响以服务于工程分析；第三，离散单元方法的计算效率较低，这也是该方法工程应用的一大瓶颈，基于GPU并行计算的开发有望为该方法的工程应用提供有力支撑.

文章来源：固体力学学报 网址: http://gtlxxb.400nongye.com/lunwen/itemid-27239.shtml

上一篇： 水下冲击波载荷作用下气背固支圆板动态毁伤实
下一篇： 农业基础科学论文_黄土高原藓结皮土壤呼吸速率对降雨量变化的响应