针对有砟道床动力特性分析的嵌入式离散元-有限元耦合方法

邵 帅， 吕泉江， 纪丹阳， 严 颖

(1.内蒙古大学 交通学院，呼和浩特 010070； 2.中国直升机设计研究所，景德镇 333001；3.大连交通大学 土木与安全工程学院，大连 116028)

1 引 言

传统的有砟铁路道床结构组成成分多样，且力学性质复杂，包括轨枕、道砟、底部道砟、路堤和地基等组成部分[1]。其中，道砟层由大量碎石构成，具有很多孔隙，是典型的散体材料[2]。枕木、路堤和地基等部分具有相对密实的材料性质，与道砟相比可视为连续介质。有砟道床中散体介质与连续体介质并存的结构特点给数值模拟造成了很多困难[3,4]。发展一种能够准确计算并预测整体有砟道床力学行为的计算模型对于有砟道床的养护、维修及安全运行均具有重要的工程意义。

目前，研究人员为预测有砟铁路道床的力学行为已经开展了大量数值仿真工作，采用的主要数值方法包括离散元方法(DEM)和有限元方法(FEM)。离散元方法最初由Cundall等[5]提出并应用于农业、化学、岩土和铁道工程等领域[6-8]。Gong等[9]采用离散元方法模拟了混合有轮胎衍生集料道砟材料的剪切强度，发现轮胎衍生集料能够减小剪切应力峰值，降低道砟材料的内聚力和摩擦角，减少直剪试验中道砟颗粒的破碎。Jing等[10]通过离散元数值模拟研究发现，道砟层对维持轨枕的稳定性起到重要作用，并通过试验进行了验证。Liu等[11]采用离散元方法模拟真实尺度下道砟材料大三轴试验，基于智能型道砟碎石开发离散元计算程序的计算精度有了很大提升。Zhang等[12]采用镶嵌单元构造道砟碎石，并通过离散元数值模拟研究了高速列车载荷频率对有砟道床累积沉降量的影响。有限元方法经过几十年的发展已经逐渐成熟，在铁道工程领域得到广泛认可[13，14]。为了能够获得不同工况下轨枕的不均匀沉降，Shih等[15]综合利用Abaqus，Python和Fortran建立了一套名为BaTrack的有限元程序，对于有砟道床力学行为分析具有较好的模拟精度。

离散元方法和有限元方法在有砟道床动力特性的模拟中各具特点，有限元方法具有足够的精度且拥有较高的计算效率，离散元方法能够在细观尺度上提供很多力学信息，为充分发挥离散元方法和有限元方法各自的优势，研究人员将离散元方法和有限元方法相互耦合并应用于铁道工程领域。道砟层采用离散元计算，路堤与地基通过有限元方法模拟，在离散元域和有限元域的接触面上依据虚功原理提出耦合接触算法，建立离散元-有限元耦合模型。此模型能够更全面地从宏细观多尺度分析整体有砟道床的力学行为[16-19]。然而，在离散元和有限元的耦合接触面上，由于球形颗粒表面相对光滑，在与有限元表面的接触过程中受外部载荷作用容易产生侧向滑动，使得数值模型不稳定，导致数值结果准确度降低。对于实际工况，底部道砟往往会嵌入到路堤中，其水平方向的运动受到很大限制，从而确保其上部结构的稳定性。

综上所述，为增加有砟道床离散元-有限元耦合模型中耦合面上的自锁能力，并考虑实际工况中底部道砟嵌入路堤中的客观事实，本文发展了一种针对有砟道床的嵌入式离散元-有限元耦合模型，有效地提高了耦合模型的数值稳定性，为计算模拟铁路工程问题提供了一种新的手段。

2 针对有砟道床结构离散元-有限元耦合模型

针对有砟道床结构的计算模拟区域及道床横截面如图1所示。根据一根枕木的作用范围，在有砟道床中截取长为4 m，宽为0.23 m的区域作为计算模拟区，如图1(a)所示。在有砟道床的横截面中考虑五个组成部分，分别为枕木、道砟、底部道砟、路堤和地基，如图1(b)所示。其中，枕木、路堤和地基相对密实，在数值模拟中可当做连续介质，其力学行为通过有限元方法进行计算。道砟和底部道砟是典型的散体材料，在数值模拟中作为离散介质，其力学行为通过离散元方法进行分析。为了实现对整体有砟道床结构的动力特性分析，需在离散元和有限元的耦合接触面上建立耦合接触算法，实现计算参数在离散元域和有限元域间的传递。在图1(b)中，利用有砟道床横截面的对称性质，在建模中只需考虑横截面的右半部分即可。

图1 有砟道床的模拟区域及道床横截面

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

道砟是一种典型的散体材料，其内部具有大量孔隙，在上部列车载荷长期作用下容易产生不均匀沉降，整体力学性质十分复杂。构成道砟层的道砟碎石具有几何形状极其不规则的结构特点，道砟颗粒间的咬合作用和自锁能力较强。为从细观尺度上分析有砟道床的动力特性以及道砟颗粒几何形状不规则的结构特点，采用镶嵌单元对道砟颗粒进行构造，如图2所示。考虑计算效率的影响，构造单个道砟颗粒单元使用的球形颗粒数目为3～11个。由于道砟形状随机性很强，采用较少数量的球形颗粒进行模拟，虽然对当前道砟形状模拟有一定误差，但能够保证整体有砟道床的模拟精度。

图2 道砟碎石以及生成的镶嵌单元

采用镶嵌单元铺设后的道砟层离散元模型如图3所示。道砟层顶部承受由枕木传递下来的列车荷载，左侧承受x方向的位移约束，右侧为自由表面，前后表面承受y方向的位移约束。构造道砟层共使用1727个镶嵌单元，10098个球形颗粒，孔隙率为0.45。颗粒间的接触模型采用赫兹模型，法向接触力计算公式为[20]

(1)

式中A为耗散系数。不考虑切向粘滞力，根据摩尔库伦摩擦定律，切向接触力计算公式为[21]

(2)

(3)

(4)

(5，6)

图3 采用镶嵌单元构造的道砟层离散单元模型

表1 道砟材料计算参数

2.2 有砟道床的有限元模型

有砟道床的有限元模型如图4所示。其中，上部的梯形和下部的矩形区域分别代表路堤和地基。采用20节点等参元对模型进行网格剖分，由于模型中枕木下方是主要的承载区域，故对网格进行了加密处理，模型共有1120个单元，5805个节点。模型左右两侧受x方向的位移约束，底面受z方向的位移约束，前后面受y方向的位移约束。上表面为自由表面，用于承受由道砟层传递下来的列车载荷。路堤和地基的材料参数列入表2。

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

表2 路堤和地基的材料参数

采用Newmark隐式时间积分计算路堤和地基的动力响应，Newmark方法的计算格式为

(7)

(8)

(9)

式中Fk + 1为时刻tk + 1时刻的等效节点载荷，δ和α为Newmark方法中的计算参数，当δ≥0.5且α≥0.25(0.5+δ)2时，Newmark方法是无条件稳定的，本文取δ=0.5，α=0.25。考虑到路堤的非线性性质，采用Newton-Raphson方法求解路堤的力学行为[22]。

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

在以往针对有砟道床结构的离散元-有限元耦合计算方法中，使用球形颗粒构造的镶嵌单元表面比较圆润，与路堤上表面很难产生较强的自锁和咬合作用，在列车载荷的作用下有砟层容易产生整体的侧向移动。此外，在实际的有砟道床中，由于路堤的塑性性质，底部道砟会嵌入路堤上表面。为使数值模型中道砟层更加稳定并考虑实际工况，本文发展了针对有砟道床结构的嵌入式离散元-有限元耦合方法。

在嵌入式离散元-有限元耦合模型中，设置一层嵌入路堤顶部的耦合过渡层，如图5所示。过渡层共有483个球形颗粒，呈规则的矩形方式排列，球形颗粒的球心位于路堤的上表面。过渡层将承受上部道砟层传递的列车载荷，并通过耦合接触算法将耦合接触力作为等效节点载荷传递至有限元

图5 离散元和有限元耦合过渡层

网格中，进而计算下部路堤和地基的变形。当下部结构的节点位移计算完成后，通过形函数插值可进一步更新过渡层中各球形颗粒的位置，从而实现离散元方法和有限元方法的耦合。在实际的应用中，过渡层中球形颗粒的粒径应参照铁路碎石道床底碴的级配要求进行选取，本文以30 mm粒径为例进行计算。

图6 嵌入式离散元-有限元耦合方法的接触模型

根据虚功原理，离散元等效接触力产生的外力虚功δW为

δW=δUTF

(10)

式中F为等效接触力矢量，U为力作用点处的位移矢量，可通过对20节点等参元节点位移进行插值求得

(i=1～20)(11)

(i=1～20)(12)

有限元节点载荷所做虚功为

(i=1～20)(13)

根据虚功原理，有限元等效节点载荷Fnodal,i的计算公式为

(i=1～20)(14)

20节点等参元的形函数为

(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)是整体坐标系下力作用点的位置坐

图7 接触力等效节点载荷

(18)

3 列车载荷作用下有砟道床结构动力特性的计算结果及分析

通过本文建立的嵌入式离散元-有限元耦合模型,对有砟道床结构进行动力特性分析，在枕木顶部施加的列车载荷如图9所示。此列车载荷由翟婉明等[23]提出的车辆-轨道耦合动力学方法获得。首先在轨道上设置一个观测点，当列车车轮通过观测点时，将产生一个力的峰值。列车的每节车厢共有2组车轮，4根轮轴，一节车厢完全通过测试点后会产生两组共4个载荷峰值，在此过程中会伴随产生小幅振动。列车载荷的最大值为20 kN，一个载荷循环周期为1.2 s。本文算例共施加了15次载荷循环。

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

图9 施加于枕木上的列车载荷

有砟道床顶部沉降量的时程曲线如图10所示。沉降量曲线的每个循环与载荷曲线形式类似，具有两组共4个峰值，伴随有小幅振动。在加载初期，由于有砟道床内部道砟相对松散，故道床顶部沉降明显。随着加载的不断进行，有砟道床内部的道砟颗粒会产生相互错动以及空间位置的重新调整，在此过程中有砟道床的密实程度不断增加，致使道床顶部沉降逐渐平缓，此时有砟道床已没有明显的累积沉降量。

有砟道床下部路堤的沉降量时程曲线如图11所示。为研究路堤顶部沿横断面方向的沉降规律，选取3个观测点进行观测，3个观测点距模型左侧边界的距离分别为0.16 m,0.78 m和1.91 m。观察发现，路堤沉降量曲线的单个循环形式与道床顶部沉降曲线形式有两点主要区别，一为在路堤的沉降量曲线中，由列车载荷引起的小幅振动已经消失，其原因是道砟层具有较好的缓冲性能，列车载荷中的小幅振动已由道砟层吸收；二为路堤的沉降量曲线中，循环载荷第二个载荷峰值产生的沉降量明显高于第一个载荷峰值产生的沉降量，主要是由道砟层散体介质的力学性质以及快速加载的载荷特性所导致，第一个载荷峰值在道砟层中产生一部分变形能，随着快速卸载，道砟层中变形能还未全部释放，第二次加载已经开始，进一步使道砟层中变形能不断增加，导致第二次载荷峰值引起的沉降量高于第一次。此外，在枕木下方(0.16和0.78观测点)，道床内部沉降略高于外侧沉降，但区别并不显著。然而在枕木施压区外侧，即1.91处观测点的位移表现出向上运动的趋势，这主要是由于路堤内侧受挤压导致外侧的膨胀。

图11 路堤顶部沉降量

有砟道床中道砟层力链分布、路堤与地基内位移和应力云图分别如图12和图13所示。发现在道砟层内部力链主要分布于枕木下方，形成承载骨架，能够将顶部的列车载荷传递至下部路堤。由于耦合模型中存在嵌入路堤中的一层道砟颗粒，耦合接触面上具有较大的摩擦系数，咬合作用明显，因此道砟层中的力链基本限制于枕木下方的矩形区域，力链没有产生向四周发散的现象。砟肩处受力较小，内部没有明显的力链分布，由于枕木的挤压作用，砟肩处有轻微隆起。路堤和地基的位移云图与道砟层中的力链分布相吻合，主要位于枕木下方。最大位移产生位置与力链分布相对应，具有一定随机性，道砟层中的主要承载骨架将传递大部分的列车载荷，进而使其下方的路堤产生较大位移。观察图13的应力云图可以发现,路堤和地基中的应力分布与道砟层中的力链分布相互吻合，应力在道床内部较大，向边缘逐渐递减。

图12 有砟道床内部力链分布及位移云图

图13 有砟道床内部力链分布及应力云图

道砟-路堤耦合交界面竖直方向与水平方向的载荷分布如图14所示。发现竖直方向载荷与水平方向载荷均与力链分布相吻合，位于枕木下方，呈现出很强的非均匀性质，主要由道砟层中道砟碎石的随机排列导致。通过耦合交界面水平载荷分布发现,水平载荷在枕木右端下方区域达到最大值，说明嵌入路堤的球形颗粒对其上端的道砟颗粒提供了明显的水平方向约束作用，有效减少道砟颗粒在上部载荷作用下向一侧滑移，对有砟道床离散元-有限元耦合模型的数值稳定性起到有益的作用。

图14 道砟-路堤耦合面上的载荷分布云图

4 结 论

本文针对有砟道床的结构特点，考虑实际铁路中底部道砟嵌入路堤的客观事实，发展了一种嵌入式离散元-有限元耦合模型,得到以下结论。

(1) 有砟道床中非规则的道砟颗粒可通过镶嵌单元模拟，然而镶嵌单元由若干球形颗粒构成，表面相对光滑，导致其与边界作用时自锁能力较差，容易产生滑动。

(2) 根据虚功原理建立了嵌入式离散元-有限元耦合模型，嵌入路堤中的球形颗粒受到路堤的约束作用很难产生侧向滑动，颗粒间具有较好的自锁能力，从而通过颗粒间咬合作用能减少上部道砟层在列车载荷作用下的侧向移动。

(3) 嵌入式离散元-有限元耦合模型能从力链分布、沉降量、位移云图及耦合交界面处力学响应等方面对有砟道床在宏细观多尺度上进行计算分析。

工程中大量问题可通过离散元-有限元耦合方法进行数值模拟，对于类似有砟道床的散体介质与连续介质共同存在的特殊结构，采用嵌入式耦合模型处理耦合交界面更加符合工程实际，具有更好的数值稳定性。