智能轮胎是未来智能化车辆控制系统的重要组成部分，主要承担着实时感知和监测轮胎与路面接触力信息的任务[1]．特别是近年来，随着传感技术和电子技术的快速发展，通过在轮胎内安装传感器监测轮胎动态数据，进而实现轮胎力实时解算的智能轮胎技术受到了国内外研究学者的广泛关注[2]．针对轮胎力学性能的解算，相关研究学者提出了多种轮胎模型[3-6]．其中，魔术公式是一种以试验数据为基础，通过对试验数据精确拟合从而计算轮胎力学性能的方法．目前，魔术公式已经发展更新到最新版本的PAC2002模型[7]，由于魔术公式的开源性，研究学者以魔术公式为基础，拓展了多种应用环境．如文献[8]将魔术公式拓展到冰雪路面的研究中，可用于计算冰雪路面轮胎驱动和制动状态的纵向力计算；文献[9]在魔术公式中引入了温度参数，并通过试验验证了仿真的准确性．这些优化和拓展进一步加深了魔术公式的工程应用价值．重载子午轮胎是重载特种车辆的关键接地和承载部件[10-11]，通常具有垂向载荷范围大、承载高的特点．为保证使用车辆越野性，与乘用车和普通卡车轮胎相比，该型轮胎扁平比较高．为了探究魔术公式在重载子午轮胎力学研究中的应用特性，本研究基于轮胎六分力测试试验，针对重载子午轮胎垂向载荷范围大、轮胎扁平比大的特点，分析了重载工况下轮胎垂向载荷对轮胎侧偏特性和纵滑特性的影响规律；并采用粒子群算法对魔术公式轮胎模型的特征参数辨识进行研究，分析了辨识结果对不同载荷作用下的轮胎力学特性匹配程度与影响因素；最后基于相似性原理对重载工况的魔术公式进行了修正．1 轮胎六分力试验及性能分析根据轮胎动力学原理结合广州孔辉科技轮胎六分力试验设备，开展了不同垂向载荷作用下纯纵滑和纯侧偏工况的轮胎六分力测试试验研究．试验轮胎型号为GL073A，试验平台为六自由度轮胎力学特性试验台，试验路面为4 m聚四氟乙烯涂层面滑台，滑台运行速度200 mm/s，纯纵滑工况试验中滑移率测试范围为-0.012~0.000，纯侧偏工况试验中侧偏角测试范围为-10°~10°．试验分别测试了在轮胎压力为810 kPa工况下，垂向载荷为24.835，34.755，44.669，54.530 kN时的轮胎力学性能，试验结果如图1~3所示，图中：θ为侧偏角；Fy为侧偏力；My为回正力矩；λ为滑移率；Fx为纵向力；Fz为垂向载荷．图1~3结果显示如下．10.13245/j.hust.211007.F001图1侧偏力变化曲线10.13245/j.hust.211007.F002图2回正力矩变化曲线10.13245/j.hust.211007.F003图3纵向力变化曲线a. 当侧偏角不超过1.5°时，侧偏力和侧偏角近似成线性关系，侧偏刚度随着垂向载荷的增大从-3.508 kN/(°)变为-5.805 kN/(°)，回正刚度随着垂向载荷的增大从97 N·m/(°)变为265 N·m/(°)；当侧偏角超过1.5°时，随着侧偏角的增加，侧偏力变化速率减慢，轮胎开始出现部分侧滑，随着侧偏角度的增加，回正力矩开始减小．无量纲侧偏刚度和无量纲回正刚度通常被用来对比不同工况下的操稳能力，计算不同垂向载荷作用下的无量纲侧偏刚度和无量纲回正刚度如图4所示，有Cy=ky/Fz;CM=kM/Fz，式中：Cy为无量纲侧偏刚度；ky为试验侧偏刚度；CM为无量纲回正刚度；kM为试验回正刚度．10.13245/j.hust.211007.F004图4不同垂向载荷作用下的无量纲侧偏特性图4结果显示：随着垂向载荷的增大，该型重载子午轮胎的无量纲侧偏刚度和无量纲回正刚度随之增大；同时，无量纲侧偏刚度和无量纲回正刚度的变化与垂向载荷呈现出明显的非线性关系．b. 试验测试的滑移率范围在-1.2%~0.0%之间，轮胎的滑移率主要是由胎面的弹性变形引起的，轮胎纵滑刚度随着垂向载荷的增大从1.944 00×105变为2.932 99×105．计算不同垂向载荷作用下的无量纲纵滑刚度如图5所示，有Cx=kx/Fz，其中：Cx为无量纲纵滑刚度；kx为试验纵滑刚度．10.13245/j.hust.211007.F005图5不同垂向载荷作用下的无量纲纵滑特性由图5可以看出：随着垂向载荷的增大，该型重载子午轮胎的无量纲纵滑刚度随之减小，轮胎的制动能力越差；无量纲纵滑刚度与垂向载荷也呈现出明显的非线性关系．图1~5结果表明：重载子午轮胎侧偏特性和纵滑特性受垂向载荷变化影响较大，并呈现出较为明显的非线性关系．2 魔术公式参数辨识及载荷匹配分析2.1　特征参数辨识纯工况下魔术公式轮胎模型的特征参数共有58个，其中侧偏力待辨识的特征参数有18个，回正力矩待辨识的特征参数有25个，纵向力待辨识的特征参数有15个[7]．针对魔术公式轮胎模型的多特征参数辨识问题，本研究采用粒子群优化算法进行参数辨识，粒子群优化算法由于参数少且易于实现，在非线性、多峰问题中具有较强的全局搜索能力，在科学研究与工程实践中均得到了广泛关注．假设在一个N维特征参数向量解搜索空间中，有M个粒子组成一个群，其中第i个粒子在搜索空间中的位置可以表示为Xi=(xi1,xi2,⋯,xiN) (i=1,2,⋯,N)．第i个粒子在搜索空间中的飞行速度可以表示为Vi=(vi1,vi2,⋯,viN)．第i个粒子目前搜索到的最优解为pbest=(pi1,pi2,⋯,piN)．整个粒子群目前搜索到的最优解为gbest=(gi1,gi2,⋯,giN)．迭代过程中粒子更新后粒子速度和位置为v(t+1)=wv(t)+c1r1[p(t)-x(t)]+c2r2[g(t)-x(t)]; (1)x(t+1)=x(t)+v(t+1)，式中：v(t)和x(t)分别为目前粒子速度和位置；p(t)和g(t)分别为目前粒子搜索结果和粒子群搜索结果；w为惯性权重；t为迭代次数；c1和c2为学习因子；r1和r2为[0，1]范围内的均匀随机数．式(1)中w表征粒子在多大程度保留原来的速度．相关实验研究发现，较大的惯性权重能够使算法全局收敛能力强，而较小的惯性权重能够使算法局部收敛能力强．因此，本研究选择文献[12]提出的线性递减动态权值策略，其表达式为w=wmax-(wmax-wmin)t/Tmax，式中：Tmax为最大进化代数；wmax和wmin分别表示最大惯性权重和最小惯性权重，在大多数的应用中，通常取wmax=0.9，wmin=0.4．从该型重载轮胎六分力测试试验数据图1~3可以看出：与侧偏力和纵向力规律相比，回正力矩的变化规律在垂向载荷较大时更为明显和全面，因此本研究选择垂向载荷为5.453 0×104 N时的回正力矩与垂向载荷为2.483 5×104 N时的侧偏力和纵向力为基础辨识数据，以此辨识出魔术公式的特征参数对其他垂向载荷工况进行推算．对辨识结果进行残差定量分析与评价，有σ=∑i=1A(ysim-ytest)2/∑i=1Aytest21/2×100%，式中：ysim为辨识结果；ytest为试验结果．辨识和推算残差结果如图6所示，图6结果表明：以基础辨识数据进行辨识的残差均小于10%，辨识结果可信．其中：2.483 5×104 N的侧偏力辨识残差为1.65%，5.453 0×104 N的回正力矩辨识残差为9.1%，2.483 5×104 N的纵向力辨识残差为9.46%，回正力矩辨识残差相对较高是由于如图4所示的回正力矩自身的强非线性，纵向力辨识残差较高的原因是试验数据量相对较少．随着施加垂向载荷与基础辨识数据所对应的垂向载荷的差值逐渐增大，辨识残差也逐渐增大，其中回正力矩辨识残差的增加相对明显，即直接采用未经修正的公式计算得到的数据与试验数据的符合度偏差相对较大．纵向力辨识残差的增加幅度相对较小，这是由于纵滑试验的滑移率范围在-1.2%~0.0%之间，轮胎的滑移率主要是由胎面的弹性变形引起，纵向力呈现出如图3所示的近似线性变化．侧偏力辨识残差的增加幅度最小，一方面是试验数据量的足够充分，另一方面是侧偏力自身的非线性属性与回正力矩相比较小．10.13245/j.hust.211007.F006图6辨识和推算残差2.2　垂向载荷匹配分析按照辨识残差分析，辨识推算结果受垂向载荷差值影响十分明显，且与垂向载荷相关的名义垂向载荷增量在不同垂向载荷作用下偏差较小，总体幅值较小．对侧偏力而言，与垂向载荷相关的系数主要有侧偏峰值因子Dy，Ky及侧偏垂向漂移Svy；纵向力与垂直载荷相关的系数为纵滑峰值因子Dx，Kx及纵滑垂向漂移Svx；回正力矩与垂向载荷相关的系数为残余力矩的峰值因子Dr、轮胎拖距b及Fy．分别对侧偏力、回正力矩和纵向力进行无量纲处理有：CFy0=Fy0/Fz；(2)CMz0=Mz0/Fz；(3)CFx0=Fx0/Fz．(4)从式(2)中无量纲侧偏力表达式可知：无量纲侧偏力消除了Dy和Svy受垂向载荷的影响，此时侧偏力的变化主要由Ky决定．从式(3)中无量纲回正力矩表达式可知：由于无量纲侧偏力趋于一致，无量纲回正力矩无法消除b和Dr受垂向载荷的影响，因此此时回正力矩的变化仍然由b和Dr共同决定．从式(4)中无量纲纵向力表达式可知：无量纲纵向力能够有限地消除Dx和Svx受垂向载荷的影响，且与滑移率呈现出相关性．3 轮胎魔术公式模型修正及验证载荷匹配分析辨识推算结果与试验结果存在误差的主要原因在以下三个方面．a. 侧偏力的误差主要表现为辨识推算的侧偏力曲率与试验侧偏力曲率相比偏小，且随着施加垂向载荷与基础辨识数据所对应的垂向载荷的差值逐渐增大，曲率差越大．这是由于侧偏力的曲率变化主要由Ky决定，魔术公式在重载工况中的Ky计算与垂向载荷呈现出的非线性关系不明显，这与图4分析的试验所得无量纲侧偏刚度不符，因此引入侧偏刚度修正公式对魔术公式进行修正，有ηy=1-py4[(Fz-Fzj)/Fzj]2，其中：py4为增加的侧偏力特征参数；Fzj为基础辨识数据对应的垂向载荷．b. 回正力矩的误差主要表现为辨识推算的回正力矩与试验回正力矩在负侧偏角处呈现出幅值较大的偏差，这主要是由于魔术公式计算的侧偏力在负侧偏角处幅值较大，轮胎拖距和残余力矩也无法消除垂向载荷的影响，因此引入回正力矩修正公式对魔术公式进行修正，有ηM=pz1{1-[(Fzj-Fz)/Fzj]2sign α}，式中pz1为增加的回正力矩特征参数．c. 纵向力的误差主要表现为辨识推算的侧偏力与试验侧偏力相比垂向偏差较大，且随着施加垂向载荷与基础辨识数据所对应的垂向载荷的差值逐渐增大，偏差越大，同时随着滑移率绝对值的增大，这种偏差也呈现出逐渐变大的趋势，但与回正力矩相比，偏差相对较小．这一方面是由于纵向力的曲率变化由Kx决定，另一方面不同滑移率时的垂向漂移Svx差值较大，因此可以从这两个方面进行修正，考虑到滑移率的影响，本研究引入垂向漂移修正公式对魔术公式进行修正，有ηx=px4(Fz-Fzj){1-exp[px5(λ+0.002)]},式中px4和px5为增加的纵向力特征参数．计算结果如图7~9和表1~3所示．10.13245/j.hust.211007.F007图7侧偏力修正结果10.13245/j.hust.211007.F008图8回正力矩修正结果10.13245/j.hust.211007.F009图9纵向力修正结果10.13245/j.hust.211007.T001表1侧偏力残差σFz/kN34.75544.66954.530修正前0.0350.0580.089修正后0.0280.0290.038%10.13245/j.hust.211007.T002表2回正力矩残差σFz/kN34.75544.66954.530修正前0.7450.6110.216修正后0.2090.1820.113%10.13245/j.hust.211007.T003表3纵向力残差σFz/kN34.75544.66954.530修正前0.1930.2200.260修正后0.0620.0400.086%图7~9和表1~3结果表明：针对重载工况修正后的魔术公式轮胎模型，能够有效降低辨识推算残差，提高模型的符合度．其中：侧偏力修正后的推算结果与修正前相比降低了20%~57%的残差；回正力矩修正前后相比降低了47%~72%的残差；纵向力修正前后相比降低了67%~82%的残差，修正后的魔术公式与试验结果符合度较高．4 结论a. 重载子午轮胎侧偏特性和纵滑特性受垂向载荷变化影响较大，并呈现出较为明显的非线性关系．b. 不同垂向载荷作用下侧偏力的变化主要由Ky决定；不同垂向载荷作用下回正力矩的变化由t和Dr共同决定；不同垂向载荷作用下的纵向力受到Dx和Svx影响，且与滑移率呈现出相关性．c. 针对重载工况修正的魔术公式轮胎模型，能够有效降低辨识推算残差．