黃通,劉志浩,高欽和,王冬,馬棟
(火箭軍工程大學(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)用提供理論支撐.
根據(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.
魔術(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
粒子群優(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.
數(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矩陣.
基于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é)果提升相對較弱.
輪胎模型特征參數(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).
在采樣數(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)確.
在采樣數(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ù).
在采樣數(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)整的便捷性.
本文以某型輪胎六分力測試試驗(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ī)律.