铁路道床有砟-无砟过渡段动力特性的DEM-FEM耦合分析

邵 帅， 吕泉江， 严 颖， 季顺迎＊

（1.内蒙古大学 交通学院，呼和浩特010070；2.大连理工大学 工业装备结构分析国家重点实验室，大连116024；3.大连交通大学 土木与安全工程学院，大连116028）

1 引 言

以混凝土或沥青砂浆组成的无砟轨道结构形式在我国高速铁路线路上广泛应用，但对于不良地质和大跨度桥梁等地段，仍为有砟轨道设计，不可避免地会遇到有砟和无砟道床的过渡段。在有砟-无砟轨道过渡段，由于材料和刚度等方面的差异，产生和存在的沉降差异较大。随着沉降差的增大，将导致列车产生垂向加速度，对轨道产生很大的竖向冲击力，加速轨道结构的劣化。因此，在过渡段采取相应的措施减少沉降差十分必要。在施工过程中可采用将道砟颗粒嵌入无砟道床的方式增加道砟颗粒之间的咬合力，降低沉降量，保证过渡段道床的稳定。

有砟道床作为一种散体材料，具有很强的非均匀、非连续性和各项异性等非线性性质［1］。道砟颗粒的非规则形状对整体道床的稳定性起到重要作用［2］。离散元方法作为一种有效的数值方法，在模拟散体材料时具有很大的优势［3－4］。Cundall等［5］提出离散单元法之后，该方法广泛地用于分析道砟材料的力学行为［6－9］。目前，对非规则道砟的形态描述主要有粘结单元和镶嵌单元。粘结单元是由球形颗粒在不同位置进行组合生成的可破碎组合单元，当颗粒间的力大于粘结强度时，粘结单元可发生破坏［10］。采用粘结单元可模拟单个道砟的破碎过程［11］，也可以模拟列车荷载下道砟道床的颗粒破碎［8，12，13］。镶嵌单元是在颗粒排列中存在一定重叠量的组合单元，采用镶嵌单元可模拟道砟颗粒［14，15］，并且在计算过程中将整个镶嵌单元视为刚体，保持单元不会破碎［16－18］。此外，采用离散单元法可分析土工格栅加固下道砟直剪和沉降的力学规律［19－21］。

有限单元法可有效地分析连续体材料的动力学特性，在铁路道床动力分析中也有较为广泛的应用［22－24］。将铁路道床视为连续体材料，并赋予材料属性，利用基于连续介质力学的有限元法对铁路道床进行分析研究，得出的参数可作为铁路设计的有益参考［25－28］。列车荷载引起的道床振动是铁路设计的重要参考因素，采用有限元法可模拟分析列车荷载下的道床振动，并通过试验对比验证计算模型的可靠性［26］。此外，采用子结构方法可分析高速列车载荷下铁轨的振动，以及列车提速时铁轨的动力响应［29］。本文将通过有限元法构造铁路过渡段无砟道床模型，分析道床在循环荷载作用下的动力响应。

铁路道床过渡段是一个复杂的结构体系，既包括散体道砟材料，又包括连续体道床材料。采用离散元-有限元耦合算法可在不同的区域分别采用离散单元法和有限单元法，充分发挥两种算法的优势。离散元和有限元之间的参数传递是耦合过程中的关键问题，可以通过引入过渡层实现两种算法的平滑过渡［30－31］，也可以采用罚函数法实现不同区域之间力学信息的传递［32］。此外，自适应耦合算法可在局部发生破坏时，将部分有限元转化成离散单元，其他区域仍采用有限元法计算，提高了计算效率［33－34］。李锡夔等［35］基 于 连 接尺度方 法，建 立了细尺度上离散颗粒模型和粗尺度上Cosserat连续体模型的离散元和有限元连接尺度模型。虽然离散元-有限元耦合算法在不断发展和完善，但在铁路道床有砟-无砟过渡区域中的应用尚未开展。

本文将建立一种分析循环载荷作用下，铁路道床有砟-无砟过渡段的三维DEM-FEM耦合模型。在道砟颗粒嵌入无砟道床的情况下，探讨有砟-无砟道床过渡段的动力特性，分析有砟和无砟道床的沉降量差异。

2 铁路道床有砟-无砟过渡段的DEM-FEM耦合模型

2.1 有砟道床的离散单元法

在采用离散元方法模拟有砟道床时，需要生成道砟颗粒集合体模型。道砟颗粒外形复杂，每个颗粒的外形均不相同，其几何形状对道床的运动有显著影响。为建立真实的有砟道床模型，本文采用镶嵌球体单元模拟道砟颗粒。镶嵌颗粒的内部球体单元之间无法计算接触力，其在运动过程中也不会破碎，可视为刚体。

图1 道砟颗粒及其镶嵌球体离散元模型Fig.1 Ballast particle and its DEM with overlapped spheres

为准确生成非规则道砟颗粒的几何外形，确定道砟颗粒的边界，采用三维激光扫描仪确定道砟颗粒的形貌，如图1所示，采用道砟颗粒边界，可对其内部进行球体颗粒填充，并通过颗粒在边界内膨胀的方法得到具有镶嵌效果的道砟颗粒单元。该方法可使模型更加光滑，并接近真实的道砟颗粒形态，如图1（d）所示。在非规则形态道砟颗粒的球体镶嵌构造中，球体单元的粒径越小，构造的道砟颗粒越真实，但所采用的球体单元也越多，从而降低计算效率。考虑以上因素，每个道砟颗粒采用10～20个非均一粒径球体单元进行随机构造，将镶嵌球体单元设为刚体颗粒，在描述道砟颗粒现象时具有一定的局限性，但仍可以有效地表征道砟颗粒的非规则几何形态和受力自锁特性。道砟颗粒间的接触力基于镶嵌球体单元，采用非线性粘弹性模型进行计算［1，18］。

2.2 无砟道床的有限单元法

无砟道床由混凝土或沥青混合料组成，可采用有限单元法进行计算。目前，有很多可以求解连续结构动力学特性的离散差分格式，如中心差分法、Wilson-θ方法和Newmark方法等。为提高计算效率，本文采用隐式积分的Newmark方法对无砟道床部分进行动力学分析，采用20节点等参单元构造无砟道床。

常用的Newmark法的形式为

式中 Fk＋1为tk＋1时刻结构的载荷，δ和α均为Newmark法的计算参数，当满足δ≥0.5和α≥0.25（0.5＋δ）2时，计算可无条件稳定，通常取δ＝0.5和α＝0.25。

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

将计算得到的xk＋1代入式（2，3），可分别得到和。

2.3 有砟-无砟过渡段的DEM-FEM耦合算法

在离散元-有限元耦合过程中，关键问题是耦合过渡区离散元和有限元间的边界条件处理，即有限元和离散元耦合区域内力学参数的合理传递，使耦合过渡区域内的力学和物理性质实现平滑过渡。由于耦合区域是离散元和有限元的重叠区域，离散元求出的道砟颗粒的合力和作用位置作为有限元的力边界条件；同时，有限单元和离散单元需要保持位移协调，即耦合区域内的离散单元位移在每个时刻应与有限单元内的节点保持一致，从而为离散元计算提供位移边界条件。

（1）耦合区域内的位移协调

有砟-无砟过渡段耦合区域内的位移分为有限元域位移uFEM和离散元域位移uDEM，其中uFEM由有限元域内实体单元节点的独立位移δi对形函数插值所得，即

形函数矩阵Ni为20个节点在每个镶嵌单元质心坐标上所取形函数的值Nαi，即

耦合区域内离散元的位移表示为

由于有砟-无砟耦合区域内离散元和有限元的位移不协调，离散元的位移可以分解为随有限元一起运动的位移uFEM和单独运动位移uSEP，即uDEM＝uFEM＋uSEP。为保证耦合区域内两种方法的位移协调，提出一种嵌入单元模型，即在耦合区域内离散元外围加入一层嵌入单元。考虑实际的力学过程，嵌入单元的单独位移uSEP＝0，因此uDEM＝uFEM，即耦合区域内的嵌入单元颗粒的位移完全由有限元的位移决定，从而保证区域内的位移协调。即

有砟-无砟过渡段离散元-有限元耦合的计算模型如图2所示。图2（a）中嵌入到有限元模型中的颗粒单元形成耦合过渡层，嵌入单元的位置关系和相互作用如图2（b）所示。采用上述道砟颗粒嵌入无砟道床板的结构形式，有助于提高有砟-无砟过渡界面上的摩擦作用，从而增强有砟-无砟过渡段刚度变化的连续性。

（2）离散元-有限元耦合区域内的参数传递

耦合区域内可以通过颗粒间的接触模型来确定离散元嵌入单元颗粒的合力，并将该合力作为有限元的力边界条件进行计算。在完成离散元每个时间步的计算后，可以确定其合力FDEMj，然后将其作为有限元的外载荷累加到有限元的节点上，得到有限元的等效节点载荷。单个嵌入单元与有限单元耦合如图3所示。依据能量守恒原理可以得到有限元等效节点力与嵌入单元所受合力的关系式，即

式中 Fnodal，i为有限单元节点上的等效外力，Ni为嵌入单元形心处的形函数值，j为嵌入单元编号，n为嵌入单元总个数，FDEMj为每个嵌入单元所受自由离散单元的合力。

（3）局部坐标求解

在耦合区域的参数传递中，需计算嵌入单元坐标的形函数。局部坐标通过求解方程（13）确定，即

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

图3 从坐标系 （x，y，z）到坐标系 （ε，η，ξ）的一一映射Fig.3 Mapping between global and local coordinate systems

式中xi，yi和zi为单元节点在总体坐标内的坐标，Ni为用局部坐标表示的插值函数。

令

假设函数f（ε，η，ξ）在点 （ε0，η0，ξ0）的某一邻域内连续且有2阶连续偏导数，（ε0＋l，η0＋m，ξ0＋n）为此邻域内任意一点，则有

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

即

同理可得

求解式（17～19）组成的方程组，当 （fξgεgξfε）（gηhε－hηgε）－（gξhε－hξgε）（fηgε－gηfε）≠0时，有

至此，求得嵌入单元的局部坐标。根据求得的嵌入单元合力求解有限元的边界力，完成离散元到有限元的力参数传递，进而得到有限单元的变形。根据有限元网格节点位移以及形函数更新嵌入单元的位置坐标，嵌入单元的运动作为位移边界条件传入离散元的计算，由此完成离散元和有限元之间的耦合过程。

3 有砟-无砟过渡段的DEM-FEM耦合模拟

3.1 有砟-无砟过渡段的DEM-FEM耦合模型

采用DEM-FEM耦合算法建立列车循环荷载作用下铁路道床有砟-无砟过渡段的计算模型。该模型包括无砟混凝土道床、枕木和有砟散体道床三部分，如图4所示。考虑高速铁路轨道的对称性，并节省计算时间，本文取铁路道床的一半进行建模分析，对于无砟道床的钢筋混凝土结构，采用线弹性模型分析其位移和变形沉降。在两种结构交界面上采用有砟颗粒嵌入无砟道床的方式，以增加有砟道床的咬合力。分别选取1＃和2＃位置作为列车荷载下的有砟道床和无砟道床沉降观测点。

在有砟-无砟过渡段的DEM-FEM耦合模型中，离散元模型采用镶嵌单元建模。有砟道床计算模型的长度x＝0.4m，高度z＝0.38m，厚度y＝0.23m，整个有砟道床由736个道砟单元颗粒、4291个球体单元组成。在离散元-有限元耦合区域内，共有46个嵌入单元，包括264个颗粒单元。有限元模型由道床和轨枕组成。模型采用20节点等参单元，划分为775个单元，总共4018个节点。采用线弹性模型建模分析，弹性模量为35GPa，泊松比为0.22。模型左侧为自由表面，右侧施加x方向的约束，前面施加y方向的约束，下侧施加z方向的约束。将列车荷载简化到轨枕上，在轨枕上施加z方向的荷载。由于离散元的时间步长远小于有限元的时间步长，为了提高计算效率，有限元的时间步长取为离散元时间步长的400倍。离散元和有限元的主要计算参数列入表1和表2。考虑道砟颗粒在往复荷载作用下因不断摩擦而趋于光滑，因此在离散元计算中，道砟颗粒间的摩擦系数随往复荷载次数的增加而不断降低。

图4 铁路道床有砟-无砟过渡段有限元-离散元耦合模型Fig.4 Coupled discrete-finite element model of transition zone

3.2 计算结果分析

为分析高速列车在250km／h的速度下通过有砟-无砟过渡段时的动力响应，在枕木上施加幅度23kN，频率0.36Hz的循环荷载。荷载的时程曲线如图5（a）所示，循环90次，时间为32.4s。一次循环中列车荷载会有四个峰值，时程曲线如图5（b）所示。

表1 道砟离散元计算参数Tab.1 Computational parameters of ballast in DEM

表2 无砟道床有限元计算参数Tab.2 Material parameters of track slab in FEM

图5 施加的列车循环荷载Fig.5 Applied cyclic load of train

在列车循环荷载的作用下，有砟-无砟过渡段1＃和2＃观测点的沉降过程如图6所示。有砟道床在循环加载的初期，道床的沉降量变化明显并且一直增大，之后随着循环次数的增加，沉降量趋于平稳并且缓慢下降，最终沉降量达到15mm左右，如图6（a）所示。无砟道床则在前几次加载后很快趋于平稳，之后保持一个相对稳定的状态，最终总体沉降量较少，约为0.005mm，如图6（b）所示。将两种道床结构对比发现，无砟道床很快达到稳定状态，而有砟道床则需要较长时间趋于稳定，同时在稳定之后两种结构会有一个明显的沉降差。这是由于在加载之前有砟道床的道砟石子排列比较松散，颗粒之间有较大的空隙，前几次的循环加载可以使有砟道床的颗粒重新排列，颗粒之间空隙变小，道床进一步压实，产生较大的沉降。对于无砟道床，是由混凝土或沥青混合料构成的轨道道床结构，在循环加载过程中，道床的杨氏模量保持不变，刚度较大，道床可以较快达到平稳。

循环荷载作用下，枕木的荷载-沉降量曲线如图7所示。由于一个时程内列车荷载具有两个峰值，并且在两个峰值之间有一次卸载过程，荷载-沉降量曲线也表现出相应的规律。在前期循环加载过程中，有砟道床沉降量明显，表明整体道床发生了明显的塑性变形；随着加载次数增加，颗粒进一步压实，整体道床平缓沉降，单次循环加卸载曲线所包围的面积也逐渐减少。有砟道床载荷-沉降量曲线的斜率逐步增大并趋于稳定，表明道砟材料的宏观弹性模量逐渐增大并趋于稳定，整体道床逐渐表现出非线性弹性材料的性质。

图6 铁路道床过渡段沉降量曲线Fig.6 Settlement of ballasted-ballastless track transition zone

图7 循环载荷作用下载荷-沉降量曲线Fig.7 Force-settlement curve under cyclic load

图8 有砟道床力链分布和无砟道床应力分布Fig.8 Force chain distribution and stress contour

图9 有砟道床力链分布和无砟道床位移云图Fig.9 Force chain distribution and displacement contour

循环加载达到峰值时，有砟道床中的力链分布和无砟道床中的应力分布如图8所示。有砟道床的力链主要分布在枕木的下方，呈梯形分布，从而将枕木上受到的载荷分散在道床上，荷载主要由道砟颗粒骨架承受。无砟道床的应力主要分布在道床的下部区域，应力分布区域与有砟道床道砟骨架相对应。无砟道床左侧边界为自由表面，仅受到有砟道床的约束，右侧边界受水平方向的位移约束，有砟道床的约束效果低于无砟道床的约束效果，因此右侧应力略大于左侧应力。从力链分布可以看出，梯形力链骨架并不是左右对称的，右面的力链分布较大，主要是由于在耦合区域嵌入单元受有限元约束，相当于增加了有砟道床和无砟道床交界面上的粗糙度，导致耦合区域的道砟颗粒具有较大的咬合力，从而限制道床过渡段附近道砟颗粒的流动，一定程度上减少两种道床的沉降差。

循环加载达到峰值时，无砟道床x方向的位移主要集中在道床的左半部分，位移变形与道砟力链分布相对应，如图9（a）所示。由于道砟颗粒荷载作用于无砟道床左侧表面，导致道床左侧表面x方向位移变形较大。无砟道床整体位移较大的区域位于枕木下方，同时道床左侧位移略大于右侧位移，如图9（b）。这也是由于无砟道床左侧位移约束略小于右侧位移约束，与道床内下部应力分布相对应。

4 结 论

针对铁路道床有砟-无砟过渡段在循环载荷下的沉降特性，发展了离散元-有限元耦合方法。对道砟颗粒和无砟道床分别采用离散元和有限元方法，并针对有砟-无砟过渡段的道砟颗粒嵌入模式，在离散元-有限元耦合区域发展了一种参数传递的耦合算法。计算结果表明，有砟道床的沉降量显著，主要是由于道砟颗粒的重新排列，在加载后期趋于稳定；无砟道床沉降较小，并在前几次加载沉降后很快达到平稳。在加载前期，有砟道床整体产生塑性变形，但随着加载的持续，道砟材料的宏观弹性模量逐渐增大并趋于稳定，整体道床逐渐表现出非线性弹性材料的性质。有砟道床力链主要分布在枕木的下方，荷载主要由道砟颗粒骨架承受。有砟-无砟过渡段嵌入单元的存在限制了道砟颗粒的运动，力链呈非对称梯形分布，在一定程度上减少了两种道床的沉降差异。无砟道床的应力和位移较大区域位于枕木下方，应力分布同样呈梯形分布，与道砟力链分布趋势相对应。