基于離散元與有限元耦合的有砟道床沉降特性研究

王 桐，王立華，陳佳明，栗先增

(昆明理工大學(xué) 機(jī)電工程學(xué)院，云南 昆明 650500)

鐵路作為我國(guó)的重要基礎(chǔ)設(shè)施，在交通運(yùn)輸體系中占據(jù)重要地位，也是目前應(yīng)用最廣泛的運(yùn)載方式之一，對(duì)我國(guó)的經(jīng)濟(jì)建設(shè)和發(fā)展有著十分重大的意義。由于我國(guó)鐵路近幾年發(fā)展迅速，且與各個(gè)領(lǐng)域之間的交叉融合越來越多，所以對(duì)于鐵路方面的研究也愈加深入。有砟鐵路作為我國(guó)鐵路運(yùn)載的重要組成部分，是鐵路發(fā)展的基礎(chǔ)，對(duì)其深入的研究就顯得尤為重要。有砟軌道道床質(zhì)量對(duì)于保障列車運(yùn)行穩(wěn)定安全方面具有重要意義。

在對(duì)有砟鐵路進(jìn)行研究時(shí)，由于組成道床的道砟屬于離散體顆粒，而道床上部軌枕屬于連續(xù)體，所以在對(duì)二者進(jìn)行分析時(shí)不能使用單一的方法。有限元法只能從宏觀上反應(yīng)軌枕與道床的整體特性。Ishikawa[1]等通過室內(nèi)試驗(yàn)研究有砟道床在定點(diǎn)載荷和移動(dòng)載荷作用下的力學(xué)行為。其結(jié)果表明，在移動(dòng)載荷作用下有砟道床產(chǎn)生的塑性變形比定點(diǎn)位移作用下的大。Mc Dowell[2-4]等基于真實(shí)顆粒外形的最小半徑離散單元數(shù)值模型的構(gòu)造方法，用所生成的組合球單元對(duì)復(fù)雜外形的道砟顆粒進(jìn)行了模擬，并對(duì)比分析了組合單元顆粒與球形單元顆粒在荷載作用下的力學(xué)特性。井國(guó)慶[5]等基于二維離散單元法顆粒軟件Particle Flow Code 2D(PFC 2D)，建立了道砟軌枕離散單元模型，研究了列車循環(huán)荷載作用下道砟顆粒破碎與軌枕沉降的關(guān)系。Shaer[6]等采用等比例縮減尺寸的鐵軌-枕木-道床-路基模型得到整體道床的沉降與枕木和道砟的振動(dòng)加速度密切相關(guān)。高亮等[7-8]利用離散單元法，分別建立球形道砟顆粒和顆粒簇形成道砟顆粒的軌枕與散體道床的離散元模型，對(duì)比分析了兩種不同模型在循環(huán)荷載下的力學(xué)性能，并利用離散元顆粒數(shù)值模型構(gòu)造道砟顆粒模型。同時(shí)，在此基礎(chǔ)上建立了循環(huán)荷載道砟箱數(shù)值模型，研究在列車高速及重載線路條件下，道砟級(jí)配對(duì)散體道床沉降力學(xué)特性的影響情況，并從細(xì)觀角度分析了散體道床的沉降原理。但若要同時(shí)從宏觀和微觀方面對(duì)有砟道床進(jìn)行分析，就需要引入有限元和離散元耦合的方法。

針對(duì)離散元與有限元耦合方法的研究，馮春[9]將對(duì)有限元模型進(jìn)行加載達(dá)到臨界狀態(tài)時(shí)，把有限元替換為離散元研究。嚴(yán)穎等[10-12]在離散元與有限元的接觸面處提出一種新的接觸算法，實(shí)現(xiàn)有砟鐵路軌道的離散元-有限元耦合法的研究分析，又通過鑲嵌單元模擬來達(dá)到耦合的方式并通過Fortran語言編程計(jì)算有限元與離散元的耦合結(jié)果。嚴(yán)穎等依據(jù)接觸面上的能量守恒原理，確定兩種介質(zhì)間的接觸力實(shí)現(xiàn)離散元-有限元單元區(qū)間力子參數(shù)傳遞。尹超[13]利用MATLAB編制程序接口，通過間接計(jì)算實(shí)現(xiàn)離散元軟件與有限元軟件的數(shù)據(jù)計(jì)算與傳遞。

上述文獻(xiàn)都是基于耦合方法研究軌枕與散體道砟顆粒的相互作用力，需要通過外部程序來進(jìn)行有限元和離散元顆粒之間的數(shù)據(jù)傳遞。為簡(jiǎn)化數(shù)據(jù)之間傳遞的繁瑣過程，本文采用有限元軟件得到軌枕網(wǎng)格坐標(biāo)，導(dǎo)入離散元中用離散元顆粒以及鍵連接構(gòu)造軌枕模型。后采用球顆粒組合模型構(gòu)造非規(guī)則道砟，得到縮小比例工程尺度的枕木-道床離散元模型。本研究通過對(duì)軌枕施加循環(huán)正弦激振，研究了軌枕和散體道砟顆粒的相互作用和軌枕的受力、變形情況。此外，本文還研究了不同激振頻率對(duì)有砟道床力學(xué)性質(zhì)的影響，得到道床的沉降規(guī)律并探究了道床累積沉降的內(nèi)在機(jī)理。

1 有砟道床的離散元模型

1.1 典型形狀道砟模型構(gòu)建

離散單元法(Discrete Element Method，DEM)主要用于離散顆粒物料研究領(lǐng)域。非規(guī)則形狀顆粒建模耗費(fèi)大量時(shí)間，且模型精度對(duì)顆粒間的動(dòng)力學(xué)響應(yīng)特性有較大影響，故本文采用離散單元法構(gòu)建道床模型，旨在提高非規(guī)則形狀顆粒建模效率及精度。在有限元軟件中建立軌枕模型并劃分網(wǎng)格，同時(shí)獲取網(wǎng)格單元坐標(biāo)信息，并在離散元軟件中設(shè)置接觸半徑及物理半徑，然后通過工程離散單元法(Engineering Discrete Element Method，EDEM)中的應(yīng)用程序編程接口實(shí)現(xiàn)多球快速聚合顆粒建模。道砟顆粒的形態(tài)會(huì)對(duì)散體道床的動(dòng)力特性產(chǎn)生顯著的影響，構(gòu)成簇顆粒的球單元越多，其對(duì)道砟形態(tài)的描述就越真實(shí)。利用以上構(gòu)造方法可以對(duì)任意形態(tài)的道砟進(jìn)行構(gòu)造?？紤]到離散元計(jì)算效率，本文采用4～10顆小球顆粒組成的簇顆粒單元對(duì)不規(guī)則道砟顆粒進(jìn)行近似構(gòu)造，模型示意圖如圖1所示。組成簇顆粒的球顆粒之間存在很大的重疊量，且保證了道砟內(nèi)部不存在孔隙。顆粒簇在離散元模擬計(jì)算中不考慮其破碎，每個(gè)顆粒簇在運(yùn)動(dòng)過程中被視為一個(gè)整體，顆粒簇之間采用線性接觸模型去模擬顆粒之間的接觸行為。

圖1 不同形狀的鐵路道砟顆粒簇模型

正常工作的有砟鐵路使用的道砟顆粒尺寸一般在37～54 mm之間，如上圖所示粒徑范圍分別在表1道砟級(jí)配中的30～45 mm、45～60 mm、>60 mm這3個(gè)范圍內(nèi)選取。

表1 道砟級(jí)配

道砟與道砟之間，道砟與構(gòu)成軌枕的離散元顆粒間的接觸參數(shù)如表2所示

表2 道砟間的接觸參數(shù)

1.2 節(jié)點(diǎn)坐標(biāo)提取與離散元軟件顆粒的替換

在離散元中將分析對(duì)象視為充分?jǐn)?shù)量的離散介質(zhì)進(jìn)行計(jì)算，其將每個(gè)離散顆粒作為一個(gè)單元，然后求解出每一時(shí)刻中顆粒的相互作用、接觸力及離散單元的運(yùn)動(dòng)狀態(tài)。隨后，通過引入時(shí)間步長(zhǎng)的方式，依次求解出離散單元體的運(yùn)動(dòng)形態(tài)。離散單元體一般使用球形顆粒，根據(jù)實(shí)際問題的不同，將顆粒模型分為軟球模型和硬球模型兩種，本文所采用的是軟球模型。

利用離散元軟件EDEM所提供的應(yīng)用程序編譯外接口(Application Programing Interface，API)編寫插件來自定義接觸模型、顆粒生成、外部耦合及顆粒模型[14]。軌枕顆粒建模仿真實(shí)驗(yàn)中利用Visual Studio編程軟件創(chuàng)建的顆粒工廠動(dòng)態(tài)鏈接庫(kù)文件，得到顆粒ID、顆粒三維坐標(biāo)、顆粒實(shí)際物理半徑，顆粒接觸半徑等所需的模板顆粒信息，并對(duì)顆粒進(jìn)行替換?？紤]到顆粒替換得到的模型與顆粒之間的粘結(jié)鍵及仿真所需要的計(jì)算時(shí)間，文中將組成軌枕的顆粒小球?qū)嶋H物理半徑設(shè)置為1.25 cm，顆粒接觸半徑設(shè)置為1.6 cm。由離散元顆粒替換構(gòu)成的軌枕模型如圖2所示。

圖2 離散元顆粒替換構(gòu)成的軌枕

1.3 軌枕-散體道床耦合模型的建立

將道砟顆粒模型建模后，根據(jù)已有相關(guān)研究[15]和道砟顆粒模型的尺寸建立道砟箱模型。當(dāng)?shù)理娜萜鞯某叽绱笥陬w粒平均粒徑的8倍時(shí)，容器的邊界效應(yīng)可以被忽略。本文建立的道砟顆粒模型的平均粒徑不超過45 mm，容器最小要求尺寸為340 mm，建立鐵路道砟箱三維模型，其中道砟箱的尺寸為1 000 mm×700 mm×600 mm(長(zhǎng)×寬×高)。考慮到仿真所面臨的時(shí)間和內(nèi)存問題，本文僅建立由一根軌枕和有砟道床組成的離散元模型。簡(jiǎn)化模型中有砟道床及軌枕的尺寸與實(shí)際道床及軌枕的尺寸比例為1∶8，軌枕由簡(jiǎn)單的長(zhǎng)方體簡(jiǎn)化代替。

在有限元軟件中對(duì)節(jié)點(diǎn)數(shù)據(jù)進(jìn)行處理后，在離散元軟件中，設(shè)置單獨(dú)母體小球，在到達(dá)需要放置的軌枕位置后，通過顆粒替換將離散元顆粒按照軌枕網(wǎng)格劃分坐標(biāo)在離散元軟件中排布組合成新軌枕。建立符合要求的有砟道床離散元模型，如圖3所示。

圖3 軌枕與道床的三維離散元模型

完成軌枕與道床的離散元模型建立后，為其添加接觸參數(shù)，確定軌枕顆粒和道砟顆粒運(yùn)動(dòng)邊界條件。由于對(duì)道砟顆粒的研究需要確定顆粒間準(zhǔn)確的接觸力變化情況，且暫不考慮顆粒破碎的情況，因此采用Hertz-Mindlin(no slip)無滑動(dòng)接觸模型。由于軌枕由散體顆粒構(gòu)成，需要設(shè)置bond鍵來連接各個(gè)離散元顆粒以達(dá)到整體作為軌枕的效果。軌枕采用Hertz-Mindlin with bonding基礎(chǔ)模型。用于模擬破碎、斷裂等問題時(shí)，則采用小顆粒粘結(jié)成大塊物料。外力作用下顆粒間粘結(jié)力會(huì)發(fā)生破壞，從而產(chǎn)生破碎及斷裂效果，但本文不考慮軌枕的劣化和破損問題，故bond鍵強(qiáng)度設(shè)置較大。最終所確定的道砟顆粒與顆粒，顆粒與軌枕間的恢復(fù)系數(shù)、靜摩擦系數(shù)和動(dòng)摩擦系數(shù)如表3所示。

表3 軌枕粘結(jié)模型物理參數(shù)

由于無法直接對(duì)軌枕施加激振載荷，因此需要將作用在軌枕上的載荷轉(zhuǎn)化為不同的形式施加在軌枕上。將軌枕所受到的激振力轉(zhuǎn)化成位移的形式輸入到軌枕上，即在仿真過程中通過在軌枕兩側(cè)施加激振源，對(duì)激振源設(shè)置位移函數(shù)，水平激振力屬于正弦激振力，存在頻率和幅值。因此將水平激振力轉(zhuǎn)化成頻率和位移幅值的形式施加在軌枕上，并通過改變激振頻率來對(duì)不同激振頻率下軌枕的沉降狀態(tài)進(jìn)行研究。對(duì)軌枕施加激振，完成穩(wěn)定過程仿真，以離散元軟件為主體來實(shí)現(xiàn)耦合過程。

2 穩(wěn)定作用下軌枕受力及道床動(dòng)力特性分析

2.1 軌枕下方不同區(qū)域的速度變化趨勢(shì)

在軌枕下方的道砟在穩(wěn)定作用下，道砟的移動(dòng)方向整體向下，其中以軌枕下方的道砟顆粒尤為明顯。為了直觀分析不同載荷頻率對(duì)道床的沖擊程度，如圖4所示，給出了載荷頻率為36 Hz時(shí)枕木受到穩(wěn)定作用，道床中間厚度1 cm垂向截面區(qū)域內(nèi)所有道砟顆粒的速度分布圖?？煽闯稣砟緝蓚?cè)道砟顆粒的速度均較小，說明穩(wěn)定載荷對(duì)軌枕下方道砟起主要作用。道砟顆粒的速度幅值呈現(xiàn)近似扇形的衰減，且軌枕下方道砟顆粒具有更加明顯的垂直向下運(yùn)動(dòng)趨勢(shì)，越靠近軌枕中心的道砟顆粒，其垂向速度響應(yīng)越大。

軌枕在受到激振時(shí)，向左右兩側(cè)產(chǎn)生微小位移，在軌枕向左向右移動(dòng)時(shí)，砟肩部分的道砟隨著軌枕產(chǎn)生同向位移。由于軌枕下方道砟顆粒始終為向下移動(dòng)，所以在激振影響中，道砟不斷向下移動(dòng)，達(dá)到重新排布效果。整個(gè)過程中，軌枕下方道砟不斷下沉，使上方軌枕下沉。

2.2 軌枕的力學(xué)特性分析

選取軌枕底部顆粒，提取該顆粒在穩(wěn)定過程中的位移量來得到軌枕整體的沉降量，如圖5所示，顆粒在穩(wěn)定過程中的沉降量為20 mm左右。由圖中可以看出，在對(duì)不同激振頻率下道床穩(wěn)定作業(yè)過程的仿真中，隨著穩(wěn)定作業(yè)的進(jìn)行道床軌枕沉降量逐漸增大，并且隨著激振頻率的增加軌枕總沉降量逐漸增大，在36 Hz時(shí)達(dá)到最大值；當(dāng)激振頻率繼續(xù)增大時(shí)，沉降量又逐漸減小。

圖5 軌枕沉降量與時(shí)間的關(guān)系

3 設(shè)計(jì)道砟箱試驗(yàn)與模型驗(yàn)證

3.1 道砟箱試驗(yàn)的設(shè)計(jì)

根據(jù)前文軌枕和道砟顆粒的選型設(shè)計(jì)搭建試驗(yàn)平臺(tái)，如圖6所示，建立與仿真同比例的道砟箱試驗(yàn)。文中所使用的箱體尺寸(長(zhǎng)×寬×高)為1 000 mm×700 mm×550 mm，所用材料為3 mm厚鋼板。通過安裝偏心振動(dòng)電機(jī)在水泥軌枕上部來模擬穩(wěn)定過程，并使用調(diào)頻器達(dá)到不同激振頻率的效果。該振動(dòng)電機(jī)可以提供水平方向的激振力，在垂直下壓力方面只考慮軌枕和振動(dòng)電機(jī)的質(zhì)量。

圖6 道砟箱試驗(yàn)

為研究試驗(yàn)中模擬穩(wěn)定作業(yè)軌枕沉降量的變化規(guī)律，通過激光測(cè)距儀可以準(zhǔn)確地得出動(dòng)力穩(wěn)定作業(yè)過程中道床的垂向位移和軌枕的水平激振幅值以及水平激振頻率與軌枕沉降的關(guān)系。將激光測(cè)距儀安裝在軌枕正上方距離3 cm處，同時(shí)避免與道砟箱接觸，防止振動(dòng)對(duì)測(cè)量產(chǎn)生影響。在進(jìn)行測(cè)量時(shí)，測(cè)距儀的激光照射到軌枕的振動(dòng)電機(jī)安裝鋼板上，當(dāng)模擬穩(wěn)定作業(yè)過程開始時(shí)，激光測(cè)距儀發(fā)出紅外線照射到軌枕上方的平滑鋼板上進(jìn)行數(shù)據(jù)采集，得到每一個(gè)時(shí)刻軌枕與激光發(fā)射口的距離大小。將初始位置設(shè)置為零點(diǎn)，即可得到軌枕的沉降量。

圖7 激光測(cè)距儀實(shí)物圖

動(dòng)力穩(wěn)定車在對(duì)劣化道床進(jìn)行穩(wěn)定作業(yè)的過程中，水平激振頻率的工作頻率范圍是0～45 Hz，劉暢[16]得到的動(dòng)力穩(wěn)定作業(yè)過程中最優(yōu)激振頻率為30 Hz；項(xiàng)永志[17]通過仿真得到動(dòng)力穩(wěn)定作業(yè)最優(yōu)激振頻率為31～33 Hz。本文建立的試驗(yàn)平臺(tái)為小型道砟箱試驗(yàn)平臺(tái)，由于激振電機(jī)的頻率限制，激振頻率在30 Hz以下的激振作用振動(dòng)比較微弱，所對(duì)應(yīng)的穩(wěn)定效果不明顯，因此最終選擇設(shè)置梯度遞增的激振頻率，分別為30 Hz、32 Hz、34 Hz、36 Hz、38 Hz和40 Hz，共6組。由于偏心振動(dòng)電機(jī)的額定頻率為50 Hz，無法直接改變其頻率，因此通過變頻器來調(diào)節(jié)電機(jī)的頻率，達(dá)到本文對(duì)于不同頻率的試驗(yàn)要求。

3.2 數(shù)據(jù)的提取與驗(yàn)證

在比例縮小模型實(shí)驗(yàn)中，使用激光測(cè)距儀測(cè)量軌枕下沉量，與上文中的軌枕下沉量進(jìn)行對(duì)比分析。為方便統(tǒng)計(jì)數(shù)據(jù)和計(jì)算，在仿真模擬中選取如圖8所示軌枕底面小球顆粒編號(hào)11 316作為分析對(duì)象，提取此處顆粒在穩(wěn)定過程中的位移來進(jìn)行數(shù)據(jù)對(duì)比分析。

圖8 道砟顆粒選取位置

提取此位置道砟顆粒在穩(wěn)定過程中的位移即軌枕的沉降量，如圖9所示得到仿真過程中軌枕顆粒的位移情況，可以發(fā)現(xiàn)在穩(wěn)定過程中，軌枕的沉降量隨著激振頻率的增加而增大，當(dāng)激振頻率為36 Hz時(shí)達(dá)到最大值，之后隨著頻率的增加逐漸減小。

由圖9可知激振頻率不變時(shí)，動(dòng)力穩(wěn)定作業(yè)下道床沉降趨勢(shì)與循環(huán)載荷作用產(chǎn)生的道床沉降趨勢(shì)基本相同。由于道砟顆粒相互接觸時(shí)具有一定的彈性，因此道床沉降過程具有一定的反復(fù)，試驗(yàn)結(jié)果和文獻(xiàn)[16]中動(dòng)力穩(wěn)定下的道床沉降趨勢(shì)基本相同。這是因?yàn)榈来渤跗诒容^松散，密實(shí)度低，穩(wěn)定作業(yè)下道砟顆粒之間迅速密實(shí)，沉降急劇增大，隨著穩(wěn)定作業(yè)的進(jìn)行，道床逐漸密實(shí)，穩(wěn)定性逐漸提高，道床沉降量逐漸變緩；當(dāng)激振頻率增加時(shí)，道床沉降量呈現(xiàn)先增大后減小的趨勢(shì)，在36 Hz時(shí)達(dá)到最大，其道床的沉降穩(wěn)定效果最好。

圖9 道床沉降趨勢(shì)圖

將實(shí)驗(yàn)得到的軌枕隨時(shí)間變化的沉降量變化情況繪制成軌枕沉降量相對(duì)于時(shí)間變化關(guān)系圖像，并繪制經(jīng)過仿真得出的軌枕沉降與時(shí)間關(guān)系圖像對(duì)比實(shí)驗(yàn)及仿真結(jié)果。如圖10所示，將離散元有限元耦合建模與試驗(yàn)獲取的數(shù)據(jù)結(jié)果進(jìn)行對(duì)比可知，試驗(yàn)獲得的軌枕沉降曲線由于在道砟箱堆放道砟時(shí)道砟顆粒之間空隙較大堆疊稀松，所以在激振在初期沉降量上升較快。試驗(yàn)與仿真的沉降趨勢(shì)基本一致，說明本文建立的有砟道床仿真分析模型是可靠的，所選取的參數(shù)正確，所確定的邊界條件合理，仿真分析結(jié)果有效。

圖10 實(shí)驗(yàn)與仿真對(duì)比圖

經(jīng)過對(duì)軌枕在36 Hz激振力作用下仿真得到的軌枕沉降量與時(shí)間變化可看出其趨勢(shì)變化規(guī)律形如指數(shù)函數(shù)。為了方便研究軌枕沉降量變化規(guī)律并得到普遍性結(jié)論，文中對(duì)仿真得到的軌枕沉降量變化趨勢(shì)進(jìn)行指數(shù)函數(shù)擬合，指數(shù)函數(shù)階數(shù)為2，擬合結(jié)果如圖11所示，則36 Hz激振力影響下軌枕沉降量與時(shí)間的演化規(guī)律大致服從圖中擬合過后得到的曲線發(fā)展規(guī)律。

分析處理得到道床累計(jì)下沉量隨時(shí)間變化的函數(shù)模型方程為

f(x)=-106×e-0.006 388x+106×e-0.000 638 7x

(1)

式中選取的常數(shù)項(xiàng)為最佳數(shù)值，即擬合度最高的沉降量計(jì)算式。當(dāng)常數(shù)項(xiàng)在符合要求的區(qū)間時(shí)，可對(duì)相應(yīng)的激振頻率進(jìn)行擬合，可靠率達(dá)95%。

圖11 擬合曲線與仿真曲線對(duì)比圖

4 結(jié)束語

文中通過有限元軟件建立了縮小尺寸軌枕模型，得到網(wǎng)格劃分后的節(jié)點(diǎn)坐標(biāo)，在離散元軟件中通過顆粒替換將軌枕單元顆粒按照節(jié)點(diǎn)坐標(biāo)排列，得到有砟軌道結(jié)構(gòu)的三維離散元模型。建立相應(yīng)道砟箱模型，分析在穩(wěn)定載荷下軌枕對(duì)道砟顆粒的作用情況、軌枕的受力情況和應(yīng)力分布情況。

本研究根據(jù)仿真模型搭建了試驗(yàn)平臺(tái)，綜合分析離散單元法在仿真道床實(shí)際狀態(tài)的優(yōu)勢(shì)和特點(diǎn)，通過將模擬仿真分析結(jié)果與實(shí)際試驗(yàn)數(shù)據(jù)對(duì)比，驗(yàn)證了此耦合方法在穩(wěn)定作業(yè)過程中仿真道床實(shí)際狀態(tài)的可行性和有效性，并利用離散單元法進(jìn)一步深入研究。本研究發(fā)現(xiàn)軌枕在鋼軌與道床之間，上部穩(wěn)定車傳遞過來的激振效果通過軌枕來影響道床，且與下部結(jié)構(gòu)散體道砟顆粒直接接觸。軌枕正下方道砟顆粒振動(dòng)速度均隨著道床深度的增加呈扇形迅速衰減，且隨著穩(wěn)定作業(yè)的進(jìn)行，軌枕下方道砟顆粒不斷向下移動(dòng)，當(dāng)激振頻率為36 Hz時(shí)，道床沉降量最大。利用EDEM軟件進(jìn)行了作業(yè)過程的模擬仿真并提取了相關(guān)數(shù)據(jù)，得到了不同作業(yè)工況下任意時(shí)刻道床下沉量等相關(guān)數(shù)值，并給出了沉量隨時(shí)間變化的函數(shù)模型方程。