图1 界面单元及耦合原理
Fig.1 Interface element and coupling principle
肖军华，张 德，王延海，郭佳奇
(同济大学道路与交通工程教育部重点实验室，上海 201804)
摘 要： 引入离散元法(DEM)-有限差分法(FDM)耦合算法对普通铁路碎石道床-土质基床的界面应力进行分析。首先，基于激光扫描以及室内三轴试验，实现碎石道砟颗粒二维精细化离散元建模，并通过设置界面单元进行道砟与土两类不同介质层间速度与力的相互传递，实现离散元法与有限差分法耦合，建立了轨枕-碎石道床-连续土体耦合模型，并通过现场实测结果验证模型的合理性。在此基础上，计算列车通过时耦合模型中基床动应力分布特征，讨论轨枕-道床接触状态对其影响，并将计算结果与轨枕-道床-基床的整体有限元模型的计算结果进行对比。结果表明，轨枕-道床未完全密贴接触状态下，道床中力链呈离散型传递，反之，道床中力链呈向下扩散发展；两种接触状态下，基床表面应力分布差别显著，但是应力峰值接近。耦合模型得到的基床表面应力峰值显著高于有限元模型的计算结果；耦合模型中因道砟散体特性引起的基床表面局部应力集中，在路基中的影响深度约为道砟最大粒径(63 mm)的8倍，即近似基床表层深度。
关键词： 有砟道床；离散元法；有限差分法；耦合算法；基床动应力
受过去经济、技术条件的限制，我国普通铁路既有线广泛采用细粒土作为路基基床填料，从而形成了碎石道床-土质基床直接接触作用的结构形式，导致线路运营中产生大量基床病害，如基床下沉、道砟陷槽以及翻浆冒泥等，给铁路运输带来了不利影响 [1] 。而随着既有线铁路实施重载化的发展战略，路基病害将进一步加剧。欧美等国家的传统铁路也存在类似的路基结构形式以及产生了相应的病害问题 [2―3] 。
道床-基床界面动应力是列车荷载通过钢轨、轨枕的分散作用并经过道床的扩散传递形成的，是基床病害形成的重要因素。目前，基于扩散角的经验法、或基于 Odemark的模量与厚度当量换算(对于道床)进而采用Boussinesq公式计算的理论解析法 [4] 、或基于弹塑性理论的轨道-路基动力有限元法等都可以对路基面应力进行量化。然而，这些方法均是从宏观角度(即把碎石道床看作一个整体)来进行应力分析的，不能反映列车荷载在道砟颗粒体系中传递和道床-基床界面接触应力形成演化的细观实质，更不能体现随机分布的道砟形状、尺寸和粒径等对界面接触行为的影响。而采用传统土压力盒的测试方法通过现场或模型试验对有砟轨道路基面应力测量得到的结果也仅是在路基面一定区域均布化后的应力水平。
因此，基于有砟道床宏观假设及传统土压力盒测量技术得到的列车荷载下路基面动应力结果可能会因为忽视碎石道床-土质路基界面局部重要的接触应力信息，从而引起基床病害研究上的某些偏差。近年来，Abadi等 [5] 、Aikawa [6] 及 McHenry等 [7] 研究者考虑道砟细观结构特征，采用压敏纸、智能传感轨枕以及薄膜传感器等测试手段研究了碎石道床层与上下结构层间的界面接触应力，并发现轨枕与碎石道床间较低的有效接触面积导致界面存在较高的压应力，与传统计算方法所得平均应力相差数倍。
与量测方法相比，数值方法可以观测颗粒在荷载作用下的滑移、旋转等运动规律，可从细观层次来揭示荷载作用下颗粒体系内力传递和界面接触行为动态演化的机制和规律。其中，离散单元法在研究铁路道砟等散粒体的力学特性方面得到了国内外广泛应用。例如，Indraratna等 [8] 较早地采用离散元法对道砟颗粒细观相互作用机理进行了研究，国内高亮等 [9] 、井国庆等 [10] 、张徐等 [11] 也通过离散元模拟有砟道床动力特性。近年来，国内外研究者将离散元-有限差分耦合算法引入到工程计算中，分析一些典型的离散-连续介质相互作用问题的机理，在计算方法上，基本采用了较成熟的离散元-有限差分(PFC-FLAC)耦合方法 [12―13] 。同样地，针对铁路有砟道床与上下部结构间的相互作用，邵帅等 [14] 、高亮和徐旸 [15] 、Ngo等 [16] 采用类似耦合分析方法进行了研究。然而，关于道砟颗粒形态定量分类和精细模拟、列车荷载作用下道床-基床界面接触应力的宏观规律及其细观机理、以及采用离散-连续耦合分析方法与常规假设道床为连续体的有限元分析方法对基床表面应力分析有何差异等问题，仍可开展深入研究。
本文在现有文献研究基础上，结合扫描技术建立精细化道砟颗粒二维离散元模型，并通过三轴数值模拟及试验验证确定道砟离散元细观参数，在此基础上，采用离散元(DEM)模拟散体道床，采用有限差分法(FDM)模拟下部路基和基床，并基于接触单元在碎石道床与基床接触层进行离散元接触力与有限差分网格节点位移信息的持续相互传递，从而实现多层介质结构力学耦合。在此基础上，研究散体道床与连续土质基床接触面在列车荷载作用下应力特征和细观机理，并与传统对道床的连续性假设获取的基床面动应力作对比。相关结论可以为普通铁路基床病害预测和整治提供科学依据。
通过在离散单元与连续介质相互作用界面上设置一系列界面单元来实现两者的耦合。界面单元将离散介质内作用在界面单元上的力作为有限差分网格耦合界面的力边界条件，在同一时步将有限差分网格节点与界面单元耦合的位移作为离散单元界面上速度边界条件。界面单元及耦合原理如图1所示 [16] 。
图1 界面单元及耦合原理
Fig.1 Interface element and coupling principle
图中界面单元尺寸(L)与连续体单元网格相等，界面单元节点为X 1 、X 2 ，离散单元颗粒与界面单元相互接触力 F i c 、 F 1 与 F 2 ， F i c 作用在颗粒形心处， F 1 、 F 2 分别作用在界面单元两节点处，且可以沿界面单元法线与切线方向进行分解：
由X 1 节点力矩平衡可得：
或者
其中 t × t =0, n × n =0且令 t × n = k ，由式(4)得
另外，将切向力按照距离函数分配到节点上得到F 2X ，如式(6)：
由力平衡可得：
可将式(7)的力平衡方程式沿界面单元切线及法线方向表示为：
综合式(8)和式(9)可得：
由式(5)～式(11)获取界面单元节点力F 1X 、F 1Y 、F 2 X 和F 2Y 。
离散-连续耦合分析主要由如下步骤实现：
1)分别建立轨枕、道床的离散元模型以及基床、路基的有限差分模型；
2)通过界面单元，将离散元模型边界墙体运动隶属于有限差分边界网格线段；
3)施加荷载，离散元与有限差分模型各求解一步；
4)离散元模型将力分配到界面单元节点上，并传递到有限差分边界上作为节点荷载；
5)有限差分模型接收界面单元传递节点荷载，计算边界节点速度，并通过界面单元传递到离散元模型墙体相应节点上；
6)离散元与有限差分模型同时计算下一加载步，并返回步骤4)，直至加载结束。
其中，步骤 4)～步骤 5)通过 Socket(I/O)接口完成离散元模型与有限差分模型每个时步的数据交换，如图2所示。
图2 DEM-FDM耦合模型数据交换过程图
Fig.2 The DEM-FDM coupled data exchange process
已有研究 [17] 指出，道砟颗粒形状对颗粒力学行为产生重要影响，对道砟颗粒外形特征以及细观力学参数的合理描述是实现散体道床数值模拟的关键。为此，本文采用Creamform手持激光扫描仪对单个道砟颗粒进行扫描，获取该道砟颗粒三维空间点云坐标，并基于点云坐标构建三角形薄片指数方程 [18] 重构道砟颗粒形状轮廓。图 3中颗粒形状由5118个扫描顶点及10232个三角网格拓扑面组成，在此基础上取粒径(d)、球度指数 [19] (Sphericity，S p )与棱角尖锐系数 [20] (Ψ)量化道砟颗粒整体形态。
图3 单颗粒激光扫描及道砟三维点云图
Fig.3 Laser scanning and 3D point cloud of single ballast
以图3中扫描道砟颗粒为例，由各表面点坐标计算颗粒最长轴(P + P - )为 Φ 1 ，沿长轴 Φ 1 取与之垂直的最长轴(AC)为Φ 2 ，交点为I，并作垂直横切面S，在垂直剖面上选取与 Φ 1 、Φ 2 同时垂直轴(BD)为 Φ 3 ，如图 4(a)所示。Φ 1 、Φ 2 、Φ 3 可由式(12)表示：
球度指数 S p 以及棱角尖锐系数 Ψ可由式(13)表示：
因颗粒在与短轴 Φ 3 轴垂直的平面方向最容易稳定，且棱角尖锐系数主要依赖长轴 Φ 1 的信息，为此，选取包含Φ 1 轴、Φ 3 轴所在平面作为颗粒二维投影面；对投影面轮廓进行分析，并采用与轮廓线相切的圆盘组合黏结成无重叠二维离散单元来描述道砟颗粒不规则形状，如图4(b)。
图4 颗粒形状指标计算及二维模型
Fig.4 Shape index of particle and 2D numerical model
扫描道砟粒径(d)满足《铁路碎石道砟》(TB/T 2140－2008) [21] 一级碎石道砟要求，粒径范围为16 mm、25 mm、35.5 mm、45 mm、56 mm及63 mm。每个粒径范围内，综合S p 与Ψ对道砟样本形状可分为 4类(类三角形、基本型 [22] 或类矩形、扁平型以及狭长型)，共形成 20种道砟颗粒形状。其中25 mm～35.5 mm粒径范围形状分类如表1(a)所示，其余各粒径颗粒二维离散单元如表1(b)所示。
表1 典型形状道砟颗粒二维离散单元
Table1 2D DEM model of typical ballast particles
离散元模拟生成的道砟颗粒级配曲线与试验值、规范值进行比较，如图 5。说明本文对道砟颗粒的形态量化及离散元模拟是合理的。
图5 试验及离散元数值模拟级配曲线
Fig.5 Gradation of testing and numerical simulation
通过大型三轴试验标定道砟细观参数，试验采用美国 GCTS振动大型三轴仪，试样高度H=600 mm，直径D=300 mm，如图6所示。对试样分层压实，密度为1533 kg/m 3 。之后对试样施加周期正弦荷载，围压保持为60 kPa，最大加载轴向应力幅值为 420 kPa，加载频率为 5 Hz。道砟数值模型采用伺服控制实现围压加载，并用墙体模拟竖向加载压头。
根据试验得到的循环加载累积应变曲线对离散单元模型细观参数进行标定。其标定的道砟颗粒细观参数见表 2，数值计算结果与试验加载结果对比如图7所示。
图6 三轴试验及离散元数值模型
Fig.6 Triaxial testing and DEM simulation model
表2 道砟颗粒细观参数
Table2 Mesoscopic parameters of ballast
图7 三轴试验与数值模拟结果对比
Fig.7 Comparison of triaxial test and numerical simulation
由于主要关注列车垂向荷载下道床-路基的相互作用应力，并假设沿铁路线路方向为无限长，因而可忽略沿线路纵向的道砟颗粒运动，而仅考虑线路横断面内道砟颗粒的平动及转动。为此，可将铁路线路简化为二维问题进行数值模拟。
在前述开展的精细化道砟二维离散单元建模基础上，取道床顶面宽度为3.4 m，厚度为0.35 m，边坡坡度为 1∶1.75。为保证边坡形状，先设置与边坡角度相同的临时倾斜墙体约束，同时采用落雨法生成道砟颗粒，并按照《铁路轨道工程施工质量验收标准》 [23] (TB 10413―2003)规定的道砟铺设作业中须采用压强不小于160 kPa的机械碾压，利用伺服机制设置顶部墙采用200 kPa幅值，5 Hz动载对道床道砟分3层压实。
待道床颗粒稳定后，撤去边坡及顶部压实加载墙体。并在道床上方生成规则排列的圆盘颗粒簇cluster模拟III型混凝土轨枕(轨枕长度为2.6 m)，轨枕-碎石道床离散元数值模型如图8所示。道砟颗粒细观参数采用本文试验标定的模型参数，轨枕颗粒细观参数则采用文献[24]由现场动力测试数据标定的参数，如表3所示。
图8 轨枕-碎石道床离散元模型
Fig.8 Sleeper-ballast DEM simulation model
表3 轨枕离散元模型细观参数 [24]
Table3 Mesoscopic parameters of sleeper DEM model [24]
建立的土质路基有限差分模型如图9所示，路基面宽度为7.9 m，路堤高度为2.5 m，路基边坡坡度为 1:1.5，地基深度为 15 m。路基材料计算参数参考普通铁路路基土工试验数据 [25] 选取，见表4。土体阻尼比取0.1。
图9 土质路基有限差分模型
Fig.9 FDM model of subgrade
表4 土质路基基本计算参数 [25]
Table4 Parameters of subgrade [25]
在道床模型底部设置界面单元，单元尺寸与路基模型的基床表面单元网格大小一致。
根据某普通铁路有砟轨道上列车以不同时速通过产生的轨枕加速度实测值与耦合模型仿真结果进行对比，验证耦合模型的正确性。这里，以DF 11 机车，轴重23 t，列车速度110 km/h为例进行对比说明。本文将轨道结构简化为连续弹性支承梁模型 [26] ，将列车荷载简化为一系列常力P，以匀速v在钢轨上移动，得到多个轮载同时作用下，一根轨枕上的反力曲线如图10所示，详细计算方法见文献[26]。
图10 DF 11 机车以110 km/h通过时轨枕反力时程曲线Fig.10 Time-history curve of sleeper counterforce under DF 11 train at speed of 110 km/h
将所得反力时程曲线经等效平面应变处理后作为输入荷载施加到模型中 [27] ，计算得到列车经过时轨枕加速度并与实测结果对比，如图11。
图11 耦合模型计算结果与实测结果对比
Fig.11 The comparison between coupling model and measured results
由图11(a)～图11(b)可以看出，因现场线路存在不平顺，导致实测结果存在明显的毛刺现象，但现场实测轨枕加速度与耦合模型计算结果的幅值及波形相近，从而验证了耦合模型的正确性。其他不同时速的列车荷载下轨枕加速度测试和计算结果相对比，均能得到同样的验证结果。
McHenry等 [7] 采用薄膜传感器对轨枕与碎石道床间的接触应力进行测试，发现道床与轨枕接触界面的道砟颗粒状态(新道砟、中度脏污道砟以及细砂含量较多的重度脏污道砟)会明显影响轨枕与道床接触应力分布特征。因此，可推断道砟颗粒的离散特征同样会影响列车荷载在道床-基床的传递，从而影响到基床表面应力。下面考虑轨枕-道床两种不同的接触状态下，基床表面应力的特征。
在实际工程中，虽然道床按照规范要求填筑压密，但由于未经列车荷载作用的磨合，初始时道床与轨枕可能未完全密贴，存在局部空隙(如图12所示，称为“工况1”，模拟初始新道砟状态)。
图 13为该耦合模型在列车荷载(图 10荷载波形)最大值时道床中力链(即颗粒接触网络)及基床应力分布情况。可以看出，由于轨枕与道床面未完全密贴，道床主力链没有在轨枕加载点正下方形成明显的扩散传递，而是呈离散型分布，且在轨枕中部下方也有主力链出现，如图14。这与McHenry等 [7] 通过薄膜传感器进行现场试验揭示的现象一致，该试验表明，对于新道砟，列车荷载下轨枕底部的最大压应力不在轨枕加载点正下方。同时，计算还表明，由轨枕通过道床传递到基床表面的力链所形成的基床表面应力也呈离散型分布，在主力链处形成局部应力集中，对于本计算工况，基床表面应力最大值为110.5 kPa，且并非在轨枕加载点正下方。
图12 模拟初始新道砟与轨枕接触状态(工况1)
Fig.12 Simulating interaction between initial new ballast and sleeper (Case 1)
图13 轨枕-道床-基床耦合模型中道床力链及基床应力图(工况1)
Fig.13 Force chain and stress of sleeper-ballast-subgrade coupling model (Case 1
图14 工况1中轨枕及道砟力链详图
Fig.14 Force chain of sleeper and ballast for Case 1
为了对比轨枕-道床不同的接触状态对基床面应力的影响，假设铁路进入运营一段时间后，由于轨枕底部和道床表面的磨合，使轨枕与道床接触达到密实(如图15，称为“工况2”，模拟运营后道砟状态)。采用在轨枕与道床间填充细颗粒来模拟该接触状态。
图15 模拟进入运营后道砟与轨枕接触状态(工况2)
Fig.15 Simulating interaction between sleeper and ballast in operation stage (Case 2)
图16 运营期轨枕-道床-基床耦合模型中道床力链及基床应力图(工况2)
Fig.16 Force chain and stress of sleeper-ballast-subgrade coupling model in operation period (Case 2)
图 16为模拟运营后道砟状态的耦合模型在列车荷载(图10荷载波形)最大值时，道床力链及基床应力分布情况。可以看出，道砟主力链由轨枕加载点下方形成竖直向下扩散分布，这同样与 McHenry等 [7] 通过薄膜传感器进行现场试验揭示的现象一致。扩散角α范围为19°～23°，如图17。且基床应力呈对称分布，对于本计算工况，基床表面应力最大值为106.8 kPa，位于轨枕加载点正下方。因此，道床与轨枕的接触状态，直接影响基床表面应力的分布。
另外，将轨枕与道床均假设为连续体，采用有限单元法建立轨枕-道床-基床整体有限元模型，计算有限元模型中基床表面及路基不同深度处的动应力，并与离散-连续耦合模型计算结果进行对比。有限元模型中轨枕及道床参数参考既有线实际的试验数据 [25,27] ，并取与上述耦合模型中等同的轨枕及道床条件：取轨枕弹性模量为3×10 4 MPa，泊松比为 0.15，重度为 21.8 kN/m 3 ；道床弹性模量为200 MPa，泊松比为0.3，重度为20 kN/m 3 。
图17 工况2轨枕及道砟力链详图
Fig.17 Force chain of sleeper and ballast for Case 2
图18为两类数值模型(分别为有限元模型、离散-连续耦合模型)计算得到的基床表面应力(列车荷载取图10荷载波形)。可见，有限元模型计算得到的基床面最大动应力(约 55 kPa)与采用《铁路路基设计规范》(TB 10001―2016) [28] 推荐的路基动应力经验公式σ d =2.6P(1+αv)，对DF 11 机车、轴重P=23 t、列车速度v=110 km/h计算得到的路基动应力最大值59.8 kPa比较接近。
此外，耦合模型工况1、工况2得到的基床表面应力峰值相近，但分布特征差异大，表明轨枕-道床接触状态对基床表面应力最大值影响不大，但是影响其分布。然而，耦合模型得到的基床面最大应力均比有限元模型的计算值大，对于本次计算工况，相差约1倍，定性上符合McHenry等 [7] 通过薄膜传感器进行现场试验揭示的现象，即轨枕-碎石道床界面实际的压应力，可比假定轨枕与道床之间的接触应力是均匀、平均分布的情况高数倍。
再者，由基床表面应力分布宽度可以看出，有限元模型中道床中应力扩散角大致为 31°，较耦合模型中应力扩散角大。
图18 基床表面应力对比曲线
Fig.18 The contrast curves for stress of subgrade surface
图 19为两类数值模型计算的基床应力沿深度变化曲线，均取列车通过时(列车荷载取图10荷载波形)产生的峰值应力随深度的衰减情况。由图 19可以看出，距基床表面0.5 m深度时，耦合模型计算的应力值与有限元模型计算的应力值接近，因此，可认为由道砟力链传递至基床表面引起的局部应力集中，在路基中的影响深度约为道砟最大粒径(63 mm)的8倍，即近似为基床表层厚度。这说明，基于道床宏观假设得到的路基应力分布，会低估基床表层范围内的实际应力水平。
图19 基床应力沿深度变化图
Fig.19 Stress variation along with subgrade depth
(1)基于激光扫描以及室内三轴试验，实现碎石道砟颗粒二维精细化建模，并通过设置界面单元进行离散介质与连续介质间速度和力的相互传递实现离散单元法(DEM)与有限差分法(FDM)耦合，建立了轨枕-碎石道床-连续土体耦合模型。通过与现场实测结果对比，验证了模型的合理性。
(2)通过DEM/FDM耦合模型讨论轨枕与道床不同接触状态对基床应力的影响。结果表明在模拟初始新道砟与轨枕接触状态下，道床中道砟力链呈离散型分布，而模拟运营后道砟与轨枕接触状态下，道砟力链由轨枕加载点下方呈竖向扩散传递到基床表面，扩散角大约为 19°～23°。两种接触状态下，基床表面应力最大值相差不大，但应力分布差异显著。
(3)由耦合模型模拟列车通过时基床表面动应力，并与传统将道床考虑为连续介质的有限元模型的计算结果进行对比，表明耦合模型得到的基床面应力峰值显著高于有限元模型的计算值。
(4)对比有限元模型和耦合模型中基床应力随深度的变化规律，当距路基表面深度约0.5 m时，两类数值模型计算得到的应力值接近，表明由道砟力链传递至基床表面引起的局部应力集中，在路基中的影响深度约为道砟最大粒径(63 mm)的 8倍，即近似为基床表层厚度。
今后，还需开展碎石道床-土质基床界面接触应力分布的试验工作，以进一步揭示道砟-基床土两种介质的接触应力特性和细观动态演化规律。
参考文献：
[1]周镜, 叶阳升.基床结构设计的探讨[J].铁道科学与工程学报, 2004, 1(1): 1―6.Zhou Jing, Ye Yangshen.Design of roadbed structure [J].Journal of Railway Science and Engineering, 2004, 1(1):1―6.(in Chinese)
[2]Duong T V, Cui Y J, Tang A M, et al.Investigating the mud pumping and interlayer creation phenomena in railway sub-structure [J].Engineering Geology, 2014,171(3): 45―58.
[3]Li D Q, James H, Ted S, et al.Railway geotechnics [M].Boca Raton: CRC Press, 2015.
[4]刘钢, 罗强, 张良, 等.高速铁路有砟轨道路基设计荷载分析[J].铁道科学与工程学报, 2015, 12(3): 475―481.Liu Gang, Luo Qiang, Zhang Liang, et al.Analysis of the design load on the high-speed railway ballasted track subgrade [J].Journal of Railway Science and Engineering, 2015, 12(3): 475―481.(in Chinese)
[5]Abadi T, Le Pen L, Zervos A, et al.Measuring the area and number of ballast particle contacts at sleeper/ballast and ballast/subgrade interfaces [J].The International Journal of Railway Technology, 2015, 4(2): 45―72.
[6]Aikawa A.Dynamic characterization of a ballast layer subject to traffic impact loads using three-dimensional sensing stones and a special sensing sleeper [J].Construction and Building Materials, 2015, 92(S1):23―30.
[7]McHenry M T, Brown M, LoPresti J, et al.Use of matrix-based tactile surface sensors to assess fine-scale ballast-tie interface pressure distribution in railroad track[J].Transportation Research Record, 2015, 2476(1):23―31.
[8]Ngo N T, Indraratna B, Rujikiatkamjorn C.A study of the geogrid–subballast interface via experimental evaluation and discrete element modelling [J].Granular Matter,2017, 19(3): 54-1―54-16.
[9]徐旸, 高亮, 井国庆, 等.脏污对道床剪切性能影响及评估指标的离散元分析[J].工程力学, 2015, 32(8):96―102.Xu Yang, Gao Liang, Jing Guoqing, et al.Shear behavior analysis of fouled railway ballast by DEM and its evaluation index [J].Engineering Mechanics, 2015,32(8): 96―102.(in Chinese)
[10]井国庆, 王子杰, 施晓毅.多围压下三轴压缩试验与不可破裂颗粒离散元法分析[J].工程力学, 2015,32(10): 82―88.Jin Guoqing, Wang Zijie, Shi Xiaoyi.Ballast triaxial tests and non-breakable particle discrete element method analysis under different confining pressures [J].Engineering Mechanics, 2015, 32(10): 82―88.(in Chinese)
[11]Zhang X, Zhao C, Zhai W.Dynamic behavior analysis of high-speed railway ballast under moving vehicle loads using discrete element method [J].International Journal of Geomechanics, 2017, 17(7): 04016157-1―04016157-14.
[12]周健, 邓益兵, 贾敏才, 等.基于颗粒单元接触的二维离散-连续耦合分析方法[J].岩土工程学报, 2010,32(10): 1479―1484.Zhou Jian, Deng Yibing, Jia Mincai, et al.Coupling method of two-dimensional discontinuum-continuum based on contact between particle and element [J].Chinese Journal of Geotechnical Engineering, 2010,32(10): 1479―1484.(in Chinese)
[13]Indraratna B, Ngo N T, Rujikiatkamjorn C, et al.Coupled discrete element–finite difference method for analyzing the load-deformation behaviour of a single stone column in soft soil [J].Computers and Geotechnics, 2015, 63(1):267―278.
[14]邵帅, 严颖, 季顺迎.铁路有砟道床动力特性的离散元-有限元耦合模型[J].计算力学学报, 2014, 31(增刊1): 184―189.Shao Shuai, Yan Ying, Ji Shunyin.A coupled discrete-finite element model for dynamic behavior analysis of railway ballasted track [J].Chinese Journal of Computational Mechanics, 2014, 31(Suppl1): 184―189.(in Chinese)
[15]高亮, 徐旸.离散元-有限差分耦合法在铁路有砟道床研究中的应用[J].北京交通大学学报, 2016, 40(4):37―44.Gao Liang, Xu Yang.Application of coupled discrete element-finite difference method in railway ballast bed research [J].Journal of Beijing Jiaotong University,2016, 40(4): 37―44.(in Chinese)
[16]Ngo N T, Indraratna B, Rujikiatkamjorn C.Simulation ballasted track behavior: numerical treatment and field application [J].International Journal of Geomechanics,2017, 17(6): 04016130-1―04016130-12.
[17]边学成, 李伟, 李公羽, 等.基于颗粒真实几何形状的铁路道砟剪切过程三维离散元分析[J].工程力学,2015, 32(5): 64―75.Bian Xuecheng, Li Wei, Li Gongyu, et al.Three-dimensional discrete element analysis of railway ballast’s shear process based on particles’ real geometry[J].Engineering Mechanics, 2015, 32(5): 64―75.(inChinese)
[18]Anochie-Boateng J K, Komba J J, Mvelase G M.Three-dimensional laser scanning techniques to quantify aggregate and ballast shape properties [J].Construction &Building Materials, 2013, 43(6): 389―398.
[19]Krumbein W C.Measurement and geological significance of shape and roundness of sedimentary particles [J].Journal of Sedimentary Petrology, 1941,11(2): 64―72.
[20]Lee J R J, Smith M L, Smith L N.A new approach to the three-dimensional quantification of angularity using image analysis of the size and form of coarse aggregates[J].Engineering Geology, 2007, 91(2): 254―264.
[21]TB/T 2140－2008, 铁路碎石道砟[S].北京：中国铁道出版社, 2008.TB/T 2140―2008, Gravel ballast of railway [S].Beijing:China Railway Publishing House, 2008.(in Chinese)
[22]Eliáš J.Simulation of railway ballast using crushable polyhedral particles [J].Powder Technology, 2014,264(264):458―465.
[23]TB 10413―2003, 铁路轨道工程施工质量验收标准[S].北京：中国铁道出版社, 2003.TB 10413―2003, Standard quality of railway track engineering [S].Beijing: China Railway Publishing House, 2003.(in Chinese)
[24]Zhang Z, Zhang X, Qiu H, et al.Dynamic characteristics of track-ballast-silty clay with irregular vibration levels generated by high-speed train based on DEM [J].Construction & Building Materials, 2016, 125(24):564―573.
[25]刘建坤, 肖军华, 杨献永.提速条件下粉土铁路路基动态稳定性研究[J].岩土力学, 2009, 30(2): 399―406.Liu Jiankun, Xiao Junhua, Yang Xianyong.Study on dynamic stability of silt subgrade subjected to train moving loading [J].Rock and Soil Mechanics, 2009,30(2): 399―406.(in Chinese)
[26]狄宏规, 冷伍明, 薛继连, 等.朔黄铁路重载扩能的路基强度评估[J].铁道学报, 2014, 36(8): 84―90.Di Honggui, Leng Wuming, Xue Jilian, et al.Assessment of subgrade strength for transport capacity enlargement of Shuo-Huang heavy-haul railway [J].Journal of the China Railway Society, 2014, 36(8): 84―90.(in Chinese)
[27]陈建国, 肖军华, 李前进, 等.提速列车荷载作用下铁路路基动力特性的研究[J].岩土力学, 2009, 30(7):1944―1950.Chen Jianguo, Xiao Junhua, Li Qianjin, et al.Dynamic responses of exiting-railway subgrade with train speed increasing [J].Rock and Soil Mechanics, 2009, 30(7):1944―1950.(in Chinese)
[28]TB 10001―2016, 铁路路基设计规范[S].北京：中国铁道出版社, 2016.TB 10001―2016, Code for design on subgrade of railway [S].Beijing: China Railway Publishing House,2016.(in Chinese)
XIAO Jun-hua , ZHANG De , WANG Yan-hai , GUO Jia-qi
(Key Laboratory of Road and Traffic Engineering of the Ministry of Education, Tongji University, Shanghai 201804, China)
Abstract: The discrete element method (DEM)and finite difference method (FDM)coupling algorithm were introduced into the analysis of interface stress between railway ballast bed and soil subgrade.The crushed ballast was simulated by the refined DEM modelling, combined with laser scanning and triaxial testing.Then, an interface element was established between interlayer structure boundaries of ballast bed and subgrade, and consequently the DEM and FDM were dynamically coupled through exchanging the data of forces and velocities between two layers within each time step.Thus, the coupling simulation model was established, which had been validated with the field measured results.On this basis, the characteristics of the dynamic stress distribution of subgrade under the load of one train passing were calculated by the coupling model, and the effect of tie-ballast bed contact state was investigated.The calculated results were compared with the ones derived from the finite element method model of tie-ballast bed-subgrade.The results show that: the force chain in ballast bed is discrete in the coupling model if tie-ballast bed not fully close contact, and on the contrary, the force chain is verticaldiffusion; the distributions of subgrade surface stress for the two contact conditions of tie-ballast bed are significantly different, but the stress peak is close for each other.Comparing the stress peak derived from a coupling model and a FEM model, the former is significantly higher than the latter.The influence depth of local stress concentration on subgrade surface caused by irregular ballast was about 8 times as the maximum particle size of ballast (i.e.63 mm), approximating the depth of a subgrade surface layer.
Key words: ballast bed; discrete element method; finite difference method; coupling algorithm; dynamic stress of subgrade
STUDY ON INTERFACE STRESS BETWEEN BALLAST AND SUBGRADE FOR TRADITIONAL RAILWAY BASED ON COUPLED DEM-FDM
郭佳奇(1993―)，男，山西人，硕士生，主要从事铁路路基动力学研究 (E-mail: class5126386@163.com).
王延海(1994―)，男，安徽人，硕士生，主要从事铁路路基动力学研究 (E-mail: 2539607905@qq.com);
作者简介：
张 德(1990―)，男，江苏人，博士生，主要从事铁路路基动力学研究 (E-mail: dz2015@tongji.edu.cn);
基金项目： 国家自然科学基金项目(51678447)；高速铁路线路工程教育部重点实验室开放研究基金项目(2015-HRE-03)
通讯作者： 肖军华(1980―)，男，湖北人，教授，博士，主要从事铁路路基动力学研究(E-mail: jhxiao@tongji.edu.cn).
文章编号： 1000-4750(2018)09-0170-10
收稿日期： 2017-06-05；修改日期：2017-11-15
文献标志码： A
doi: 10.6052/j.issn.1000-4750.2017.06.0422
中图分类号： U213.7+1
0.6＜S 
0＜S 
0＜S 
(b)各粒径道砟颗粒二维离散单元56 mm～63 mm 
45 mm～56 mm 
35.5 mm～45 mm 
16 mm～25 mm 
