孫珊珊， 嚴 穎， 趙春發(fā)， 季順迎

(1. 大連理工大學 工業(yè)裝備結構分析國家重點實驗室，遼寧 大連 116023； 2. 大連交通大學 土木與安全工程學院，遼寧 大連 116028； 3. 西南交通大學 牽引動力國家重點實驗室，四川 成都 610031)

有砟軌道是當前鐵路軌道的主要結構形式，在鐵路運輸向高速和重載方向發(fā)展過程中，其力學性能不斷改進[1-2]。鐵路道床在列車往復荷載作用下因道砟顆粒的重新排列和破碎粉化發(fā)生沉降，導致整體彈性的降低，影響其使用性能[3]。

相關學者對材料的力學性能開展了系統(tǒng)研究，主要集中在累積沉降量和有效彈性等方面。道砟材料的三軸試驗結果表明，在密實過程中道砟材料減震和排水性能均有所降低[4-5]。在高頻重載條件下，道砟材料的沉降和道砟顆粒破壞現(xiàn)象尤為突出[6]。在往復荷載作用下，道砟顆粒的破碎是引起道床變形的重要原因[7]。將鐵路道床假設為連續(xù)體，用有限元方法可分析其動力特性和變形規(guī)律[8-9]，但僅限小變形行為，且不能描述道砟顆粒的細觀受力特征。顆粒的幾何形狀及級配對其宏觀力學行為的影響顯著[10-11]。

為真實模擬道砟材料的幾何形態(tài)和力學行為，基于非連續(xù)材料的離散單元模型迅速發(fā)展[12-14]。在顆粒材料的離散元數(shù)值模擬中，采用最普遍的是球形顆粒單元。對于非規(guī)則形態(tài)的道砟顆粒，球形顆粒很難有效地模擬道砟材料的力學行為[15]，可采用球形顆粒的不同黏接或鑲嵌組合方式構造相應的道砟顆粒[14,16-17]，模擬道砟顆粒在外荷載作用下的破碎現(xiàn)象，對有砟道床的沉降特性進行離散元分析[7,17-18]。在球形顆粒的組合中，通過增加球形顆粒單元的數(shù)量可更精確地模擬非規(guī)則顆粒，甚至可采用上千個球形顆粒單元構造一個道砟塊石[16]。但顆粒數(shù)目的增加急劇降低離散元的計算效率，不利于開展大規(guī)模計算。

為模擬非球形顆粒的幾何形態(tài)，發(fā)展了不同的離散單元模型。采用二次或超二次曲面方程可定義不同形態(tài)的顆粒單元[19-20]。對類似道砟材料的多棱邊顆粒，多面體單元能夠合理描述顆粒的幾何形態(tài)以及顆粒間的相互作用，準確模擬其宏觀力學特征[21]。多面體單元的作用力計算需要繁瑣的接觸判斷和鄰居搜索，嚴重影響其計算效率和工程應用性。采用Minkowski Sum方法構造擴展多面體顆粒模型，準確快速地模擬非規(guī)則顆粒間的相互作用。在離散單元法的數(shù)值研究中，這種非規(guī)則顆粒模型已成功用于二維的擴展多邊形[22-23]和三維的擴展多面體計算[24-26]。擴展多面體單元的主要特點是在多面體單元的棱邊及外表面設置若干個球形顆粒，通過改變球形顆粒的直徑靈活控制多面體棱邊及棱角的尖銳度[27]。該模型可有效簡化多面體單元的接觸判斷過程，在滿足計算精度的條件下有效地提高離散元計算效率。

本文采用擴展多面體模型對有砟道床在往復荷載作用下的動力過程進行離散元分析，確定其累積沉降量和有效模量隨荷載往復次數(shù)的演化規(guī)律。

1 道砟顆粒的擴展多面體離散元模型

基于Minkowski Sum 方法，道砟顆粒的擴展多面體單元由基本多面體單元和擴展球體單元疊加而成。Minkowski Sum方法是處理復雜圖像和構造復雜顆粒模型的一種有效方法[25-26]。將2個空間體A、B疊加，即在歐幾里得幾何空間中將一個空間點集掃過另一個點集，得到新的集合，其定義式為

A⊕B={x+y|x∈A，y∈B}

( 1 )

式中:A、B分別為2個空間體；x、y分別為其對應的空間坐標。

擴展單元為球體時，它將空間上的一個點擴展為球體，將一條線段擴展為球形圓柱體，將任一平面擴展為由圓柱體包絡的具有一定厚度的板，將任一多面體擴展為具有光滑棱邊及棱角的多面體[28-29]。

本文構造的擴展多面體單元是將多面體作為基本單元，將球體作為擴展單元，將2個空間體疊加，形成具有光滑棱邊和角點的非規(guī)則形態(tài)單元體，見圖1。

在離散元數(shù)值模擬中，通過改變擴展球體的直徑得到不同光滑度的非規(guī)則顆粒單元。設基本立方體單元和四面體單元邊長均為a，球體顆粒半徑r分別為a/10、a/5時構造的擴展多面體單元見圖2。從圖2可以看出球體顆粒半徑對擴展多面體單元棱邊和棱角尖銳度的影響。

針對道砟顆粒的非規(guī)則幾何形態(tài)，本文采用擴展多面體模型構造不同形態(tài)的顆粒單元，見圖3。由于擴展多面體單元具有非規(guī)則性，用微元劃分法確定其體積，繼而確定其質量、形心及對形心的轉動慣量。在顆粒單元的作用力和運動計算中，采用全局坐標和局部坐標2種坐標。全局坐標固定在數(shù)值模擬的總體計算局域上，局部坐標固定在每個單元的形心上并隨顆粒的轉動而改變。采用歐拉角表征每個單元局部坐標的方位，通過四元數(shù)方法計算局部坐標與全局坐標之間的轉換矩陣，對顆粒間作用力和轉動分量進行坐標轉化[30]。

2 鐵路道砟單元接觸判斷

擴展多面體道砟顆粒主要由平面、棱邊和棱角3部分組成，它們在道砟材料動力過程中發(fā)生相互接觸，主要有3類接觸形式：(1) 棱角-平面、棱角-棱邊及棱角-棱角接觸；(2) 棱邊-棱邊、棱邊-平面接觸；(3)平面-平面接觸。本文選取有代表性的棱角-平面、棱邊-棱邊及平面-平面接觸進行接觸判斷。

2.1 棱角-平面接觸

道砟材料的棱角與平面間的作用狀態(tài)見圖4(a)。擴展多面體單元中道砟棱角為球形顆粒，其與平面的接觸判斷轉化為球體與平面的接觸判斷，見圖4(b)。已知球體球心O1和平面中心O2的位置矢量，由幾何關系可確定球體與平面點的重疊量。

確定球心O1到平面的距離矢量

Δ=n(O21·n)

( 2 )

式中：O21為O2到O1的位置矢量；n為平面的法向單位矢量。

O21在平面上的分量

P21=O21-Δ

( 3 )

式中：P21為O2到球心O1在平面上的投影P的位置矢量。

判斷P點是否在平面內。P點不在平面內，棱角球體與平面不發(fā)生接觸；P點在平面內時，球體到平面的最短距離

δ=|Δ|-2R

( 4 )

式中：R為擴展多面體單元選用的球體顆粒半徑；δ為發(fā)生接觸時2個顆粒單元間的法向重疊量。

δ<0時發(fā)生接觸，否則不接觸。發(fā)生接觸時，接觸力法線與Δ平行，接觸點分別為Pc1、Pc2。對棱角-平面的接觸位移進行計算，確定接觸作用力。

2.2 棱邊-棱邊接觸

2.3 平面-平面接觸

擴展多面體單元中平面-平面的接觸判斷見圖6。為簡化計算，將平面-平面接觸轉化為接觸面上的多點接觸。圖6中，將2個平面間的接觸力通過4個點接觸進行計算，其中接觸點A、C為棱角-平面接觸，B、D為棱邊-棱邊接觸。4個接觸點的接觸力確定后，通過各接觸力對顆粒質心的力矩確定接觸面間的力矩。該方法具有模型簡單、計算效率高的特點。

3 道砟單元間的接觸模型

擴展多面體單元間接觸模式不同，接觸面積差異大，導致單元間的接觸剛度和阻尼系數(shù)不同。計算擴展多面體單元的接觸力時需考慮接觸點處單元的曲率半徑和接觸面積。

3.1 道砟單元間的接觸力模型

擴展多面體道砟單元間的接觸力主要取決于2個接觸單元間的重疊量和相對速度。在不同接觸模式下，接觸單元間的法、切向力分別為

( 5 )

( 6 )

式( 5 )、式( 6 )中各參數(shù)與擴展多面體單元的接觸模式密切相關。

在接觸的切線方向，依據(jù)Mohr-coulomb摩擦定律，有

( 7 )

式中:μ為道砟顆粒材料的摩擦系數(shù)。

擴展多面體單元的棱角-棱角接觸中，接觸方式為擴展球體的接觸。采用Hertz非線性接觸模型計算其接觸力[31-32]，其法向接觸剛度為

( 8 )

式中：G、ν分別為道砟材料的剪切模量和泊松比；R*為接觸球體的有效半徑，R*=R1R2/(R1+R2)，其中R1、R2為2個擴展球體的曲率半徑。

切向剛度與法向剛度成線性關系，Kt=βKn。系數(shù)β由道砟材料的泊松比確定，β=3(1-ν)/(2-ν)。

計算接觸單元間法向和切向黏滯力時，法向、切向黏滯力系數(shù)分別為

( 9 )

(10)

3.2 時間步長的確定

對于2個球體單元，計算時間步長采用簡諧振動模型確定[33]。忽略法向黏滯力的影響，基于簡諧振動方程得最小振動周期

(11)

式中:Mmin為最小顆粒質量。

按2個球體單元的接觸模式確定時間步長

Δt=Tmin/40

(12)

4 道砟材料動力特性的擴展多面體離散元分析

4.1 道砟箱試驗的離散元模型

采用道砟箱試驗對道砟材料在往復荷載下的動力特性進行離散元分析。道砟顆粒材料的隨機生成過程中，調整道砟顆粒的長寬比，生成多種顆粒單元。道砟顆粒最大粒徑在36～63 mm均勻分布，所有擴展多面體道砟顆粒的擴展球體半徑相同。主要計算參數(shù)見表1。

表1 道砟動力特性離散元模擬中的主要計算參數(shù)

考慮有砟道床的真實結構和軌枕距離，參考文獻[15]中采用球體組合顆粒對道砟材料動力特性的離散元分析，取計算域長L=700 mm，寬B=300 mm；軌枕長Ls=300 mm，寬Bs=300 mm，質量Ms=34 kg。為模擬有砟道床在無限長度的動力特性，在長度方向采用周期邊界條件。有砟道床結構的道砟厚度一般不小于350 mm，取Hs=370 mm。

將道砟顆粒自由下落到道砟箱內，在蓋板上施加頻率3 Hz、大小0～3 kN的循環(huán)荷載，持續(xù)3 s以壓實。將軌枕放到道床上，對軌枕以低速率加載至3 kN，完成道砟顆粒的初始排列。在軌枕上施加正弦往復荷載，最大、最小值分別為40、3 kN，荷載頻率3 Hz。只考慮軌枕在外荷載和道砟反力的作用下的豎直位移，忽略其水平位移和扭轉運動。

4.2 往復荷載下道砟材料的沉降量及變形模量

對道砟材料在往復荷載作用下的動力過程進行數(shù)值分析，獲得90 s內軌枕位移和相應的往復荷載隨時間的變化，見圖8。從軌枕位移u的變化趨勢可以發(fā)現(xiàn)，軌枕的豎直沉降量隨往復荷載呈周期性變化，其累計沉降量隨加載次數(shù)的增加不斷增大，增長趨勢趨于平緩。對78～80 s的計算結果進行放大，可清晰地看出軌枕位移和荷載的對應關系，見圖8(b)。軌枕的上下震動與往復荷載在頻率和相位上均保持一致。

在90 s內的270個加載周期內，外加荷載與軌枕沉降量的對應關系見圖9。從圖9可以看出，最初幾次加載時，軌枕的位移明顯較大，因為道砟材料的初始排列相對疏松，在外荷載的作用力下達到密實排列。隨著加載次數(shù)的增加，道砟不斷密實，位移變形曲線越來越密集。在往復荷載作用下，道砟材料的累計沉降量趨于穩(wěn)定。在一個加載周期內，加載過程中位移的增長速率比卸載時的位移降低速率快。圖9體現(xiàn)了道砟材料變形的2個特性：(1)彈性變形，即每次加卸載過程中軌枕的豎向位移變化量，體現(xiàn)有砟道床的彈性性能，主要由道砟顆粒的彈性變形引起；(2)塑性變形，即在多次往復加卸載過程中軌枕的累計沉降量，主要由道砟顆粒的重新排列引起。

為分析道床在往復荷載下彈性性能的演化趨勢，通過每次往復荷載中軌枕的豎直位移變化量確定道床的有效剛度Ke

Ke=ΔF/Δu

(13)

式中：ΔF為往復荷載的幅度；Δu為一次加卸載過程中軌枕的豎向位移。

往復荷載過程中道砟材料有效剛度和軌枕累計沉降量隨加載次數(shù)的變化見圖10。從圖10(a)可以看出，隨著加載次數(shù)的增加，道床有效剛度不斷增加。在加載的前幾個周期，由于道床剛度的幾何狀態(tài)尚不穩(wěn)定，其剛度的波動性較大。特別是第一個加載周期，有效剛度較低，為7 kN/mm。隨著加載次數(shù)的增加，道床結構趨于穩(wěn)定，其有效剛度增加到11 kN/mm。圖10(b)中，umin和umax分別為每次往復荷載達到最小和最大值時所對應的軌枕沉降量。道床的累計沉降量有一個先快速沉降后趨于穩(wěn)定的過程。有效剛度和累計沉降量在一定程度上反映了有砟道床在往復荷載作用下的密實過程，分別表征了道床彈性和永久變形的變化規(guī)律。

5 結語

為描述鐵路有砟道床中道砟顆粒的非規(guī)則幾何形態(tài)，本文基于Minkowski Sum方法，將多面體模型與球體單元疊加，形成具有非規(guī)則形態(tài)的光滑顆粒單元。將塊體接觸轉化為球體單元的接觸，有效簡化了多面體單元間的接觸判斷和接觸力計算。通過改變擴展球體粒徑，構造具有不同尖銳度的擴展多面體單元，合理模擬非規(guī)則顆粒間相互作用過程。

本文采用擴展多面體顆粒單元對道砟材料在往復荷載下的動力過程進行了離散元分析，確定了其有效剛度和累計沉降量隨加載次數(shù)的演化規(guī)律。結果表明，軌枕位移在隨往復荷載上下震蕩過程中，有效剛度和累計沉降量在初始階段波動較大，隨道砟材料的密實度增加而趨于穩(wěn)定。該擴展多面體單元能夠合理地模擬道砟材料的動力行為，由微觀作用過程揭示相應的宏觀演化規(guī)律。在此工作基礎上，今后將對該模型進一步改進，系統(tǒng)地分析道砟顆粒的級配、摩擦系數(shù)、彈性模量等參數(shù)及道砟在往復荷載下的破碎和粉化過程對道床力學性質的影響。

參考文獻：

[1] 練松良. 軌道工程[M]. 上海: 同濟大學出版社, 2006.

[2] 趙國堂. 高速鐵路無砟軌道結構[M]. 北京: 中國鐵道出版社, 2006.

[3] SUN Q D, INDRARATNA B, NIMBALKAR S. Effect of Cyclic Loading Frequency on the Permanent Deformation and Degradation of Railway Ballast[J]. Geotechnique, 2014，64(9): 746-751.

[4] SUIKER A S J, SELIG E T, FRENKEL R. Static and Cyclic Triaxial Testing of Ballast and Subballast[J].Journal of Geotechnical and Geoenvironmental Engineering, 2005, 131(6): 771-782.

[5] ASCE B I F, THAKUR P K, VINOD J S,et al. Semiempirical Cyclic Densification Model for Ballast Incorporating Particle Breakage[J].International Journal of Geomechanics, 2012, 12(3): 260-271.

[6] LACKENBY J, INDRARATNA B, MCDOWELL G, et al. Effect of Confining Pressure on Ballast Degradation and Deformation under Triaxial Loading[J].Geotechnique,2007,57(6):527-536.

[7] HOSSAIN Z, INDRARATNA B, DARVE F,et al. DEM Analysis of Angular Ballast Breakage under Cyclic Loading[J].Geomechanics and Geoengineering: An International Journal, 2007, 2(3): 175-181.

[8] ZHAI W M, WANG K Y, LIN J H. Modeling and Experiment of Railway Ballast Vibrations[J]. Journal of Sound and Vibration，2004，270(4-5): 673-683.

[9] 董亮, 蔡德鉤, 葉陽升, 等. 列車循環(huán)荷載作用下高速鐵路路基累積變形預測方法[J].土木工程學報，2010, 43(6): 100-108.

DONG Liang, CAI De-gou, YE Yang-sheng, et al. A Method for Predicting the Deformation of High-speed Railway Subgrades under Cyclic Train Loads[J]. China Civil engineering journal, 2010, 43(6): 100-108.

[10] 井國慶, 封坤, 高亮, 等. 循環(huán)荷載作用下道砟破碎老化的離散元仿真[J]. 西南交通大學學報, 2012, 47(2): 187-191.

JING Gou-qing, FENG Kun, GAO Liang, et al. DEM Simulation of Ballast Degradation and Breakage under Cyclic Loading[J]. Journal of Southwest Jiaotong University, 2012, 47(2): 187-191.

[11] 肖宏, 高亮, 侯博文. 鐵路道床振動特性的三維離散元分析[J]. 鐵道工程學報, 2009, (9): 14-17.

XIAO Hong, GAO Liang, HOU Bo-wen. Analysis of Ballast Dynamic Behavior with Three-dimensional Discrete Element Method[J]. 2009, (9): 14-17.

[12] MCDOWELL G R, LIM W L, COLLOP A C, et al. Comparison of Ballast Index Tests for Railway Trackbeds[J]. Proceeding of the Institution of Civil Engineers Geotechnical Engineering, 2004, 157(3): 151-161.

[13] 姜衛(wèi)利, 范俊杰. 散粒體道床離散單元法分析[J]. 鐵道學報, 2001, 23(4): 98-101.

JIANG Wei-li, FAN Jun-jie. Using of Distinct Element Method to Analyze Granular Ballast Bed[J]. Journal of the China Railway Society, 2001,23(4):98-101.

[14] 嚴穎, 狄少丞, 蘇勇, 等. 風沙影響下鐵路道砟變形模量的離散元法數(shù)值分析[J]. 計算力學學報, 2012, 29(3): 339-445.

YAN Ying, DI Shao-cheng, SU Yong, et al.Discrete Element Analysis of Elastic Modulus of Railway Ballasts in Wind with Different Sand Contents[J]. Chinese Journal of Computational Mechanics, 2012, 29(3): 339-445.

[15] LIM W L, MCDOWELL G R. Discrete Element Modeling of Railway Ballast[J]. Granular Matter, 2005, 7(1): 19-29.

[16] LU M, MCDOWELL G R.The Importance of Modeling Ballast Particle Shape in the Discrete Element Method[J]. Granular Matter,2007, 9(1-2): 69-80.

[17] LOBO-GUERRERO S, VALLEJO L. Discrete Element Method Analysis of Rail Track Ballast Degradation During Cyclic Loading[J].Granular Matter,2006, 8(3-4): 195-204.

[18] INDRARATNA B, THAKUR P, VINOD J. Experimental and Numerical Study of Railway Ballast Behavior under Cyclic Loading[J]. International Journal of Geomechanics, 2010, 10(4): 136-144.

[19] WILLIAMS J R, PENTL A P. Superquadrics and Modal Dynamics for Discrete Elements in Interactive Design [J]. Engineering Computations, 1992, 9(2): 115-127.

[20] CLEARY P W. Large Scale Industrial DEM Modeling[J]. Engineering Computations, 2004, 21(2/3/4): 169-204.

[21] HOHNER D, WIRTZ S, SCHERER V. A Numerical Study on the Influence of Particle Shape on Hopper Discharge within the Polyhedral and Muti-sphere Discrete Element Method[J]. Powder Technology, 2012, (226): 16-28.

[22] PENA A A, LIZCANO A, ALONSO-MARROQUIN F, et al. Biaxial Test Simulations Using a packing of Polygonal Particles[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2008, 32(2):143-160.

[23] HOSSEININIA E S. Discrete Element Modeling of Inherently Anisotropic Granular Assemblies with Polygonal Particles[J]. Particuology, 2012,(10):542-552.

[24] GALINDO-TORRES S A, MUNOZ J D. Minkowski - Voronoi Diagrams as a Method to Generate Random Packings of Spheropolygons for the Simulation of Soils[J]. Physical Review E, 2010, 82(5): 1-12.

[25] ALONSO-MARROQUIN F, WANG Y C. An Efficient Algorithm for Granular Dynamics Simulations with Complex-shaped Objects [J]. Granular Matter, 2009, 11(5), 317-329.

[26] GALINDO-TORRES S A, ALONSO-MARROQUIN F, WANG Y C, et al. Molecular Dynamics Simulation of Complex Particles in Three Dimensions and the Study of Friction Due to Nonconvexity[J]. Physical Review E, 2009, 79(6): 060301.

[27] HOPKINS M A. Discrete Element Modeling with Dilated Particles[J]. Engineering Computations, 2003, 21(2/3/4): 422-430.

[28] POURNIN L, WEBER M, TSUKAHARA M,et al. Three-dimensional Distinct Element Simulation of Spherocylinder Crystallization[J]. Granular Matter, 2005, 7(2): 119-126.

[29] GALINDO-TORRES S A, PEDROSO D M, WILLIAMS D J,et al. Breaking Processes in Three-dimensional Bonded Granular Materials with General Shapes [J]. Computer Physics Communications, 2012, 183(2):266-277.

[30] YAN Y, JI S. Discrete Element Modeling of Direct Shear Tests for a Granular Material[J]. International Journal for Numerical and Analytical Methods in Geomechanics,2010,34(9):978-990.

[31] JI S, SHEN H H. Effect of Contact Force Models on Granular Flow Dynamics[J]. Journal of Engineering Mechanics, 2006, 132(11): 1252-1259.

[32] JOHNSON K L. Contact Mechanics[M].Cambridge: Cambridge University Press，1987.

[33] 孫其誠, 王光謙. 顆粒物質力學導論[M]. 北京: 科學出版社. 2009.