陳智鋒， 李穎雄， 賴遠明， 蘇 成，4

（1.華南理工大學土木與交通學院，廣東廣州510640； 2.西南交通大學智慧城市與交通學院，四川成都611756；3.中國科學院西北生態(tài)環(huán)境資源研究院凍土工程國家重點實驗室，甘肅蘭州730000； 4.華南理工大學亞熱帶建筑科學國家重點實驗室，廣東廣州510640）

0 引言

我國為世界第三凍土大國，多年凍土面積占國土總面積的21.5%，其中高原凍土面積更是居于世界之最［1］。隨著“一帶一路”倡議深入推進實施，凍土區(qū)的鐵路建設與維護得到高度重視。同時，國民經(jīng)濟的發(fā)展使得列車行駛頻次有所提升，加之受全球變暖［2］等氣候因素的影響，凍土區(qū)鐵路病害問題越來越突出。為保障凍土區(qū)鐵路的安全運營，亟須發(fā)展高效準確的計算方法，分析列車行駛對凍土路基動力響應的影響。

凍土路基在列車荷載作用下的動力響應是當前的研究熱點之一。馬巍等［3］和朱占元等［4］開展了列車荷載作用下凍土路基動力響應的現(xiàn)場監(jiān)測研究，為列車通過凍土路段時路基穩(wěn)定性的預測分析提供了實測資料依據(jù)。Stevens［5］、孔祥兵等［6］和李雙洋等［7-9］采用數(shù)值分析方法，研究了一次行車荷載作用下路基的位移和應力等動力響應問題。此外，還有部分學者［10-13］進行了列車荷載作用下凍土路基的穩(wěn)定性、蠕變、軌枕響應和車速對路基響應影響等方面的數(shù)值模擬研究。

上述研究中，都采用了確定性的列車荷載，選取某一次列車駛過產(chǎn)生的荷載樣本計算路基響應，其結果難以反映鐵路實際運行中由軌道不平順及列車類型和軸重差異導致的荷載隨機性［14］。對于路基的特定橫斷面，列車駛入和駛出的過程導致其頂部荷載的統(tǒng)計特征具有顯著的時變特性。因此，在研究路基動力響應時應將其頂部列車荷載按非平穩(wěn)隨機荷載考慮［7，15］。由于凍土獨特的物理、熱力學性質以及非平穩(wěn)隨機振動分析的復雜性，目前關于隨機列車荷載作用下凍土路基動力響應的研究較少。

目前，解決結構非平穩(wěn)隨機振動問題的主要方法有蒙特卡羅模擬法、虛擬激勵法［16］、概率密度演化法［17］以及時域顯式法［18］等。其中，時域顯式法通過構建系統(tǒng)動力響應的時域顯式表達式，在隨機分析中可以針對任意關鍵響應進行降維，具有理想的計算效率，已成功應用于高層建筑［19］、大跨度橋梁［20］和車橋耦合系統(tǒng)［21］等大型復雜結構系統(tǒng)的隨機振動分析。在隨機振動分析中，若已知隨機激勵的均值和互相關函數(shù)等統(tǒng)計特征，則可按統(tǒng)計矩的運算規(guī)律根據(jù)時域顯式表達式直接求得響應的前兩階統(tǒng)計矩，這一方法稱為時域顯式直接法。若通過數(shù)值模擬手段生成大量隨機激勵樣本，則可結合時域顯式表達式進行高效蒙特卡羅模擬，這一方法稱為時域顯式蒙特卡羅模擬法。

為了解決隨機列車荷載作用下凍土路基動力響應分析這一工程問題，本文首先提出一種隨機列車荷載的簡化模擬方法，以獲取隨機分析所需的荷載樣本。然后，通過有限元脈沖響應分析構建路基結構動力響應的時域顯式表達式。在此基礎上，進一步利用時域顯式蒙特卡羅模擬法，對隨機列車荷載作用下凍土路基的關鍵響應開展統(tǒng)計分析，可大幅提高凍土路基結構隨機振動的計算效率。最后，通過數(shù)值算例展示時域顯式蒙特卡羅模擬法在凍土路基非平穩(wěn)隨機振動分析中的計算精度與效率。

1 隨機列車荷載模擬方法

在凍土路基隨機振動問題中，需要列車駛過引起的隨機荷載的統(tǒng)計特性或大量樣本作為輸入。本文采用常見的二維路基模型［7］進行計算，此類模型需要列車行駛引起的隨機軌枕豎向作用力作為荷載輸入［15］。目前，相關實測數(shù)據(jù)較少，難以滿足凍土路基隨機振動分析的需求。為了獲取分析所需輸入，可以對確定性列車荷載的數(shù)值模擬方法［22-24］加以改進，在其中引入隨機參數(shù)，進而生成大量隨機列車荷載樣本。

列車在軌道上行駛的簡圖如圖1 所示，其中列車行駛速度為v。為了描述列車在沿軌道方向上的位置，建立如圖1 所示的直線坐標系Ox。列車前端的初始坐標假定為-L′0。

圖1 列車行駛簡圖Fig. 1 Moving train schematic

考慮如圖1 所示的列車模型，其車廂的數(shù)量為N，包括機車和普通兩類車廂。在列車行駛過程中，每個輪對的荷載通過軌道及軌枕傳遞到路基頂面。對于某一關注軌枕，列車對路基的豎向作用力F(t)的表達式為

式中：fi(t)表示第i節(jié)車廂的荷載傳遞到關注軌枕下方的部分，可進一步表示為該車廂所有輪對產(chǎn)生的豎向軌枕力之和。若該車廂為機車車廂，則其表達式為

式中：w為車廂輪對數(shù)量，其中機車和普通兩類車廂的輪對數(shù)分別為6 和4；Pij(t)為第i節(jié)車廂第j組輪對在軌道上產(chǎn)生的荷載；K(xij)為關注軌枕對輪對作用力的分擔系數(shù)，按日本規(guī)范的經(jīng)驗假設輪對荷載. 附 近5 根. 枕 按0.1∶0.2∶0.4∶0.2∶0.1 分擔［25］，可進一步表示為

式中：Ls為軌枕間距，如圖2 所示；xij為輪對在Ox坐標系中的坐標。當?shù)趇節(jié)車廂為機車車廂時，xij具體表示為

圖2 道床頂面荷載沿線路方向分布圖Fig. 2 Distribution of loading on top surface of road bed along the road

式中：ai第i節(jié)車廂同一個轉向架內(nèi)相鄰兩個輪對的軸距；bi為第i節(jié)車廂不同轉向架相鄰兩個輪對的軸距；Ld，k為第k節(jié)車廂長度；L0為車廂后端與該車廂最后一組輪對水平方向上的距離；-L′0為列車前端在該坐標系內(nèi)的初始坐標。當?shù)趇節(jié)車廂為普通車廂時，也可以參考圖1 采用類似的方式計算輪對坐標。

列車第i節(jié)車廂第j組輪對產(chǎn)生的作用力Pij(t)可以采用反映行車平穩(wěn)性、附加動載和軌面波形磨耗效應三個控制條件的激勵力來模擬［24］，其表達式為

軌道幾何不平順曲線波長及相應的正矢和列車行駛速度都具有一定的不確定性［14］。可以引入相應的隨機變量，以量化這些參數(shù)的不確定性。對描述相關參數(shù)的隨機變量進行抽樣，將樣本代入式（5）計算輪對在軌道上產(chǎn)生的作用力，進而根據(jù)式（1）可以計算所關注軌枕的豎向作用力，即可獲得后續(xù)分析所需的隨機列車荷載樣本。

2 凍土路基隨機振動分析的時域顯式蒙特卡羅模擬法

作為一種多孔結構的物質，凍土中包含了土、水、冰、氣四種成分，其動態(tài)力學行為較為復雜。為了實現(xiàn)對凍土路基這一復雜工程結構的隨機振動分析，需要在反映凍土動力響應本質特征的前提下，對問題進行簡化以降低建模難度和提升分析計算效率。凍土力學性質對季節(jié)性溫度變化十分敏感，然而本文關注的列車荷載作用下凍土路基動力響應分析問題最多僅考慮數(shù)分鐘的時程，因此可將列車荷載作用期間凍土溫度及力學性質視為不變［26-27］。此外，列車荷載引起的凍土應變相對較小，特別是在凍結狀態(tài)下，一般不會造成路基下方凍土的塑性變形［28］。因此，為降低后續(xù)非平穩(wěn)隨機振動分析的難度，可將凍土路基中各土層假定為各向同性線彈性材料［29-30］。

2.1 動力響應顯式表達式

根據(jù)上述假定，可將列車荷載作用下凍土路基的動力響應分析問題簡化為線性非平穩(wěn)隨機振動問題。時域顯式蒙特卡羅模擬法是解決這類問題的高效手段。在該方法中，首先需要建立結構響應的時域顯式表達式。將凍土路基離散為具有nd個自由度的有限元模型，在非平穩(wěn)隨機激勵作用下該結構的運動微分方程可表示為

式中：M，C，K分別表示凍土路基結構的質量矩陣、阻尼矩陣和剛度矩陣；U，U?，U?分別為位移、速度和加速度向量；L為定位隨機激勵的nd× 1 階常向量；F(t)為具有非平穩(wěn)性的隨機列車荷載。

將時間步長記為Δt，時程積分步數(shù)記為n，并把列車荷載F(t)離散為各時刻t1，t2，…，tn處的荷載F1，F(xiàn)2，…，F(xiàn)n，定義狀態(tài)向量為V=[UTU?T]T。采用任意一種數(shù)值積分方法求解式（7），均可得到路基結構動力響應關于各時刻列車荷載的顯式表達式。不失一般性，可假設結構初始狀態(tài)V0=V(0) = 0以及荷載F0=F(0) = 0，則ti=iΔt時刻結構的狀態(tài)向量Vi=V(ti)的顯式表達式為

式中：Ai，1，Ai，2，…，Ai，i為只與路基結構參數(shù)有關的系數(shù)向量，反映結構參數(shù)對結構動力響應的影響。

參考文獻［18］和［31］推導了系數(shù)向量的閉合公式，如下式所示：

式中：T，Q1和Q2的表達式由求解式（7）時所采用的積分格式?jīng)Q定，當采用Newmark-β積分格式［31］，可進一步表示為

式中：γ= 0.5，β= 0.25 為Newmark-β算法中的計算參數(shù)。

式（9）展示了各系數(shù)向量間的內(nèi)在聯(lián)系。根據(jù)該關系，可將各時刻結構響應顯式表達式所需的系數(shù)向量排列成如表1 所示。由表1 可以看出，僅需計算第一列系數(shù)向量Ai，1(i= 1，2，…，n)即可獲得所有時刻響應的顯式表達式，其計算量與1 次確定性時程分析的計算量相當。該列系數(shù)向量的物理含義如圖3 所示，可方便地通過在結構上施加單位三角脈沖激勵進行時程分析求得。在存儲量方面，由于結構隨機響應分析中并不需要關注結構的全部響應，因此Ai，1(i= 1，2，…，n)中的元素不需要全部儲存。假設所關心的結構響應量個數(shù)為m，則需要存儲的元素個數(shù)僅為2mn，其中n為時程分析步數(shù)?？梢?，需儲存的系數(shù)向量元素個數(shù)與結構自由度數(shù)無關，即使對于大型復雜路基結構隨機振動問題，在矩陣元素的存儲方面也不會出現(xiàn)問題。從本質上看，這是由于結構的物理演變機制已經(jīng)在式（8）中得到全面反映，結構響應已經(jīng)全部解耦，可以僅對所關心的響應量進行后續(xù)隨機分析而不涉及其他響應量，這是采用時域顯式表達式的主要優(yōu)勢之一。

表1 各時刻結構響應對應的系數(shù)向量Table 1 Coefficient vectors for structural responses at each time instant

圖3 在時刻t1施加全三角單位脈沖激勵Fig. 3 Unit impulse excitation at time t1

2.2 關鍵響應均值、標準差和平均峰值的計算

在工程結構分析中，通常只關注少量的結構關鍵響應。對于結構某一關鍵響應r，由式（8）可以得到關于r的一維顯式表達式為

其中

式中：q為由狀態(tài)向量V轉換得到關鍵響應r的轉換行向量，當r為某一位移或速度響應時，q由元素0和1 組成；當r為某一應力響應時，q依賴于相應的本構關系。

通過本文第1節(jié)所述方法生成大量隨機列車荷載樣本，假設一共得到M個列車荷載F(t)的樣本。在第k個向量樣本Fk(t)(k= 1，2，…，M)作用下，由式（11）可以得到關鍵響應r在各時刻的值為

則關鍵響應r的均值、標準差和平均峰值分別為

可見，利用所建立的時域顯式表達式，結合蒙特卡羅模擬，可以高效地計算關鍵響應的統(tǒng)計值，這種方法稱為時域顯式蒙特卡羅模擬法。而在傳統(tǒng)蒙特卡羅模擬法中，對于每一個列車荷載樣本，均需要根據(jù)時程分析法求解式（7）所示的運動微分方程。因此，時域顯式蒙特卡羅模擬法的計算量遠小于傳統(tǒng)蒙特卡羅模擬法的計算量，從而可以使蒙特卡羅模擬法直接應用于凍土路基的隨機振動分析。

3 數(shù)值算例

3.1 隨機列車荷載

本文主要考慮軌道幾何不平順曲線波長、相應的正矢以及列車行駛速度的不確定性。引入服從高斯分布的隨機變量［24，32-33］，以量化軌道不平順波長、正矢和列車行駛速度的不確定性，如表2 所示。列車荷載模擬過程中采用的時間步長為Δt=0.005 s，考慮列車編組的方式為2 節(jié)NJ2 機車+4 節(jié)YZ25T 客車，即式（1）中N= 6，列車模型參數(shù)如表3所示。圖1 中列車前端的初始坐標可以假定為L′0=0，圖2中軌枕間距為Ls= 0.556 m。對表2所列的隨機變量進行抽樣，將一組樣本代入式（5）和式（6）中，就可以得到輪對在軌道上產(chǎn)生的作用力，然后根據(jù)式（1）和式（2）即可得到所關注軌枕的豎向作用力，獲得一個荷載樣本。

表2 波長、正矢和車速隨機變量Table 2 Random variables of wavelength，versine and vehicle speed

表3 青藏鐵路列車模型參數(shù)［34］Table 3 Vehicle model parameters of Qinghai-Tibet Railway［34］

為驗證本文提出的隨機列車荷載模擬方法的準確性，將生成的一個荷載樣本與文獻［15］中采用車輛-軌道耦合動力學計算所得的軌枕豎向作用力時程進行對比，如圖4 所示。由圖4 可知，本文方法得到的荷載樣本在幅值和波形上與文獻［15］的方法得到的荷載時程結果吻合，說明本文所提的隨機列車荷載模擬方法是可行的。

圖4 列車荷載時程樣本Fig. 4 Time history sample of train load

3.2 路基結構參數(shù)

以青藏鐵路某傳統(tǒng)道砟路基為研究對象，計算模型如圖5所示。計算區(qū)域中的土層依次為道碴層、路基填土層、砂土層、粉質黏土層及弱風化巖層，如圖5（a）所示。凍土的力學參數(shù)不僅依賴于各土層材料，還取決于土層的溫度。本文選取青藏鐵路路基修建完成后第10 年夏冬兩個典型季節(jié)氣溫最高時路基內(nèi)部的溫度分布，根據(jù)參考文獻［7］中的擬合公式計算路基各土層的力學參數(shù)，結果見表4～5。

圖5 二維凍土路基計算模型Fig. 5 Two-dimensional permafrost embankment model：geometry of embankment［7］（a），finite-element model of permafrost embankment（b）

表4 路基各土層的參考溫度和力學參數(shù)（夏季）Table 4 Mechanical parameters of the soils in embankment model（Summer）

表5 路基各土層的參考溫度和力學參數(shù)（冬季）Table 5 Mechanical parameters of the soils in embankment model（Winter）

根據(jù)表4～5 所得的各土層的力學參數(shù)，在通用有限元軟件ANSYS 中建立二維凍土路基有限元計算模型。建模時按平面應變問題考慮［34-35］，采用PLANE42 單元，模型單元數(shù)為3 500，節(jié)點數(shù)為3 631，總自由度數(shù)為7 262，如圖5（b）所示?？紤]在路基兩側及底部布置黏彈性人工邊界，在ANSYS中采用COMBIN14單元進行模擬，單元的剛度和阻尼參數(shù)根據(jù)文獻［36］所述的方法確定。分析過程中選用瑞利阻尼［7］，阻尼參數(shù)為α=β= 0.03，時程積分步長取Δt= 0.005 s。在路基設計和監(jiān)測過程中，路基填土層頂面和底面的豎向位移、豎向速度和豎向正應力是設計人員較為關注的關鍵響應，因此選擇路基填土層頂面和底面的中點作為關注位置，分別記為1#點和2#點，如圖5 所示。由于時域顯式表達式獨特的降維計算優(yōu)勢，在后續(xù)計算過程中可以只針對這兩個關注位置進行降維計算，從而大幅度提升蒙特卡羅模擬的效率。

3.3 計算結果討論

為驗證所建立的時域顯式表達式的計算精度，本文分別采用時域顯式表達式和ANSYS 直接時程分析法計算凍土路基在一個列車荷載樣本作用下的動力響應，結果如圖6 所示。由圖6 可知，采用兩種計算方法得到的結果一致，從而驗證了時域顯式表達式在凍土路基動力響應分析中的計算精度。

圖6 路基1#點豎向動力響應時程（夏季）Fig. 6 Time histories of dynamic responses at point 1#（Summer）：time history of vertical displacement（a），time history of vertical velocity（b），time history of vertical normal stress（c）

為展示時域顯式蒙特卡羅模擬法相比于傳統(tǒng)蒙特卡羅模擬法的效率優(yōu)勢，分別采用了不同數(shù)量的樣本作為輸入進行隨機分析，兩種方法所需計算時間如表6 所示。由第2.1 節(jié)中的分析可知，時域顯式蒙特卡羅模擬法的計算時間包括求取系數(shù)向量的計算時間和實施蒙特卡羅模擬的計算時間兩部分。由于該方法在實施蒙特卡羅模擬時無需反復求解凍土路基結構運動微分方程，第二部分耗時很短。而傳統(tǒng)蒙特卡羅模擬法在每次樣本分析中均需求解路基結構運動微分方程，計算非常耗時，其效率遠低于時域顯式蒙特卡羅模擬法的計算效率。由表6還可以看出，隨著樣本數(shù)的增大，傳統(tǒng)蒙特卡羅模擬法計算時間呈線性增長，而時域顯式蒙特卡羅模擬法的計算時間增長不多，說明當樣本數(shù)量規(guī)模較大時，時域顯式蒙特卡羅模擬法的計算優(yōu)勢更加明顯。

表6 兩種不同方法的計算時間比較Table 6 Comparison of time costs between two methods

當隨機激勵樣本數(shù)分別取500、1 000 和2 000時采用時域顯式蒙特卡羅模擬法獲得的夏季凍土路基豎向位移、豎向速度和豎向正應力計算結果同時顯示于圖7 中。由圖7 可以發(fā)現(xiàn)，當樣本量達到1 000時，凍土路基各隨機響應的均值和標準差結果已經(jīng)收斂。由表6 可知，采用傳統(tǒng)蒙特卡羅模擬法處理1 000 個樣本耗時約4 d，而時域顯式蒙特卡羅模擬法耗時約6 min。

圖7 不同樣本數(shù)下路基1#點2#點豎向動力響應統(tǒng)計矩時程（夏季）Fig. 7 Statistical moment time histories of dynamic responses at point 1#and point 2#with different numbers of samples（Summer）：mean value time histories of vertical displacements at point 1#and point 2#（a），mean value time histories of vertical velocities at point 1#and point 2#（b），mean value time histories of vertical normal stresses at point 1#and point 2#（c），stadnard deviation time histories of vertical displacements at point 1#and point 2#（d），standard deviation time histories of vertical velocities at point 1#and point 2#（e），standard deviation time histories of vertical normal stresses at point 1#and point 2#（f）

由時域顯式蒙特卡羅模擬法還可以得到夏季和冬季凍土路基1#點和2#點動力響應平均峰值，結果如表7 所示。由表7 可知，列車荷載作用下凍土路基的動力響應不僅與路基深度有關，還受季節(jié)變化影響。總的來說，夏季時路基淺層的位移和速度響應最為劇烈，這是因為凍土是一種對溫度敏感的土體。夏季溫度相對冬季較高，導致凍土的彈性模量等力學參數(shù)降低，從而導致了路基中更為劇烈的位移和速度響應。

表7 夏季和冬季凍土路基響應平均峰值Table 7 Mean peak value of permafrost embankment in Summer and Winter

4 結論

凍土路基隨機振動問題十分復雜，目前對凍土路基隨機振動的理論研究較少。本文對上述問題展開了探索性的分析與研究，取得了一定的研究成果?；谝陨戏治鲅芯浚玫降闹饕Y論如下：

（1）采用本文提出的隨機列車荷載模擬方法，可以快速生成凍土路基隨機振動分析所需的大量列車荷載樣本用于蒙特卡羅模擬。

（2）本文采用時域顯式蒙特卡羅模擬法計算了凍土路基在兩個典型季節(jié)下的動力響應平均峰值，發(fā)現(xiàn)路基的位移和速度響應在夏季最為劇烈。

（3）對比了時域顯式蒙特卡羅模擬法和傳統(tǒng)蒙特卡羅模擬法在處理不同數(shù)量的樣本時所需的計算時間，展示了時域顯式蒙特卡羅模擬法在處理凍土路基隨機振動問題中的計算精度和效率。

（4）本文中的時域顯式表達式僅適用于線性系統(tǒng)，因此將凍土理想化為線彈性材料，這一理想化帶來的誤差尚需進一步分析。

由于目前尚無隨機列車荷載的現(xiàn)場實測數(shù)據(jù)，因此本文僅從計算理論上進行初步探索，以期能為凍土路基的設計、維護及研究提供理論依據(jù)和參考。凍土是一種多孔多相材料，特別是在夏季高溫時其力學行為十分復雜。為反映這一特性，后續(xù)研究將結合等效線性化方法，考慮非線性凍土本構關系，計算凍土路基永久變形和可靠度。