高速鐵路車-軌-橋耦合系統(tǒng)實時混合試驗的高效計算方法

古泉 張德宇 國巍 武際亮 袁彬 周惠蒙 王濤

(1.廈門大學 建筑與土木工程學院，福建 廈門 361005；2.中南大學 土木工程學院，湖南 長沙 410075；3.中國地震局工程力學研究所 地震工程與工程振動重點實驗室，黑龍江 哈爾濱 150080)

實時混合試驗是在擬動力試驗的基礎上發(fā)展而來的一種試驗方法，又稱實時子結(jié)構(gòu)擬動力試驗[1]。與擬動力試驗相同，實時混合試驗把試驗對象分為兩部分，將結(jié)構(gòu)中難以進行準確數(shù)值模擬的部分作為試驗子結(jié)構(gòu)，其他具有明確力學特性的部分作為數(shù)值子結(jié)構(gòu)[2- 4]。對于實時混合試驗中數(shù)值子結(jié)構(gòu)的計算，國內(nèi)外學者提出了許多積分算法來提高積分的速度、準確性和穩(wěn)定性，例如Newmark族積分算法[5]、CR[6]和KR-α[7]積分算法。而對于試驗中出現(xiàn)的時滯問題，學者們也提出了許多經(jīng)典的時滯補償方法[8- 11]。在高速鐵路系統(tǒng)的實時混合試驗中，可以將列車作為試驗子結(jié)構(gòu)用振動臺加載，橋梁作為數(shù)值子結(jié)構(gòu)用有限元計算。但由于列車的運行速度往往很快，可以達到350 km/h，軌道梁的數(shù)值模型必須足夠長，且列車試驗子結(jié)構(gòu)的控制要求精確的計算，必須采用足夠精細化的車和橋梁模型，這就導致模型規(guī)模非常大，計算極其耗時[12]；同時，實時混合試驗要求實時性，每一步都需要在極短時間內(nèi)完成，對數(shù)值計算速度的要求很高[13]。因此，已有的積分算法方法不能滿足高速鐵路系統(tǒng)實時混合試驗對計算速度和精度的要求。為解決這一難題，文中提出一種高速鐵路車-軌-橋耦合系統(tǒng)實時混合試驗的高效計算方法——移動荷載積分法。在進行實時混合試驗之前將橋梁數(shù)值子結(jié)構(gòu)計算所需要的數(shù)據(jù)進行存儲，在實際計算時只需要將每一時步對應的數(shù)據(jù)進行抽取和計算，即可得到當前時步橋梁數(shù)值子結(jié)構(gòu)的位移響應，從而可以在保證計算精度的同時獲得較高的計算效率。

1 實時混合試驗的移動荷載積分法

1.1 基于線彈性系統(tǒng)的移動荷載積分法

對于線彈性橋梁系統(tǒng)，隨時間連續(xù)移動的輪軌作用力可以按照以下方式疊加，作用于橋梁結(jié)構(gòu)上：

(1)

式中：x為橋梁位置；t為時間；n為計算的總時間步數(shù)；下標k為計算時步；tk表示第k時間步，xk表示當前時刻輪軌力的作用位置，對于以速度v運行的列車，tk與xk的關系為xk=vtk；F(xk，tk)為tk時刻作用于橋梁xk位置處的輪軌作用力，由混合試驗測得；Lk(x，t)為tk時刻作用于橋梁xk位置處的三角形單位脈沖函數(shù)。f(x，t)與Lk(x，t)的函數(shù)圖像如圖1所示。

圖1 f(x，t)與Lk(x，t)的函數(shù)圖像

在線彈性結(jié)構(gòu)中，橋梁在任意(x，t)處的位移可根據(jù)線性疊加原理計算：

(2)

式中：uk(x，t)為三角形單位脈沖Lk(x，t)所引起的位移響應，是位置x與時間t的二元函數(shù)，當k=1時，其函數(shù)圖像如圖2所示；u(x，t)也是位置x與時間t的二元函數(shù)，其函數(shù)圖像如圖3所示。

由式(2)可知，u(x，t)的疊加計算需要先計

圖2 uk=1(x，t)的函數(shù)圖像

圖3 u(x，t)的函數(shù)圖像

算出uk(x，t)并進行存儲，計算時再進行讀取和疊加，所需要的總數(shù)據(jù)量是由計算的總時間步數(shù)決定的。在高速鐵路車-軌-橋耦合系統(tǒng)的實時混合試驗中，橋梁模型具有較大的總長度且時間步長極小，由式(2)計算u(x，t)的過程中需要存儲和輸入極大的數(shù)據(jù)量，每一步的計算耗時過長，無法滿足實時混合試驗的控制要求。

車橋?qū)崟r混合試驗模擬橋上行車的過程，在確定的車速以及時間步長下，每一時刻對應的橋上的輪軌作用力位置都是唯一的。由于在實時混合試驗中每一時間步僅需要計算當前時刻輪軌作用力作用位置處的橋梁響應，而不需要計算橋梁的完整響應，因此對于各個時刻作用于對應位置上的單位脈沖來說，無需存儲完整的結(jié)構(gòu)反應函數(shù)，只需要存儲單位脈沖作用時間后各時刻輪軌力在橋梁作用點處的位移響應。從大量數(shù)據(jù)中抽取這些特定點的數(shù)據(jù)，只存儲最少的有效數(shù)據(jù)(如圖4所示)，并在此基礎上進行疊加計算，其他數(shù)據(jù)即為無效數(shù)據(jù)。

對于線性系統(tǒng)，在i時刻輪軌力作用位置處施加單位脈沖力，在j時刻對應位置處產(chǎn)生的位移響應記為uij(t)。對于動力分析來說，uij為一個矩陣，整個結(jié)構(gòu)儲存位移響應的數(shù)據(jù)結(jié)構(gòu)如式(3)所示。

圖中圓點表示有效數(shù)據(jù)點

在同時存在多個移動荷載的情況下，每個荷載分別采用這個數(shù)據(jù)結(jié)構(gòu)，最終計算時對所有的數(shù)據(jù)進行線性疊加。在進行疊加計算時，只需要提取當前時刻之前作用的單位脈沖的響應數(shù)據(jù)中對應當前時刻的那部分進行疊加，而不需要對所有存儲的響應數(shù)據(jù)進行疊加，可極大提高計算效率。

(3)

以i時刻為例，使用優(yōu)化后的數(shù)據(jù)存儲矩陣，特定時刻橋梁特定位置處的位移響應的計算過程為

(4)

式中，fn為tn時刻作用于橋梁對應位置處的真實輪軌作用力，即混合試驗中所測得的實時輪軌力的大小。

由式(4)可知，為了計算每一時步輪軌力作用位置處的位移，從單位移動脈沖荷載作用下記錄的整體響應中提取所需要位置的響應，形成數(shù)據(jù)矩陣，計算時僅需要提取對應數(shù)據(jù)進行線性疊加。因此，對于n時刻的響應計算，文中算法只需對n個數(shù)據(jù)進行疊加，極大地提升了計算速度。

由上文可知，移動荷載積分法要求在計算前事先存儲好并在計算過程中實時讀取一個維度為n×n的數(shù)據(jù)矩陣(n為計算的總時間步數(shù))。即使這樣，在實時混合試驗中n的取值也可能比較大，因此這一過程會耗時較多。而事實上輪軌力的影響范圍有限，足夠遠處(比如橋梁3～5跨以外)的影響可忽略不計，因此在實際計算中，矩陣U的列數(shù)n取有限值。考慮到移動荷載的有效影響范圍確定，而行數(shù)會一直增加直到橋梁最末端，值得說明的是，該方法適用于橋梁跨度較小的連續(xù)梁橋或簡支梁橋的情況，對于大跨度橋梁并不適用。在混合試驗中，可以分批次讀取數(shù)據(jù)矩陣，每次讀入適當行(比如n行)，減少計算過程中單次讀取的數(shù)據(jù)量，從而減少內(nèi)存占用量并提高計算效率。

上述算法是基于線彈性系統(tǒng)和位移疊加法進行的，下面基于彈塑性數(shù)值子結(jié)構(gòu)方法[14- 15]將移動荷載積分法拓展到局部非線性系統(tǒng)。

1.2 基于彈塑性數(shù)值子結(jié)構(gòu)的移動荷載積分法

對于橋梁局部某些橋墩處可能出現(xiàn)非線性的情況，基于彈塑性數(shù)值子結(jié)構(gòu)方法，橋梁整體仍采用最初的線性模型進行分析，而非線性區(qū)域構(gòu)件則被隔離出來采用非線性模型計算[14]。計算線性和非線性構(gòu)件的內(nèi)力差值，將其作為非線性修正力加到橋梁節(jié)點上，非線性修正力作為外加荷載所引起的位移疊加在線性橋梁位移上，得出真實的非線性橋梁的位移響應，相當于在移動荷載引起的橋梁響應基礎上再疊加一個非移動的荷載導致的響應。值得說明的是，非線性修正力引起的位移響應也可作為移動荷載疊加法的特例，基于式(4)計算得到，只是式(3)中矩陣取法不同，需保持荷載位置不變。

下面簡單介紹數(shù)值子結(jié)構(gòu)方法中非線性修正力的計算方法。

基于有限元的整體結(jié)構(gòu)運動方程為

(5)

將式(5)按Newmark-β法對時間t進行離散:

(6)

式中，α、β為引入?yún)?shù)，F(xiàn)n+1為等效外荷載。

在式(6)的左右兩邊同時加上Kun+1項并整理得

(7)

(8)

文中將非線性材料的本構(gòu)關系直接編寫進程序中，計算非線性構(gòu)件的內(nèi)力，并使用式(8)得到非線性修正力。

2 實時混合試驗中的移動荷載積分法驗證

文中試驗的數(shù)值子結(jié)構(gòu)模型采用如圖5所示的軌道-橋梁系統(tǒng)，總長度為330 m，共11跨，軌道與橋梁的參數(shù)取值見表1。

圖5 軌道-橋梁系統(tǒng)

表1 軌道-橋梁系統(tǒng)的參數(shù)表

試驗中采用如圖6(a)所示的振動臺系統(tǒng)，該振動臺為電液式振動臺，臺面尺寸為1 m×1 m，可承載的最大荷載質(zhì)量為10 t，可加載的有效頻率為0.1～50.0 Hz，最大輸出加速度為1.2g。選取四分之一CRH380A高鐵機車車體(1個輪對)作為試驗子結(jié)構(gòu)。試驗的時間步長為1/128 s，行車速度為230 km/h。為簡化研究問題，對列車1個轉(zhuǎn)向架中的1個輪對單獨劃分進行研究，如圖6(b)所示。試驗子結(jié)構(gòu)的各項參數(shù)見表2。

(a)振動臺系統(tǒng)

(b)簡化后的四分之一車體模型

表2 試驗子結(jié)構(gòu)的參數(shù)1)

1)表中，mc、mt、mw分別為車體、轉(zhuǎn)向架和輪對的質(zhì)量，Ks、Kp分別為車體與轉(zhuǎn)向架、轉(zhuǎn)向架與輪對之間的彈簧剛度，Cs、Cp分別為車體與轉(zhuǎn)向架、轉(zhuǎn)向架與輪對之間的阻尼，msum為總質(zhì)量(含輪對)。

在實時混合試驗中，將軌道-橋梁系統(tǒng)作為數(shù)值子結(jié)構(gòu)，通過移動荷載積分法計算出當前時步的位移響應，車體作為試驗子結(jié)構(gòu)，通過振動臺加載并根據(jù)設備上的傳感器實際測得的加速度計算出當前時步下的輪軌作用力。采用自適應時間序列(Adaptive Time Series)的時滯補償方法[16](簡稱ATS方法)以減小橋梁數(shù)值子結(jié)構(gòu)計算得出的位移命令信號與車輛子結(jié)構(gòu)接收的輸入信號之間的幅值及相位誤差。

文中的實時混合試驗使用Pulsar迭代控制方法(ICS)來進行系統(tǒng)響應的復現(xiàn)，通過Matlab/Simulink平臺實現(xiàn)實時仿真，利用Simulink Real-time實時系統(tǒng)工具箱將Simulink編譯的模型運行在另一臺目標主機(xPC-Target)上，Simulink Real-time支持SCRAMNet接口，使用SCRAMNet將Pulsar系統(tǒng)中的各類信號實時向外發(fā)送，同時將外部的控制指令讀入Pulsar。

2.1 線彈性系統(tǒng)下的移動荷載積分法驗證

本節(jié)驗證移動荷載積分法應用于實時混合試驗中的線彈性系統(tǒng)的可行性與效果。在應用移動荷載積分法進行計算前，需要先計算出車輛行駛前自重荷載所引起的軌道-橋梁系統(tǒng)的初始位移，并將單位脈沖荷載按照實際計算時的時間步依次加載到軌道-橋梁模型上的相應位置，記錄當前加載情況下橋軌所有節(jié)點的全時程位移響應。根據(jù)所記錄的全部加載情況對應的位移響應，選取其中單位脈沖作用時間之后各時刻輪軌力在橋梁作用點處的位移響應，即圖4中的有效數(shù)據(jù)點，形成式(3)中的移動荷載積分法的基本數(shù)據(jù)矩陣。然后根據(jù)式(4)即可計算得到每一時步下軌道-橋梁系統(tǒng)的位移響應。

試驗流程如圖7所示。將第i步移動荷載積分法計算出的位移加載到試驗子結(jié)構(gòu)上，通過振動臺測出的加速度計算出第i+1步的輪軌力，將其加載到數(shù)值子結(jié)構(gòu)，再計算出第i+1步的位移，將其加載到試驗子結(jié)構(gòu)上，不斷重復這個過程直至試驗完成。

圖7 線彈性系統(tǒng)下的實時混合試驗流程

圖8是試驗過程中輪軌力隨時間變化的情況。在OpenSees中建立與試驗相同的軌道-橋梁數(shù)值模型，將試驗過程中記錄下的輪軌力按對應時間步進行加載，對比數(shù)值求解的結(jié)果與試驗子結(jié)構(gòu)記錄的命令信號，如圖9所示。使用ATS方法進行時滯補償后，本次試驗的平均時滯為9.9 ms。將試驗中的命令信號時滯處理后，與數(shù)值解進行對比，結(jié)果如圖10所示。所得結(jié)果驗證了移動荷載積分法應用于線彈性系統(tǒng)下的實時混合試驗的可靠性。

圖8 線彈性系統(tǒng)下的試驗輪軌力

圖9 線彈性系統(tǒng)下數(shù)值信號與命令信號的對比

圖10 線彈性系統(tǒng)下消除時滯后數(shù)值信號與命令信號的對比

2.2 非線性系統(tǒng)下的移動荷載積分法驗證

本節(jié)驗證基于彈塑性子結(jié)構(gòu)方法的移動荷載積分法應用于實時混合試驗中的非線性系統(tǒng)時的可行性。在2.1節(jié)算例的基礎上將跨中處的支座替換為非線性支座來模擬橋梁局部非線性問題。非線性支座采用Duffing材料，材料本構(gòu)關系如下：

σ=E(ε+γε3)

(9)

式中，E=13.89 MPa，γ=-3 000，σ、ε分別為應力和應變。

與2.1節(jié)相同，選取四分之一CRH380A高鐵車體機車(1個輪對)作為試驗子結(jié)構(gòu)，試驗的時間步長為1/128 s，行車速度為230 km/h。

該算法應用于非線性系統(tǒng)時的計算過程比較復雜，下面結(jié)合算例作具體介紹。首先將本算例中的非線性結(jié)構(gòu)用2.1節(jié)算例中的線性結(jié)構(gòu)來代替，此時整個橋軌系統(tǒng)均為線彈性(即轉(zhuǎn)化為線彈性問題)，將單位荷載按照時間步依次作用于軌道-橋梁模型上的相應位置，記錄當前加載情況下橋軌所有節(jié)點的全時程位移響應。根據(jù)記錄的全部加載情況對應的位移響應，形成移動荷載積分法的基本數(shù)據(jù)矩陣?？紤]非線性系統(tǒng)與線性系統(tǒng)之間的內(nèi)力區(qū)別，在基本數(shù)據(jù)矩陣之外，需要編制其他兩類數(shù)據(jù)矩陣。第1類矩陣(矩陣A)是單位荷載在當前時步作用下非線性節(jié)點處的位移響應。本算例中的非線性支座的兩端均為非線性節(jié)點，故需要記錄單位荷載按照時間步依次作用于軌道-橋梁模型上相應位置時兩個節(jié)點的全時程位移響應。根據(jù)記錄的非線性節(jié)點的位移響應，形成對應的數(shù)據(jù)矩陣，用于后續(xù)計算非線性修正力。第2類矩陣(矩陣B)是單位荷載作用于非線性節(jié)點時軌道-橋梁模型各節(jié)點的位移響應。在本算例中，將單位荷載按照時間步依次作用于非線性節(jié)點上，記錄下當前加載情況下軌道-橋梁模型各節(jié)點的全時程位移響應，然后形成對應的數(shù)據(jù)矩陣。

在編制完所需的數(shù)據(jù)矩陣進行實際計算時，先根據(jù)移動荷載積分法的基本矩陣計算出線性情況下當前時步對應位置處的位移響應，再根據(jù)矩陣A計算出當前時步下非線性節(jié)點處的位移，使用該位移并結(jié)合程序中非線性材料的本構(gòu)關系公式計算當前時步下非線性節(jié)點處的非線性修正力，根據(jù)矩陣B計算出非線性修正力作用于非線性節(jié)點時荷載作用位置處的位移響應，疊加之前得到的線性情況下的位移響應，最終得到結(jié)構(gòu)真實的位移響應。

具體的實時混合試驗流程如圖11所示。將第i步移動荷載積分法計算出的位移加載到試驗子結(jié)構(gòu)上，通過振動臺測出的加速度計算出第i+1步的輪軌力，將其加載到數(shù)值子結(jié)構(gòu)。數(shù)值子結(jié)構(gòu)首先根據(jù)輪軌力計算出線彈性系統(tǒng)下當前時步輪軌力作用位置處的位移響應，然后計算出當前時步的非線性修正力，由非線性修正力對線彈性系統(tǒng)下的位移響應進行修正，最終計算出第i+1步的位移，將其加載到試驗子結(jié)構(gòu)上，不斷重復這個過程直至試驗完成。

圖11 非線性系統(tǒng)下的實時混合試驗流程

圖12是非線性系統(tǒng)下試驗過程中輪軌力隨時間變化的情況。同樣將試驗過程中記錄下的輪軌力按對應時間步加載到OpenSees中的數(shù)值模型進行計算，對比數(shù)值求解的結(jié)果與試驗子結(jié)構(gòu)記錄的命令信號，如圖13所示。使用ATS方法進行時滯補償后，本次試驗的平均時滯為8.6 ms，將命令信號根據(jù)試驗時滯處理后，與數(shù)值解進行對比，結(jié)果如圖14所示。所得結(jié)果驗證了基于彈塑性數(shù)值子結(jié)構(gòu)的移動荷載積分法應用于非線性系統(tǒng)下的實時混合試驗的可靠性。

圖12 非線性系統(tǒng)下的試驗輪軌力

圖13 非線性系統(tǒng)下數(shù)值信號與命令信號的對比

圖14 非線性系統(tǒng)下消除時滯后數(shù)值信號與命令信號的對比

3 結(jié)語

為高效、精確地進行高速鐵路車-軌-橋耦合系統(tǒng)實時混合試驗的計算和模擬，提出了一種新型計算方法——移動荷載積分法，并得到如下結(jié)論。

1)將移動荷載積分法應用于實際的高速鐵路系統(tǒng)實時混合試驗中，試驗中每一時步的操作(數(shù)值計算、數(shù)據(jù)交換、時滯補償和作動器加載等環(huán)節(jié))均可在1/128 s內(nèi)完成，滿足實時混合試驗對于計算速度的要求。需要注意的是移動荷載積分法僅適用于線彈性橋梁結(jié)構(gòu)；當出現(xiàn)局部非線性時，結(jié)合彈塑性數(shù)值子結(jié)構(gòu)方法可以較好地解決。

2)將消除時滯后數(shù)值求解的結(jié)果與試驗子結(jié)構(gòu)記錄的命令信號進行對比，發(fā)現(xiàn)試驗結(jié)果擁有可靠的精度——在線彈性系統(tǒng)下，誤差為1×10-10量級；在非線性系統(tǒng)下，誤差為1×10-5量級。

文中選取四分之一CRH380A高鐵機車車體(1個輪對)作為試驗子結(jié)構(gòu)，未考慮輪對間的動力疊加效應，在實際應用中，應考慮其影響。