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

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

0 引言

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

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

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

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

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

1 隨機列車荷載模擬方法

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

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

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

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

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

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

式中：Ls為軌枕間距，如圖2 所示；xij為輪對在Ox坐標系中的坐標。當第i節車廂為機車車廂時，xij具體表示為

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

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

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

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

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

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

2.1 動力響應顯式表達式

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

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

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

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

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

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

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

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

表1 各時刻結構響應對應的系數向量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為由狀態向量V轉換得到關鍵響應r的轉換行向量，當r為某一位移或速度響應時，q由元素0和1 組成；當r為某一應力響應時，q依賴于相應的本構關系。

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

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

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

3 數值算例

3.1 隨機列車荷載

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

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

表3 青藏鐵路列車模型參數［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 路基結構參數

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

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

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

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

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

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）

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

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

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

圖7 不同樣本數下路基1#點2#點豎向動力響應統計矩時程（夏季）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 可知，列車荷載作用下凍土路基的動力響應不僅與路基深度有關，還受季節變化影響。總的來說，夏季時路基淺層的位移和速度響應最為劇烈，這是因為凍土是一種對溫度敏感的土體。夏季溫度相對冬季較高，導致凍土的彈性模量等力學參數降低，從而導致了路基中更為劇烈的位移和速度響應。

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

4 結論

凍土路基隨機振動問題十分復雜，目前對凍土路基隨機振動的理論研究較少。本文對上述問題展開了探索性的分析與研究，取得了一定的研究成果。基于以上分析研究，得到的主要結論如下：

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

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

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

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

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