牛留斌,劉金朝,肖炳環(huán),徐曉迪
(中國鐵道科學研究院集團有限公司基礎設施檢測研究所,北京 100081)
軌道的高平順性是車輛安全運行的前提,軌下基礎良好的支承性能是保持軌道空間形位的關鍵。線路上不可避免的軌道不平順是引起輪軌振動的主要激勵源,其波長和幅值是決定輪軌振動頻率和烈度的主要因素。軌道不平順對車輛的影響特征與不平順波長有關,如短波軌道不平順引起輪軌高頻振動和環(huán)境噪聲,長波軌道不平順引起車體低頻振動[1],特別是引起車輛振動頻率接近其固有頻率1 Hz 左右的波長被定義為敏感波長[2]。軌道不平順引起輪軌間附加動態(tài)力對軌道的影響不低于靜輪重[3]。過大的鋼軌動態(tài)位移會加速軌道結(jié)構累積變形不利于保持軌道幾何空間形位,不但增加運營線路的養(yǎng)護維修工作量,還可能會引發(fā)軌道局部結(jié)構破壞,增大車輛運行風險。垂向動載荷作用下鋼軌撓曲變形是軌下基礎服役狀態(tài)的綜合體現(xiàn),鋼軌撓曲位移為線路的質(zhì)量評價、維修管理等提供了參考依據(jù),如我國新建線路驗收規(guī)范[4]中規(guī)定了鋼軌垂向位移的最大允許值和基準值;國外采用高頻數(shù)字相機[5]、位移傳感器[6]等技術手段測量關注區(qū)的鋼軌動態(tài)位移,根據(jù)其變化特征分析軌下基礎支承的劣化狀態(tài)。
鋼軌撓曲位移除了與其承受的垂向載荷有關外,還受到軌道垂向位移導納特性[7]、軌下基礎剛度及車輛運行速度等因素影響,是一個復雜的車輛-軌道-基礎系統(tǒng)動力學問題。Chen Y H 等[8]利用勻速運動坐標系下無限長鐵木辛柯(Timoshen?ko)梁與歐拉-伯努利(Bernouli-Euler)梁在黏彈性地基上引起的鋼軌計算位移差異;Basu D 等[9]基于溫克爾(Winkler)彈性地基無限長連續(xù)Ber?nouli-Euler梁理論的穩(wěn)態(tài)響應解析解,研究移動載荷速度對鋼軌撓曲位移的影響。除了理論解析求解穩(wěn)態(tài)鋼軌撓曲位移外,數(shù)值仿真技術也應用于軌道基礎變化時鋼軌撓曲求解中。如Giannakos 等[10]建立路基-道床-軌枕-鋼軌模型研究無砟線路上因軌道板出現(xiàn)裂紋時引起軌道整體剛度劣化規(guī)律及鋼軌撓曲位移突變特征,理論結(jié)果在現(xiàn)場測試中得到了驗證;Chung W 等[11]建立的有限元模型中將軌道板連結(jié)方式等效為剪切彈簧,研究了移動載荷作用下浮置板軌道的撓度特征;Auersch 等[12]采用輪軌耦合動力學理論建模,研究了不同頻段的道床、支承剛度不平順對鋼軌垂向振動加速度、鋼軌撓曲位移等時頻域變化特點的影響;Kacimi A E等[13]建立三維軌道-車輛有限元模型,數(shù)值模擬了高速列車單車通過鋼軌時輪軌間復雜的振動響應,以及列車通過速度對軌道動態(tài)撓曲的影響;Ang K K 等[14]利用移動單元分析方法建立非均勻黏彈性地基三維多剛體列車模型,研究了移動荷載在速度和軌道不平順的組合工況下軌道剛度突變惡化導致輪跳現(xiàn)象,以及輪跳時動力放大系數(shù)和車體加速度變化情況;Shan Y 等[15]采用多體力學和有限元相結(jié)合的方法構建單自由度模型,研究了25 m·s-1速度條件下低頻軌道幾何不平順激勵鋼軌撓度及車體垂向加速度變化特性;宋小林等[16]通過車輛-軌道耦合動力學模型,仿真計算高速條件下的動態(tài)輪軌力作為輸入條件時,板式無砟軌道鋼軌垂向位移在時空頻域上的分布情況。以上研究中,較多關注了固定載荷或者實測動態(tài)載荷作用下求解鋼軌撓曲變化特性,而將軌道設置成理想平順狀態(tài),軌道不平順激勵的輪軌力折算成固定頻率的簡諧函數(shù)[17],將其引起的撓曲位移耦合在鋼軌總撓曲位移之中。開展軌道高低不平順激勵的鋼軌撓曲位移變化規(guī)律的研究,有利于為我國高鐵線路重點關注波長的精準選擇,為完善軌道不平順的波長管理提供科學依據(jù)。
本文基于我國典型高鐵線路及服役車輛主要參數(shù),建立車輛-軌道耦合動力學仿真模型,在模型中輸入我國高鐵無砟軌道譜反演的軌道高低不平順,研究速度350 km·h-1條件下軌道不平順幅值和波長對鋼軌撓曲位移在時域和頻域上的變化特征,利用非線性最小二乘法[18]擬合幅值和波長引起的鋼軌撓曲最大值之間的函數(shù)表達式;結(jié)合研究結(jié)果及中國高速鐵路無砟軌道譜曲線,給出高鐵波長管理建議。
車輛運行過程中,鋼軌在輪軌垂向力作用下產(chǎn)生彈性撓曲變形,軌道靜態(tài)不平順發(fā)展成軌道動態(tài)不平順,輪軌垂向力FP可分解為靜輪重和動態(tài)附加力2部分[17],即
式中:FP0為靜輪重,kN;FPdyn為動態(tài)附加輪軌垂向力,kN。
與FP0和FPdyn對應,鋼軌的總位移y由2 部分構成,即
式中:ys和yd分別為靜輪重FP0和動態(tài)附加力FPdyn引起的鋼軌撓曲位移,mm。
構建輪軌耦合動力學模型研究yd與軌道靜態(tài)高低不平順之間的關聯(lián)特性,耦合模型由車輛子模型和軌道結(jié)構子模型組成,如圖1所示。圖中:M和Mb分別為車體和轉(zhuǎn)向架的質(zhì)量;J和Jb分別為車體和轉(zhuǎn)向架點頭轉(zhuǎn)動慣量;z,zbi(i=1,2)及zwj(j=1,2,3,4)分別為車體、前后轉(zhuǎn)向架及4 個輪對垂向沉浮位移;θ和θbi分別為車體和前后轉(zhuǎn)向架的俯仰角;Pj為輪軌接觸力;k1,c1,k2及c2分別為車輛一系、二系懸掛剛度和阻尼;L為轉(zhuǎn)向架中心到車體中心的距離;Lb為同一轉(zhuǎn)向架上2 個輪對軸距的一半;ks和cs分別為鋼軌支承扣件系統(tǒng)的等效剛度和阻尼;kH為輪軌接觸剛度。
圖1 輪軌耦合動力學模型示意圖
車輛剛體模型選用CRH380A 型動車組參數(shù)建模,包括1 個車體、2 個構架及4 個輪對,車輛一系、二系懸掛系統(tǒng)用彈簧阻尼元件模擬。鋼軌撓曲位移主要與垂向激勵有關,車輛子模型中包含了車體俯仰角θ,前后構架俯仰角θbi,車體沉浮位移z,前后構架沉浮位移zbi及4 個輪對沉浮位移zwj等共計10個自由度。車輛子模型的動力學方程為
其中,
式中:Mu,Cu及Ku分別為車輛子模型的質(zhì)量、阻尼及剛度矩陣;uu和P分別為車輛子模型的位移向量和載荷向量;g為重力加速度。
軌道子模型包括鋼軌、軌道板及混凝土底座等部件,采用CRTS Ⅲ型軌道結(jié)構和CN60廓形鋼軌參數(shù)建模,有限元網(wǎng)格均采用三維實體單元進行精細化建模,扣件系統(tǒng)等效為支承彈簧,建立軌道結(jié)構子模型簡圖如圖2所示。
圖2 軌道結(jié)構子模型
軌道結(jié)構子模型的動力學方程可表示為
其中,
Pr=(Pr1Pr2Pr3Pr4)
式中:Mr,Cr,Kr及ur分別為軌道結(jié)構子模型質(zhì)量矩陣、阻尼矩陣、剛度矩陣及位移矩陣;Pr為輪軌接觸位置鋼軌受到的輪軌接觸力。
車輛子模型與軌道結(jié)構子模型通過式(3)和式(4)中輪軌接觸力相聯(lián)系,在輪軌接觸位置,車輛與軌道受到的輪軌力是1 對相互作用力,即Prj=-Pj。采用Newmark-β 法顯式求解振動方程式(3),得到車輛子模型位移向量uu;采用Abaqus 軟件和軌道子模型仿真得到式(4)中的位移向量ur。輪軌接觸力Pj由Hertz 非線性彈性接觸理論[19]得到,即
式中:δj為采用罰函數(shù)法計算輪軌接觸點間的相對位移。
某時刻輪對j穿透鋼軌表面的距離為δj,輪軌接觸示意圖如圖3所示。
圖3 輪軌接觸示意圖
采用修改鋼軌橫斷面上全部網(wǎng)格節(jié)點坐標的方式在軌道模型上施加軌道高低不平順,如圖4所示。圖中,鋼軌網(wǎng)格節(jié)點A處的軌道不平順幅值為dzA時,將原節(jié)點坐標A(xA,yA,zA)修正為新的節(jié)點坐標A’(XA’,YA’,ZA’),換算式為
圖4 輪軌動力學軌道模型網(wǎng)格調(diào)整示意圖
式中:θ1為軌道坡角度;h為CN60 廓形鋼軌高度,176 mm。
鋼軌撓曲位移及軌下基礎頻響范圍主要集中在0~40 Hz[20],速度350 km·h-1條件下對應的軌道不平順波長大于2.4 m,因此主要研究波長在1.5~120.0 m 范圍內(nèi)的軌道靜態(tài)高低不平順引起的鋼軌撓曲變化特性(如不特殊說明,下文模型中輸入的軌道高低不平順均為軌道靜態(tài)高低不平順)。軌道高低不平順采用我國高速鐵路無砟軌道譜反演樣本和余弦型模擬波形。余弦型高低不平順作為輸入時,計算區(qū)域包含10個完整周期的波形。我國高速鐵路無砟軌道高低不平順譜計算式[21]為
式中:f為空間頻率,m-1;Sv(f)為軌道高低不平順譜,mm2·m;ψ和B均為分段擬合系數(shù),rad·m-1。
式(7)中擬合系數(shù)見表1。
表1 軌道高低不平順譜擬合系數(shù)
綜上,輪軌耦合動力學模型的數(shù)據(jù)處理步驟如下。
(1)建立軌道結(jié)構子模型,施加位移邊界條件,并在鋼軌上設置軌道高低不平順。
(2)運用Abaqus 軟件隱式模塊計算重力作用下位移場,得到輪軌耦合動力學模型初始時刻t0時的位移向量。
(3)在t1時刻,由式(5)求出軌道不平順條件下的輪軌垂向力向量Pj,時間步長為10-5s。
(4)采用Newmark-β法和式(3),計算車輛子模型在軌道不平順條件下的位移向量uu,1。
(5)利用MATLAB 修改inp 文件中輪軌接觸力位置參數(shù),并提取輪軌接觸處節(jié)點的位移ur,1。
(6)由uu,1,ur,1和式(5)計算得到輪軌垂向力向量Pru,1,重復上述步驟(4)和(5),得到uu,2和ur,2;如果計算結(jié)果滿足式(8),則uu,2,ur,2及prw,2為t1時刻車輛軌道振動位移和接觸力的結(jié)果;如果不滿足式(8),則繼續(xù)上述步驟(4)—步驟(6),直到前后2次計算結(jié)果滿足式(8)。
式中:ε為極小數(shù),取10-6。
(7)重復上述步驟(3)—步驟(6),計算t2,t3,…,結(jié)束時刻tend車輛、輪對、軌道等位移和輪軌接觸力結(jié)果。
(8)依次輸出鋼軌上所有節(jié)點的撓曲位移,統(tǒng)計該工況下的鋼軌撓曲位移等變化規(guī)律。
(9)修改步驟(1)中鋼軌節(jié)點坐標,仿真計算新軌道高低不平順工況下的鋼軌撓曲位移曲線。
(10)匯總各計算工況下鋼軌撓曲位移曲線,提取波形曲線時頻域特征,進行趨勢擬合等分析處理。
本文所建輪軌耦合動力學模型的仿真參數(shù)見表2和表3。
表2 車輛子模型主要參數(shù)
表3 軌道子模型主要參數(shù)
軌道結(jié)構的振動頻段較低,支承彈簧等部件的阻尼對鋼軌撓曲位移影響較?。?2]?;赪inkler彈性地基理論,鋼軌在輪軌垂向力FP作用下?lián)锨灰苰z0解析式[9,23]為其中,
式中:z0為輪軌接觸位置縱坐標,mm;k為軌下基礎剛度,N·m-1;β為鋼軌的剛比系數(shù);E為鋼軌材料的彈性模量;J0為鋼軌截面慣性矩。
通過疊加式(9)的單輪結(jié)果,得到多輪對作用下的鋼軌撓曲位移ys為
式中:I為輪對個數(shù);當z-zI大于5 m 時輪對之間撓曲影響非常小,我國車輛同一轉(zhuǎn)向架輪對間軸距小于3 m,式(10)中輪對個數(shù)n取2。
通過放大靜態(tài)撓曲位移,得到移動載荷作用下鋼軌撓曲位移yd為
式中:αv為速度放大系數(shù),在車輛運行速度為350 km·h-1時α350取值1.5[24]。
車輛以速度350 km·h-1通過未施加高低不平順的軌道時,鋼軌撓曲波形如圖5所示。由圖5可以看出:因扣件的不連續(xù)支承作用引起鋼軌撓曲位移在靜輪重FP0引起的ys附近周期性波動,波動幅度為0.005 mm,波動的長度等于支承彈簧的間距。
圖5 未施加軌道高低不平順時鋼軌撓曲位移圖
圖5中鋼軌某點D垂向撓曲位移時程曲線及由式(11)計算得到的解析解對比如圖6所示。
圖6 鋼軌撓曲位移理論與仿真計算結(jié)果
由圖6可以看出:前輪對接近D點時,同一轉(zhuǎn)向架前后輪對引起的鋼軌撓曲位移分量逐漸增大并疊加在一起;前輪在D1時刻到達D點,此時鋼軌撓曲位移的最大撓度值為-0.76 mm;前輪對通過D點后,前輪對引起的鋼軌撓曲位移逐步減小,后輪對接近D點,引起的鋼軌撓曲位移在逐步增大,在D2時刻前輪對遠離D點而減小的撓曲位移量分量與后輪對靠近D點而增大的撓曲位移分量相等,此時鋼軌撓曲位移為-0.66 mm;后輪引起的鋼軌撓曲分量逐步增大,在D3時刻后輪對到達D點時,鋼軌撓曲位移達到最大值為-0.96 mm;此后,后輪對遠離D點,前后輪對引起的鋼軌撓曲位移分量均逐漸減小,直至D4時刻鋼軌D點近似回到平衡位置;圖1輪軌耦合動力學模型輸出的鋼軌撓曲位移仿真結(jié)果與理論解析解變化趨勢一致,波形吻合良好。
在鋼軌上施加波長1.5 m 幅值1 mm 軌道高低不平順時仿真得到的輪軌垂向力變化曲線,與文獻[25]在相同參數(shù)下的仿真結(jié)果進行對比,如圖7所示。
圖7 輪軌垂向力仿真波形與文獻結(jié)果對比
由圖7可以看出:2 種仿真計算得到的輪軌垂向力波形吻合較好,峰值相近,驗證了模型仿真結(jié)果的正確性;0~20 m 里程內(nèi)是未施加軌道高低不平順區(qū)段,該區(qū)段輪軌垂向力的小幅度波動是由扣件不連續(xù)支承引起的。
圖6中鋼軌撓曲位移隨低通截止頻率的變化曲線如圖8所示。由圖8可以看出:在截止頻率為500 Hz 時,鋼軌撓曲位移趨于的穩(wěn)定值。因此,輸出結(jié)果濾波均采用500 Hz低通濾波。
圖8 鋼軌垂向位移隨截止頻率變化曲線
我國高鐵無砟軌道譜反演得到的軌道高低不平順樣本如圖9所示。
圖9 軌道高低不平順樣本
為了研究不同幅值軌道高低不平順激勵下鋼軌撓曲位移變化規(guī)律,對其幅值進行放大,放大系數(shù)AMF(Amplitude Multiple Factor) 取2,4,5,8。利用式(6)在鋼軌模型中輸入軌道高低不平順樣本,輪軌垂向力及采用不同放大系數(shù)進行幅值放大后的鋼軌撓曲位移曲線仿真結(jié)果分別如圖10 和圖11所示。
圖10 不同AMF下輪軌垂向力波形曲線
由圖10 和圖11 可以看出:隨著AMF的增大,鋼軌撓曲位移yd的變化趨勢相似,即相位未發(fā)生變化,但波動范圍在逐步增大;AMF為8 時,部分輪軌垂向力為0,說明輪軌間出現(xiàn)脫離,但對鋼軌撓曲位移極大值影響較小。
圖11 不同AMF下鋼軌撓曲位移曲線
不同放大系數(shù)AMF條件下?lián)锨灰婆cAMF=1時的撓曲位移之差反映了幅值變化量對鋼軌撓曲位移的影響,定義幅值增大倍數(shù)AIF(Amplitude In?crease Factor)的值為AMF-1,不同AIF下的鋼軌撓曲變化曲線如圖12所示。
圖12 不同AIF下鋼軌撓曲位移增量
鋼軌撓曲位移曲線波谷位置對應的縱坐標值為鋼軌撓曲變化量,采用AIF=7 鋼軌撓曲位移曲線2階導數(shù)為正值的方法得到12 組曲線波谷位置,如圖中圈選位置所示。過曲線波谷位置做平行于縱坐標軸的點狀豎線,依次定義線1、線2、…、線11。每條線與鋼軌撓曲位移曲線的交點,自上而下依次為AIF=1,AIF=3,AIF=4,AIF=7 時的鋼軌撓曲位置局部最小值,統(tǒng)計線1—線11上波谷位置對應的縱坐標值,結(jié)果見表4。由表4可知:每條線上4個交點(波谷)處對應撓曲位移值的比例關系約為1∶3∶4∶7,與AIF值在對應位置取值的比例關系基本一致。這說明隨著軌道高低不平順幅值增大,鋼軌撓曲位移yd的波形和相位不變,幅值與其最大值近似成正比。
表4 圖12中鋼軌撓曲位移局部小值mm
因此,波長為λ、幅值為Aλ的軌道高低不平順引起的鋼軌撓曲位移yd(Aλ,λ)可記為
波長λ=1.5 m、幅值Aλ=1 mm 余弦型軌道高低不平順條件下鋼軌輪軌垂向力、撓曲位移仿真結(jié)果如圖13所示。
圖13 幅值為1 mm、波長為1.5 m軌道高低不平順時模型仿真結(jié)果
由圖13 可以看出:在軌道高低不平順激勵下輪軌垂向力呈現(xiàn)余弦型變化特征,兩者周期一致,而鋼軌在輪軌垂向力作用下的yd呈現(xiàn)有規(guī)律的波動,不再具有單頻率正余弦變化特征,其頻譜曲線如圖14所示。
圖14 yd頻譜曲線及諧波成分波形
由圖14 可以看出:yd頻譜曲線中包含多個等間隔峰值,yd是由多階周期性諧波成分疊加而成的,基頻為0.666 7 m-1,對應的空間波長為1.5 m;4 階以上諧波的峰值與基頻峰值之比小于10%,1~4次諧波成分所占比重較大。
因此,幅值為1 mm、波長為λ的軌道高低不平順引起的yd(1,λ,z)的函數(shù)方程可用前4 階傅里葉級數(shù)展開式表示為
式中:α0為靜輪重FP0為引起的鋼軌撓曲位移,mm;ak,βk分別為k階諧波成分振幅、諧波與軌道高低不平順之間延遲相位;ωλ為軌道高低不平順激勵角頻率,
由式(13)可知,yd的各諧波分量的振幅a與波長λ有關,波長λ分別為1.5,2,3及4 m條件下的輪軌垂向力如圖15 所示。由圖15 可知:不同波長條件下的輪軌垂向力在靜輪重FP0附近波動范圍依次為±40.9,±22.2,±9.6和±5.5 kN,呈現(xiàn)出輪軌垂向力隨著波長λ增大而減小的趨勢。
圖15 不同波長條件下輪軌垂向力
輪軌垂向力引起的各諧波分量振幅a與其頻率有關,軌道結(jié)構子模型中鋼軌垂向?qū)Ъ{位移特性曲線如圖16所示。圖中:速度350 km·h-1條件1.5~120.0 m 波長的軌道高低不平順激勵角頻率為0.8~65.0 Hz,遠離共振頻率[26]frail且不存在共振突變峰曲線,如深色區(qū)域所示。
圖16 鋼軌垂向?qū)Ъ{位移曲線
由圖16可以看出:隨著軌道高低不平順波長λ增大其引起的激勵角頻率變小,1 kN 的輪軌垂向力引起的鋼軌垂向撓曲位移減小。
由以上分析可知,隨著波長λ的增大和軌道高低不平順激勵的輪軌垂向力減小,輪軌垂向力頻率也變小,引起的鋼軌垂向撓曲位移變小。將式(13)中yd(1,λ,z)的最大值記為α(λ),即α(λ)=max(yd(1,λ,z)),波長λ在1.5~12.0 m、幅值1 mm 時鋼軌撓曲位移最大值α由0.38 mm 下降到0.02 mm,如圖17所示。
圖17 鋼軌撓曲位移最大值
在輪軌耦合動力學模型中輸入波長λ為1.5~120.0 m 幅值1 mm 的軌道高低不平順,仿真計算對應的鋼軌撓曲位移α(λ),并采用非線性最小二乘法[18]對波長λ與α(λ)之間進行擬合,方程為
式中:p1,p2,q1及q2均為待定擬合系數(shù)。
在擬合優(yōu)度大于99.5%時,α(λ)與波長λ之間的擬合曲線及系數(shù)如圖18所示。
圖18 鋼軌撓曲最大值隨波長變化散點及擬合曲線
由式(2)和式(14)可知,波長λ幅值Aλ的軌道高低不平順激勵下鋼軌最大撓曲位移ymax的計算式為
當波長λ大于20 m 時,式(15)中的值較小,可簡化為
式中:比例系數(shù)p1取0.151(見圖18)。
由式(16)可知,軌道高低不平順幅值Aλ與波長λ比值近似與鋼軌撓曲位移最大值成正比。
現(xiàn)行的軌道高低不平順局部峰值[24,27]管理模式未能充分兼顧波長對線路的動態(tài)影響,利用式(15)可以從軌道動態(tài)撓曲位移控制的角度維護線路,精準找出軌道管理重點關注的波長范圍。某速度350 km·h-1高鐵線路軌道高低不平順實測波形,其幅值處于±2 mm 之間,對其進行傅立葉變換,再按式(16)得到的鋼軌最大撓曲位移曲線yd,max,結(jié)果如圖19所示。
由圖19 可以看出:該線路的軌道高低不平順中20~35 m 波長成分的幅值較大,而波長在1.5~10.0 m 范圍內(nèi)的軌道高低不平順成分能夠引起較大的鋼軌撓曲位移。
圖19 某線路實測軌道高低不平順數(shù)據(jù)
(1)構建輪軌耦合動力學模型研究速度350 km·h-1條件下軌道高低不平順引起的鋼軌撓曲位移特性,通過對比鋼軌撓曲位移仿真結(jié)果與理論解析的方式解驗證了模型輸出的準確性。我國高鐵無砟線路軌道譜反演的高低不平順樣本數(shù)據(jù)仿真結(jié)果表明:軌道高低不平順的幅值不影響鋼軌動態(tài)撓曲位移yd響應特征且yd的變化量與軌道高低不平順的幅值變化量近似成正比。
(2)鋼軌垂向位移導納特性隨著波長在1.5~120.0 m 范圍內(nèi)的高低不平順激勵頻率增大而增大。隨波長的增大軌道高低不平順激勵的輪軌垂向力值及其頻率均減小,而鋼軌撓曲位移隨著垂向力值及激勵頻率減小而減小,并由軌道高低不平順引起的多次諧波疊加而成,基波頻率是由波長和車輛運行速度決定,4 階及以上諧波成分所占比重小于10%。
(3)利用非線性最小二乘法及有理式方程擬合了速度350 km·h-1條件下鋼軌撓曲位移與軌道高低不平順幅值、波長之間關系曲線。在波長大于20 m 時,幅值與波長比值近似和最大鋼軌撓曲位移量成正比,比例系數(shù)為0.151。