劉付山,曾志平,2,王衛(wèi)東,2
(1.中南大學(xué) 土木工程學(xué)院,湖南 長沙 410075;2.中南大學(xué) 重載鐵路工程結(jié)構(gòu)教育部重點(diǎn)實(shí)驗(yàn)室,湖南 長沙 410075)
軌道隨機(jī)不平順是由不同波長、不同幅值、不同相位的不同不平順波疊加而成的隨線路里程變化的復(fù)雜隨機(jī)過程[1]。軌道隨機(jī)不平順常通過線路實(shí)地測(cè)量獲得,并采用功率譜密度函數(shù)表示其統(tǒng)計(jì)特征,我國及英國、日本、德國等都測(cè)定了各自的軌道不平順譜密度和相關(guān)函數(shù)。然而,在非線性車輛-軌道耦合系統(tǒng)中,功率譜密度函數(shù)無法直接輸入動(dòng)力系統(tǒng),最有效的辦法是以軌道隨機(jī)不平順的時(shí)域樣本作為輸入。翟婉明等[2]提出一種根據(jù)軌道隨機(jī)不平順功率譜求出頻譜的幅值和隨機(jī)相位,再通過傅里葉逆變換得到軌道隨機(jī)不平順時(shí)域樣本的數(shù)值模擬方法。除此以外,三角級(jí)數(shù)法也是采用較多的一種軌道隨機(jī)不平順時(shí)域樣本數(shù)值模擬方法[3]。眾所周知,隨機(jī)分析計(jì)算效率隨著所需樣本數(shù)增加而增加,為了更加高效的研究軌道隨機(jī)不平順激勵(lì)下的車輛-軌道耦合系統(tǒng)隨機(jī)振動(dòng)特性,需要用更少的隨機(jī)樣本來反映軌道隨機(jī)不平順的統(tǒng)計(jì)特性。無論是常用的三角級(jí)數(shù)法還是傅里葉逆變換法,在模擬軌道隨機(jī)不平順時(shí)都需要大量的隨機(jī)變量,而分析所需的隨機(jī)樣本數(shù)與隨機(jī)變量數(shù)正相關(guān)。因而用更少的隨機(jī)變量來描述軌道隨機(jī)不平順,可極大提高的車輛-軌道耦合系統(tǒng)隨機(jī)振動(dòng)分析的效率及難度。
文獻(xiàn)[4]中,通過隨機(jī)諧和函數(shù)法生成軌道隨機(jī)不平順時(shí)域樣本,并通過頻率與相位分別隨機(jī)以減少所需的隨機(jī)變量數(shù),其中隨機(jī)頻率數(shù)為50個(gè),考慮軌道高低、方向、水平、軌距不平順情況下,各種不平順相互獨(dú)立,隨機(jī)頻率和隨機(jī)相位共計(jì)需400個(gè)隨機(jī)變量,所需隨機(jī)變量依然很多,且單個(gè)隨機(jī)樣本中頻率分量減少必然會(huì)導(dǎo)致模擬的軌道不平順與實(shí)際偏差增大,并不能很好地反映不平順的時(shí)域特征,該方法可用于線性車輛-軌道耦合系統(tǒng)中軌道隨機(jī)不平順的模擬,但不適用于非線性車輛-軌道耦合動(dòng)力學(xué)系統(tǒng)。
本文基于譜表示方法,將譜表示方法中的標(biāo)準(zhǔn)正交隨機(jī)變量表達(dá)為基本隨機(jī)變量的正交函數(shù)形式,從而實(shí)現(xiàn)用較少隨機(jī)變量來描述軌道隨機(jī)不平順的目的。減少隨機(jī)變量數(shù)的同時(shí),保留更多的頻率及相位分量,可更準(zhǔn)確反映軌道隨機(jī)不平順統(tǒng)計(jì)特性以及時(shí)域特征,提高非線性車輛-軌道耦合系統(tǒng)隨機(jī)振動(dòng)計(jì)算效率。
一維、單變量、功率譜密度函數(shù)為S(ω)的隨機(jī)過程η(t)可以表示為[5]
(1)
式中:Ut(ω)和Vt(ω)為隨機(jī)過程η(t)的譜過程,其增量dUt(ω)和dVt(ω)滿足隨機(jī)過程譜表示的基本條件,即
其中,E[·]表示數(shù)學(xué)期望;ω、ω′為隨機(jī)過程中的頻率;S(·)為雙邊功率譜密度函數(shù)。
對(duì)式(1)進(jìn)行離散,并引入一組標(biāo)準(zhǔn)正交隨機(jī)變量{Xk,Yk},當(dāng)頻率間隔Δω足夠小時(shí),式(1)可表述為
(2)
式中:dUt(ω)≈ΔUt(ω)=σkXk; dVt(ω)≈ΔVt(ω)=σkYk;σk=[2S(ωk)Δω]1/2,S為軌道不平順功率譜密度函數(shù);Nω為頻率離散點(diǎn)數(shù);頻率分量均勻離散情況下,ωk=ωl+k·Δω,Δω=(ωu-ωl)/Nω,ωu、ωl為頻率離散點(diǎn)上、下限值。
(3)
式中:cas(x)=cos(x)+sin(x)為Hartley正交基函數(shù);隨機(jī)變量Θ1、Θ2相互獨(dú)立且服從區(qū)間[0,2π)上均勻分布。直接生成的標(biāo)準(zhǔn)正交變量并不能直接代入式(2),還需要構(gòu)造隨機(jī)映射將標(biāo)準(zhǔn)正交變量映射到標(biāo)準(zhǔn)正交隨機(jī)變量{Xk,Yk},由此將描述隨機(jī)過程的隨機(jī)變量數(shù)降低為2。
以TB/T 3352—2014《高速鐵路無砟軌道不平順譜》[6]為例,軌道不平順譜包括高低、方向、軌距和水平四種軌道不平順譜,計(jì)算不平順譜時(shí),將實(shí)測(cè)得到的軌道不平順按1 024 m劃分單元,剔除異常值和趨勢(shì)項(xiàng),加窗處理后進(jìn)行傅里葉變換(FFT)計(jì)算,對(duì)所有單元的軌道不平順譜進(jìn)行擬合,得到軌道不平順功率譜密度函數(shù)。我國的軌道不平順譜擬合公式為
(4)
式中:f為空間頻率;A、k為擬合系數(shù)。
統(tǒng)計(jì)誤差可表示為
(5)
需要指出的是,軌道不平順功率譜密度函數(shù)是在大量實(shí)測(cè)數(shù)據(jù)基礎(chǔ)上,按一定單元長度組成軌道不平順樣本集,由樣本集計(jì)算得到的均值譜。單個(gè)軌道不平順分段的隨機(jī)誤差達(dá)到了100%,隨著分段次數(shù)的增加隨機(jī)誤差逐漸減小,如果分段100次則隨機(jī)誤差可以降到10%[6]。因而采用單個(gè)軌道隨機(jī)不平順樣本或單元計(jì)算得到的不平順功率譜與現(xiàn)有軌道不平順功率譜對(duì)比,并不能說明數(shù)值模擬方法的準(zhǔn)確性,需采用多個(gè)不平順樣本或分段所組成的樣本集不平順功率譜與現(xiàn)有軌道不平順功率譜對(duì)比驗(yàn)證。
對(duì)于多樣本模擬,常采用蒙特卡羅法選擇隨機(jī)變量,但其計(jì)算效率偏低。擬蒙特卡羅法是通過數(shù)論方法產(chǎn)生高度均勻化的擬隨機(jī)數(shù)序列代替蒙特卡羅隨機(jī)數(shù)模擬,以減少隨機(jī)樣本數(shù)從而提高計(jì)算效率,但當(dāng)隨機(jī)變量較多,即高維度情況下,數(shù)論選點(diǎn)法選取的隨機(jī)數(shù)列不可避免的會(huì)產(chǎn)生叢聚現(xiàn)象[7]。采用譜表示-隨機(jī)函數(shù)法模擬生成軌道不平順時(shí)域樣本,可減少所需的隨機(jī)變量數(shù),避免高維數(shù)論選點(diǎn)中叢聚現(xiàn)象對(duì)計(jì)算結(jié)果的影響。
綜上所述,采用譜表示-隨機(jī)函數(shù)法生成軌道隨機(jī)不平順樣本集過程可表示為
Step1根據(jù)生成的軌道隨機(jī)不平順類型確定隨機(jī)變量數(shù)N,當(dāng)需同時(shí)生成軌道高低、方向、軌距、水平不平順時(shí),由于各類型不平順間相互獨(dú)立,因而需要8個(gè)隨機(jī)變量(各類型不平順均為2個(gè))。
Step2根據(jù)所需隨機(jī)變量數(shù),采用擬蒙特卡羅法生成隨機(jī)變量的低偏差序列Θ,常用的低偏差序列有Corput序列、Halton序列、Sobol序列等[7]。本文采用Sobol序列生成隨機(jī)變量[5],采用Matlab直接計(jì)算Sobol的命令。
Θ={θ1j,θ2j,…,θNj}j=1,2,…,Npt
式中:N為隨機(jī)變量數(shù);Npt為樣本集合的樣本數(shù)。此時(shí)生成的隨機(jī)變量序列服從[0,1)上均勻分布, 需將其變換為[0,2π)上均勻分布隨機(jī)變量Θ′
Θ′=2πΘ
(6)
高低不平順ηV(t)為
方向不平順ηA(t)為
水平不平順ηC(t)為
軌距不平順ηG(t)為
式中:SV、SA、SC、SG分別為高低、方向、水平、軌距不平順功率譜密度:ωk為軌道隨機(jī)不平順中不同的頻率分量,其表達(dá)式為
(7)
采用譜表示-隨機(jī)函數(shù)模擬方法生成樣本數(shù)為300的軌道隨機(jī)不平順樣本集,其中軌道不平順功率譜計(jì)算方法中單元長1 024 m,采樣間隔0.25 m,單個(gè)單元共計(jì)4 096個(gè)點(diǎn),因而模擬軌道隨機(jī)不平順時(shí)頻率離散點(diǎn)數(shù)Nω取4 096。
對(duì)生成的軌道不平順樣本集中各樣本,參考文獻(xiàn)[8]中高速鐵路無砟軌道不平順譜的生成方法及軌道不平順譜頻率分布規(guī)律檢驗(yàn)方法,采用周期圖法計(jì)算得到各樣本的不平順譜值,并求得樣本集不平順譜均值,對(duì)各頻點(diǎn)不平順譜值Sk,j按下式進(jìn)行變換
(8)
采用Kolmogorov-Smirnov檢驗(yàn)方法,對(duì)λk,j進(jìn)行顯著性水平α=0.05、自由度為2的χ2分布檢驗(yàn)。軌道不平順單樣本情況、單樣本功率譜、樣本均值、樣本標(biāo)準(zhǔn)差、部分頻點(diǎn)功率譜概率分布CDF等評(píng)價(jià)指標(biāo)以及樣本集總體功率譜PSD與目標(biāo)譜對(duì)比情況如圖1~圖4所示。
圖1 譜表示-隨機(jī)函數(shù)法高低不平順樣本評(píng)價(jià)指標(biāo)對(duì)比
圖2 譜表示-隨機(jī)函數(shù)法方向不平順樣本評(píng)價(jià)指標(biāo)對(duì)比
圖3 譜表示-隨機(jī)函數(shù)法水平不平順樣本評(píng)價(jià)指標(biāo)對(duì)比
圖4 譜表示-隨機(jī)函數(shù)法軌距不平順樣本評(píng)價(jià)指標(biāo)對(duì)比
從圖1~圖4中可以看出,計(jì)算得到的軌道高低、方向、水平、軌距不平順單個(gè)樣本功率譜密度值與軌道不平順功率譜存在較大的誤差,而軌道不平順樣本集計(jì)算得到的功率譜均值與目標(biāo)譜吻合較好,其中各種不平順樣本總體均值與目標(biāo)譜最大相差分別為0.18、0.10、0.07、0.08 mm,標(biāo)準(zhǔn)差最大偏差分別為4.68%、5.77%、6.97%、6.99%,樣本集功率譜最大偏差分別為14.54%、10.23%、10.07%、9.87%,各不平順98.00%、96.33%、98.00%、98.33%的頻點(diǎn)功率譜概率分布滿足自由度為2的χ2分布??梢钥闯鰳颖炯髦笜?biāo)與目標(biāo)譜吻合良好,且僅需8個(gè)隨機(jī)變量便可同時(shí)模擬生成軌道高低、方向、水平、軌道不平順時(shí)域樣本,同時(shí)還具有大量頻率分量。
目前在工程中應(yīng)用較多的軌道隨機(jī)不平順數(shù)值模擬方法有三角級(jí)數(shù)和傅里葉逆變換法。三角級(jí)數(shù)法模擬軌道隨機(jī)不平順可表示為[2]
(9)
式中:ak為相互獨(dú)立的均值為0、標(biāo)準(zhǔn)差為σk的高斯隨機(jī)變數(shù);φk為相互獨(dú)立的隨機(jī)相位。
傅里葉逆變換法模擬軌道隨機(jī)不平順可表示為[1]
η=ifft(X(k))
(10)
式中:X(k)為時(shí)間序列,其為復(fù)序列,實(shí)部關(guān)于Nr/2偶對(duì)稱,虛部關(guān)于Nr/2奇對(duì)稱。設(shè)模擬總時(shí)間為Ts,則時(shí)域采樣點(diǎn)數(shù)為Ts/ΔT,一般需要在末尾添0以保證采樣點(diǎn)數(shù)為2的整數(shù)次冪,即為Nr。
X(k)表達(dá)式為
k=0,1,…,Nr-1
(11)
式中:ξn為獨(dú)立隨機(jī)相位序列,ξn=cosφn+isinφn=exp(iφn),φn服從[0,2π]上的均勻分布;Δf=1/(NrΔT)。
三角級(jí)數(shù)法中,ak為服從高斯分布的隨機(jī)變量,φk為服從均勻分布的隨機(jī)變量,因而三角級(jí)數(shù)模擬軌道隨機(jī)不平順時(shí)所需隨機(jī)變量數(shù)為2Nω,且存在高斯分布與均勻分布的兩種類型隨機(jī)變量。傅里葉逆變化法中,僅有φn服從正態(tài)分布且相互獨(dú)立,而所需隨機(jī)變量數(shù)與時(shí)域采樣點(diǎn)數(shù)相關(guān),一般為2的整數(shù)次冪,當(dāng)采樣間隔0.25 m,模擬1 km軌道隨機(jī)不平順時(shí)所需隨機(jī)變量數(shù)為4 096個(gè)。三種方法模擬軌道隨機(jī)不平順時(shí)域樣本時(shí),不平順頻率點(diǎn)計(jì)算方法基本一致,采用相同頻率數(shù)時(shí),具有相同的不平順波長成分,因而文中僅對(duì)不平順時(shí)域樣本集的幅值進(jìn)行對(duì)比分析。
采用擬蒙特卡羅法生成低偏差隨機(jī)向量序列,通過三角級(jí)數(shù)法、傅里葉逆變換法模擬得到300個(gè)軌道隨機(jī)不平順時(shí)域樣本。計(jì)算樣本集均值、標(biāo)準(zhǔn)差、功率譜,并檢驗(yàn)樣本集各頻點(diǎn)軌道譜分布,各項(xiàng)對(duì)比情況如圖5和圖6所示。各種方法生成的軌道不平順樣本集不平順幅值范圍情況見表1。
圖5 三角級(jí)數(shù)法高低不平順樣本評(píng)價(jià)指標(biāo)對(duì)比
圖6 傅里葉逆變換法高低不平順樣本評(píng)價(jià)指標(biāo)對(duì)比
表1 樣本集不平順幅值范圍 mm
可以看出,除了傅里葉逆變換法計(jì)算得到高低不平順幅值范圍偏小外,三種方法計(jì)算得到的軌道隨機(jī)不平順幅值范圍大體相當(dāng),最大相差5.51%。三角級(jí)數(shù)法計(jì)算得到的軌道不平順單個(gè)樣本功率譜與目標(biāo)功率譜存在較大誤差,各方向軌道不平順樣本集與目標(biāo)譜均值最大相差分別為0.22、0.09、0.11、0.08 mm,標(biāo)準(zhǔn)差最大偏差分別為5.14%、6.78%、4.95%、7.25%,功率譜密度最大偏差分別為15.82%、13.93%、10.78%、9.65%,各不平順97.33%、96.67%、97.33%、98.67%的頻點(diǎn)功率譜概率分布均滿足自由度為2的χ2分布。傅里葉逆變換法計(jì)算得到的軌道隨機(jī)不平順,無論是單樣本還是樣本集,其功率譜值均能很好滿足目標(biāo)功率譜,最大偏差僅分別為0.30%、0.24%、0.35%、0.13%,但各頻點(diǎn)功率譜概率分布均不滿足文獻(xiàn)[8]中關(guān)于軌道不平順譜分布的論述。除此以外,傅里葉逆變換法各方向軌道不平順樣本集與目標(biāo)譜均值最大相差分別為0.28、0.11、0.13、0.07 mm,標(biāo)準(zhǔn)差最大偏差分別為11.38%、8.53%、8.73%、6.57%,與其他兩種方法模擬結(jié)果基本一致。綜上可以看出,譜表示-隨機(jī)函數(shù)法在模擬軌道不平順時(shí)域樣本時(shí)其幅值范圍基本一致。
采用譜表示-隨機(jī)函數(shù)方法僅需8個(gè)隨機(jī)變量就可以描述四種軌道不平順,這為非線性車輛-軌道-橋梁耦合系統(tǒng)隨機(jī)振動(dòng)提供了極大的方便。文獻(xiàn)[9]中指出,由于涉及輪軌空間動(dòng)態(tài)接觸幾何關(guān)系、輪軌切向蠕滑力和輪軌法向力等多種復(fù)雜非線性求解過程,車輛-軌道耦合系統(tǒng)最有效的隨機(jī)振動(dòng)研究方法是將軌道不平順譜轉(zhuǎn)換為時(shí)域樣本,采用一段足夠長的時(shí)域樣本來代替整個(gè)隨機(jī)過程,即采用蒙特卡洛法選取足夠長的時(shí)域樣本。
分別采用三角級(jí)數(shù)法、傅里葉逆變換法、譜表示-隨機(jī)函數(shù)法模擬得到的軌道隨機(jī)不平順樣本,利用三維非線性車輛-軌道-橋梁耦合動(dòng)力學(xué)模型,計(jì)算得到系統(tǒng)動(dòng)力響應(yīng),對(duì)不同距離上的動(dòng)力響應(yīng)結(jié)果進(jìn)行統(tǒng)計(jì)分析,以輪重減載率、脫軌系數(shù)為例,樣本計(jì)算結(jié)果截口概率分布PDF見圖7。
圖7 不同計(jì)算法的響應(yīng)截口概率分布
從圖7中可以看出,在300個(gè)樣本集的情況下,三種方法計(jì)算得到的動(dòng)力響應(yīng)其截口概率分布基本相同,都具有較大的離散性,速度300 km/h、不平順波長2~120 m時(shí),輪重減載率變化范圍為±0.25,脫軌系數(shù)變化范圍為±0.12。由于文中采用中國高速鐵路無砟軌道譜,平順性較其他軌道譜好,因而輪重減載率與脫軌系數(shù)的計(jì)算結(jié)果偏小。
對(duì)于車輛-軌道耦合系統(tǒng)動(dòng)力響應(yīng),采用三種軌道不平順譜數(shù)值模擬方法計(jì)算得到結(jié)果基本一致,但所需的隨機(jī)變量數(shù)卻存在巨大差異。隨機(jī)變量數(shù)的減少方便了樣本空間概率賦予,同時(shí)也可在一定程度上減少所需的樣本數(shù),為采用廣義概率密度演化理論對(duì)車輛運(yùn)行安全性進(jìn)行評(píng)價(jià)提供了方便。
對(duì)于非線性車輛-軌道-橋梁耦合系統(tǒng)的任一動(dòng)力響應(yīng)ξ(t),車輛運(yùn)行過程中不考慮其他系統(tǒng)隨機(jī)的情況下,其隨機(jī)性來源于軌道隨機(jī)不平順中隨機(jī)參數(shù)Θ的隨機(jī)性,該系統(tǒng)為一個(gè)概率保守隨機(jī)系統(tǒng)。根據(jù)隨機(jī)系統(tǒng)的概率守恒原理,其應(yīng)滿足廣義概率密度演化方程[10]
(12)
研究中常通過與大量蒙特卡洛法(MCM)計(jì)算結(jié)果對(duì)比,對(duì)廣義概率密度演化方法(PDEM)計(jì)算結(jié)果進(jìn)行驗(yàn)證[4,12]。采用三角級(jí)數(shù)法計(jì)算模擬軌道隨機(jī)不平順樣本,結(jié)合MCM法計(jì)算得到3 000個(gè)樣本動(dòng)力響應(yīng)統(tǒng)計(jì)情況,并采用譜表述隨機(jī)函數(shù)法,結(jié)合廣義概率密度演化理論,計(jì)算得到300個(gè)樣本下的動(dòng)力響應(yīng)概率密度演化結(jié)果,兩種方法車輛運(yùn)行距離中間點(diǎn)的截口概率密度對(duì)比情況見圖8,由于三角級(jí)數(shù)法所需的隨機(jī)變量數(shù)較多,因而在3 000個(gè)樣本情況下,動(dòng)力響應(yīng)結(jié)果依然存在一定離散性,但兩種方法計(jì)算得到的截口概率密度基本一致,因而采用譜表示-隨機(jī)函數(shù)法模擬軌道隨機(jī)不平順,結(jié)合廣義概率密度演化理論可極大提高非線性車輛-軌道-橋梁耦合系統(tǒng)隨機(jī)振動(dòng)分析效率,相較于采用三角級(jí)數(shù)等方法模擬軌道隨機(jī)不平順并結(jié)合MCM法進(jìn)行車輛-軌道-橋梁耦合系統(tǒng)隨機(jī)振動(dòng)分析,其計(jì)算效率至少可提高一個(gè)量級(jí)。參考“3σ”準(zhǔn)則,取E±3σ對(duì)應(yīng)值為動(dòng)力響應(yīng)指標(biāo)最值。車輛運(yùn)行過程中輪重減載率及脫軌系數(shù)最大值、等概率密度演化曲線見圖9。由此,譜表示-隨機(jī)函數(shù)法軌道不平順數(shù)值模擬方法結(jié)合廣義概率密度演化理論可對(duì)非線性車輛-軌道-橋梁耦合系統(tǒng)車輛運(yùn)行安全性進(jìn)行評(píng)價(jià)。
圖8 MCM與PDEM截口概率密度對(duì)比
圖9 等概率密度演化曲線
軌道不平順譜是描述軌道隨機(jī)不平順最常用的統(tǒng)計(jì)函數(shù),本文以軌道隨機(jī)不平順時(shí)域樣本數(shù)值模擬方法為研究對(duì)象,采用譜表示-隨機(jī)函數(shù)法分別由2個(gè)隨機(jī)變量模擬生成軌道高低、水平、方向、軌距不平順集。
對(duì)比樣本集與目標(biāo)譜均值、標(biāo)準(zhǔn)差、軌道譜以及各頻點(diǎn)軌道譜概率分布,驗(yàn)證了譜表示-隨機(jī)函數(shù)法模擬軌道隨機(jī)不平順時(shí)域樣本準(zhǔn)確性;與三角級(jí)數(shù)法、傅里葉逆變化法相比,模擬生成軌道不平順幅值范圍基本一致,誤差在5.51%以內(nèi),其中傅里葉逆變化法生成的不平順樣本功率譜能很好地吻合目標(biāo)譜,但各頻點(diǎn)功率譜分布并不滿足與文獻(xiàn)[8]所述頻點(diǎn)功率譜分布;采用譜表示-隨機(jī)函數(shù)法模擬生成軌道隨機(jī)不平順時(shí)域樣本集,結(jié)合廣義概率密度演化理論,可用相對(duì)較少樣本計(jì)算得到車輛運(yùn)行過程中概率演化情況,相較于采用三角級(jí)數(shù)法模擬軌道隨機(jī)不平順,并結(jié)合MCM法進(jìn)行隨機(jī)振動(dòng)分析,計(jì)算效率至少提高了1個(gè)量級(jí),為非線性車輛-軌道-基礎(chǔ)結(jié)構(gòu)耦合系統(tǒng)隨機(jī)振動(dòng)分析中軌道隨機(jī)不平順的模擬提供了新思路。