孙珊珊， 严 颖， 赵春发， 季顺迎

(1. 大连理工大学 工业装备结构分析国家重点实验室，辽宁 大连 116023； 2. 大连交通大学 土木与安全工程学院，辽宁 大连 116028； 3. 西南交通大学 牵引动力国家重点实验室，四川 成都 610031)

有砟轨道是当前铁路轨道的主要结构形式，在铁路运输向高速和重载方向发展过程中，其力学性能不断改进[1-2]。铁路道床在列车往复荷载作用下因道砟颗粒的重新排列和破碎粉化发生沉降，导致整体弹性的降低，影响其使用性能[3]。

相关学者对材料的力学性能开展了系统研究，主要集中在累积沉降量和有效弹性等方面。道砟材料的三轴试验结果表明，在密实过程中道砟材料减震和排水性能均有所降低[4-5]。在高频重载条件下，道砟材料的沉降和道砟颗粒破坏现象尤为突出[6]。在往复荷载作用下，道砟颗粒的破碎是引起道床变形的重要原因[7]。将铁路道床假设为连续体，用有限元方法可分析其动力特性和变形规律[8-9]，但仅限小变形行为，且不能描述道砟颗粒的细观受力特征。颗粒的几何形状及级配对其宏观力学行为的影响显著[10-11]。

为真实模拟道砟材料的几何形态和力学行为，基于非连续材料的离散单元模型迅速发展[12-14]。在颗粒材料的离散元数值模拟中，采用最普遍的是球形颗粒单元。对于非规则形态的道砟颗粒，球形颗粒很难有效地模拟道砟材料的力学行为[15]，可采用球形颗粒的不同黏接或镶嵌组合方式构造相应的道砟颗粒[14,16-17]，模拟道砟颗粒在外荷载作用下的破碎现象，对有砟道床的沉降特性进行离散元分析[7,17-18]。在球形颗粒的组合中，通过增加球形颗粒单元的数量可更精确地模拟非规则颗粒，甚至可采用上千个球形颗粒单元构造一个道砟块石[16]。但颗粒数目的增加急剧降低离散元的计算效率，不利于开展大规模计算。

为模拟非球形颗粒的几何形态，发展了不同的离散单元模型。采用二次或超二次曲面方程可定义不同形态的颗粒单元[19-20]。对类似道砟材料的多棱边颗粒，多面体单元能够合理描述颗粒的几何形态以及颗粒间的相互作用，准确模拟其宏观力学特征[21]。多面体单元的作用力计算需要繁琐的接触判断和邻居搜索，严重影响其计算效率和工程应用性。采用Minkowski Sum方法构造扩展多面体颗粒模型，准确快速地模拟非规则颗粒间的相互作用。在离散单元法的数值研究中，这种非规则颗粒模型已成功用于二维的扩展多边形[22-23]和三维的扩展多面体计算[24-26]。扩展多面体单元的主要特点是在多面体单元的棱边及外表面设置若干个球形颗粒，通过改变球形颗粒的直径灵活控制多面体棱边及棱角的尖锐度[27]。该模型可有效简化多面体单元的接触判断过程，在满足计算精度的条件下有效地提高离散元计算效率。

本文采用扩展多面体模型对有砟道床在往复荷载作用下的动力过程进行离散元分析，确定其累积沉降量和有效模量随荷载往复次数的演化规律。

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个空间体叠加，形成具有光滑棱边和角点的非规则形态单元体，见图1。

在离散元数值模拟中，通过改变扩展球体的直径得到不同光滑度的非规则颗粒单元。设基本立方体单元和四面体单元边长均为a，球体颗粒半径r分别为a/10、a/5时构造的扩展多面体单元见图2。从图2可以看出球体颗粒半径对扩展多面体单元棱边和棱角尖锐度的影响。

针对道砟颗粒的非规则几何形态，本文采用扩展多面体模型构造不同形态的颗粒单元，见图3。由于扩展多面体单元具有非规则性，用微元划分法确定其体积，继而确定其质量、形心及对形心的转动惯量。在颗粒单元的作用力和运动计算中，采用全局坐标和局部坐标2种坐标。全局坐标固定在数值模拟的总体计算局域上，局部坐标固定在每个单元的形心上并随颗粒的转动而改变。采用欧拉角表征每个单元局部坐标的方位，通过四元数方法计算局部坐标与全局坐标之间的转换矩阵，对颗粒间作用力和转动分量进行坐标转化[30]。

2 铁路道砟单元接触判断

扩展多面体道砟颗粒主要由平面、棱边和棱角3部分组成，它们在道砟材料动力过程中发生相互接触，主要有3类接触形式：(1) 棱角-平面、棱角-棱边及棱角-棱角接触；(2) 棱边-棱边、棱边-平面接触；(3)平面-平面接触。本文选取有代表性的棱角-平面、棱边-棱边及平面-平面接触进行接触判断。

2.1 棱角-平面接触

道砟材料的棱角与平面间的作用状态见图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点不在平面内，棱角球体与平面不发生接触；P点在平面内时，球体到平面的最短距离

δ=|Δ|-2R

( 4 )

式中：R为扩展多面体单元选用的球体颗粒半径；δ为发生接触时2个颗粒单元间的法向重叠量。

δ<0时发生接触，否则不接触。发生接触时，接触力法线与Δ平行，接触点分别为Pc1、Pc2。对棱角-平面的接触位移进行计算，确定接触作用力。

2.2 棱边-棱边接触

2.3 平面-平面接触

扩展多面体单元中平面-平面的接触判断见图6。为简化计算，将平面-平面接触转化为接触面上的多点接触。图6中，将2个平面间的接触力通过4个点接触进行计算，其中接触点A、C为棱角-平面接触，B、D为棱边-棱边接触。4个接触点的接触力确定后，通过各接触力对颗粒质心的力矩确定接触面间的力矩。该方法具有模型简单、计算效率高的特点。

3 道砟单元间的接触模型

扩展多面体单元间接触模式不同，接触面积差异大，导致单元间的接触刚度和阻尼系数不同。计算扩展多面体单元的接触力时需考虑接触点处单元的曲率半径和接触面积。

3.1 道砟单元间的接触力模型

扩展多面体道砟单元间的接触力主要取决于2个接触单元间的重叠量和相对速度。在不同接触模式下，接触单元间的法、切向力分别为

( 5 )

( 6 )

式( 5 )、式( 6 )中各参数与扩展多面体单元的接触模式密切相关。

在接触的切线方向，依据Mohr-coulomb摩擦定律，有

( 7 )

式中:μ为道砟颗粒材料的摩擦系数。

扩展多面体单元的棱角-棱角接触中，接触方式为扩展球体的接触。采用Hertz非线性接触模型计算其接触力[31-32]，其法向接触刚度为

( 8 )

式中：G、ν分别为道砟材料的剪切模量和泊松比；R*为接触球体的有效半径，R*=R1R2/(R1+R2)，其中R1、R2为2个扩展球体的曲率半径。

切向刚度与法向刚度成线性关系，Kt=βKn。系数β由道砟材料的泊松比确定，β=3(1-ν)/(2-ν)。

计算接触单元间法向和切向黏滞力时，法向、切向黏滞力系数分别为

( 9 )

(10)

3.2 时间步长的确定

对于2个球体单元，计算时间步长采用简谐振动模型确定[33]。忽略法向黏滞力的影响，基于简谐振动方程得最小振动周期

(11)

式中:Mmin为最小颗粒质量。

按2个球体单元的接触模式确定时间步长

Δt=Tmin/40

(12)

4 道砟材料动力特性的扩展多面体离散元分析

4.1 道砟箱试验的离散元模型

采用道砟箱试验对道砟材料在往复荷载下的动力特性进行离散元分析。道砟颗粒材料的随机生成过程中，调整道砟颗粒的长宽比，生成多种颗粒单元。道砟颗粒最大粒径在36～63 mm均匀分布，所有扩展多面体道砟颗粒的扩展球体半径相同。主要计算参数见表1。

表1 道砟动力特性离散元模拟中的主要计算参数

考虑有砟道床的真实结构和轨枕距离，参考文献[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的循环荷载，持续3 s以压实。将轨枕放到道床上，对轨枕以低速率加载至3 kN，完成道砟颗粒的初始排列。在轨枕上施加正弦往复荷载，最大、最小值分别为40、3 kN，荷载频率3 Hz。只考虑轨枕在外荷载和道砟反力的作用下的竖直位移，忽略其水平位移和扭转运动。

4.2 往复荷载下道砟材料的沉降量及变形模量

对道砟材料在往复荷载作用下的动力过程进行数值分析，获得90 s内轨枕位移和相应的往复荷载随时间的变化，见图8。从轨枕位移u的变化趋势可以发现，轨枕的竖直沉降量随往复荷载呈周期性变化，其累计沉降量随加载次数的增加不断增大，增长趋势趋于平缓。对78～80 s的计算结果进行放大，可清晰地看出轨枕位移和荷载的对应关系，见图8(b)。轨枕的上下震动与往复荷载在频率和相位上均保持一致。

在90 s内的270个加载周期内，外加荷载与轨枕沉降量的对应关系见图9。从图9可以看出，最初几次加载时，轨枕的位移明显较大，因为道砟材料的初始排列相对疏松，在外荷载的作用力下达到密实排列。随着加载次数的增加，道砟不断密实，位移变形曲线越来越密集。在往复荷载作用下，道砟材料的累计沉降量趋于稳定。在一个加载周期内，加载过程中位移的增长速率比卸载时的位移降低速率快。图9体现了道砟材料变形的2个特性：(1)弹性变形，即每次加卸载过程中轨枕的竖向位移变化量，体现有砟道床的弹性性能，主要由道砟颗粒的弹性变形引起；(2)塑性变形，即在多次往复加卸载过程中轨枕的累计沉降量，主要由道砟颗粒的重新排列引起。

为分析道床在往复荷载下弹性性能的演化趋势，通过每次往复荷载中轨枕的竖直位移变化量确定道床的有效刚度Ke

Ke=ΔF/Δu

(13)

式中：ΔF为往复荷载的幅度；Δu为一次加卸载过程中轨枕的竖向位移。

往复荷载过程中道砟材料有效刚度和轨枕累计沉降量随加载次数的变化见图10。从图10(a)可以看出，随着加载次数的增加，道床有效刚度不断增加。在加载的前几个周期，由于道床刚度的几何状态尚不稳定，其刚度的波动性较大。特别是第一个加载周期，有效刚度较低，为7 kN/mm。随着加载次数的增加，道床结构趋于稳定，其有效刚度增加到11 kN/mm。图10(b)中，umin和umax分别为每次往复荷载达到最小和最大值时所对应的轨枕沉降量。道床的累计沉降量有一个先快速沉降后趋于稳定的过程。有效刚度和累计沉降量在一定程度上反映了有砟道床在往复荷载作用下的密实过程，分别表征了道床弹性和永久变形的变化规律。

5 结语

为描述铁路有砟道床中道砟颗粒的非规则几何形态，本文基于Minkowski Sum方法，将多面体模型与球体单元叠加，形成具有非规则形态的光滑颗粒单元。将块体接触转化为球体单元的接触，有效简化了多面体单元间的接触判断和接触力计算。通过改变扩展球体粒径，构造具有不同尖锐度的扩展多面体单元，合理模拟非规则颗粒间相互作用过程。

本文采用扩展多面体颗粒单元对道砟材料在往复荷载下的动力过程进行了离散元分析，确定了其有效刚度和累计沉降量随加载次数的演化规律。结果表明，轨枕位移在随往复荷载上下震荡过程中，有效刚度和累计沉降量在初始阶段波动较大，随道砟材料的密实度增加而趋于稳定。该扩展多面体单元能够合理地模拟道砟材料的动力行为，由微观作用过程揭示相应的宏观演化规律。在此工作基础上，今后将对该模型进一步改进，系统地分析道砟颗粒的级配、摩擦系数、弹性模量等参数及道砟在往复荷载下的破碎和粉化过程对道床力学性质的影响。

参考文献：

[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] 董亮, 蔡德钩, 叶阳升, 等. 列车循环荷载作用下高速铁路路基累积变形预测方法[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] 井国庆, 封坤, 高亮, 等. 循环荷载作用下道砟破碎老化的离散元仿真[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] 姜卫利, 范俊杰. 散粒体道床离散单元法分析[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] 严颖, 狄少丞, 苏勇, 等. 风沙影响下铁路道砟变形模量的离散元法数值分析[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.