基于DEM-MFBD方法的有砟轨道路基不均匀沉降影响分析

余翠英 ,雷红博 ,罗文俊 ,马斌

（1.华东交通大学 a.理学院；b.土木建筑学院，南昌 330013；2.中铁二局集团有限公司，成都 610013）

在铁路运输高速化、重载化背景下，有砟轨道因其经济性及适应性较强，在世界范围内广泛铺设。由于地质条件复杂及其基床动应力较大[1]，容易发生路基不均匀沉降及轨枕空吊等病害[2]。路基不均匀沉降已成为影响有砟轨道结构性能和耐久性的关键问题[3]。

传统试验方法成本高、可重复性差，难以观测道床细观力学行为[4]，当今大多数学者采用数值模拟方法开展轨道沉降受力变形及动力响应研究[5]，尤其是离散元法，作为一种反映颗粒受力与运动状态的数值方法，可准确表达散粒体道床的力学特性，被广泛应用于有砟轨道相关研究[6]。但是，单纯离散元方法难以准确描述轨道结构的连续体特征，部分学者采用离散元—有限元等耦合方法尝试开展有砟轨道性能研究[7]。离散元—多柔性体动力学(DEM-MFBD)耦合方法作为一种可以有效反映离散体和连续体之间的接触及其各自的物理力学特性的仿真方法，目前鲜有应用于有砟轨道的相关研究。

针对有砟轨道上部结构连续性与道床散粒体特征，笔者采用离散元—多柔性体动力学耦合方法，建立有砟轨道结构沉降损伤模型，研究路基不均匀沉降对有砟轨道结构整体受力变形的影响，以期为路基沉降控制和有砟轨道结构综合养护、维修提供理论依据和技术支持。

1 离散元与多柔性体动力学耦合方法

表1 为不同数值模拟方法的特性及适用范围的对比。通过对比表1 中多种耦合方法的特性，采用离散元—多柔性体动力学耦合方法探讨有砟轨道路基沉降大变形问题。

表1 不同数值模拟方法的对比Table 1 Comparison of numerical simulation methods

离散元接触模型采用Hertz-Mindlin 非滑动模型[8]，该模型适于有砟轨道道砟非黏性颗粒的计算。由图1 可见，道砟颗粒簇由若干个球单元镶嵌而成，道砟颗粒簇之间的接触即转化为若干个球单元间的两两接触，接触单元间的接触力由单元间的相互重叠量及相对速度决定。

图1 Hertz-Mindlin 非滑动接触模型Fig.1 Hertz-Mindlin (no slip) contact model

离散元—多柔性体动力学耦合求解中，墙单元网格面上的接触力往往不在网格节点上，且每个三角形网格上可能有多个接触点，需采用形函数插值法，结合虚功原理，将颗粒作用在耦合面上的接触力转换为等效节点力，其整体坐标与局部坐标转换公式[9]为

式中：Ttrans3×3为转换矩阵，则耦合边界上等效接触力计算公式[9]为

2 模型的建立

图2 为DEM-MFBD 耦合过程流程图，有砟轨道结构建模过程分为离散元道床建模和多柔性体动力学上部结构建模，通过轨枕墙单元作为耦合过渡边界，进行轨道上部结构和离散元道床耦合分析。

图2 离散元−多柔性体动力学模型耦合流程图Fig.2 Flowchart of coupling DEM-MFBD model

2.1 有砟轨道道床离散元模型

为准确模拟道床内摩擦力和抗剪强度等力学行为，应考虑道砟的真实几何特征与粒径级配[10-11]。选取华东交通大学轨道交通实验基地有砟线路的特级道砟为样本，利用三维扫描技术获取的不规则道砟外形轮廓作为道砟颗粒外形边界，填充外形边界，建立道砟颗粒簇。如图3 所示，根据《铁路碎石道砟》(TB/T 2140—2008)[11]规范要求，首先剔除道砟针状、片状等异常样本；其次，结合文献[10,12]的经验并多次优化计算，最终选取10 种道砟颗粒样本，需要注意的是，样本需借助NX 修改粒径比例，以满足颗粒级配要求，故选取样本仅考虑其外形特征影响。其中，每个道砟颗粒平均由9.23 个大小不同的球单元构成，用以平衡计算精度及计算效率。

图3 道砟颗粒样本建模过程Fig.3 The modeling process of ballast-particle sample

基于离散元的方法目前存在计算耗时高、占用内存大等局限性，加上路基不均匀沉降问题对模型纵向长度的要求[13]，建立全尺度轨道模型难以同时满足计算效率与精度的双重要求，为此，需进行轨道横断面上的简化。参考既有研究[6]并综合考虑计算精度及计算效率，最终确定轨枕和道床横断面均缩尺为0.2 m 宽的2.5 维简化模型，如图4 所示。离散元道床模型长为6 m、宽为0.2 m、厚度为0.33 m，模型中共9 558 个道砟颗粒簇，考虑为特级道砟级配，如图5 所示。

图4 离散元道床模型Fig.4 The discrete element model of ballast bed

图5 颗粒级配曲线Fig.5 Particle gradation curve

2.2 有砟轨道上部结构多柔性体动力学模型

有砟轨道上部结构采用多柔性体动力学建模（见图6），将钢轨、轨枕考虑为连续体，钢轨视为柔性体，轨枕视为刚体。选用60 kg/m 钢轨，长12.5 m，为实体单元；Ⅲ型弹条扣件，刚度为6×107N/m，采用bushing 单元模拟；Ⅲ型混凝土轨枕，轨枕间距为0.6 m，将轨枕作为耦合边界导出为墙单元。为保证计算效率，仅将路基沉降区段6 m 范围的道床进行离散元建模。考虑到有砟道床对上部结构的支承作用及不承受拉力的特点，且模型试验中离散元区域外道床受力较小，考虑为无张力弹性，忽略其塑性变形。

图6 离散元−多柔性体动力学耦合模型Fig.6 The coupling model by DEM-MFBD

模型简化关键在于能否准确反映道床应力状态。假定列车荷载及上部结构自重由轨枕传递至轨枕底面压力为均匀分布，基于轨枕宽度不变，轨枕长度简化为0.2 m，为保证道床顶面所受平均应力不变，将轨枕底面压力等比例缩减。相对于列车荷载，钢轨和扣件自重影响较小，故仅考虑轨枕底面压力受列车荷载单一因素影响，对输入荷载进行等效折减。

式中：F为输入总荷载；P为静荷载；Fmax为动荷载幅值；k为缩减系数；f为输入荷载频率。

钢轨外形尺寸与实际情况一致，由材料力学挠曲线近似微分方程可知，钢轨抗弯刚度与弹性模量成正比，外荷载等效折减后应对钢轨进行刚度折减，以准确反映轨道的抗弯特性。

式中：E为钢轨输入弹性模量；E′为设计弹性模量。

3 参数标定及模型验证

3.1 模型参数及标定

校准离散元模型参数的方法通常可分为两种：1）基于微观视角的直接测量法，直接在颗粒或接触层面测量参数值；2）基于宏观视角的方法，参数值在宏观角度上反映颗粒堆积体的整体力学行为。显然，直接测量法难以反映所有相关的物理效应，在计算效率及精度的限制下，很难完全恢复颗粒的不规则形状，且颗粒间的相互咬合也会被简化。在宏观方法中，不同参数组合可以得到相同的模拟结果，但参数的真实物理意义不明确[14]。为应对参数组合的模糊性，需通过多种标准化试验进行验证。结合文献[10,12]的参数选择，在多次试算后确定模型参数，见表2。

表2 模型参数Table 2 Model parameters

为验证计算参数反映道砟宏观剪切特性的正确性，首先，参考文献[10]的方法进行三轴压缩试验的数值模拟。如图7 所示，三轴压缩试验样本尺寸为0.4 m×0.2 m×0.2 m，以上述离散元道砟模型建立样本，密度为1 410 kg/m3，孔隙率为45.8%。对模型施加50 kPa 的围压，将盖板以缓慢的匀速运动下压至应变为16.8%，记录盖板与道砟间的接触力，以计算偏应力及应力比；进一步地，参照文献[12]模拟多级剪切试验，模型包括上下两个箱体，在箱体内部装满上述离散元道砟模型，样本尺寸为0.3 m×0.3 m×0.2 m，密度为1 448 kg/m3，孔隙率为46.3%，对上方箱体向下阶段性施加16.5、31.5、56.5、81.5、106.5 kPa 的压强，并以0.01 m/s 的速度匀速平移下方箱体，记录上方箱体侧壁压强，以计算切应力。

图7 数值样本Fig.7 Numerical sample

水平应力σ2=σ3，参考文献[10]中的方法，通过平均有效应力p'及差应力q计算内摩擦角ϕ。

式中：σ1和σ2分别为轴向应力和水平应力。

图8 所示为三轴压缩试验的差应力及应力比与竖向应变的关系，图9 所示为多级剪切试验中的剪应力变化，数值模拟计算抗压试验中内摩擦角约为41°，多级剪切试验中内摩擦角为50.5°。可见，数值模拟结果与相关试验[10,12]拟合较好，本文结果略大于文献[10,12]结果，考虑到样本级配的差异及本文道砟颗粒不可破碎，误差在可接受范围内，说明选取的本构模型可以还原道床的宏观特性。

图8 三轴压缩试验数值样本的偏应力及应力比Fig.8 The differential stress and stress ratio of the numerical sample by triaxial compression test

图9 多级剪切试验数值样本的应力—应变行为Fig.9 The stress-displacement behavior of numerical sample by multi-stage shear test

3.2 模型验证

为验证模型表达道床应力状态的正确性，基于提出的荷载折减方法，并考虑轮轨力为枕上压力的3 倍[15]，对道砟箱试验[16]工况进行模拟，将模型计算结果与试验结果[6,16]进行比较。图10 显示了列车荷载沿轨道纵向5 根轨枕分散传递至道床（红线表示接触力走向，红色越明显表示接触力越大），中间3根轨枕承担75%以上的列车荷载，与既有研究结论一致[6]；模型测得最大接触力为4 287.3 N，文献[16]中最大接触力为4 233 N，误差约为1.3%，表明模型能够反映真实道床应力状态。

图10 道床接触力链Fig.10 Contact force chain of ballast bed

为进一步验证模型反映轨道结构变形特性的正确性，参照同济大学开展的有砟轨道路基不均匀沉降实尺模型试验[4]，对钢轨中心相距1.8 m 的两点施加正弦荷载以模拟单一转向架的动力作用。荷载分为静载与动载两部分，其中静载为45 kN，动载幅值为49 kN。路基不均匀沉降为余弦型，路基沉降波长为2.4 m，沉降幅值为25 mm。通过试算并参考文献[17]，最终确定荷载频率为5 Hz，加载40 个完整周期，可得到钢轨中心点随加载过程的沉降曲线及基本稳定状态下的钢轨沉降变形曲线。

由图11 可以看出，钢轨顶面中心点的竖直沉降量严格地随着周期性循环荷载的作用而产生周期性变化，并产生弹性变形；由于钢轨累计沉降量不断增大，加载初期钢轨沉降量较大，其增长趋势随循环加载次数的增加逐渐趋于平缓，与相关研究中对道床初期沉降与道床塑性变形的结论一致[18]。

图11 轨面中心点随循环加载过程的沉降曲线Fig.11 Settlement curve of central point of rail with loading process

如图12 所示，将数值模拟得到钢轨沉降变形与文献[4]对比，仅两数据点相对误差较大，为11%和11.8%，其余数据点相对误差范围为0.2%～6.9%。数值模拟中钢轨最大沉降量为7.47 mm，略小于文献[4]中的结果，沉降曲线趋势性拟合较好。考虑到试验的离散性和趋势性，模拟结果合理，数值模型可行，能够反映道床真实应力状态、路基不均匀沉降下轨道结构变形特性及沉降产生的历程及机理。

图12 钢轨沉降变形曲线与文献[4]的对比Fig.12 Comparison between the rail settlement deformation curve and literature [4]

4 计算结果分析

4.1 沉降波长对轨道受力变形的影响

为研究路基不均匀沉降波长下轨道的受力变形规律，鉴于文献[4]足尺试验工况中沉降波长均小于3 m，工况设定考虑波长范围为5 根轨枕间距，故设定路基沉降幅值为10 mm，沉降波长分别为0.6、1.2、1.8、2.4 m 等4 种工况，加载50 个周期，荷载频率为5 Hz，分析各路基沉降波长下轨道的结构受力变形特征。

如图13 所示，轨枕沉降与钢轨沉降变形大致相同。路基沉降波长经道床、轨枕传递到钢轨表面会有一定程度的扩大。测得路基—钢轨沉降波长扩散角θ为30.3°～34.4°，略大于文献[4]试验结果，结果表明，路基—轨面沉降波长扩散角受路基沉降波长、幅值变化影响较小，由于道床的塑性形变，路基沉降幅值及沉降面积传递至钢轨表面均有减小。

图13 沉降幅值为10 mm，波长为2.4 m 时轨道的整体变形Fig.13 The overall deformation of track structure with settlement amplitude10 mm and wavelength 2.4 m

如图14 所示，路基沉降幅值为10 mm，波长由0.6 m 增大至2.4 m，钢轨不平顺幅值由1.82 mm 增大至6.34 mm，增大了约2.5 倍。同时，钢轨两端的上拱减少，钢轨整体变形呈下沉趋势。当路基沉降波长为0.6 m 时，钢轨最大沉降量为−1.79 mm，而钢轨上拱最高达到1.81 mm，由上拱导致的不平顺较为严重。当路基产生小波长沉降时，钢轨上拱导致的不平顺将影响列车运营的舒适性。

图14 不同路基沉降波长下钢轨变形Fig.14 Rail deformation under different subgrade settlement wavelength

由图15 可知，在路基沉降区相邻轨枕处，道床应力集中现象明显，随着路基沉降波长的增加，应力集中位置由7 号轨枕逐渐外扩至8 号轨枕处，最大接触力整体呈上升趋势。正常情况下，最大接触力为5 841.5 N，位于7 号轨枕下方；沉降波长为0.6、2.4 m 时，最大接触力分别为7 227.8、10 405.0 N，分别上升了23.7%及78.1%。图15 中数字为该位置轨枕与道床间接触力，红色突出显示为其中最大的接触力。

图15 不同路基沉降波长下道床接触力链Fig.15 Contact force chain of ballast bed under different subgrade settlement wavelengths

4.2 沉降幅值对轨道受力变形的影响

为研究路基不均匀沉降幅值对轨道受力变形的影响规律，基于工务维修规则和易操作性，同时借鉴文献[4]的统计和现场调研，设定路基沉降波长为1.8 m，沉降幅值为10、15、20 mm 三种工况，加载情况同上，分析不同路基沉降幅值下轨道结构受力变形特征。

如图16 所示，随着路基沉降幅值由10 mm 增大至15、20 mm，钢轨不平顺幅值分别增大了16%和29%。钢轨沉降波长略有扩大，波长扩散角θ随之减小，而端部上拱差别不大。

图16 不同路基沉降幅值下的钢轨变形Fig.16 Rail deformation under different subgrade settlement amplitudes

如图17所示，随着路基沉降幅值的增大，7号轨枕与道床的接触力从7 679.3 N 逐渐降低至3 693.3 N，同时，8 号轨枕与道床接触力从8 202.5 N 升高至10 827.5 N，相较于正常线路，上升了85.4%，应力集中位置逐渐外扩（由7 号轨枕向8 号轨枕迁移）。

图17 不同路基沉降幅值下的道床接触力链Fig.17 Contact force chain of ballast bed under different subgrade settlement amplitude

4.3 路基沉降对轨枕空吊的影响规律

不同沉降波长及幅值的组合数值模拟结果如表3 所示。由表3 可知，轨枕空吊情况可由钢轨沉降面积与路基沉降面积之比(S1/S0)划分，当S1/S0>1.03 时，均没有发现轨枕空吊；而当S1/S0≤1.03时，均产生不同程度的轨枕空吊。说明通过轨面沉降面积与路基沉降面积的比值反映轨枕空吊情况具有合理性和可行性，但对于该指标反映发生轨枕空吊的临界值及影响机理仍需进一步研究。

表3 不同工况计算结果Table 3 The calculation results of different working condition

路基沉降波长为1.8 m 时，随着沉降幅值由10 mm 增大至15、20 mm，S1/S0由1.22 下降至0.95和0.83，并产生了不同程度的轨枕空吊。沉降幅值为15 mm 时，轨枕在列车荷载作用下与道床接触，而随着列车荷载的减小，逐渐与道床分离，即为非完全空吊；沉降幅值继续增大至20 mm 时，在列车荷载作用下，轨枕始终不与道床接触，即为完全空吊。当路基沉降幅值为15 mm 时，沉降波长由0.6 mm 增大至2.4 mm，轨枕与道床间隙由0 mm 增大至4.3 mm，轨枕空吊加剧。

在路基沉降波长为0.6～2.4 m 的条件下，沉降幅值为10 mm 时，均未发生轨枕空吊；沉降幅值为15 mm 时，则发生了轨枕空吊。可见，路基沉降波长介于0.6～2.4 m 之间时，轨枕空吊临界值的路基沉降幅值应介于10～15 mm 之间，建议将路基沉降幅值控制在10 mm 以内。

5 结论

1）基于DEM-MFBD 耦合方法，提出一种2.5维简化模型，该模型可以正确反映有砟道床的宏观力学行为，可应用于有砟轨道沉降变形机理分析。

2）钢轨沉降面积与路基沉降面积的比值可以反映轨枕空吊情况；当路基沉降波长范围小于4 根轨枕间距时，轨枕空吊的沉降幅值临界值介于10～15 mm 之间；建议沉降幅值控制在10 mm 以内，避免轨枕空吊。

3）路基沉降波长和幅值的增加均导致轨道不平顺明显增大，并加剧轨枕空吊现象，建议加强路基沉降监测。

4）路基不均匀沉降促使道床产生应力集中现象，道床应力增大；随着沉降波长和幅值的增加，进一步导致道床应力集中范围外延，长期运营将导致道砟破碎，道床板结；期望为工务维修“科学修，精确修”提供理论依据。