赵春来，臧孟炎
（华南理工大学机械与汽车工程学院 广州）
摘 要：采用离散元与有限元耦合方法（FEM/DEM）模拟轮胎沙石路面相互作用，提出 “路面更替法”，以控制仿真模型规模，提高计算效率。建立刚性有限元轮胎在离散元沙石路面上行驶的三维仿真模型，使用自主开发软件CDFP对刚性车轮砂石路面行驶行为进行仿真分析。获取了轮胎挂钩牵引力、下陷量、沙石颗粒流动趋势等行驶行为参数。通过与相关实验结果对比，验证了该方法的可行性。
关键词：越野车辆；通过性能；沙石路面；FEM/DEM；路面更替
中图分类号：U461.5+1 文献标志码：A
深入研究轮胎-沙粒间的相互作用对高水平越野车的设计、性能评价等具有重要的指导意义。然而，单纯使用有限元方法（FEM）虽然能够很好地描述轮胎特性[1-5]（包括轮胎的复杂花纹形状、充气轮胎的变形等），但在研究越野车辆沙地行驶行为评价方面遇到了挑战：由于沙石路面具有典型的非连续介质特性，用传统的连续介质理论（如：FEM）不能合理描述其运动行为。离散元方法（DEM）在描述散体介质的非连续特性方面具有明显优势[6-9]。
为了兼顾FEM和DEM的优点，Nakashima[10-12]等采用二维离散元与有限元耦合方法（FEM/DEM）研究轮胎在沙石路面上的行驶行为，并开发了相应的计算程序。但是，二维方法在研究轮胎行驶行为时存在很多局限性，例如：无法描述轮胎花纹、沙石颗粒的侧向流动及轮胎侧向力等。本文介绍了一种我们自主开发的、适用于越野车轮胎行驶行为仿真分析的三维FEM/DEM计算方法，然后针对越野车辆沙地行驶过程的特点，提出了“路面更替法”以控制仿真模型规模，并基于FORTRAN 95语言开发了相应的专用分析软件。
1 方法介绍
1.1 FEM/DEM算法介绍
本文所述FEM/DEM方法在车轮-沙地相互作用分析的本质是使用FEM来模拟轮胎，用DEM方法来描述沙粒。离散单元及与沙石接触轮胎有限单元节点的运动受牛顿第二定律控制。算法中包含图1（a）、（b）所示的两种接触类型：（1）离散单元间的接触；（2）离散单元与有限单元间的接触。这两种接触类型分别用来描述沙石颗粒间的相互作用和轮胎与沙石颗粒间的相互作用。其中前者采用Williams等提出的‘C-grid’方法[13]进行接触判断，后者采用作者开发的适用于球形离散单元和六面体有限单元的NTS算法[14]。
（a） 离散单元间相互作用 （b） 离散元与有限元相互作用 （c） 单元间作用力
图1 单元接触模型
图1中，hij为单元间的重叠量，vi、vj、ωi、ωj分别为单元i和j的速度和角速度。Fn和Fs分别为单元间的法向力和切向力，如式（1）、（2）所示，其中切向力Fs满足库伦摩擦定律。
（1）
（2）
式中，Fn,e、Fn,v、Fs,e、Fs,v分别为单元间的法向弹簧力、法向阻尼力、切向弹簧力和切向阻尼力，通过Hertz-Mindlin理论求解[15]。μ为单元间的摩擦系数。求解有限单元与离散单元间作用力时，将有限单元视为半径无穷大的离散单元[16]。
轮胎与沙石颗粒相互作用的二维示意图如图2所示，车轮用有限元描述，沙石路面用离散元描述。轮胎沿
正方向行驶。轮胎的挂钩牵引力N、轮胎所受的路面竖直方向反作用力P及轮胎滑转率s分别为
（3）
（4）
（5）
其中，G=Σfx+；R=Σfx-；f为轮胎与沙石颗粒间的相互作用力；V为轮胎质心水平速度； r为轮胎的滚动半径，对刚性轮胎，滚动半径等于轮胎自然半径；
为轮胎的角速度。
图2 轮胎与路面间的相互作用示意图
2.2 路面置换法
离散元方法在描述非连续的沙石路面时具有优势，但该法计算效率较低，而影响计算效率的主要因素是参与计算的离散单元的数量。实际当中，车轮后方一定范围外，被碾压过的沙石颗粒对车轮行驶行为的影响可以忽略，同时，车轮前方一定距范围外，沙石颗粒对车轮行驶行为影响也较小。基于此，我们提出一种“路面更替法”，以控制参与计算的离散单元数目，提高计算效率。下面通过二维示意图描述该方法的主要步骤：
（1）首先，在给定区域内生成离散单元集，作为初始沙石路面样本，记为DES A0，长度为a，如图3（a）所示；
（2）第二，将上一步中的离散单元集DES A0复制两份，分别记为DES A1和DES A2，并按顺序排列，组成长度为2a的初始路面，如图3（b）所示；
（3）第三，将车轮放置在路面DES A1中间位置，车轮在自重及外载荷作用下自由下沉，直到轮胎所受的路面法向反作用力等于车轮自重及外载荷之和。然后给车轮施加恒定的角速度和相应的水平速度，使车轮在某一滑转率下行驶，如图3（c）所示；
（4）第四，当车轮行驶到DES A2的中间位置时，DES A1对车轮行驶行为影响可以忽略，且DES A2前端没有发生明显破坏，此时开始进行第一次路面更替：将DES A1删除，同时重新备份路面样本，记为DES A3，放置在DES A2前端，与DES A2组成新的路面，如图3（d）所示。
反复执行步骤（4），参与计算的沙石路面长度始终等于2a，可以实现较少数量的离散单元来模拟任意长度的沙石路面，大幅提高计算效率的目的。对于由多种结构路面构成的越野车辆沙地通过性能评价，只需对每种结构路面生成样本路面，然后根据车轮即将行驶的路面结构调入对应结构样本路面即可。
（a） 沙石区域样本 （b） 组成初始沙石区域
（c） 放置轮胎 （d） 路面更替
图3 路面置换方法的二维示意图
同种路面结构相应的程序流程图如图4所示，其中Tter为计算终止时间，Δt为显式算法的时间步长，T为当前计算时间。
图4 执行流程图
3 数值算例
为验证方法的有效性，基于室内土槽实验[17]，建立了轮胎在沙石路面行驶的三维仿真模型，分析了不同滑转率下轮胎的行驶行为。
3.1 离散元沙石路面
首先，采用作者之前开发的分层生成方法[18]，在指定区域内，随机生成给定半径范围的球形离散单元集，区域边界采用刚性墙约束。然后，让初始生成的离散单元在自重作用下达到稳定状态，以模拟自然沙地。最终生成的离散元在X: （0,735）、Y: （0,480）、Z: （0,280）范围，离散单元集孔隙率约为0.32，以此作为沙石路面初始样本，如图5所示。本文侧重于对开发方法的验证，因此，离散单元半径范围取6~7mm，大于实际沙粒半径。自重压实过程离散单元总势能的时间历程如图6所示，在1.1s以后，其趋势基本平稳，此时可以认为压实过程结束。
图5 初始沙石路面样本 图6 离散元沙石颗粒重力势能时间历程
杨氏模量：75000 MPa，泊松比：0.3，密度：2400 kg/m3，单元数目：45551，摩擦系数：0.3
3.2 初始路面及车轮模型
初始沙石路面区域由两段经自重压实后的“路面样本”组成，单元总数为91102个，长度等于路面样本长度二倍的1470mm。有限元车轮采用如图7所示的简化模型，暂不考虑轮胎变形。轮胎与沙路间的摩擦系数为0.4，轮胎放置在如图所示位置。给轮胎施加自重及外载荷（计1295N）同时加载恒定的角速度0.5rad/s和相应的水平速度，使轮胎在给定滑转率下沿X轴正方向行驶。
（a） 主视图 （b） 左视图
图7 初始路面和有限元轮胎模型
轮胎参数：杨氏模量：2 MPa，泊松比：0.49，密度：1800 kg/m3，单元数目：1344
3.3 仿真结果
当车轮行驶距离为k×735 mm（k为路面更替次数）时执行路面更替。图8给出轮胎在30%滑转率下，离散单元Z方向的位移云图：图8（a）为初始状态；图8（b）为轮胎行驶到第二段样本路面中间位置时的车辙；图（c）为执行第一次沙石路面更替后的车辙；图8（d）给出了第二路面更替前的车辙状况；图8（e）执行了第二次路面更替；图8（f）为计算终止状态，车轮行驶距离为1560mm。
（a） （b）
（c） （d）
（e） （f）
图8 车轮滑转率为30%时沙石颗粒Z方向位移云图
图9为轮胎在滑转率为30%时的轮胎挂钩牵引力时间历程。挂钩牵引力在初始阶段有一个较大峰值，主要原因是初始阶段，轮胎与路面间竖直方向有相对冲击速度，产生一个冲击力，引起路面法向反作用力出现峰值，与此相应，挂钩牵引力初始阶段也出现峰值。随后挂钩牵引力快速减小，但仍有较大波动，主要原因是离散单元颗粒半径取值较大，对接触力计算结果影响所致。
图9 滑转率为30%时车轮行驶过程 图10 滑转率为30%时车轮行驶过程中的
中的轮胎挂钩牵引力时间历程 轮胎下陷量时间历程
图10为轮胎行驶过程中的下陷量时间历程，由图可知：在初始阶段，轮胎下陷量快速增加，主要原因是轮胎竖直方向载荷未达到平衡状态。在行驶200mm之后，下陷量基本上趋于一个稳定的状态，均值约为45mm，且在路面更替时未出现明显波动。
图11为轮胎在30%滑转率下行驶时，沙石颗粒的流动趋势，图中箭头表示离散单元在X-Z平面内的速度矢量。由图可知：沙石颗粒的流动趋势可分为两个区域：前区，绕顺时针方向；后区，绕逆时针方向。该现象与庄继德[19]等人的实验研究结果基本一致。
（a） 仿真结果 （b） 实验结果
图11 车轮行驶过程中沙石颗粒的运动趋势
滑转率是影响轮胎行驶行为的主要因素，为研究不同滑转率下轮胎的行驶行为，图12和图13分别给出了轮胎下陷量和挂钩牵引力随滑转率的变化关系，其中下陷量和挂钩牵引力值取各工况稳定状态下的均值。由图可知：轮胎的下陷量随滑转率的增加而逐渐增加，而且其增加趋势逐渐“陡峭”，该趋势与文献[17]研究结果基本一致。
挂钩牵引力随滑转率的增加逐渐增加，当滑转率小于25%时，增加趋势较陡峭，当滑转率大于25%时，增加趋势缓慢，该现象也与文献[17]实验结果整体趋势一致。
图12 轮胎下陷量随滑转率变化关系 图13 挂钩牵引力随滑转率变化关系
4 结 论
将FEM/DEM方法应用于车轮在沙石路面的行驶行为仿真分析，并提出了一种“路面更替法”。路面更替法可以实现一定数量的离散单元来模拟任意长度的沙石路面，达到控制模型规模，提高计算效率的目的。
轮胎行驶行为，包括路面竖直方向反作用力，挂钩牵引力，下陷量及轮胎行驶过程中沙石颗粒的流动状况等都可很方便地通过计算求得，数值试验结果与实验结果对比验证了该方法的有效性。
参考文献
[1] Xia K. Finite element modeling of tire/terrain interaction: Application to predicting soil compaction and tire mobility[J]. Journal of Terramechanics. 2011;48:113-23.
[2] Xia K, Yang YM. Three-dimensional finite element modeling of tire/ground interaction[J]. International Journal For Numerical and Analytical Methods in Geomechanics. 2012;36:498-516.
[3] Li H, Schindler C. Analysis of soil compaction and tire mobility with finite element method[J]. Proceedings of the Institution of Mechanical Engineers, Part K: Journal of Multi-body Dynamics. 2013;227:275-91.
[4] González Cueto O, Iglesias Coronel CE, Recarey Morfa CA, Urriolagoitia Sosa G, Hernández Gómez LH, Urriolagoitia Calderón G, et al. Three dimensional finite element model of soil compaction caused by agricultural tire traffic[J]. Computers and Electronics in Agriculture. 2013;99:146-52.
[5] Li H, Schindler C. Three-dimensional finite element and analytical modelling of tyre–soil interaction[J]. Proceedings of the Institution of Mechanical Engineers, Part K: Journal of Multi-body Dynamics. 2013;227:42-60.
[6] Smith W, Peng H. Modeling of wheel-soil interaction over rough terrain using the discrete element method[J]. Journal of Terramechanics. 2013;50:277-87.
[7] Li W, Huang Y, Cui Y, Dong S, Wang J. Trafficability analysis of lunar mare terrain by means of the discrete element method for wheeled rover locomotion[J]. Journal of Terramechanics. 2010;47:161-72.
[8] 张锐, 刘芳, 曾桂银, 李建桥. 非规则车轮 散体模拟月壤相互作用仿真系统研究[C]. 颗粒材料计算力学研究进展2012.
[9] Jiang M, Liu F, Shen Z, Zheng M. Distinct element simulation of lugged wheel performance under extraterrestrial environmental effects[J]. Acta Astronautica. 2014;99:37-51.
[10] 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:423-34.
[11] 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:55-65.
[12] 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-37.
[13] Williams J, Perkins E, Cook B. A contact algorithm for partitioning N arbitrary sized objects[J]. Engineering Computations. 2004;21:235-48.
[14] Zang MY, Gao W, Lei Z. A contact algorithm for 3D discrete and finite element contact problems based on penalty function method[J]. Computational Mechanics. 2011;48:541-50.
[15] Balevi?ius R, D?iugys A, Ka?ianauskas R. Discrete element method and its application to the analysis of penetration into granular media[J]. Journal of Civil Engineering and Management. 2004;1:3-14.
[16] Han K, Peric D, Crook A, Owen D. A combined finite/discrete element simulation of shot peening processes–Part I: studies on 2D interaction laws[J]. Engineering Computations. 2000;17:593-620.
[17] 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:61-6.
[18] Zang M, Zhao C. Numerical Simulation of Rigid Wheel Running Behavior on Sand Terrain[C]. 2013.
[19] 庄继德. 计算汽车地面力学[M]: 机械工业出版社; 2002.
