邵 帥， 呂泉江， 嚴(yán) 穎， 季順迎＊

（1.內(nèi)蒙古大學(xué) 交通學(xué)院，呼和浩特010070；2.大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室，大連116024；3.大連交通大學(xué) 土木與安全工程學(xué)院，大連116028）

1 引 言

以混凝土或?yàn)r青砂漿組成的無砟軌道結(jié)構(gòu)形式在我國(guó)高速鐵路線路上廣泛應(yīng)用，但對(duì)于不良地質(zhì)和大跨度橋梁等地段，仍為有砟軌道設(shè)計(jì)，不可避免地會(huì)遇到有砟和無砟道床的過渡段。在有砟-無砟軌道過渡段，由于材料和剛度等方面的差異，產(chǎn)生和存在的沉降差異較大。隨著沉降差的增大，將導(dǎo)致列車產(chǎn)生垂向加速度，對(duì)軌道產(chǎn)生很大的豎向沖擊力，加速軌道結(jié)構(gòu)的劣化。因此，在過渡段采取相應(yīng)的措施減少沉降差十分必要。在施工過程中可采用將道砟顆粒嵌入無砟道床的方式增加道砟顆粒之間的咬合力，降低沉降量，保證過渡段道床的穩(wěn)定。

有砟道床作為一種散體材料，具有很強(qiáng)的非均勻、非連續(xù)性和各項(xiàng)異性等非線性性質(zhì)［1］。道砟顆粒的非規(guī)則形狀對(duì)整體道床的穩(wěn)定性起到重要作用［2］。離散元方法作為一種有效的數(shù)值方法，在模擬散體材料時(shí)具有很大的優(yōu)勢(shì)［3－4］。Cundall等［5］提出離散單元法之后，該方法廣泛地用于分析道砟材料的力學(xué)行為［6－9］。目前，對(duì)非規(guī)則道砟的形態(tài)描述主要有粘結(jié)單元和鑲嵌單元。粘結(jié)單元是由球形顆粒在不同位置進(jìn)行組合生成的可破碎組合單元，當(dāng)顆粒間的力大于粘結(jié)強(qiáng)度時(shí)，粘結(jié)單元可發(fā)生破壞［10］。采用粘結(jié)單元可模擬單個(gè)道砟的破碎過程［11］，也可以模擬列車荷載下道砟道床的顆粒破碎［8，12，13］。鑲嵌單元是在顆粒排列中存在一定重疊量的組合單元，采用鑲嵌單元可模擬道砟顆粒［14，15］，并且在計(jì)算過程中將整個(gè)鑲嵌單元視為剛體，保持單元不會(huì)破碎［16－18］。此外，采用離散單元法可分析土工格柵加固下道砟直剪和沉降的力學(xué)規(guī)律［19－21］。

有限單元法可有效地分析連續(xù)體材料的動(dòng)力學(xué)特性，在鐵路道床動(dòng)力分析中也有較為廣泛的應(yīng)用［22－24］。將鐵路道床視為連續(xù)體材料，并賦予材料屬性，利用基于連續(xù)介質(zhì)力學(xué)的有限元法對(duì)鐵路道床進(jìn)行分析研究，得出的參數(shù)可作為鐵路設(shè)計(jì)的有益參考［25－28］。列車荷載引起的道床振動(dòng)是鐵路設(shè)計(jì)的重要參考因素，采用有限元法可模擬分析列車荷載下的道床振動(dòng)，并通過試驗(yàn)對(duì)比驗(yàn)證計(jì)算模型的可靠性［26］。此外，采用子結(jié)構(gòu)方法可分析高速列車載荷下鐵軌的振動(dòng)，以及列車提速時(shí)鐵軌的動(dòng)力響應(yīng)［29］。本文將通過有限元法構(gòu)造鐵路過渡段無砟道床模型，分析道床在循環(huán)荷載作用下的動(dòng)力響應(yīng)。

鐵路道床過渡段是一個(gè)復(fù)雜的結(jié)構(gòu)體系，既包括散體道砟材料，又包括連續(xù)體道床材料。采用離散元-有限元耦合算法可在不同的區(qū)域分別采用離散單元法和有限單元法，充分發(fā)揮兩種算法的優(yōu)勢(shì)。離散元和有限元之間的參數(shù)傳遞是耦合過程中的關(guān)鍵問題，可以通過引入過渡層實(shí)現(xiàn)兩種算法的平滑過渡［30－31］，也可以采用罰函數(shù)法實(shí)現(xiàn)不同區(qū)域之間力學(xué)信息的傳遞［32］。此外，自適應(yīng)耦合算法可在局部發(fā)生破壞時(shí)，將部分有限元轉(zhuǎn)化成離散單元，其他區(qū)域仍采用有限元法計(jì)算，提高了計(jì)算效率［33－34］。李錫夔等［35］基 于 連 接尺度方 法，建 立了細(xì)尺度上離散顆粒模型和粗尺度上Cosserat連續(xù)體模型的離散元和有限元連接尺度模型。雖然離散元-有限元耦合算法在不斷發(fā)展和完善，但在鐵路道床有砟-無砟過渡區(qū)域中的應(yīng)用尚未開展。

本文將建立一種分析循環(huán)載荷作用下，鐵路道床有砟-無砟過渡段的三維DEM-FEM耦合模型。在道砟顆粒嵌入無砟道床的情況下，探討有砟-無砟道床過渡段的動(dòng)力特性，分析有砟和無砟道床的沉降量差異。

2 鐵路道床有砟-無砟過渡段的DEM-FEM耦合模型

2.1 有砟道床的離散單元法

在采用離散元方法模擬有砟道床時(shí)，需要生成道砟顆粒集合體模型。道砟顆粒外形復(fù)雜，每個(gè)顆粒的外形均不相同，其幾何形狀對(duì)道床的運(yùn)動(dòng)有顯著影響。為建立真實(shí)的有砟道床模型，本文采用鑲嵌球體單元模擬道砟顆粒。鑲嵌顆粒的內(nèi)部球體單元之間無法計(jì)算接觸力，其在運(yùn)動(dòng)過程中也不會(huì)破碎，可視為剛體。

圖1 道砟顆粒及其鑲嵌球體離散元模型Fig.1 Ballast particle and its DEM with overlapped spheres

為準(zhǔn)確生成非規(guī)則道砟顆粒的幾何外形，確定道砟顆粒的邊界，采用三維激光掃描儀確定道砟顆粒的形貌，如圖1所示，采用道砟顆粒邊界，可對(duì)其內(nèi)部進(jìn)行球體顆粒填充，并通過顆粒在邊界內(nèi)膨脹的方法得到具有鑲嵌效果的道砟顆粒單元。該方法可使模型更加光滑，并接近真實(shí)的道砟顆粒形態(tài)，如圖1（d）所示。在非規(guī)則形態(tài)道砟顆粒的球體鑲嵌構(gòu)造中，球體單元的粒徑越小，構(gòu)造的道砟顆粒越真實(shí)，但所采用的球體單元也越多，從而降低計(jì)算效率。考慮以上因素，每個(gè)道砟顆粒采用10～20個(gè)非均一粒徑球體單元進(jìn)行隨機(jī)構(gòu)造，將鑲嵌球體單元設(shè)為剛體顆粒，在描述道砟顆?，F(xiàn)象時(shí)具有一定的局限性，但仍可以有效地表征道砟顆粒的非規(guī)則幾何形態(tài)和受力自鎖特性。道砟顆粒間的接觸力基于鑲嵌球體單元，采用非線性粘彈性模型進(jìn)行計(jì)算［1，18］。

2.2 無砟道床的有限單元法

無砟道床由混凝土或?yàn)r青混合料組成，可采用有限單元法進(jìn)行計(jì)算。目前，有很多可以求解連續(xù)結(jié)構(gòu)動(dòng)力學(xué)特性的離散差分格式，如中心差分法、Wilson-θ方法和Newmark方法等。為提高計(jì)算效率，本文采用隱式積分的Newmark方法對(duì)無砟道床部分進(jìn)行動(dòng)力學(xué)分析，采用20節(jié)點(diǎn)等參單元構(gòu)造無砟道床。

常用的Newmark法的形式為

式中 Fk＋1為tk＋1時(shí)刻結(jié)構(gòu)的載荷，δ和α均為Newmark法的計(jì)算參數(shù)，當(dāng)滿足δ≥0.5和α≥0.25（0.5＋δ）2時(shí)，計(jì)算可無條件穩(wěn)定，通常取δ＝0.5和α＝0.25。

將式（2，3）代入式（1）可得

將計(jì)算得到的xk＋1代入式（2，3），可分別得到和。

2.3 有砟-無砟過渡段的DEM-FEM耦合算法

在離散元-有限元耦合過程中，關(guān)鍵問題是耦合過渡區(qū)離散元和有限元間的邊界條件處理，即有限元和離散元耦合區(qū)域內(nèi)力學(xué)參數(shù)的合理傳遞，使耦合過渡區(qū)域內(nèi)的力學(xué)和物理性質(zhì)實(shí)現(xiàn)平滑過渡。由于耦合區(qū)域是離散元和有限元的重疊區(qū)域，離散元求出的道砟顆粒的合力和作用位置作為有限元的力邊界條件；同時(shí)，有限單元和離散單元需要保持位移協(xié)調(diào)，即耦合區(qū)域內(nèi)的離散單元位移在每個(gè)時(shí)刻應(yīng)與有限單元內(nèi)的節(jié)點(diǎn)保持一致，從而為離散元計(jì)算提供位移邊界條件。

（1）耦合區(qū)域內(nèi)的位移協(xié)調(diào)

有砟-無砟過渡段耦合區(qū)域內(nèi)的位移分為有限元域位移uFEM和離散元域位移uDEM，其中uFEM由有限元域內(nèi)實(shí)體單元節(jié)點(diǎn)的獨(dú)立位移δi對(duì)形函數(shù)插值所得，即

形函數(shù)矩陣Ni為20個(gè)節(jié)點(diǎn)在每個(gè)鑲嵌單元質(zhì)心坐標(biāo)上所取形函數(shù)的值Nαi，即

耦合區(qū)域內(nèi)離散元的位移表示為

由于有砟-無砟耦合區(qū)域內(nèi)離散元和有限元的位移不協(xié)調(diào)，離散元的位移可以分解為隨有限元一起運(yùn)動(dòng)的位移uFEM和單獨(dú)運(yùn)動(dòng)位移uSEP，即uDEM＝uFEM＋uSEP。為保證耦合區(qū)域內(nèi)兩種方法的位移協(xié)調(diào)，提出一種嵌入單元模型，即在耦合區(qū)域內(nèi)離散元外圍加入一層嵌入單元?？紤]實(shí)際的力學(xué)過程，嵌入單元的單獨(dú)位移uSEP＝0，因此uDEM＝uFEM，即耦合區(qū)域內(nèi)的嵌入單元顆粒的位移完全由有限元的位移決定，從而保證區(qū)域內(nèi)的位移協(xié)調(diào)。即

有砟-無砟過渡段離散元-有限元耦合的計(jì)算模型如圖2所示。圖2（a）中嵌入到有限元模型中的顆粒單元形成耦合過渡層，嵌入單元的位置關(guān)系和相互作用如圖2（b）所示。采用上述道砟顆粒嵌入無砟道床板的結(jié)構(gòu)形式，有助于提高有砟-無砟過渡界面上的摩擦作用，從而增強(qiáng)有砟-無砟過渡段剛度變化的連續(xù)性。

（2）離散元-有限元耦合區(qū)域內(nèi)的參數(shù)傳遞

耦合區(qū)域內(nèi)可以通過顆粒間的接觸模型來確定離散元嵌入單元顆粒的合力，并將該合力作為有限元的力邊界條件進(jìn)行計(jì)算。在完成離散元每個(gè)時(shí)間步的計(jì)算后，可以確定其合力FDEMj，然后將其作為有限元的外載荷累加到有限元的節(jié)點(diǎn)上，得到有限元的等效節(jié)點(diǎn)載荷。單個(gè)嵌入單元與有限單元耦合如圖3所示。依據(jù)能量守恒原理可以得到有限元等效節(jié)點(diǎn)力與嵌入單元所受合力的關(guān)系式，即

式中 Fnodal，i為有限單元節(jié)點(diǎn)上的等效外力，Ni為嵌入單元形心處的形函數(shù)值，j為嵌入單元編號(hào)，n為嵌入單元總個(gè)數(shù)，F(xiàn)DEMj為每個(gè)嵌入單元所受自由離散單元的合力。

（3）局部坐標(biāo)求解

在耦合區(qū)域的參數(shù)傳遞中，需計(jì)算嵌入單元坐標(biāo)的形函數(shù)。局部坐標(biāo)通過求解方程（13）確定，即

圖2 離散元-有限元耦合模型示意圖Fig.2 Sketch of coupled DEM-FEM method

圖3 從坐標(biāo)系 （x，y，z）到坐標(biāo)系 （ε，η，ξ）的一一映射Fig.3 Mapping between global and local coordinate systems

式中xi，yi和zi為單元節(jié)點(diǎn)在總體坐標(biāo)內(nèi)的坐標(biāo)，Ni為用局部坐標(biāo)表示的插值函數(shù)。

令

假設(shè)函數(shù)f（ε，η，ξ）在點(diǎn) （ε0，η0，ξ0）的某一鄰域內(nèi)連續(xù)且有2階連續(xù)偏導(dǎo)數(shù)，（ε0＋l，η0＋m，ξ0＋n）為此鄰域內(nèi)任意一點(diǎn)，則有

于是式（14）可近似表示為

即

同理可得

求解式（17～19）組成的方程組，當(dāng) （fξgεgξfε）（gηhε－h(huán)ηgε）－（gξhε－h(huán)ξgε）（fηgε－gηfε）≠0時(shí)，有

至此，求得嵌入單元的局部坐標(biāo)。根據(jù)求得的嵌入單元合力求解有限元的邊界力，完成離散元到有限元的力參數(shù)傳遞，進(jìn)而得到有限單元的變形。根據(jù)有限元網(wǎng)格節(jié)點(diǎn)位移以及形函數(shù)更新嵌入單元的位置坐標(biāo)，嵌入單元的運(yùn)動(dòng)作為位移邊界條件傳入離散元的計(jì)算，由此完成離散元和有限元之間的耦合過程。

3 有砟-無砟過渡段的DEM-FEM耦合模擬

3.1 有砟-無砟過渡段的DEM-FEM耦合模型

采用DEM-FEM耦合算法建立列車循環(huán)荷載作用下鐵路道床有砟-無砟過渡段的計(jì)算模型。該模型包括無砟混凝土道床、枕木和有砟散體道床三部分，如圖4所示?？紤]高速鐵路軌道的對(duì)稱性，并節(jié)省計(jì)算時(shí)間，本文取鐵路道床的一半進(jìn)行建模分析，對(duì)于無砟道床的鋼筋混凝土結(jié)構(gòu)，采用線彈性模型分析其位移和變形沉降。在兩種結(jié)構(gòu)交界面上采用有砟顆粒嵌入無砟道床的方式，以增加有砟道床的咬合力。分別選取1＃和2＃位置作為列車荷載下的有砟道床和無砟道床沉降觀測(cè)點(diǎn)。

在有砟-無砟過渡段的DEM-FEM耦合模型中，離散元模型采用鑲嵌單元建模。有砟道床計(jì)算模型的長(zhǎng)度x＝0.4m，高度z＝0.38m，厚度y＝0.23m，整個(gè)有砟道床由736個(gè)道砟單元顆粒、4291個(gè)球體單元組成。在離散元-有限元耦合區(qū)域內(nèi)，共有46個(gè)嵌入單元，包括264個(gè)顆粒單元。有限元模型由道床和軌枕組成。模型采用20節(jié)點(diǎn)等參單元，劃分為775個(gè)單元，總共4018個(gè)節(jié)點(diǎn)。采用線彈性模型建模分析，彈性模量為35GPa，泊松比為0.22。模型左側(cè)為自由表面，右側(cè)施加x方向的約束，前面施加y方向的約束，下側(cè)施加z方向的約束。將列車荷載簡(jiǎn)化到軌枕上，在軌枕上施加z方向的荷載。由于離散元的時(shí)間步長(zhǎng)遠(yuǎn)小于有限元的時(shí)間步長(zhǎng)，為了提高計(jì)算效率，有限元的時(shí)間步長(zhǎng)取為離散元時(shí)間步長(zhǎng)的400倍。離散元和有限元的主要計(jì)算參數(shù)列入表1和表2?？紤]道砟顆粒在往復(fù)荷載作用下因不斷摩擦而趨于光滑，因此在離散元計(jì)算中，道砟顆粒間的摩擦系數(shù)隨往復(fù)荷載次數(shù)的增加而不斷降低。

圖4 鐵路道床有砟-無砟過渡段有限元-離散元耦合模型Fig.4 Coupled discrete-finite element model of transition zone

3.2 計(jì)算結(jié)果分析

為分析高速列車在250km／h的速度下通過有砟-無砟過渡段時(shí)的動(dòng)力響應(yīng)，在枕木上施加幅度23kN，頻率0.36Hz的循環(huán)荷載。荷載的時(shí)程曲線如圖5（a）所示，循環(huán)90次，時(shí)間為32.4s。一次循環(huán)中列車荷載會(huì)有四個(gè)峰值，時(shí)程曲線如圖5（b）所示。

表1 道砟離散元計(jì)算參數(shù)Tab.1 Computational parameters of ballast in DEM

表2 無砟道床有限元計(jì)算參數(shù)Tab.2 Material parameters of track slab in FEM

圖5 施加的列車循環(huán)荷載Fig.5 Applied cyclic load of train

在列車循環(huán)荷載的作用下，有砟-無砟過渡段1＃和2＃觀測(cè)點(diǎn)的沉降過程如圖6所示。有砟道床在循環(huán)加載的初期，道床的沉降量變化明顯并且一直增大，之后隨著循環(huán)次數(shù)的增加，沉降量趨于平穩(wěn)并且緩慢下降，最終沉降量達(dá)到15mm左右，如圖6（a）所示。無砟道床則在前幾次加載后很快趨于平穩(wěn)，之后保持一個(gè)相對(duì)穩(wěn)定的狀態(tài)，最終總體沉降量較少，約為0.005mm，如圖6（b）所示。將兩種道床結(jié)構(gòu)對(duì)比發(fā)現(xiàn)，無砟道床很快達(dá)到穩(wěn)定狀態(tài)，而有砟道床則需要較長(zhǎng)時(shí)間趨于穩(wěn)定，同時(shí)在穩(wěn)定之后兩種結(jié)構(gòu)會(huì)有一個(gè)明顯的沉降差。這是由于在加載之前有砟道床的道砟石子排列比較松散，顆粒之間有較大的空隙，前幾次的循環(huán)加載可以使有砟道床的顆粒重新排列，顆粒之間空隙變小，道床進(jìn)一步壓實(shí)，產(chǎn)生較大的沉降。對(duì)于無砟道床，是由混凝土或?yàn)r青混合料構(gòu)成的軌道道床結(jié)構(gòu)，在循環(huán)加載過程中，道床的楊氏模量保持不變，剛度較大，道床可以較快達(dá)到平穩(wěn)。

循環(huán)荷載作用下，枕木的荷載-沉降量曲線如圖7所示。由于一個(gè)時(shí)程內(nèi)列車荷載具有兩個(gè)峰值，并且在兩個(gè)峰值之間有一次卸載過程，荷載-沉降量曲線也表現(xiàn)出相應(yīng)的規(guī)律。在前期循環(huán)加載過程中，有砟道床沉降量明顯，表明整體道床發(fā)生了明顯的塑性變形；隨著加載次數(shù)增加，顆粒進(jìn)一步壓實(shí)，整體道床平緩沉降，單次循環(huán)加卸載曲線所包圍的面積也逐漸減少。有砟道床載荷-沉降量曲線的斜率逐步增大并趨于穩(wěn)定，表明道砟材料的宏觀彈性模量逐漸增大并趨于穩(wěn)定，整體道床逐漸表現(xiàn)出非線性彈性材料的性質(zhì)。

圖6 鐵路道床過渡段沉降量曲線Fig.6 Settlement of ballasted-ballastless track transition zone

圖7 循環(huán)載荷作用下載荷-沉降量曲線Fig.7 Force-settlement curve under cyclic load

圖8 有砟道床力鏈分布和無砟道床應(yīng)力分布Fig.8 Force chain distribution and stress contour

圖9 有砟道床力鏈分布和無砟道床位移云圖Fig.9 Force chain distribution and displacement contour

循環(huán)加載達(dá)到峰值時(shí)，有砟道床中的力鏈分布和無砟道床中的應(yīng)力分布如圖8所示。有砟道床的力鏈主要分布在枕木的下方，呈梯形分布，從而將枕木上受到的載荷分散在道床上，荷載主要由道砟顆粒骨架承受。無砟道床的應(yīng)力主要分布在道床的下部區(qū)域，應(yīng)力分布區(qū)域與有砟道床道砟骨架相對(duì)應(yīng)。無砟道床左側(cè)邊界為自由表面，僅受到有砟道床的約束，右側(cè)邊界受水平方向的位移約束，有砟道床的約束效果低于無砟道床的約束效果，因此右側(cè)應(yīng)力略大于左側(cè)應(yīng)力。從力鏈分布可以看出，梯形力鏈骨架并不是左右對(duì)稱的，右面的力鏈分布較大，主要是由于在耦合區(qū)域嵌入單元受有限元約束，相當(dāng)于增加了有砟道床和無砟道床交界面上的粗糙度，導(dǎo)致耦合區(qū)域的道砟顆粒具有較大的咬合力，從而限制道床過渡段附近道砟顆粒的流動(dòng)，一定程度上減少兩種道床的沉降差。

循環(huán)加載達(dá)到峰值時(shí)，無砟道床x方向的位移主要集中在道床的左半部分，位移變形與道砟力鏈分布相對(duì)應(yīng)，如圖9（a）所示。由于道砟顆粒荷載作用于無砟道床左側(cè)表面，導(dǎo)致道床左側(cè)表面x方向位移變形較大。無砟道床整體位移較大的區(qū)域位于枕木下方，同時(shí)道床左側(cè)位移略大于右側(cè)位移，如圖9（b）。這也是由于無砟道床左側(cè)位移約束略小于右側(cè)位移約束，與道床內(nèi)下部應(yīng)力分布相對(duì)應(yīng)。

4 結(jié) 論

針對(duì)鐵路道床有砟-無砟過渡段在循環(huán)載荷下的沉降特性，發(fā)展了離散元-有限元耦合方法。對(duì)道砟顆粒和無砟道床分別采用離散元和有限元方法，并針對(duì)有砟-無砟過渡段的道砟顆粒嵌入模式，在離散元-有限元耦合區(qū)域發(fā)展了一種參數(shù)傳遞的耦合算法。計(jì)算結(jié)果表明，有砟道床的沉降量顯著，主要是由于道砟顆粒的重新排列，在加載后期趨于穩(wěn)定；無砟道床沉降較小，并在前幾次加載沉降后很快達(dá)到平穩(wěn)。在加載前期，有砟道床整體產(chǎn)生塑性變形，但隨著加載的持續(xù)，道砟材料的宏觀彈性模量逐漸增大并趨于穩(wěn)定，整體道床逐漸表現(xiàn)出非線性彈性材料的性質(zhì)。有砟道床力鏈主要分布在枕木的下方，荷載主要由道砟顆粒骨架承受。有砟-無砟過渡段嵌入單元的存在限制了道砟顆粒的運(yùn)動(dòng)，力鏈呈非對(duì)稱梯形分布，在一定程度上減少了兩種道床的沉降差異。無砟道床的應(yīng)力和位移較大區(qū)域位于枕木下方，應(yīng)力分布同樣呈梯形分布，與道砟力鏈分布趨勢(shì)相對(duì)應(yīng)。