基于多体单元组合法的柔性履带系统动力学建模方法研究与行走振动试验分析

李沐桐，何 林，张 轩，陈 华

（广东省现代农业装备研究所，广东 广州 510630）

0 引言

近年来，履带底盘在农业机械领域内的应用逐渐增加，尤其在提出发展丘陵山区农业之后，其较优的抓地力、减震效果、包络特性及越障过沟能力，保证了其在不规则地形和非结构环境下的作业性能相比轮式有着显著的优势[1,2]。随着装备的多样化升级与发展，履带在结构和形式上也逐渐发生了多样化的衍生，但目前在农机装备和智能农业机器人的研制和应用中，多数以主流的橡胶履带为主。在丘陵山区农机装备的示范与应用过程中，常由于缺乏精确的模型和可靠的仿真手段，导致目前多数学者仍不能准确的研究分析履带底盘的行走特性，严重阻碍了履带底盘向多元化、高适应性发展的进程[3]。故此，针对橡胶类柔性履带系统，探索出能用于精确、快速模拟计算的履带建模方法，是为未来农机向智能化、特色化、小型集成化发展的重要科技支撑。

当前多数履带底盘系统主要包括驱动轮、支重轮、导向轮、张紧装置、托带轮等，其底盘也多为刚性非变形式梁架，少数工程师在军事、勘探、航天领域中通过设置弹性元件或多自由度联动机构来实现支重轮的浮动仿形功能[4]。近几年，部分学者针对丘陵山区农用履带底盘设计研制多种可变形浮动轮系底盘，如李沐桐[5]针对果园坡地路面研制一种用于除草机器人行走的旋点式双侧浮动行走底盘、谢铌[6]研制的小型山地履带底盘和聂建军[7]设计的新型弓腰式移动底盘等，均具有较好的坡地仿形能力且增加了爬坡越障能力。然而对于上述新型底盘的试验研究仅针对部分转弯、爬坡和姿态等表象性能的测试，预深入进行履带内部力学行为（如振动滤波特性、变形机理或与土壤内部剪切驱动测试）时却不能有效的做出测定，需结合虚拟仿真试验来进行深入探究。目前国内外对于柔性履带的建模方式，多采用多刚体铰接方[8]、柔性化大变形[9]建模，多数学者利用多体动力学软件ADAMS_ATV和RecurDyn_Track(LM)等内置工具包实现履带建模[10-12]，然而目前主要适合刚性铁板履带，模型相对过于简单固化，不能灵活的满足市场应用和实际产品设计仿真需求，也无法准确地得到近真实的试验结论。

本文将丘陵山地农机装备常用的橡胶履带系统作为研究对象，提出一种多体单元组合法的柔性履带系统建模方式。选取3SYLZ 型果园遥控除草机作为案例，对履带单元、连接形式、驱动系统、接触特性等方面进行了分析，绘制了履带单元三维结构，进一步创建了橡胶履带的整体动力学模型，并对接触类型的优化与参数试验标定。最后结合数值模拟与实地试验对比，验证了目标样机履带系统模型的可靠性，为未来丘陵山区农业机器人行走底盘的设计研究提供技术支撑和模型参考。

1 履带系统基本组成与建模

通过切割真实履带局部样块得到履带内部的基本机构组成，主要包括芯铁、钢丝帘线（纤维线）和橡胶结构，如图1 所示，其组成结构和连接方式对履带单元建模后的连接力学特性设定有着重要的参考依据。

图1 橡胶履带内部结构组成示意图

1.1 履带单元模型创建

首先需要说明的是，本文研究对象主要针对小型农业机器人和特色农机装备，其多数履带尺寸较小，且主要功能多针对越障、爬坡，对地面的仿形能力具有一定的需求，故一般多采用无托带轮的轮系结构，具有易维修、成本低、结构简单的优点。本文以此类履带底盘系统作为案例进行建模研究，不作为通用履带系统的标准参考。

在履带单元建模时，需通过单位节距内的断面尺寸图逐一绘制每个履带单元三维模型，本文所选用的履带来自广东弘科农业机械研究开发有限公司研制的9GC750 型果园智能除草机器人所用，行走底盘的主要参数如表1 所示。其履带单元的断面尺寸由履带生产厂家杭州中策橡胶有限公司提供，规格型号为150×72×30（宽度×节距×单元个数），外表面选取“L”型[S3]类花纹样式，其断面尺寸示意图与三维模型如图2 所示。

表1 整机履带底盘主要参数

图2 履带单元尺寸及其3D模型示意图

1.2 履带缠绕约束设定

截止目前，国内外学者在对履带建模时，主要根据底盘轮系分布位置和节圆直径等参数，逐一绘制连续相切的闭环曲线作为履带中性层轮廓线。本文在前期建模中发现，即使是实际使用中的橡胶履带，其履带单元的弯曲变形能力均较差，若按上述理论轮廓曲线进行缠绕建模，履带单元顺序缠绕到最后的封闭点可能无法得到定位，具体因为理论单元的弯曲度通常远大于实际弯曲度，故将产生较大的周长误差，造成后续的虚拟仿真试验中出现单元间不均匀裂缝，进而发生接触出错，造成了无效模拟计算。

传统刚性铰链式履带单元连接模型中多采用添加铰接副来实现，橡胶履带的连接形式在本质上与其类似，然而差异主要集中在制作工艺过程中的橡胶填充成型环节，相当于在之前履带单元间的铰接副内添加了弹性轴套力，来替代高强度弹性橡胶的柔弹物理特性，而所添加的铰接副回转中心均位于履带单元内部的中性层截面线上。因此，在上述履带单元模型建成后，定位到履带单元中心剖面处两侧的中性层端点位置，以此来设定履带每个单元间的铰接点，如图3 所示。按照底盘轮系前提确定分布位置和默认履带节距（72 mm）及单元个数（30 个），将各单元添加铰接约束，并按同方向（逆时针或顺时针）逐个调节各单元的安装角度，保持履带单元不嵌入轮系三维特征内，避免导入动力学软件后出现无法设置接触错误。

图3 履带各单元连接设定与缠绕初始状态

此处值得说明的是，上述建模与约束设定环境均可以通用于多数机械工业类三维建模软件，需要注意的是，在设置顺序缠绕履带时，需将张紧装置或张紧轮调节到最小位置处，保证履带单元总体缠绕后成松弛状态，以方便后续在动力学软件中计算最佳张紧长度。

2 运动系统建模

2.1 基本运动副添加

在上述工作的前提下，将底盘轮系模型与缠绕后的各履带单元总成的几何体以step 格式导入到RecurDyn 动力学分析软件中。首先分别设置驱动轮转动副（DW_rev）、支重轮转动副（SW_rev）和导向轮转动副（Idle_rev），如图4 所示；其次按图4（a）中所示，以每个履带前期设置的单元铰接点为中心，以履带宽度的中间面为基准，逐次添加履带单元间的平面副（Inplane_(i)）(其中i 表示第i 个履带单元)。此处需要说明的是，在实际运动中，横向力变化主要表征为与轮系侧壁的接触压力，对行走动力的结果影响可忽略不计，因此在履带单元间设置了平面约束，以保证每个履带单元不受横向接触力干扰，同时也减少后续不必要的计算量。最后，如图4（b）所示，同以上方式类似的将履带单元逐次的添加转动副（seg_rev（i）），直至完成闭环约束。至此，完成履带的基本运动副添加。

图4 刚性履带单元的初始运动副添加

2.2 衬套力的添加

据上文所述，由于橡胶履带为模具式成型，与铁板履带不同，在每个履带单元间配装有高弹性的纤维线，并粘附着高密度橡胶结构，通常在强度和韧性上具有特殊的优势。因此，在优化橡胶类柔性履带的建模中，本文选取衬套力作为代替刚性单元间的复杂弹性力，在上述每个履带单元间所添加的转动副上进一步增加衬套力(Bushing(i))，其中包括三轴方向在平移与转动自由度内的刚度值和阻尼值，本文根据目前橡胶履带的参考文献默认设定刚度（Stiffness）为100，阻尼（Damping）为0.1，如图5 所示。需要说明的是，由于不同履带的泊松比、密度、弹性模量等物理特性差异较大，故在此处设定时要根据基础力学测试依据恰当的选择其值，这与仿真计算结果和真实差异性具有直接关系。本文所设数值主要以默认案例提供参考方法，以仿真橡胶履带在外力过程中的抗弯、回弹等特性。

图5 履带单元间的衬套力添加示意图

2.3 内部接触设定

履带的接触对象主要为两种，一种为内部接触，即履带内壁面与底盘系统内的轮系接触作用，另一种为外部接触，即履带外表面或履带花纹面与路面表层的接触、碰撞或摩擦作用。首先需要确定橡胶履带内壁面结构与底盘轮系的接触形式。在部分动力学软件的履带工具插件中主要借助命令流形式实现每个轮系对各履带单元的批量接触设置，但多为实体-实体形式的接触（solid-solid），在计算中大幅增加了不必要的计算量，且常附加着刚性轮系的实体变形计算。这样不仅在很大程度上降低了仿真计算速度，其结果也将与实际产生巨大的差异。

在本文的前期接触设定中，为提高计算速度，起初对各接触轮进行了中心环线提取，以支重轮的中心环线（SW_curve）为例，与履带单元铁芯处设置的截面线（seg(i)_inner_curve）间添加2D 曲线接触（cur-cur），但仿真中发现，由于曲线线条的接触刚度和阻尼需设定较大值才能保证不出现穿透现象，故此选用面-面形式的接触（sur-sur），如图6 所示。同样以支重轮为例，在支重轮接触表面所设置的外环面（SW_Surface）与履带单元铁芯内所设置的内表面（seg(i)_inner_surface）之间，添加面接触，其主要接触参数如表2 所示。

图6 履带单元-轮系接触示意图

表2 履带单元与轮系接触参数

表2 所示的参数为默认参数，在特殊案例中，其参数需根据对比实际试验进行虚拟标定确定，需要说明的是，在本案例中履带单元与接触轮系的位置较近，一般选取最大穿透值为1 mm，若起始的直线距离相差较远，仿真计算时由于高速冲击导致的惯性碰撞易发生溢出现象，需改变最大穿透值为10 mm，以保证碰撞过程不发生失真问题。

3 初始状态计算

根据上文所述，已完成了履带几何模型的创建、运动副添加、特殊力的赋予以及内部接触的设定，然而此时张紧杆处于缩回状态，履带整体还处于松弛状态，需要张紧滑杆（Tension Rod）前端添加压簧进而计算出添加合理预紧力后的平衡位置，并通过弹性张紧计算得到张紧滑杆在装配时的最优伸长量，为实际固定张紧的安装提供理论依据。如图7所示，因此需在仿真模型中的张紧滑杆与底盘间添加平动副（模型文件名称为：TR_tra，tra 为translate简写），同时在张紧滑杆的起始端空间内添加压簧力（Tension_Spring）。参考文献[13-15]可知，缓冲弹簧预紧力不能太小,该预紧力保证缓冲弹簧不会在外来微小冲击时产生附加变形而引起履带经常振动，同时又必须保证导向轮在各种正常工作情况下(如倒退、转弯等)不会后移，以避免脱带。但预紧力也不宜过大，否则就失去了缓冲作用。故一般选择预紧力P0范围为

图7 平动副与张紧压簧的添加过程

缓冲弹簧最大变形时的弹簧压力（被限位时）为

式中 G为整机质量，kg。

在此前提下根据整机参数和部件结构设置弹簧的自然长度、丝径、圈数和刚度等参数，进一步开始仿真计算，仿真时间设定为0.5 s，步长为1000，在后处理计算结果中以张紧滑杆作为标记对象，得到在x 轴方向上的位置变化图表。如图8 所示，张紧滑杆相对于坐标系中心点的x 轴向初始位置约为799.20 mm 处，当压簧自然释放后张紧滑杆在x 轴方向的位移瞬间增大，在与橡胶履带的非弹性碰撞后其变化量逐渐收敛在807.50 mm 处，故仿真计算的张紧滑杆伸长量约为8.30 mm，其仿真计算结果作为与实体张紧试验的对比依据。

图8 张紧滑杆x方向的位移变化曲线

4 驱动添加与设定

在完成初始化预紧后的履带已呈张紧状态，需要进一步添加驱动轮的动力参数，如图9 所示，在上述已添加的驱动轮转动副（DW_rev）上添加驱动参数（DW_rev_motion），其中本文设定的参数单位为（°）/s。按本案例计算，当底盘以1 km/h的速度前进时，转换为驱动轮节圆角速度为138.9（°）/s，即数值设定为138.9 d*time。需要注意的是，本文规定在仿真预紧力计算0.5 s后再激活驱动参数生效，故在驱动设定中采用了IF条件函数，即表达式为：“IF（time-0.5:0,0，-138.9d*time）”，进而设定总仿真时间为5 s，步长为1000，至此，完成了橡胶柔性履带底盘的空载仿真运行一次。

图9 驱动参数的添加

5 履带单元柔性化处理

根据上文梳理的动力学模型的设定和仿真计算可知，目前默认的刚性履带单元设置了材料物理参数和基本的运动副，并配合添加了衬套力等来模拟橡胶履带单元的变形及回弹特性，该方法能适用于多数普通橡胶履带底盘，而针对特殊大变形履带材料，上述方式不能较好的计算每个单元结构的细化变形。当履带单元个数或轮系个数增多时，计算结果将与实际产生较大差异，因此，需要在上述基础上，将每个履带单元进行柔性化处理，如图10 所示，其中图（a）为无其他单元拉力情况时的未变形状态，图（b）为当履带单元处于大变形处位置时，同时受两侧履带的切向拉力FT和内表面轮系的支撑力Fns作用，进而产生了履带单元的内部独立变形，并可根据实际情况，调节履带网格的密度、尺寸和极限应力的数值。其可作为特殊计算需求时的扩展手段。

图10 驱动参数的添加

需要说明的是，由于对履带单元进行了有限元网格化处理，在针对大量履带单元的联合计算时，一般需要较长的时间才可完成计算，因此在本文的案例中未采用该方法进行拓展计算，其常规橡胶履带材料多数高密度橡胶制品，普遍视为局部的刚性单元，不考虑内部微小变形情况。

6 整机行走虚拟试验

6.1 环境设定

考虑到路面环境的多样性和复杂性，本文初步选用不可压缩干性硬质土路面作为行走载体进行初步试验验证，以保证其仿真计算结果与实际试验的对比具有实践价值。首先，在模型中创建1500 mm×5000 mm×200 mm的实体地块，并转换为step格式导入到RecurDyn中，并参考文献[16-19]中对路面土壤的标定参数，赋予硬土基本参数值。本文模型中选用道路密度为2.6 g/cm3、弹性模量为1×107MPa、泊松比为0.25，其中橡胶履带单元与土壤表层的接触模型为面-面（sur-sur）形式，静摩擦系数设定为0.2、动摩擦系数设定为0.1、接触刚度系数设定为100000、阻尼系数为100，仿真试验初步采用直线行走计算，速度为1 km/h，仿真时间为10 s，步长设定为1000。

图11 整机仿真环境示意图

6.2 指标选取与试验设计

根据文献[20-23]参考可知，本文对于仿真行走结果的评价主要以底盘不同点位置处在垂直方向的振动变化为参考指标，具体为Z 轴方向的加速度变化规律和变化范围作为分析对象。考虑到履带底盘内部的前端支重轮为接地的初始轮，在行走过程中的运动状态变化具有较显著的特点，故在上述仿真结果的前提下，初选履带底盘系统中的前端支重轮为观测点，在实际试验中选取DYTRAN-4400B3 型三轴加速度传感器作为记录装置，并同时在模型内前支重轮回转中心的侧板表面固装传感器模型，模型与实物安装如图12 所示。

图12 试验方法与现场测试

6.3 数据后处理与对比分析

通过后处理功能提取出传感器安置点处在垂直方向（虚拟环境方向为Z 轴，实物为Y 轴）的振动加速度变化曲线，并通过低通四阶低通滤波器处理后得到加速度变化曲线对比图表，如图13 所示。其中红色为仿真计算后得到的Z 轴方向加速度变化值，蓝色为实地测试后记录在传感器内的垂直方向加速度变化值，通过曲线数据分析可发现，在0～1 s 内，由于初始启动产生的加速度突变在仿真与实测中的极限值可视为一致，在1 s 后，驱动轮进去匀速转动，并带动整机直线前进，其仿真数据的最大加速度变化范围在0～21 m/s2内，整体呈均匀收敛状态，而实测值的变化范围在25～28 m/s2范围内，产生误差的主要问题可能由于仿真计算中的动力装置不进行振动计算，而实际中的汽油机、皮带和液压系统内部的运转振动较为复杂，在驱动过程中产生交叉影响的强化作用，因此，导致了加速度变化范围过大的误差。总体而言，在1～8 s 范围内，仿真结果与实测结果的加速度变化趋势几乎接近一致，由于路面-履带-轮系的接触力学作用非常复杂，对于产生的振动范围误差在本文的研究中初步视为可接受[24-28]。

图13 前支重轮中心处垂直方向的振动加速度变化曲线

7 结语

本文主要针对丘陵山区农机装备与农业机器人常用的履带底盘行走机构，提出了一种多刚体单元组合形式的柔性履带建模方法，并借助当前热点产品——果园智能除草机器人作为研究对象，通过整机建模、运动副创建、衬套力添加、接触模型分析等逐步对橡胶履带系统完成了数字化动力模型的建设。根据实际需求，研究了刚-柔模型转化的方法，通过数值模拟与实物测试试验，验证了该建模方法的可行性，分析探讨了产生误差的原因和改进方向，为未来智能农机行走可靠性的研究提供重要手段和参考依据。