重載車輛輪胎模型參數(shù)辨識與靈敏度分析

黃通，劉志浩，高欽和，王冬，馬棟

（火箭軍工程大學(xué)導(dǎo)彈工程學(xué)院，陜西西安 710025）

重載子午輪胎作為重載車輛與地面直接接觸的部件，其力學(xué)性能直接影響車輛動力學(xué)和行駛平順特性.與常見轎車輪胎相比，重載子午輪胎具有氣壓高、阻尼低、花紋粗大、扁平率大的特點(diǎn).研究分析重載子午輪胎的力學(xué)性能對設(shè)計(jì)和分析先進(jìn)的底盤系統(tǒng)和優(yōu)化先進(jìn)車輛系統(tǒng)結(jié)構(gòu)至關(guān)重要［1］.為了研究輪胎力學(xué)性能，國內(nèi)外研究者提出了多種輪胎模型［2-3］.

魔術(shù)公式（Magic Formula，MF）輪胎模型以試驗(yàn)數(shù)據(jù)為基礎(chǔ)，通過精確擬合試驗(yàn)數(shù)據(jù)，描述輪胎的力學(xué)性能.基于MF 模型，相關(guān)研究者已發(fā)展了許多適應(yīng)于各種工況的不同改進(jìn)形MF 模型［4-5］.PAC2002模型是Pacejka［6］不斷更新和完善的一種廣泛應(yīng)用于車輛動力學(xué)仿真和分析的MF模型，能夠研究不同載荷、傾角和胎壓下的輪胎力學(xué)性能，具有較強(qiáng)的工程應(yīng)用背景.但特征參數(shù)繁多，而且高度非線性，給特征參數(shù)的辨識帶來很大困難.

目前，輪胎模型的特征參數(shù)辨識基本采用數(shù)值優(yōu)化算法和智能搜索算法.遺傳算法因?yàn)榫哂休^強(qiáng)的魯棒特性，最先開始被應(yīng)用于輪胎模型的特征參數(shù)辨識［7-8］.遺傳算法雖然能夠在全局范圍內(nèi)逼近最優(yōu)解，但存在局部搜索能力較差，收斂速度慢的問題.為此，相關(guān)研究者將全局范圍逼近較強(qiáng)的智能搜索算法和局部搜索較強(qiáng)的數(shù)值優(yōu)化算法結(jié)合起來，產(chǎn)生了多種輪胎參數(shù)辨識的混合優(yōu)化算法，這些混合算法普遍提升了輪胎模型的辨識精度［9-11］.

基于此，本文選擇適合在動態(tài)和多目標(biāo)優(yōu)化環(huán)境中尋優(yōu)的粒子群算法和Levenberg-Marquardt（LM）算法進(jìn)行純縱滑和純側(cè)偏工況下的輪胎模型的特征參數(shù)辨識，并對辨識后的輪胎模型進(jìn)行殘差分析和結(jié)果驗(yàn)證，最后基于Sobol 靈敏度分析對該型輪胎進(jìn)行參數(shù)分析和優(yōu)化，研究結(jié)果可以為該型輪胎的應(yīng)用提供理論支撐.

1 輪胎六分力測試試驗(yàn)

根據(jù)輪胎動力學(xué)原理，結(jié)合現(xiàn)有的試驗(yàn)設(shè)備，進(jìn)行不同垂向載荷作用下純縱滑和純側(cè)偏工況的輪胎六分力測試試驗(yàn).試驗(yàn)輪胎型號為GL073A，試驗(yàn)平臺為六自由度輪胎力學(xué)特性試驗(yàn)臺，試驗(yàn)路面為4 m水泥路面滑臺，滑臺運(yùn)行速度為200 mm∕s.純縱滑工況試驗(yàn)中滑移率為-0.018～0；在純側(cè)偏工況試驗(yàn)中，側(cè)偏角為-10°～13°.輪胎六分力測試試驗(yàn)如圖1所示.

圖1 輪胎六分力測試試驗(yàn)Fig.1 Six component force test for tire

在輪胎壓力為810 kPa 下，當(dāng)垂向載荷分別為24 800 N、34 700 N 和44 600 N 時，分別測試輪胎力學(xué)性能.測試試驗(yàn)結(jié)果分別如圖2～圖4所示.

圖2 側(cè)偏力與側(cè)偏角關(guān)系曲線Fig.2 Curve of sideforce and sideslip angle

圖3 回正力矩與側(cè)偏角關(guān)系曲線Fig.3 Curve of aligning torque and sideslip angle

圖4 縱向力與滑移率關(guān)系曲線Fig.4 Curve of longitudinal force and slip ratio

根據(jù)輪胎力學(xué)，側(cè)偏剛度是側(cè)偏角度為0 時的側(cè)偏力曲線斜率；回正剛度是側(cè)偏角度為0 時的回正力矩曲線斜率；縱滑剛度是滑移率為0 時的縱向力曲線斜率.由圖2～圖4可知：

1）當(dāng)側(cè)偏角不超過6°時，側(cè)偏力與側(cè)偏角呈線性關(guān)系，側(cè)偏剛度從3 828 N∕（°）增加到5 320 N∕（°）；

2）側(cè)偏角為6°時，回正力矩達(dá)到最大值，回正剛度從151 N·m∕（°）增加到287 N·m∕（°）；

3）縱滑剛度隨垂向載荷增大，從301 799 增加到375 493.

2 魔術(shù)公式輪胎模型

魔術(shù)公式輪胎模型是基于試驗(yàn)數(shù)據(jù)通過一個公式即可表示輪胎的縱向力、側(cè)偏力和回正力矩，公式中的特征參數(shù)物理意義明確.

基本的魔術(shù)公式為純側(cè)偏和純縱滑工況下的側(cè)偏力和縱向力：

式中：y為純縱滑工況的縱向力或純側(cè)偏工況的側(cè)偏力；x為對應(yīng)的滑移率或側(cè)偏角；B為剛度因子；C為形狀因子；D為峰值因子；E為曲率因子.

對基本魔術(shù)公式進(jìn)行修正，可以得到純側(cè)偏工況下的側(cè)偏力計(jì)算公式：

式中：γ為外傾角；Fwz為輪胎垂向載荷；Fwz0為輪胎名義垂向載荷；dfwz為名義垂向載荷增量；Cy、Dy、Ey、Ky、By為側(cè)偏工況側(cè)偏力計(jì)算因子；SHy、SVy分別為側(cè)偏工況橫向和垂向漂移量；α為側(cè)偏角；αy為修正側(cè)偏角；Fy0為側(cè)偏力；py為側(cè)偏力特征參數(shù).

在純側(cè)偏工況中，側(cè)偏力計(jì)算公式中待辨識的特征參數(shù)共18個［12］，如表1所示.

表1 側(cè)偏力的特征參數(shù)Tab.1 Characteristic parameters of sideforce

純側(cè)偏工況下的回正力矩可以用輪胎拖距和殘余力矩來計(jì)算：

式中：t為輪胎拖距；Mzr為殘余力矩.均可采用余弦函數(shù)表達(dá).

式中：Ct、Dt、Et、Kt、Bt為側(cè)偏工況回正力矩計(jì)算因子；SHt、SHr分別為側(cè)偏工況拖距和殘余力矩的橫向漂移量；αt為拖距修正側(cè)偏角；αr為殘余力矩修正側(cè)偏角；R0為輪胎半徑；qz為特征參數(shù).

式（6）中各因子均能由類似式（3）的計(jì)算公式表達(dá)出來，在純側(cè)偏工況中，回正力矩計(jì)算公式中待辨識的特征參數(shù)共25個［12］，如表2所示.

表2 回正力矩的特征參數(shù)Tab.2 Characteristic parameters of aligning torque

對基本魔術(shù)公式進(jìn)行修正，可以得到純縱滑工況下的縱向力計(jì)算公式：

式中：Cx、Dx、Ex、Kx、Bx為側(cè)偏工況回正力矩計(jì)算因子；SHx、SVx分別為縱滑工況橫向和垂向漂移量；λ為滑移率；λx為修正滑移率；Fx0為縱向力；px為特征參數(shù).

在純縱滑工況中，縱向力計(jì)算公式中待辨識的特征參數(shù)共15個［12］，如表3所示.

表3 縱向力的特征參數(shù)Tab.3 Characteristic parameters of longitudinal force

3 特征參數(shù)辨識

3.1 粒子群優(yōu)化算法

粒子群優(yōu)化算法的提出來源于鳥類群體行為的規(guī)律性，通過模擬鳥類的覓食行為，將尋求優(yōu)解問題的搜索空間比作鳥類的飛行空間，每只鳥抽象成一個沒有質(zhì)量和體積的粒子，用該粒子表征問題的一個可能解，從搜索空間中的隨機(jī)解出發(fā)，通過迭代運(yùn)算尋找最優(yōu)解.粒子群優(yōu)化算法由于參數(shù)少且易于實(shí)現(xiàn)，在非線性、多峰問題中具有較強(qiáng)的全局搜索能力，在科學(xué)研究與工程實(shí)踐中均得到廣泛關(guān)注.

對于魔術(shù)公式輪胎模型的特征參數(shù)辨識問題，假設(shè)在一個D維特征參數(shù)向量解搜索空間中，有N個粒子組成一個群，其中第i個粒子在搜索空間中的位置可以表示為：

第i個粒子在搜索空間中的“飛行速度”可以表示為：

第i個粒子目前搜索到的最優(yōu)解為：

整個粒子群目前搜索到的最優(yōu)解為：

迭代過程中粒子更新后的粒子速度和位置分別為：

式中：v(t)、x(t)分別為目前粒子速度和位置；w為慣性權(quán)重；t為迭代次數(shù)；c1、c2為學(xué)習(xí)因子；r1、r2為［0，1］內(nèi)的均勻隨機(jī)數(shù).

式（14）中慣性權(quán)重w表示粒子在多大程度保留原來的速度.相關(guān)試驗(yàn)研究發(fā)現(xiàn)，較大的慣性權(quán)重能夠使算法全局收斂能力強(qiáng)，而較小的慣性權(quán)重能夠使算法局部收斂能力強(qiáng).因此，本文選擇由Shi 等［13］提出的線性遞減動態(tài)權(quán)值策略，其表達(dá)式如下：

式中：Tmax為最大進(jìn)化代數(shù)；wmax、wmin分別為最大慣性權(quán)重和最小慣性權(quán)重.在大多數(shù)應(yīng)用中，通常取wmax=0.9，wmin=0.4.

3.2 Levenberg-Marquardt優(yōu)化算法

數(shù)值優(yōu)化算法具有相對較強(qiáng)的局部搜索能力.本文將粒子群優(yōu)化算法辨識的結(jié)果作為數(shù)值優(yōu)化算法的初始值，選擇收斂精度較高、非奇異性較好的Levenberg-Marquardt（LM）數(shù)值優(yōu)化算法進(jìn)行局部最優(yōu)搜索.

LM 算法通常被視為Gauss-Newton 算法和最速下降法的結(jié)合，其基本原理公式為：

式中：x為特征參數(shù)向量；g(t)為殘差；μ為參數(shù)因子；J為g(t)的Jacobian矩陣.

3.3 辨識結(jié)果

基于MATLAB 編程實(shí)現(xiàn)粒子群優(yōu)化算法和LM優(yōu)化算法.首先，利用粒子群算法的全局搜索能力辨識出魔術(shù)公式輪胎模型特征參數(shù)的近似解，并將其作為LM 算法的初始值進(jìn)行局部精確搜索，得到特征參數(shù)的最終解.然后，采用式（18）對辨識結(jié)果進(jìn)行殘差定量分析與評價.

輪胎側(cè)偏力、回正力矩和縱向力的辨識結(jié)果分別如圖5～圖7所示.

圖5 側(cè)偏力辨識結(jié)果Fig.5 Identification results of sideforce

圖6 回正力矩辨識結(jié)果Fig.6 Identification results of aligning torque

圖7 縱向力辨識結(jié)果Fig.7 Identification results of longitudinal force

由圖5～圖7 可以看出，通過粒子群優(yōu)化算法的辨識結(jié)果與試驗(yàn)數(shù)據(jù)大致吻合，說明該算法能夠有效獲得魔術(shù)公式輪胎模型特征參數(shù)的近似解；通過混合優(yōu)化算法的辨識結(jié)果與試驗(yàn)數(shù)據(jù)吻合度較高，辨識結(jié)果殘差相對減小，說明混合優(yōu)化算法能夠獲得更為精確的解.

受試驗(yàn)數(shù)據(jù)采集量多的影響，垂直載荷為44 100 N 時的辨識殘差相對于垂直載荷分別為24 500 N 和34 300 N 時的辨識殘差要低，辨識結(jié)果與試驗(yàn)數(shù)據(jù)的吻合度相對較高.由圖7 可知，在縱向力辨識結(jié)果中，由于試驗(yàn)數(shù)據(jù)的類線性特征，混合優(yōu)化算法的辨識結(jié)果對粒子群算法的辨識結(jié)果提升相對較弱.

4 多特征參數(shù)敏感性分析

輪胎模型特征參數(shù)靈敏度能夠反映特征參數(shù)變化對輪胎模型響應(yīng)的影響程度，靈敏度較大的特征參數(shù)對輪胎模型擁有更強(qiáng)的調(diào)整能力，通?？梢援?dāng)作主導(dǎo)特征參數(shù)以反映輪胎模型的變化情況，特別是魔術(shù)公式這類特征參數(shù)較多的輪胎模型，確定主導(dǎo)特征參數(shù)能夠有效提高模型調(diào)整的便捷性.

Sobol 靈敏度分析方法是基于多重積分的分解方法，可以對非線性的復(fù)雜模型進(jìn)行分析，研究各參數(shù)同時變化下各參數(shù)之間的耦合作用.本文選擇Sobol 靈敏度分析方法對魔術(shù)公式輪胎模型特征參數(shù)的影響進(jìn)行分析，并以各特征參數(shù)的一階靈敏度和總階靈敏度作為評價標(biāo)準(zhǔn).

4.1 側(cè)偏力特征參數(shù)

在采樣數(shù)為5 000時，側(cè)偏力特征參數(shù)的一階和總階靈敏度分析結(jié)果如圖8所示.

圖8 側(cè)偏力特征參數(shù)靈敏度Fig.8 Sensitivity of characteristic parameters for sideforce

由圖8 可知，垂直漂移因子特征參數(shù)pvy1的一階靈敏度最大，約為0.35，總階靈敏度稍大，約為0.18，說明pvy1不僅對模型響應(yīng)影響較大，而且與其他參數(shù)存在耦合關(guān)系，這是因?yàn)樵谑剑?）和式（3）中，pvy1作為一個疊加項(xiàng)獨(dú)立存在，對模型影響比較直接.重點(diǎn)分析剩余一階靈敏度較大的剛度因子特征參數(shù)pky1和pky2和水平漂移因子特征參數(shù)phy2.

在pky1和pky2以及phy2不同值下，各主要特征參數(shù)總階靈敏度如圖9所示.

圖9 側(cè)偏力主要特征參數(shù)的總階靈敏度Fig.9 Distribution of total order sensitivity for sideforce main characteristic parameters

由圖9可知：

1）峰值因子特征參數(shù)pdy1受pky1、pky2和phy2的影響較大，這與圖8 所示的pdy1一階靈敏度較小，總階靈敏度較大的結(jié)果是一致的.其靈敏度隨剛度因子特征參數(shù)pky1和pky2的增大而減小，隨水平漂移因子特征參數(shù)phy2的增大而增大.

2）垂直漂移因子特征參數(shù)pvy1靈敏度隨剛度因子特征參數(shù)pky1和pky2的增大而增大，隨水平漂移因子特征參數(shù)phy2的增大而減小.

3）剛度因子特征參數(shù)pky1和pky2以及水平漂移因子特征參數(shù)phy2之間的耦合靈敏度較高，pky1和pky2相互耦合靈敏度呈正比關(guān)系.當(dāng)phy2為14 時，pky1和pky2的靈敏度均達(dá)到最高.

在側(cè)偏力18 個特征參數(shù)中，選擇剛度因子特征參數(shù)pky1和pky2、水平漂移因子特征參數(shù)phy2和垂直漂移因子特征參數(shù)pvy1為主導(dǎo)特征參數(shù).通過調(diào)整這4個參數(shù)來反映魔術(shù)公式輪胎模型側(cè)偏力的變化規(guī)律，并選擇峰值因子特征參數(shù)pdy1為主要特征參數(shù)，可以通過調(diào)整峰值因子特征參數(shù)，使主導(dǎo)特征參數(shù)更為準(zhǔn)確.

4.2 回正力矩特征參數(shù)

在采樣數(shù)為5 000時，回正力矩特征參數(shù)的一階和總階靈敏度分析結(jié)果如圖10所示.

圖10 回正力矩特征參數(shù)靈敏度Fig.10 Sensitivity of characteristic parameters for aligning torque

由圖10 可知，一階靈敏度大于0.10 的特征參數(shù)有qez1、qez4、qbz1和qhz1.其中，曲率因子特征參數(shù)qez1的一階靈敏度最大，約為0.38，但總階靈敏度相對較小，說明qez1對模型影響較大，但與其他參數(shù)耦合相對較??；曲率因子特征參數(shù)qez4和水平漂移因子特征參數(shù)qhz1也存在相同的規(guī)律.而峰值因子特征參數(shù)qdz1和形狀因子特征參數(shù)qcz1雖然一階靈敏度較小，但總階靈敏度相對較大，qdz1的總階靈敏度最大，約為3.75.因此，重點(diǎn)分析總階靈敏度較大的qdz1和qcz1，如圖11（a）和圖11（b）所示.

圖11 回正力矩主要特征參數(shù)的總階靈敏度Fig.11 Distribution of total order sensitivity for aligning torque main characteristic parameters

由圖11可知，峰值因子特征參數(shù)qdz1主要與形狀因子特征參數(shù)qcz1耦合相關(guān)，二者相互耦合靈敏度較高，隨著qdz1的增加，qcz1靈敏度增加幅度較小，基本控制在0.4～0.45，其他特征參數(shù)靈敏度也相對穩(wěn)定.隨著qcz1的增加，qdz1靈敏度增長幅度較大，當(dāng)qcz1為-0.5 時，qdz1靈敏度達(dá)到最大，約為0.95；同時隨著qcz1的增大，其他特征參數(shù)的靈敏度均逐漸減小.

在回正力矩的25 個特征參數(shù)中，選擇曲率因子特征參數(shù)qez1和qez4、剛度因子特征參數(shù)qbz1、水平漂移因子特征參數(shù)qhz1以及峰值因子特征參數(shù)qdz1為主導(dǎo)特征參數(shù)；選擇形狀因子特征參數(shù)qcz1為主要特征參數(shù).

4.3 縱向力特征參數(shù)

在采樣數(shù)為5 000時，縱向力特征參數(shù)的一階和總階靈敏度分析結(jié)果如圖12所示.

圖12 縱向力特征參數(shù)靈敏度Fig.12 Sensitivity of characteristic parameters for longitudinal force

由圖12 可知，剛度因子特征參數(shù)pkx1的一階靈敏度和總階靈敏度最大，分別約為0.37 和0.25，說明pkx1對模型響應(yīng)和其他參數(shù)的耦合關(guān)系均有較大影響.剩余一階靈敏度較大的有：水平漂移因子特征參數(shù)phx1、垂直漂移因子特征參數(shù)pvx1和峰值因子特征參數(shù)pdx1.其中pvx1為模型疊加項(xiàng)且總階靈敏度較小，pdx1和phx1的總階靈敏度都相對較大.因此，重點(diǎn)分析峰值因子特征參數(shù)pdx1和水平漂移因子特征參數(shù)phx1變化下的總階靈敏度分布規(guī)律.

pdx1和phx1變化下的總階靈敏度如圖13 所示.由圖13 可知，隨著pdx1的增大，剛度因子特征參數(shù)pkx1和水平漂移因子特征參數(shù)phx1均表現(xiàn)出較高的靈敏度，并隨之增加；其他特征參數(shù)靈敏度隨著pdx1的增大而逐漸減小.隨著水平漂移因子phx1的逐漸增大，pdx1和pkx1靈敏度雖然占據(jù)較大分量，但靈敏度總體呈減小趨勢.

圖13 縱向力主要特征參數(shù)的總階靈敏度Fig.13 Distribution of total order sensitivity for longitudinal force main characteristic parameters

因此，在縱向力15 個特征參數(shù)中，選擇剛度因子特征參數(shù)pkx1、水平漂移因子特征參數(shù)phx1、垂直漂移因子特征參數(shù)pvx1以及峰值因子特征參數(shù)pdx1為主導(dǎo)特征參數(shù).

將通過Sobol 靈敏度分析選擇出的主導(dǎo)特征參數(shù)代入輪胎模型進(jìn)行計(jì)算，結(jié)果如表4所示.

表4 辨識結(jié)果對比Tab.4 Comparison of identification results

由表4 可知，采用Sobol 靈敏度分析所得的主導(dǎo)參數(shù)進(jìn)行辨識的結(jié)果與直接采用混合優(yōu)化進(jìn)行辨識的結(jié)果相比，辨識殘差均相對增大，但增大的幅值較小.其中垂向載荷為34 300 kN 的回正力矩殘差增幅相對較大，殘差增幅為0.138%，這是由于主導(dǎo)參數(shù)辨識省略了其他特征參數(shù)影響，使殘差增大，但仍小于單獨(dú)粒子群算法的辨識殘差.而穩(wěn)定迭代次數(shù)相對于直接采用混合優(yōu)化進(jìn)行辨識的結(jié)果減小明顯，模型收斂速度增加，最大增幅為30.4%，對輪胎模型擁有較強(qiáng)的調(diào)整能力，采用Sobol 靈敏度分析所得出的主導(dǎo)參數(shù)能夠有效提高模型調(diào)整的便捷性.

5 結(jié)論

本文以某型輪胎六分力測試試驗(yàn)為基礎(chǔ)，對該型重載子午輪胎的魔術(shù)公式輪胎模型進(jìn)行參數(shù)辨識和靈敏度分析.得出以下結(jié)論：

1）基于粒子群算法能夠?qū)崿F(xiàn)對魔術(shù)公式輪胎模型特征參數(shù)的辨識，且辨識結(jié)果殘差分布在5%左右.混合優(yōu)化算法能夠提高特征參數(shù)辨識的精確性，將辨識結(jié)果殘差控制在5%以內(nèi).

2）基于Sobol 靈敏度分析方法可以從魔術(shù)公式輪胎模型的58 個特征參數(shù)中選擇出13 個特征參數(shù)作為主導(dǎo)特征參數(shù)，可以優(yōu)先調(diào)整這13 個主導(dǎo)參數(shù)來控制魔術(shù)公式輪胎模型的變化規(guī)律.