邵 帥， 呂泉江， 紀(jì)丹陽， 嚴(yán) 穎

(1.內(nèi)蒙古大學(xué) 交通學(xué)院，呼和浩特 010070； 2.中國直升機(jī)設(shè)計研究所，景德鎮(zhèn) 333001；3.大連交通大學(xué) 土木與安全工程學(xué)院，大連 116028)

1 引 言

傳統(tǒng)的有砟鐵路道床結(jié)構(gòu)組成成分多樣，且力學(xué)性質(zhì)復(fù)雜，包括軌枕、道砟、底部道砟、路堤和地基等組成部分[1]。其中，道砟層由大量碎石構(gòu)成，具有很多孔隙，是典型的散體材料[2]。枕木、路堤和地基等部分具有相對密實的材料性質(zhì)，與道砟相比可視為連續(xù)介質(zhì)。有砟道床中散體介質(zhì)與連續(xù)體介質(zhì)并存的結(jié)構(gòu)特點給數(shù)值模擬造成了很多困難[3,4]。發(fā)展一種能夠準(zhǔn)確計算并預(yù)測整體有砟道床力學(xué)行為的計算模型對于有砟道床的養(yǎng)護(hù)、維修及安全運行均具有重要的工程意義。

目前，研究人員為預(yù)測有砟鐵路道床的力學(xué)行為已經(jīng)開展了大量數(shù)值仿真工作，采用的主要數(shù)值方法包括離散元方法(DEM)和有限元方法(FEM)。離散元方法最初由Cundall等[5]提出并應(yīng)用于農(nóng)業(yè)、化學(xué)、巖土和鐵道工程等領(lǐng)域[6-8]。Gong等[9]采用離散元方法模擬了混合有輪胎衍生集料道砟材料的剪切強(qiáng)度，發(fā)現(xiàn)輪胎衍生集料能夠減小剪切應(yīng)力峰值，降低道砟材料的內(nèi)聚力和摩擦角，減少直剪試驗中道砟顆粒的破碎。Jing等[10]通過離散元數(shù)值模擬研究發(fā)現(xiàn)，道砟層對維持軌枕的穩(wěn)定性起到重要作用，并通過試驗進(jìn)行了驗證。Liu等[11]采用離散元方法模擬真實尺度下道砟材料大三軸試驗，基于智能型道砟碎石開發(fā)離散元計算程序的計算精度有了很大提升。Zhang等[12]采用鑲嵌單元構(gòu)造道砟碎石，并通過離散元數(shù)值模擬研究了高速列車載荷頻率對有砟道床累積沉降量的影響。有限元方法經(jīng)過幾十年的發(fā)展已經(jīng)逐漸成熟，在鐵道工程領(lǐng)域得到廣泛認(rèn)可[13，14]。為了能夠獲得不同工況下軌枕的不均勻沉降，Shih等[15]綜合利用Abaqus，Python和Fortran建立了一套名為BaTrack的有限元程序，對于有砟道床力學(xué)行為分析具有較好的模擬精度。

離散元方法和有限元方法在有砟道床動力特性的模擬中各具特點，有限元方法具有足夠的精度且擁有較高的計算效率，離散元方法能夠在細(xì)觀尺度上提供很多力學(xué)信息，為充分發(fā)揮離散元方法和有限元方法各自的優(yōu)勢，研究人員將離散元方法和有限元方法相互耦合并應(yīng)用于鐵道工程領(lǐng)域。道砟層采用離散元計算，路堤與地基通過有限元方法模擬，在離散元域和有限元域的接觸面上依據(jù)虛功原理提出耦合接觸算法，建立離散元-有限元耦合模型。此模型能夠更全面地從宏細(xì)觀多尺度分析整體有砟道床的力學(xué)行為[16-19]。然而，在離散元和有限元的耦合接觸面上，由于球形顆粒表面相對光滑，在與有限元表面的接觸過程中受外部載荷作用容易產(chǎn)生側(cè)向滑動，使得數(shù)值模型不穩(wěn)定，導(dǎo)致數(shù)值結(jié)果準(zhǔn)確度降低。對于實際工況，底部道砟往往會嵌入到路堤中，其水平方向的運動受到很大限制，從而確保其上部結(jié)構(gòu)的穩(wěn)定性。

綜上所述，為增加有砟道床離散元-有限元耦合模型中耦合面上的自鎖能力，并考慮實際工況中底部道砟嵌入路堤中的客觀事實，本文發(fā)展了一種針對有砟道床的嵌入式離散元-有限元耦合模型，有效地提高了耦合模型的數(shù)值穩(wěn)定性，為計算模擬鐵路工程問題提供了一種新的手段。

2 針對有砟道床結(jié)構(gòu)離散元-有限元耦合模型

針對有砟道床結(jié)構(gòu)的計算模擬區(qū)域及道床橫截面如圖1所示。根據(jù)一根枕木的作用范圍，在有砟道床中截取長為4 m，寬為0.23 m的區(qū)域作為計算模擬區(qū)，如圖1(a)所示。在有砟道床的橫截面中考慮五個組成部分，分別為枕木、道砟、底部道砟、路堤和地基，如圖1(b)所示。其中，枕木、路堤和地基相對密實，在數(shù)值模擬中可當(dāng)做連續(xù)介質(zhì)，其力學(xué)行為通過有限元方法進(jìn)行計算。道砟和底部道砟是典型的散體材料，在數(shù)值模擬中作為離散介質(zhì)，其力學(xué)行為通過離散元方法進(jìn)行分析。為了實現(xiàn)對整體有砟道床結(jié)構(gòu)的動力特性分析，需在離散元和有限元的耦合接觸面上建立耦合接觸算法，實現(xiàn)計算參數(shù)在離散元域和有限元域間的傳遞。在圖1(b)中，利用有砟道床橫截面的對稱性質(zhì)，在建模中只需考慮橫截面的右半部分即可。

圖1 有砟道床的模擬區(qū)域及道床橫截面

2.1 有砟道床中的離散元模型

道砟是一種典型的散體材料，其內(nèi)部具有大量孔隙，在上部列車載荷長期作用下容易產(chǎn)生不均勻沉降，整體力學(xué)性質(zhì)十分復(fù)雜。構(gòu)成道砟層的道砟碎石具有幾何形狀極其不規(guī)則的結(jié)構(gòu)特點，道砟顆粒間的咬合作用和自鎖能力較強(qiáng)。為從細(xì)觀尺度上分析有砟道床的動力特性以及道砟顆粒幾何形狀不規(guī)則的結(jié)構(gòu)特點，采用鑲嵌單元對道砟顆粒進(jìn)行構(gòu)造，如圖2所示?？紤]計算效率的影響，構(gòu)造單個道砟顆粒單元使用的球形顆粒數(shù)目為3～11個。由于道砟形狀隨機(jī)性很強(qiáng)，采用較少數(shù)量的球形顆粒進(jìn)行模擬，雖然對當(dāng)前道砟形狀模擬有一定誤差，但能夠保證整體有砟道床的模擬精度。

圖2 道砟碎石以及生成的鑲嵌單元

采用鑲嵌單元鋪設(shè)后的道砟層離散元模型如圖3所示。道砟層頂部承受由枕木傳遞下來的列車荷載，左側(cè)承受x方向的位移約束，右側(cè)為自由表面，前后表面承受y方向的位移約束。構(gòu)造道砟層共使用1727個鑲嵌單元，10098個球形顆粒，孔隙率為0.45。顆粒間的接觸模型采用赫茲模型，法向接觸力計算公式為[20]

(1)

式中A為耗散系數(shù)。不考慮切向粘滯力，根據(jù)摩爾庫倫摩擦定律，切向接觸力計算公式為[21]

(2)

(3)

(4)

(5，6)

圖3 采用鑲嵌單元構(gòu)造的道砟層離散單元模型

表1 道砟材料計算參數(shù)

2.2 有砟道床的有限元模型

有砟道床的有限元模型如圖4所示。其中，上部的梯形和下部的矩形區(qū)域分別代表路堤和地基。采用20節(jié)點等參元對模型進(jìn)行網(wǎng)格剖分，由于模型中枕木下方是主要的承載區(qū)域，故對網(wǎng)格進(jìn)行了加密處理，模型共有1120個單元，5805個節(jié)點。模型左右兩側(cè)受x方向的位移約束，底面受z方向的位移約束，前后面受y方向的位移約束。上表面為自由表面，用于承受由道砟層傳遞下來的列車載荷。路堤和地基的材料參數(shù)列入表2。

圖4 路堤和地基的有限元模型

表2 路堤和地基的材料參數(shù)

采用Newmark隱式時間積分計算路堤和地基的動力響應(yīng)，Newmark方法的計算格式為

(7)

(8)

(9)

式中Fk + 1為時刻tk + 1時刻的等效節(jié)點載荷，δ和α為Newmark方法中的計算參數(shù)，當(dāng)δ≥0.5且α≥0.25(0.5+δ)2時，Newmark方法是無條件穩(wěn)定的，本文取δ=0.5，α=0.25。考慮到路堤的非線性性質(zhì)，采用Newton-Raphson方法求解路堤的力學(xué)行為[22]。

2.3 道砟與路堤交界面的嵌入式離散元-有限元耦合模型

在以往針對有砟道床結(jié)構(gòu)的離散元-有限元耦合計算方法中，使用球形顆粒構(gòu)造的鑲嵌單元表面比較圓潤，與路堤上表面很難產(chǎn)生較強(qiáng)的自鎖和咬合作用，在列車載荷的作用下有砟層容易產(chǎn)生整體的側(cè)向移動。此外，在實際的有砟道床中，由于路堤的塑性性質(zhì)，底部道砟會嵌入路堤上表面。為使數(shù)值模型中道砟層更加穩(wěn)定并考慮實際工況，本文發(fā)展了針對有砟道床結(jié)構(gòu)的嵌入式離散元-有限元耦合方法。

在嵌入式離散元-有限元耦合模型中，設(shè)置一層嵌入路堤頂部的耦合過渡層，如圖5所示。過渡層共有483個球形顆粒，呈規(guī)則的矩形方式排列，球形顆粒的球心位于路堤的上表面。過渡層將承受上部道砟層傳遞的列車載荷，并通過耦合接觸算法將耦合接觸力作為等效節(jié)點載荷傳遞至有限元

圖5 離散元和有限元耦合過渡層

網(wǎng)格中，進(jìn)而計算下部路堤和地基的變形。當(dāng)下部結(jié)構(gòu)的節(jié)點位移計算完成后，通過形函數(shù)插值可進(jìn)一步更新過渡層中各球形顆粒的位置，從而實現(xiàn)離散元方法和有限元方法的耦合。在實際的應(yīng)用中，過渡層中球形顆粒的粒徑應(yīng)參照鐵路碎石道床底碴的級配要求進(jìn)行選取，本文以30 mm粒徑為例進(jìn)行計算。

圖6 嵌入式離散元-有限元耦合方法的接觸模型

根據(jù)虛功原理，離散元等效接觸力產(chǎn)生的外力虛功δW為

δW=δUTF

(10)

式中F為等效接觸力矢量，U為力作用點處的位移矢量，可通過對20節(jié)點等參元節(jié)點位移進(jìn)行插值求得

(i=1～20)(11)

(i=1～20)(12)

有限元節(jié)點載荷所做虛功為

(i=1～20)(13)

根據(jù)虛功原理，有限元等效節(jié)點載荷Fnodal,i的計算公式為

(i=1～20)(14)

20節(jié)點等參元的形函數(shù)為

(i=1～8)(15a)

(i=17～20)(15b)

(i=9,11,13,15)(15c)

(i=10,12,14,16)(15d)

(16)

(i=3,4,7,8)(17a)

(i=19,20)(17b)

(i=11,15)(17c)

式中(x,y)是整體坐標(biāo)系下力作用點的位置坐

圖7 接觸力等效節(jié)點載荷

(18)

3 列車載荷作用下有砟道床結(jié)構(gòu)動力特性的計算結(jié)果及分析

通過本文建立的嵌入式離散元-有限元耦合模型,對有砟道床結(jié)構(gòu)進(jìn)行動力特性分析，在枕木頂部施加的列車載荷如圖9所示。此列車載荷由翟婉明等[23]提出的車輛-軌道耦合動力學(xué)方法獲得。首先在軌道上設(shè)置一個觀測點，當(dāng)列車車輪通過觀測點時，將產(chǎn)生一個力的峰值。列車的每節(jié)車廂共有2組車輪，4根輪軸，一節(jié)車廂完全通過測試點后會產(chǎn)生兩組共4個載荷峰值，在此過程中會伴隨產(chǎn)生小幅振動。列車載荷的最大值為20 kN，一個載荷循環(huán)周期為1.2 s。本文算例共施加了15次載荷循環(huán)。

圖8 有砟道床嵌入式離散元-有限元耦合模型

圖9 施加于枕木上的列車載荷

有砟道床頂部沉降量的時程曲線如圖10所示。沉降量曲線的每個循環(huán)與載荷曲線形式類似，具有兩組共4個峰值，伴隨有小幅振動。在加載初期，由于有砟道床內(nèi)部道砟相對松散，故道床頂部沉降明顯。隨著加載的不斷進(jìn)行，有砟道床內(nèi)部的道砟顆粒會產(chǎn)生相互錯動以及空間位置的重新調(diào)整，在此過程中有砟道床的密實程度不斷增加，致使道床頂部沉降逐漸平緩，此時有砟道床已沒有明顯的累積沉降量。

有砟道床下部路堤的沉降量時程曲線如圖11所示。為研究路堤頂部沿橫斷面方向的沉降規(guī)律，選取3個觀測點進(jìn)行觀測，3個觀測點距模型左側(cè)邊界的距離分別為0.16 m,0.78 m和1.91 m。觀察發(fā)現(xiàn)，路堤沉降量曲線的單個循環(huán)形式與道床頂部沉降曲線形式有兩點主要區(qū)別，一為在路堤的沉降量曲線中，由列車載荷引起的小幅振動已經(jīng)消失，其原因是道砟層具有較好的緩沖性能，列車載荷中的小幅振動已由道砟層吸收；二為路堤的沉降量曲線中，循環(huán)載荷第二個載荷峰值產(chǎn)生的沉降量明顯高于第一個載荷峰值產(chǎn)生的沉降量，主要是由道砟層散體介質(zhì)的力學(xué)性質(zhì)以及快速加載的載荷特性所導(dǎo)致，第一個載荷峰值在道砟層中產(chǎn)生一部分變形能，隨著快速卸載，道砟層中變形能還未全部釋放，第二次加載已經(jīng)開始，進(jìn)一步使道砟層中變形能不斷增加，導(dǎo)致第二次載荷峰值引起的沉降量高于第一次。此外，在枕木下方(0.16和0.78觀測點)，道床內(nèi)部沉降略高于外側(cè)沉降，但區(qū)別并不顯著。然而在枕木施壓區(qū)外側(cè)，即1.91處觀測點的位移表現(xiàn)出向上運動的趨勢，這主要是由于路堤內(nèi)側(cè)受擠壓導(dǎo)致外側(cè)的膨脹。

圖11 路堤頂部沉降量

有砟道床中道砟層力鏈分布、路堤與地基內(nèi)位移和應(yīng)力云圖分別如圖12和圖13所示。發(fā)現(xiàn)在道砟層內(nèi)部力鏈主要分布于枕木下方，形成承載骨架，能夠?qū)㈨敳康牧熊囕d荷傳遞至下部路堤。由于耦合模型中存在嵌入路堤中的一層道砟顆粒，耦合接觸面上具有較大的摩擦系數(shù)，咬合作用明顯，因此道砟層中的力鏈基本限制于枕木下方的矩形區(qū)域，力鏈沒有產(chǎn)生向四周發(fā)散的現(xiàn)象。砟肩處受力較小，內(nèi)部沒有明顯的力鏈分布，由于枕木的擠壓作用，砟肩處有輕微隆起。路堤和地基的位移云圖與道砟層中的力鏈分布相吻合，主要位于枕木下方。最大位移產(chǎn)生位置與力鏈分布相對應(yīng)，具有一定隨機(jī)性，道砟層中的主要承載骨架將傳遞大部分的列車載荷，進(jìn)而使其下方的路堤產(chǎn)生較大位移。觀察圖13的應(yīng)力云圖可以發(fā)現(xiàn),路堤和地基中的應(yīng)力分布與道砟層中的力鏈分布相互吻合，應(yīng)力在道床內(nèi)部較大，向邊緣逐漸遞減。

圖12 有砟道床內(nèi)部力鏈分布及位移云圖

圖13 有砟道床內(nèi)部力鏈分布及應(yīng)力云圖

道砟-路堤耦合交界面豎直方向與水平方向的載荷分布如圖14所示。發(fā)現(xiàn)豎直方向載荷與水平方向載荷均與力鏈分布相吻合，位于枕木下方，呈現(xiàn)出很強(qiáng)的非均勻性質(zhì)，主要由道砟層中道砟碎石的隨機(jī)排列導(dǎo)致。通過耦合交界面水平載荷分布發(fā)現(xiàn),水平載荷在枕木右端下方區(qū)域達(dá)到最大值，說明嵌入路堤的球形顆粒對其上端的道砟顆粒提供了明顯的水平方向約束作用，有效減少道砟顆粒在上部載荷作用下向一側(cè)滑移，對有砟道床離散元-有限元耦合模型的數(shù)值穩(wěn)定性起到有益的作用。

圖14 道砟-路堤耦合面上的載荷分布云圖

4 結(jié) 論

本文針對有砟道床的結(jié)構(gòu)特點，考慮實際鐵路中底部道砟嵌入路堤的客觀事實，發(fā)展了一種嵌入式離散元-有限元耦合模型,得到以下結(jié)論。

(1) 有砟道床中非規(guī)則的道砟顆?？赏ㄟ^鑲嵌單元模擬，然而鑲嵌單元由若干球形顆粒構(gòu)成，表面相對光滑，導(dǎo)致其與邊界作用時自鎖能力較差，容易產(chǎn)生滑動。

(2) 根據(jù)虛功原理建立了嵌入式離散元-有限元耦合模型，嵌入路堤中的球形顆粒受到路堤的約束作用很難產(chǎn)生側(cè)向滑動，顆粒間具有較好的自鎖能力，從而通過顆粒間咬合作用能減少上部道砟層在列車載荷作用下的側(cè)向移動。

(3) 嵌入式離散元-有限元耦合模型能從力鏈分布、沉降量、位移云圖及耦合交界面處力學(xué)響應(yīng)等方面對有砟道床在宏細(xì)觀多尺度上進(jìn)行計算分析。

工程中大量問題可通過離散元-有限元耦合方法進(jìn)行數(shù)值模擬，對于類似有砟道床的散體介質(zhì)與連續(xù)介質(zhì)共同存在的特殊結(jié)構(gòu)，采用嵌入式耦合模型處理耦合交界面更加符合工程實際，具有更好的數(shù)值穩(wěn)定性。