張 宇 李韶華 任劍瑩

* (石家莊鐵道大學(xué)省部共建交通工程結(jié)構(gòu)力學(xué)行為與系統(tǒng)安全國家重點實驗室,石家莊 050043)

? (石家莊鐵道大學(xué)工程力學(xué)系,石家莊 050043)

** (石家莊鐵道大學(xué)河北省工程力學(xué)基礎(chǔ)學(xué)科研究中心,石家莊 050043)

引言

在橋梁系統(tǒng)眾多力學(xué)問題中,車橋相互作用問題是一個長期的熱點課題[1-6].其核心在于如何準(zhǔn)確地預(yù)測車橋耦合振動行為.當(dāng)車輛在橋梁上運(yùn)行時,車輛和橋梁之間是相互作用和互相影響的,這稱之為車橋耦合振動.車橋耦合振動行為的研究可為橋梁結(jié)構(gòu)的設(shè)計、運(yùn)營、維護(hù)和管理提供必要的理論基礎(chǔ)、分析方法以及評估手段,具有非常重要的工程應(yīng)用價值[3].關(guān)于車橋耦合系統(tǒng)動力學(xué)問題的研究,早期研究將移動的車輛對橋梁施加的力建模為移動載荷[7-8].在這種模型中,移動的車輛和橋梁之間不存在動態(tài)耦合關(guān)系,車輛載荷只是移動在橋面上的外部激勵.然而,當(dāng)車輛模型中考慮慣性力時,橋梁和車輛在接觸點處的響應(yīng)是耦合的,即構(gòu)成車橋耦合系統(tǒng)[9].且由于車輛位置隨時間的不斷改變,描述車橋耦合系統(tǒng)的動力學(xué)方程組具有時變特性[10-11].通常情況下,即使是線性時變系統(tǒng)也很難得到精確解,因此各種近似方法被應(yīng)用于獲得車橋耦合系統(tǒng)的動力學(xué)響應(yīng).

車橋耦合系統(tǒng)通常采用有限元法[12-16]或模態(tài)疊加法[17-19]進(jìn)行建模和分析.通過相應(yīng)方法得到時變常微分動力學(xué)方程以后,進(jìn)一步采用直接積分法進(jìn)行求解,主要有線性加速度法、Wilson 法和Newmark法等[20].Zhu 等[21]采用精細(xì)積分法計算了移動載荷通過連續(xù)梁橋時的動力學(xué)響應(yīng).Zhai[22]基于Newmark法發(fā)展了兩種高效的顯式數(shù)值積分方法,提高了大規(guī)模計算時的效率和精度.Stoura 等[23]提出了一種求解車橋耦合動力學(xué)模型的動態(tài)劃分方法,該方法如結(jié)合適當(dāng)?shù)臄?shù)值格式可顯著降低計算量.處理各種時變系統(tǒng)動力學(xué)問題時,Newmark 法[10]和Runge-Kutta 法[24]是應(yīng)用較多的數(shù)值方法.此外,控制問題中經(jīng)典的時間參數(shù)凍結(jié)技術(shù)也被應(yīng)用于慢時變系統(tǒng)動力學(xué)問題[25].最近,Ge 等[26]針對時變動力學(xué)系統(tǒng)的計算問題,提出了求解車橋耦合系統(tǒng)受迫振動的時間參數(shù)凍結(jié)方法,該方法可較好地處理線性時變車橋耦合問題.

在車橋耦合問題的計算中,為了數(shù)值計算的準(zhǔn)確性,在某些工況下需考慮非線性因素,如橋梁跨度大,當(dāng)載荷作用下產(chǎn)生大變形時的幾何非線性[27].當(dāng)考慮非線性因素時,描述車橋相互作用問題的動力學(xué)方程組具有典型的時變、非線性特性,對其高效、精確的求解具有重要的意義.

1994年,鐘萬勰[28]針對線性常系數(shù)微分方程問題的計算,提出了精細(xì)積分法,該方法在求解指數(shù)矩陣時用矩陣加法代替乘法,巧妙地避免了舍入誤差.由于其優(yōu)異的數(shù)值特性,近年來精細(xì)積分法得到了廣泛的研究和應(yīng)用[29-32],并發(fā)展出了增維精細(xì)積分法[33]、廣義精細(xì)積分法[34]、辛精細(xì)積分法[35-36]和快速精細(xì)積分法[37]等方法.針對半線性微分方程問題,近年來發(fā)展了多種指數(shù)積分方法及其相應(yīng)的計算格式[38],在處理高振蕩、剛性問題中得到了良好的應(yīng)用[39-40].鄧子辰等[41]結(jié)合精細(xì)積分法和指數(shù)積分法,提出了精細(xì)指數(shù)積分法,并將其應(yīng)用于衛(wèi)星編隊飛行動力學(xué)問題中弱非線性方程的求解,但并未考慮時變問題.

車橋耦合動力學(xué)問題解析求解難度大,通常采用數(shù)值積分方法進(jìn)行求解,如經(jīng)典的Newmarkβ法、Runge-Kutta 法.但是這兩種方法受穩(wěn)定性影響較大,數(shù)值積分步長需選取較小,這會增加計算時間,降低計算效率.且這些算法本身具有一定的數(shù)值耗散性,在長時間數(shù)值積分后可能會導(dǎo)致計算結(jié)果的失效.因此,本文針對上述問題,將精細(xì)指數(shù)積分法與時間參數(shù)凍結(jié)技術(shù)相結(jié)合,提出一種新的計算方法——參數(shù)凍結(jié)精細(xì)指數(shù)積分法,并將其拓展應(yīng)用于車橋耦合系統(tǒng)動力學(xué)問題的求解中.充分發(fā)揮不同方法的特點和優(yōu)勢,使其能夠高效、精確地求解具有時變、非線性特性的車橋耦合動力學(xué)方程組,預(yù)期可為車橋耦合問題的動力學(xué)仿真計算提供一種新的方法.

1 考慮橋面鋪裝層的車橋耦合模型

考慮汽車通過跨度為L的橋梁,建模時橋梁采用包含鋪裝層和橋體的層合梁模型,汽車采用車橋耦合振動問題建模中常用的兩自由度1/4 汽車懸架模型[10,24,42].車橋耦合模型示意圖如圖1 所示.

圖1 車?橋面鋪裝層?橋耦合模型Fig.1 Vehicle-pavement-bridge coupled model

考慮汽車自橋梁左端向右以勻速v行駛,1/4 車模型有兩個描述垂向振動的動力學(xué)方程[10]

其中m1和m2分別為非簧載質(zhì)量和簧載質(zhì)量,z1和z2分別為非簧載質(zhì)量位移和簧載質(zhì)量位移,k1和k2分別為輪胎和懸架剛度系數(shù),c1和c2分別為輪胎和懸架阻尼系數(shù),w1為車橋耦合振動引起的車輪與橋面接觸點處橋面的二次位移激勵值.

車輪與橋面接觸的動態(tài)輪胎力可表示為

基于達(dá)朗貝爾原理,考慮圖1 中的具體層合梁模型,梁的振動方程可給出如下

其中 ρ為梁材料等效密度,ρ1和 ρ2分別為橋面鋪裝層和橋體的密度,ht,hm和hb分別表示鋪裝層上表面、鋪裝層與橋體連接面和橋體下表面的位置坐標(biāo),A為梁橫截面積,w表示梁的橫向位移場,M為彎矩,σx為梁軸向應(yīng)力,F(x,t)為梁橫向外載荷,符號(.)表示對時間t求導(dǎo),符號(')表示對空間x求導(dǎo).

當(dāng)式(6)做積分運(yùn)算時,z值0 點在層合梁中性軸處,即圖1 中x軸位置.因此,需根據(jù)橋面鋪裝層和橋體具體參數(shù)確定ht,hm和hb的坐標(biāo).假設(shè)鋪裝層厚度為h1,橋體厚度為h2,橋面鋪裝層和橋體的彈性模量比值為

其中E1和E2分別表示橋面鋪裝層和橋體的彈性模量.中性軸距離上邊緣距離hpt和下邊緣的距離hpb以及ht,hm,hb的坐標(biāo)分別為

橋梁結(jié)構(gòu)自身具有一定的對外界能量耗散和減振能力,其對振動能量的耗散能力可用黏性阻尼系數(shù)來體現(xiàn)[43].本文橋梁建?？紤]了Kelvin 黏彈性理論以及幾何非線性特性,梁軸向應(yīng)力?應(yīng)變關(guān)系和應(yīng)變?位移關(guān)系可給出如下

其中 εx為梁軸向應(yīng)變,η為黏性阻尼系數(shù).

將式(3)、式(6)、式(11)和式(12)代入式(4)可以得到

其中 δ為狄拉克函數(shù),滿足關(guān)系

b為梁寬度,η1,η2分別表示鋪裝層和橋體的黏性阻尼系數(shù).式(13)中如不考慮梁幾何非線性,則 ξ2和 ξ3等于0,如不考慮梁的黏彈性,則 ξ4等于0.

采用伽遼金法可將偏微分方程(13)進(jìn)行簡化,本文考慮梁兩端為簡支邊界條件,則梁橫向位移函數(shù)可表示為

其中 sin(iπx/L)為梁簡支邊界條件下的第i階模態(tài)函數(shù),qi(t)為梁振動的第i階廣義位移.對于受移動載荷作用的簡支梁橋,其一階模態(tài)函數(shù)對橫向振動起決定性作用[8].因此,本文選取了橋梁一階模態(tài)進(jìn)行研究,即N= 1.應(yīng)用伽遼金法,將式(16)代入式(13),乘以第1 階模態(tài)函數(shù)并沿梁長積分,則可給出描述梁橫向振動的非線性常微分方程如下

整理式(1)、式(2)和式(17),即為本文中描述車橋耦合振動的動力學(xué)方程組,其中廣義位移向量為X=[z1,z2,q1]T

通過觀察式(19)可以發(fā)現(xiàn),由于系數(shù)中顯含時間t,且由于 ψ3和 ψ4的存在,該動力學(xué)方程組為時變和非線性的,很難得到解析解,通常采用數(shù)值方法進(jìn)行求解.式(19)可整理成矩陣的形式,如下

系數(shù)矩陣M,C,K,F分別為質(zhì)量矩陣、阻尼矩陣、剛度矩陣和載荷列陣.為便于編程計算,式(19)中的非線性項寫入式(20)載荷列陣F中,在每一數(shù)值積分步,可采用上一積分步的計算值近似求解非線性項.式(20)即可采用常規(guī)數(shù)值積分方法進(jìn)行求解,如經(jīng)典的Newmark-β法.

2 參數(shù)凍結(jié)精細(xì)指數(shù)積分方法的應(yīng)用

將方程組(19)降階為一階動力學(xué)方程組,可表示為如下形式

如能精確求解矩陣A,即可由式(25)計算出每一離散時間步的數(shù)值解uk.對于矩陣A的求解,可采用精細(xì)積分法[28],其思路是利用指數(shù)函數(shù)的加法原理,將時間步長分為 2n等分,如令 Δt=τ/2n,則有[41]

由于 Δt=τ/2n是一個極小的時間段,如對exp(LCΔt)進(jìn)行泰勒展開,取其前5 項就能達(dá)到足夠的精度,即

其中,I為單位矩陣,A1是一個很小的矩陣.為避免計算過程中的舍入誤差,由式(26)和式(27)有

如定義矩陣序列

則由式(28)可得遞推關(guān)系如下

由式(29)求得矩陣序列An+1的值,代入式(30),與I陣相加即可求得矩陣A的值.此過程中,利用了指數(shù)函數(shù)的性質(zhì),巧妙地運(yùn)用加法代替乘法,避免了嚴(yán)重的舍入誤差,使算法具有很高的精度,通常取n=20即可保證良好的數(shù)值性能.且由于在計算中的每一時間步 τ內(nèi)插入了 2n個點,即使步長 τ值取的較大,也不影響計算的精確性[44].

精細(xì)積分法在計算一般的常系數(shù)齊次線性微分方程時,能夠達(dá)到極高的精度.但車橋耦合問題的動力學(xué)方程組系數(shù)矩陣具有時變性,即式(22)中L(t)為時變系數(shù)矩陣.此時可以采用時間參數(shù)凍結(jié)的思想,基于時間參數(shù)凍結(jié)技術(shù)的基本概念,在每一個離散步長計算間隔內(nèi),假設(shè)L(t)為常數(shù)矩陣,則系統(tǒng)定常,即可應(yīng)用精細(xì)積分法求得每一時間區(qū)間內(nèi)的響應(yīng).如式(25),在計算uk+1的數(shù)值結(jié)果時,依據(jù)平均的思想,可假設(shè)

式(31)中由于在每一離散步內(nèi)引入了近似,因此,當(dāng)時間間隔 τ越小時,近似解越接近于精確解.

進(jìn)一步考慮非線性部分的處理,式(22)中對時間離散后,解的形式可表示為

其中,右端第1 項為線性部分,即采用上述參數(shù)凍結(jié)后的精細(xì)積分法進(jìn)行求解,右端第2 部分積分項為杜哈梅積分.對式(32)采用不同的方法求解,可得到不同的指數(shù)積分法,應(yīng)用較多的是指數(shù)Runge-Kutta法,本文即采用指數(shù)Runge-Kutta 法進(jìn)行計算.s級的指數(shù)Runge-Kutta 法可表示為[41]

其中,aij(LCτ)和bi(LCτ)為待定矩陣,ci為待定常數(shù).aij(LCτ)和bi(LCτ)可以根據(jù)積分因子法按照如下的規(guī)律給出[41]

其中,系數(shù)aij,bi,ci與經(jīng)典Runge-Kutta 法中系數(shù)相同,即每一種經(jīng)典Runge-Kutta 法都對應(yīng)有一種指數(shù)積分法.式(33)中關(guān)于指數(shù)矩陣的求解均可由精細(xì)積分法進(jìn)行精確計算,這種將精細(xì)積分法和指數(shù)積分法結(jié)合得到的方法稱為精細(xì)指數(shù)積分法[41].本文基于精細(xì)指數(shù)積分法與時間參數(shù)凍結(jié)技術(shù),創(chuàng)新性地將二者結(jié)合到一起,提出參數(shù)凍結(jié)精細(xì)指數(shù)積分法,并將其應(yīng)用于車橋耦合動力學(xué)問題的求解中.參數(shù)凍結(jié)精細(xì)指數(shù)積分法的計算重點在于每一迭代步中對時變矩陣L(t)的近似處理和對指數(shù)矩陣A的精細(xì)積分求解,其數(shù)值計算格式容易構(gòu)造,且可以方便地處理時變、非線性問題.

3 算例分析

如前文所述,每一種Runge-Kutta 法都可構(gòu)造相對應(yīng)的指數(shù)積分法,本文基于簡單的四級四階顯式Runge-Kutta 法(RK4-4)[41]構(gòu)造了相應(yīng)的參數(shù)凍結(jié)精細(xì)指數(shù)積分格式(FPEI4-4).

首先,為驗證本文車橋耦合模型及理論計算的有效性,基于車橋耦合實驗平臺,進(jìn)行了實驗測試.用于測試的實驗裝置參數(shù)為: 跨度 3 m的兩端簡支實驗橋,采用空心鋁板制作,其密度為 2700 kg/m3,橫截面積為 9.24×10?4m2,彈性模量為 7 0 GPa,慣性矩為 2.67×10?7m4; 實驗小車質(zhì)量為 2 kg,輪間距為9.5×10?2m,車速為 0.2 m/s.實驗時在橋跨中布置位移傳感器,利用ARTDAQ 信號采集顯示軟件進(jìn)行數(shù)據(jù)采集,采樣頻率為1000 Hz.實驗測試模型及測試現(xiàn)場如圖2 和圖3 所示.圖2 中,1 為西門子V90 伺服電機(jī),2 為傳動齒輪,3 為傳動鏈條,4 為擋板,5 為小車,6 為實驗橋支座,7 為橋板,A 和C 為輔助梁,B 為實驗梁.

圖2 實驗測試模型示意圖Fig.2 Schematic diagram of experimental test model

圖3 實驗測試現(xiàn)場圖Fig.3 Diagram of experimental test site

實驗測試結(jié)果(test result,TR)與基于1/4 車?橋耦合模型的FPEI4-4 格式數(shù)值計算結(jié)果進(jìn)行了對比,FPEI4-4 格式選取了時間步長為 τ=1.0×10?3s,結(jié)果如圖4 所示.

圖4 數(shù)值計算和實驗測試的跨中撓度結(jié)果Fig.4 The numerical and test results of mid-span deflection

從圖4 中可以看出,實驗測試和數(shù)值計算結(jié)果的演化趨勢基本吻合,撓度峰值誤差很小.圖4 中的結(jié)果對比能夠說明本文理論建模和數(shù)值計算的有效性,其具有實際應(yīng)用價值.此外,圖4 中數(shù)值計算與實驗測試結(jié)果誤差主要來源于車輛簡化模型的誤差,以及實驗裝置的制作誤差、外界環(huán)境因素等.如采用更加精細(xì)化的車輛模型,可減小理論計算和實驗結(jié)果的誤差,但即使是選用簡單的1/4 車輛模型進(jìn)行理論分析,其依然能夠較為精確地預(yù)測撓度峰值結(jié)果.

理論結(jié)果計算的精確性主要取決于所建立的車橋耦合動力學(xué)模型的精細(xì)化程度,本文所提算法主要關(guān)注于對已建立的動力學(xué)微分方程的高效求解.后文應(yīng)用所構(gòu)造的FPEI4-4 格式與經(jīng)典的Newmark-β算法進(jìn)行了數(shù)值計算,對比了2 種算法在求解車橋耦合動力學(xué)方程組(19)時的數(shù)值特性.

數(shù)值模擬中采用了重型載貨汽車模型,具體車輛參數(shù)見表1[45],表2 給出了橋的參數(shù),包括橋面瀝青鋪裝層參數(shù)和橋體參數(shù).

表1 汽車參數(shù)Table 1 Vehicle parameters

表2 橋面瀝青鋪裝層和橋體參數(shù)Table 2 Pavement and bridge parameters

由于文中車橋耦合動力學(xué)方程組具有時變和非線性的特性,難以得到解析解,因此選用了小步長下高精度的二級4 階隱式辛Runge-Kutta 算法(SRK2-4)的計算結(jié)果作為參考.

s級隱式Runge-Kutta 法的一般格式為[46]

式(35)中當(dāng)其系數(shù)滿足以下條件時為辛Runge-Kutta 法

當(dāng)系數(shù)bij和ci取不同值時,可得到不同的辛Runge-Kutta 格式.比較常用的SRK2-4 積分格式,其系數(shù)表示如下

同時,為了能夠說明針對動力學(xué)微分方程模型數(shù)值計算結(jié)果的有效性,構(gòu)造了本文模型中橋梁位移的近似解析解.基于文獻(xiàn)[47]中的思路,式(19)中第3 個方程可簡化為

由于本文工況中廣義位移q1值和廣義速度1值較小,非線性項對振動的影響較弱,若略去兩項非線性項,并引入?yún)?shù)

則可得到橋梁廣義位移q1的近似解析解(approximate analytical solution,AAS)為

圖5 給出了車速分別為10 和50 m/s 時,近似解析解、FPEI4-4 格式和SRK2-4 格式計算的橋梁跨中撓度的結(jié)果,其中計算時FPEI4-4 格式選取了大時間步長 τ=1.0×10?2s、SRK2-4 格式選取了小時間步長 τ=1.0×10?4s.

圖5 近似解析解、時間步長為 τ=1.0×10?2 s 的FPEI4-4 格式和時間步長為 τ=1.0×10?4 s 的SRK2-4格式計算跨中撓度結(jié)果Fig.5 The results of mid-span deflection calculated by analytical approximate solution,FPEI4-4 scheme with τ=1.0×10?2 s and SRK2-4 scheme with τ=1.0×10?4 s

從圖5 中可以看出,每種工況下3 條曲線的演化趨勢相同,車速較小時(v= 10 m/s)曲線吻合得更好,這充分體現(xiàn)了數(shù)值計算結(jié)果的有效性.此外,從局部放大圖中可以看出,FPEI4-4 格式和SRK2-4 格式的計算結(jié)果始終幾乎重合,而與近似解析解有一定的偏差.這是因為近似解析解忽略了耦合因素和非線性因素,而數(shù)值格式是直接基于原微分方程(19)計算的原因.

從圖5 中的對比也可以體現(xiàn)出,當(dāng)車速較大時,車橋耦合和非線性因素對橋梁的振動特性影響更為明顯.由于所構(gòu)造的近似解析解忽略了部分因素,不能嚴(yán)格地反應(yīng)振動響應(yīng)的演化趨勢,因此后文為更好地體現(xiàn)所提出方法的數(shù)值特性,選用了SRK2-4格式的數(shù)值計算結(jié)果進(jìn)行參照對比.

假設(shè)車速分別為10 和20 m/s,利用SRK2-4 格式計算了文中車橋耦合動力學(xué)方程組.圖6 中給出了當(dāng)選取時間步長分別為 τ=1.0×10?4和 1.0×10?5時,SRK2-4 格式計算的橋梁跨中撓度結(jié)果之差.

圖6 時間步長為 τ=1.0×10?4 s 和 τ=1.0×10?5 s 時SRK2-4 格式計算跨中撓度結(jié)果之差Fig.6 The difference between the results of mid-span deflection calculated by SRK2-4 scheme under time step τ=1.0×10?4 s and τ=1.0×10?5 s

在數(shù)值積分運(yùn)算中,當(dāng)選取數(shù)值積分步長越小時,數(shù)值計算結(jié)果越收斂于解析解.從圖6 中可以看出,選取較小時間步長為 τ=1.0×10?4和 1.0×10?5時,兩種車速工況下的數(shù)值結(jié)果差值已經(jīng)非常小,達(dá)到了10?17m量級,此時已非常接近解析解.因此,可將選取小步長下的SRK2-4 格式計算結(jié)果作為參照解.由于SRK2-4 格式為隱式格式,為提高計算效率,本文選用了積分時間步長為 τ=1.0×10?4s時的SRK2-4 格式進(jìn)行數(shù)值計算.后文中的誤差均表示不同時間步長下相應(yīng)算法的計算結(jié)果與時間步長為τ=1.0×10?4s時SRK2-4 格式計算結(jié)果的差值.

3.1 不同數(shù)值積分時間步長時的計算結(jié)果對比

對比了選取不同時間步長時,應(yīng)用FPEI4-4 格式和Newmark-β算法計算文中車橋耦合模型的精度.選取車速為相對較快的20 m/s 過橋,以時間步長為 τ=1.0×10?4s時的SRK2-4 格式計算結(jié)果作為參照解,應(yīng)用本文構(gòu)造的FPEI4-4 格式和經(jīng)典的Newmark-β算法得到了橋梁跨中撓度計算結(jié)果與參照解的誤差,如圖7 所示.

圖7 FPEI4-4 格式和Newmark-β 算法求解跨中撓度計算誤差Fig.7 The numerical errors of mid-span deflection calculated by FPEI4-4 scheme and Newmark-β algorithm

從圖7 可以看出,在選取時間步長為 τ=1.0×10?3s時,FPEI4-4 格式與參照解的誤差在10?11m量級,在選取時間步長為 τ=1.0×10?4s時,與參照解的誤差達(dá)到了10?13m 量級.作為一種顯式計算格式,相比于相同數(shù)值積分步長時的經(jīng)典Newmarkβ算法,本文所構(gòu)造的FPEI4-4 格式具有明顯更小的誤差、更高的計算精度.

當(dāng)進(jìn)一步增大時間步長,取時間步長為τ=5.0×10?3s時,如圖8 所示,Newmark-β算法計算結(jié)果已出現(xiàn)明顯的數(shù)值發(fā)散現(xiàn)象,導(dǎo)致計算失敗.如繼續(xù)增大時間步長,Newmark-β算法計算結(jié)果會直接發(fā)散.而從圖8 中還可以看出,即使將時間步長增加到很大的 τ=1.0×10?1s時,FPEI4-4 格式同樣能給出較好的計算結(jié)果,這體現(xiàn)了該計算格式優(yōu)異的數(shù)值特性.由于精細(xì)積分法自身的特性,FPEI4-4 格式即使選取大時間步長,也不會出現(xiàn)數(shù)值發(fā)散現(xiàn)象,而大步長下的數(shù)值誤差主要來源于對系數(shù)矩陣L(t)參數(shù)凍結(jié)時引入的誤差.

圖8 FPEI4-4 格式與Newmark-β 算法得到的橋梁跨中撓度結(jié)果Fig.8 The results of mid-span deflection calculated by FPEI4-4 scheme and Newmark-β algorithm

表3 進(jìn)一步給出了2 種算法在選取不同時間步長時的計算時間和最大誤差.從表3 中可以看出,Newmark-β算法作為一種經(jīng)典的動力學(xué)方程數(shù)值求解方法,在相同的小時間步長下,其比FPEI4-4 格式的計算速度更快.而FPEI4-4 格式則具有明顯的精度優(yōu)勢.此外可以看出FPEI4-4 格式具有更好的數(shù)值穩(wěn)定性,能夠選取較大的積分時間步長來加快計算速度,且在大步長下仍能保持較好的計算精度,而在大步長情況下Newmark-β算法將直接發(fā)散使得計算失敗.如選取 τ=1.0×10?2s的FPEI4-4 格式與τ=1.0×10?4s的Newmark-β算法對比,FPEI4-4 格式在計算時間(0.031098 s)遠(yuǎn)小于Newmark-β算法(0.869262 s)的情況下,其最大誤差僅約為Newmarkβ算法的1/50,具有更好的效率和精度優(yōu)勢.

表3 Newmark-β 算法和FPEI4-4 格式計算跨中撓度對比Table 3 Comparison of mid-span deflection calculated by Newmark-β algorithm and FPEI4-4 scheme

后文中為了說明本文FPEI4-4 格式的優(yōu)勢,直接選取了大時間步長進(jìn)行計算,即取 τ=1.0×10?2s,在該時間步長下傳統(tǒng)的Newmark-β算法將直接發(fā)散.

3.2 數(shù)值積分的長時間穩(wěn)定性

為了更好地體現(xiàn)本文所構(gòu)造格式的長時間數(shù)值特性,進(jìn)行了長時間的數(shù)值積分運(yùn)算結(jié)果對比.

首先關(guān)注長時間振動解的穩(wěn)定性質(zhì).因車在橋上時間較短,當(dāng)車下橋以后,由式(19)描述的車橋耦合作用將消失,其后橋?qū)⒆鲇谐跷灰坪统跛俣鹊淖杂烧駝?其初位移和初速度由車下橋時橋振動位移和速度給出,其值均不大,靠近相圖(0,0)點.

關(guān)注于橋自由振動階段,此時動力學(xué)方程為

如令x=q1,y=,可得一階微分方程組為

原點(0,0)為滿足式(43)的奇點,是式(43)的零解.

如考慮橋梁黏性阻尼,式(43)的線性近似方程組特征方程為

代入數(shù)值可得特征方程兩個根均具有負(fù)實部,可知零解x=y=0為漸近穩(wěn)定的,此時奇點(0,0)為穩(wěn)定焦點,相圖中的軌線表現(xiàn)為盤旋趨近原點(0,0).

如忽略橋梁黏性阻尼,則線性近似方程組特征方程為

代入數(shù)值可得特征方程的2 個根均為純虛根,此時奇點(0,0)為中心型奇點,零解為穩(wěn)定但非漸近穩(wěn)定的,其相圖中的軌線表現(xiàn)為以(0,0)為中心的圓.

綜上所述,由于初始條件趨近(0,0)點,且(0,0)點為穩(wěn)定奇點,因此長時間積分計算的解是穩(wěn)定的.這里忽略了橋梁黏彈性特性,則汽車在過橋以后,橋梁將做幅值相等的無耗散往復(fù)自由振動.

首先給出了Newmark-β算法在選取數(shù)值積分步長 τ=1.0×10?4s時計算的跨中撓度時程曲線,以及跨中撓度?速度相圖,如圖9 所示.算例中選取車速為20 m/s,汽車下橋后繼續(xù)計算了100 s 的橋梁自由振動,計算耗時為35.948393 s.

圖9 Newmark-β 算法的計算結(jié)果Fig.9 The calculation results of Newmark-β algorithm

從圖9 中可以看出,即使選擇了較小的時間步長( τ=1.0×10?4s),Newmark-β算法在長時間積分時,仍具有明顯的數(shù)值耗散性.圖9(a) 中顯示,Newmark-β算法計算的汽車下橋后橋梁自由振動數(shù)值解具有發(fā)散的特性,長時間數(shù)值積分后會導(dǎo)致計算結(jié)果的失效.從圖9(b)中也可以看出,在自由振動的周期階段,跨中撓度?速度相圖表現(xiàn)出了明顯的發(fā)散現(xiàn)象.

選取相同的工況,圖10 進(jìn)一步給出了顯式FPEI4-4 格式計算的橋梁跨中撓度以及撓度?速度相圖的數(shù)值結(jié)果.數(shù)值積分計算時,選用了大時間步長τ=1.0×10?2s.

圖10 FPEI4-4 格式的計算結(jié)果Fig.10 The calculation results of FPEI4-4 scheme

從圖10 中可以看出,即使選擇了較大的時間步長,本文所提出的顯式FPEI4-4 格式仍具有良好的數(shù)值特性,數(shù)值計算結(jié)果保持了長時間的數(shù)值穩(wěn)定性,無發(fā)散跡象.由于此時零解為穩(wěn)定但非漸近穩(wěn)定的,其相圖中的軌線表現(xiàn)為以(0,0)為中心的圓,如圖10(b)中右側(cè)所示,圖10(b)左側(cè)為車橋耦合初始階段的軌線.由于精細(xì)積分法自身的特點,可以選取大時間步長進(jìn)行數(shù)值積分運(yùn)算,本算例中步長取為τ=1.0×10?2s,計算耗時僅為1.370017 s.文中雖然對時變矩陣L(t)在每一離散步進(jìn)行了參數(shù)凍結(jié)近似,但結(jié)果顯示在選取較大步長時依然具有良好的數(shù)值精度.通過與時間步長 τ=1.0×10?4s時的隱式SRK2-4格式計算結(jié)果對比,自由振動階段的最大誤差穩(wěn)定維持在 1.38×10?6m量級.本文選取的計算模型矩陣維數(shù)較小,如當(dāng)處理高維時變車橋耦合動力學(xué)問題時,選取大步長所帶來的計算效率優(yōu)勢將更加明顯.

3.3 不同時變程度的計算結(jié)果對比

由于汽車為移動載荷,導(dǎo)致車橋耦合動力學(xué)方程組具有時變性,由式(22)可以看出,當(dāng)汽車和橋梁參數(shù)確定的情況下,系數(shù)矩陣時變的快慢和車速v直接相關(guān).本文數(shù)值格式構(gòu)造時,在每一離散步長時間間隔內(nèi)對時變矩陣L(t)進(jìn)行了參數(shù)凍結(jié)處理,將時變系數(shù)矩陣近似為每一離散步內(nèi)的常數(shù)矩陣.這種近似帶來的誤差和矩陣時變的快慢以及步長的選取有關(guān),因此進(jìn)一步分析了系數(shù)矩陣具有不同時變程度時(不同車速)所構(gòu)造的FPEI4-4 格式的數(shù)值特性.

圖11 給出了選取大時間步長為 τ=1.0×10?2s時,不同車速下FPEI4-4 格式求解的橋梁跨中撓度與參照解的計算誤差.從圖11 中可以明顯看出,車速越快,FPEI4-4 格式計算結(jié)果與參照解的計算結(jié)果差別越大.當(dāng)車速為10 m/s 時,最大差值為 7.754 9×10?10m; 車速為20 m/s 時,最大差值為4.409 3×10?9m; 車速為50 m/s 時,最大差值為 1.432 6×10?8m.由于車速與時變矩陣變化快慢直接相關(guān),車速越大,參數(shù)凍結(jié)引起的誤差就會越大,這符合預(yù)期.但即使在車速為50 m/s 時(時速為180 km/h),誤差仍可保持在10?8m 量級,而此車速下跨中撓度的最大值為1.743 2×10?3m,跨中計算結(jié)果仍是可接受的.因此本文所構(gòu)造的FPEI4-4 格式同樣適于計算快車速、快時變問題,具有實際應(yīng)用性.

圖11 τ=1.0×10?2 s,不同車速時FPEI4-4 格式求解跨中撓度計算誤差Fig.11 The numerical errors of mid-span deflection calculated by FPEI4-4 scheme under different vehicle speeds with time step τ=1.0×10?2 s

3.4 時變系數(shù)矩陣 L(t)凍結(jié)方式的對比

本文應(yīng)用的數(shù)值格式在處理時變問題時,采用了參數(shù)凍結(jié)技術(shù),對每一積分步的時變系數(shù)矩陣進(jìn)行了近似處理,本節(jié)對比分析選取不同的近似方式時對數(shù)值計算結(jié)果的影響.下式列舉了每一離散步對時變矩陣L(t)進(jìn)行參數(shù)凍結(jié)的3 種方式,本文基于平均的思想,應(yīng)用的近似矩陣LC1.在此,對比分析分別選取LC1,LC2和LC3時的數(shù)值特性

圖12 給出了選取不同近似矩陣時,應(yīng)用顯式FPEI4-4 格式計算的橋梁跨中撓度和撓度?速度相圖的數(shù)值結(jié)果,計算時選取車速為20 m/s,積分步長為τ=1.0×10?2s,考慮了橋梁的黏彈性特性.

圖12 τ=1.0×10?2 s,FPEI4-4 格式選取不同時變矩陣凍結(jié)方式的計算結(jié)果Fig.12 The numerical results calculated by FPEI4-4 scheme under different parameter freezing forms with time step τ=1.0×10?2 s

從圖12 中可以看出,由于考慮了橋梁的黏彈性,在汽車下橋以后,跨中位移振動迅速衰減為0,在相圖中曲線逐漸趨于(0,0)點,這和實際情況相符.由于此時零解為漸近穩(wěn)定的,奇點(0,0)為穩(wěn)定焦點,如圖12(b)中右側(cè)所示,相圖中的軌線表現(xiàn)為盤旋趨近原點(0,0),圖12(b)左側(cè)為車橋耦合初始階段的軌線.同時從圖中可以看出,應(yīng)用本文構(gòu)造的顯式FPEI4-4 格式進(jìn)行數(shù)值積分計算時,選取3 種不同參數(shù)凍結(jié)方式時的圖線幾乎吻合,均具有良好的數(shù)值計算效果.雖然在對時變矩陣L(t)進(jìn)行參數(shù)凍結(jié)時引入了誤差,但數(shù)值結(jié)果顯示,即使選擇了較大的步長,依然可以得到較為滿意的結(jié)果,這體現(xiàn)了本文所提出算法在計算車橋耦合振動問題時的優(yōu)勢.

進(jìn)一步分析了采用不同參數(shù)凍結(jié)形式時的計算結(jié)果和積分步長 τ=1.0×10?4s時SRK2-4 格式計算參照解的誤差,結(jié)果如圖13 所示.從圖13 中可以看出,采用LC2和LC3參數(shù)凍結(jié)形式時產(chǎn)生的誤差峰值相近,而采用LC1的時變參數(shù)凍結(jié)方式產(chǎn)生的誤差要明顯小于LC2和LC3,具有更好的數(shù)值計算效果.表4進(jìn)一步給出了FPEI4-4 格式選取不同積分時間步長時,不同參數(shù)凍結(jié)形式與參照解的最大誤差.從表4中同樣可以看出,選取LC1的時變參數(shù)凍結(jié)方式具有更好的數(shù)值精度,而選用LC2和LC3參數(shù)凍結(jié)形式時精度稍差且峰值誤差接近,而此車速下跨中撓度的最大值為 1.395 9×10?3m,不同凍結(jié)方式的跨中計算結(jié)果均是可接受的.基于平均化的思想,在處理不同的車橋耦合動力學(xué)問題時,建議采用LC1的時變參數(shù)凍結(jié)方式.

表4 不同冰凍近似方式的對比Table 4 Comparison of different parameter freezing forms

圖13 τ=1.0×10?2 s,FPEI4-4 格式不同時變矩陣凍結(jié)方式求解跨中撓度的計算誤差Fig.13 The numerical errors of mid-span deflection calculated by FPEI4-4 scheme under different parameter freezing forms with time step τ=1.0×10?2 s

4 結(jié)論

本文基于車橋耦合振動系統(tǒng)中時變、非線性動力學(xué)方程組的求解問題,提出了一種新的數(shù)值計算方法,即時間參數(shù)凍結(jié)的精細(xì)指數(shù)積分方法.基于簡單的四級四階顯式Runge-Kutta 格式(RK4-4)構(gòu)造了相應(yīng)的時間參數(shù)凍結(jié)精細(xì)指數(shù)積分格式(FPEI4-4 格式).以考慮了橋面鋪裝層、橋梁結(jié)構(gòu)幾何非線性和黏彈性特性的車橋耦合問題為例,推導(dǎo)了描述該問題的時變、非線性動力學(xué)方程組,應(yīng)用本文提出的FPEI4-4 格式進(jìn)行了求解.首先,通過與實驗測試結(jié)果進(jìn)行對比,體現(xiàn)了本文理論建模和數(shù)值計算的有效性和實用價值.進(jìn)一步,通過與近似解析解、二級4 階隱式辛Runge-Kutta 格式(SRK2-4)計算的參照解以及經(jīng)典的Newmark-β數(shù)值積分格式計算結(jié)果進(jìn)行對比,體現(xiàn)了本文所提算法的良好數(shù)值特性.文中主要結(jié)論如下:

(1) 所提出的時間參數(shù)凍結(jié)精細(xì)指數(shù)積分方法可方便處理時變、非線性動力學(xué)問題,即使選取較大的積分步長也能得到滿意的結(jié)果,有效地提高了計算效率;

(2) 本文針對車橋耦合動力學(xué)問題提出的參數(shù)凍結(jié)精細(xì)指數(shù)積分法計算速度快,具有長時間的數(shù)值積分穩(wěn)定性,能夠有效地保持動力學(xué)問題的本質(zhì)特性;

(3) 針對快時變問題,對時變系數(shù)矩陣參數(shù)凍結(jié)引起的誤差會增大,但本文格式數(shù)值積分結(jié)果的誤差仍可保持在較小量級,計算結(jié)果仍是可接受的;

(4) 在每一離散步,對于時變系數(shù)矩陣的參數(shù)凍結(jié)方式,建議選用上一積分時刻和本積分時刻矩陣的均值,但即使選擇上一時刻矩陣值或當(dāng)前時刻矩陣值進(jìn)行近似,也能取得較滿意的數(shù)值積分結(jié)果.

本文所提出的參數(shù)凍結(jié)精細(xì)指數(shù)積分法預(yù)期可為時變、非線性的車橋耦合動力學(xué)問題研究提供一種新的計算途徑.