針對(duì)有砟道床動(dòng)力特性分析的嵌入式離散元-有限元耦合方法

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

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

1 引 言

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

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

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

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

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

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

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

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

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

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

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

(1)

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

(2)

(3)

(4)

(5，6)

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

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

2.2 有砟道床的有限元模型

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

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

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

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

(7)

(8)

(9)

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

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

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

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

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

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

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

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

δW=δUTF

(10)

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

(i=1～20)(11)

(i=1～20)(12)

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

(i=1～20)(13)

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

(i=1～20)(14)

20節(jié)點(diǎn)等參元的形函數(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)系下力作用點(diǎn)的位置坐

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

(18)

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

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

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

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

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

有砟道床下部路堤的沉降量時(shí)程曲線如圖11所示。為研究路堤頂部沿橫斷面方向的沉降規(guī)律，選取3個(gè)觀測(cè)點(diǎn)進(jìn)行觀測(cè)，3個(gè)觀測(cè)點(diǎn)距模型左側(cè)邊界的距離分別為0.16 m,0.78 m和1.91 m。觀察發(fā)現(xiàn)，路堤沉降量曲線的單個(gè)循環(huán)形式與道床頂部沉降曲線形式有兩點(diǎn)主要區(qū)別，一為在路堤的沉降量曲線中，由列車載荷引起的小幅振動(dòng)已經(jīng)消失，其原因是道砟層具有較好的緩沖性能，列車載荷中的小幅振動(dòng)已由道砟層吸收；二為路堤的沉降量曲線中，循環(huán)載荷第二個(gè)載荷峰值產(chǎn)生的沉降量明顯高于第一個(gè)載荷峰值產(chǎn)生的沉降量，主要是由道砟層散體介質(zhì)的力學(xué)性質(zhì)以及快速加載的載荷特性所導(dǎo)致，第一個(gè)載荷峰值在道砟層中產(chǎn)生一部分變形能，隨著快速卸載，道砟層中變形能還未全部釋放，第二次加載已經(jīng)開始，進(jìn)一步使道砟層中變形能不斷增加，導(dǎo)致第二次載荷峰值引起的沉降量高于第一次。此外，在枕木下方(0.16和0.78觀測(cè)點(diǎn))，道床內(nèi)部沉降略高于外側(cè)沉降，但區(qū)別并不顯著。然而在枕木施壓區(qū)外側(cè)，即1.91處觀測(cè)點(diǎn)的位移表現(xiàn)出向上運(yùn)動(dòng)的趨勢(shì)，這主要是由于路堤內(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)生位置與力鏈分布相對(duì)應(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á)到最大值，說(shuō)明嵌入路堤的球形顆粒對(duì)其上端的道砟顆粒提供了明顯的水平方向約束作用，有效減少道砟顆粒在上部載荷作用下向一側(cè)滑移，對(duì)有砟道床離散元-有限元耦合模型的數(shù)值穩(wěn)定性起到有益的作用。

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

4 結(jié) 論

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

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

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

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

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