高填路基分层填筑过程非线性有限元分析

李 宁 张海凤 周川滨

(1.中铁二院工程集团有限责任公司， 成都 610031;2.鲁南高速铁路有限公司， 济南 250101)

高填方路基工程在铁路建设中屡见不鲜，在高低起伏的地形条件下，采用高填路堤能有效满足线路的线形要求，同时消耗弃渣。高填方路基在计算分析中，应当充分考虑以下问题：首先，土体的变形模量受其应力状态变化的影响而变化，在高填方荷载作用下，将本体视为刚体或采用固定变形模量的线弹性模型计算填方的变形及稳定性显然已不适用[1]。其次，为有效减小工后沉降，高填方路基通常采用分层摊铺、碾压的施工工艺。在非线性计算分析条件下，将路基填筑加载过程等效为一次加载将无法正确的反映填筑层的应力和变形情况。为解决以上问题，更加真实、有效地对高填方路基的应力、变形进行分析计算，本文对邓肯-张(E-B)计算模型及其修正方法进行研究，以非线性有限元理论为基础，模拟了填方路基分层填筑的全过程。

1 修正邓肯-张(E-B)计算模型

土是一种复杂的多孔介质各向异性材料，作为路基填筑的粗粒料，在受到外界荷载作用的情况下，其应力与变形具有明显的非线性特性[2-3]。因此，对于高填方路基本体的应力、应变分析，所选本构关系能否正确模拟给定材料的应力、变形特性，决定了对路基本体结构变形特性预测的准确性和可靠性。

国内外有关研究工作表明[4-6]，邓肯-张(E-B)双曲线模型能够很好地反映填筑体的非线性变形特性，且与之配套的模型参数也有较为成熟的测定方法。该模型是Duncan[7]等人以σ1=σr=常量的三轴剪切试验为基础，将偏应力q=σ1-σ3和轴应变拟合为双曲线关系式，并在假定填料抗剪强度符合Mohr-Coulomb破坏准则的条件下，推导出切线模量Et和切线体积模量Bt用于弹性增量的分析。

填筑体当前应力水平S可按式(1)计算。它反映了材料强度的发挥程度，是判断填筑体内任意点是否处于极限平衡状态的一项重要指标[8]。

(1)

式中：(σ1-σ3)f——破坏时的偏应力。

将历史最大应力水平与土体当前应力水平S进行比较，可判断填筑体中任意点的加载与卸荷状态。

在逐级加载的工况下，切线弹性模量可表达为：

(2)

式中：k——切线模量基数；

pa——单位大气压力；

n——切线模量指数；

Rf——材料破坏比。

若土体处在卸载状态下，可根据三轴试验结果，采用相应的卸载回弹模量Eur替代Et。

(3)

Kur、nur为试验确定参数，与初始切线模量中的K、n确定方法相似。

切线体积模量为：

(4)

式中：Kb——体积模量系数；

m——体积模量指数。

此外，填筑模型同时还应考虑粗粒料内摩擦角φ随围压σ3的变化。

(5)

式中：φ0——σ3等于单位大气压力时的φ值；

Δφ——反映φ值随σ3降低的参数。

采用邓肯-张(E-B)进行非线性增量分析时，应力增量和应变增量的关系应适当变换。由弹性力学原理可知，B与μ、E之间存在下列关系：

(6)

由于邓肯-张(E-B)是基于σ2=σ3的常规三轴试验建立的，因此无法考虑σ2≠σ3的情况，即无法考虑中主应力σ2对强度以及变形的影响，这将造成计算分析的误差，使土体变形模量偏低，计算所得的土体变形偏大。因此，需要对该模型进行修正。Skermer提出将原邓肯-张(E-B)中的侧限压力σ3用(σ2+σ3)/2来代替，偏应力(σ1-σ3)用σ1-(σ2+σ3)/2来代替，Mohr-Coulomb准则不变[9]。本文即按此修正方法建立修正剪切模量和泊松比的计算模块。

2 路基分层填筑模拟方法

2.1 分层填筑荷载及应力场处理

路基填方一般采用分层填筑、分层碾压密实的施工方法。填方路基典型分层填筑施工过程如图1所示，(a)、(b)两种模拟方式均采用逐级加载。

图1 路基分层填筑、逐级加荷示意图

比较图1中(a)、(b)两种模拟方式[10]，在(a)模拟方式中，网格单元按填筑顺序依次激活，当级分层填筑完成后，土体的自重荷载由当级及前序已填筑部分路基本体来承担，尚未填筑部分的网格单元不承担该部分荷载。在(b)模拟方式中，网格单元在加载前一并激活，当级填筑产生的荷载由路基整体共同承担，因此，后续的分层填筑体依然受前序填筑土体自重荷载的影响，这与分层填筑过程中的实际情况显然不相符。

基于以上情况，本文采用了生死单元法来模拟分层加载过程。在分析计算中，对于尚未填筑网格单元的刚度矩阵乘以一个极小值，使其不对载荷向量生效，将其质量也设置为0。路基填筑体单元在初始状态下全部杀死，并在此状态下计算初始应力场。随后，各层网格单元在相应的分步计算中依次被激活。

2.2 分层填筑初始位移场处理

分层填筑过程中，当层填筑单元被激活后，网格单元具有零应力状态，但位移场仍受前序荷载步计算结果的影响，因此还应对其位移场进行修正，使其吻合实际施工过程中逐级加载的变形效果。路基填筑体第i层竣工时，第i+1层和第i+2层的应力结果为0，但位移结果并不为0。因此，本文处理位移场的思路是，在计算填筑第i层的位移场后，将i-1荷载步的计算结果作为初始位移场进行扣减，获得当层填筑荷载所产生的附加位移。全部荷载步计算完毕后，根据单元所处的填筑分级j，分别叠加第j步及后续荷载步的位移计算结果，合计后的位移场即为分层填筑过程中该单元的变形计算结果。

此外，在实际填筑过程中，当层的填筑厚度及变形还是会受到前序填筑顶面变形的影响。忽略土体固结等时间因素的影响，在当层填筑完成瞬时(本层填筑荷载产生的附加变形尚未发生)，填筑顶面位移应为0，底面变形应等于前序分级填筑顶面的变形。因此，本文在计算位移场叠加前序荷载步计算结果的过程中，引入了一个三角形修正系数，对其进行修正，进一步使计算结果贴近真实情况。

2.3 分层填筑计算流程

非线性有限元法可考虑路基填料的非线性影响，结合分层填筑施工过程，采用增量法分级计算，更加吻合工程实际情况，计算结果较准确。邓肯-张(E-B)模型流程如图2所示。

图2 邓肯-张(E-B)模型流程图

具体操作流程如下：

(1)建立高填方整体有限元模型，按照施工填筑方式对模型分层；将Duncan-Chang E-B模型参数赋予指定材料；杀死填方本体单元，给定初始边界条件计算初始应力场，并将结果作为分层填筑计算的初始应力条件和边界条件。

(2)填筑第一层分层模拟。读入初始应力场，复活第一层单元，由于在填筑逐级加载计算中，每级新填筑单元的初始应力为0，故将新填土层作为重液体处理，令σz=γh，σx=σy=K0γh(K0为土的静止土压力系数)，估算当层土体应力状态，并以此修正模型刚度矩阵。

(3)将当层材料修正得到初始材料参数后的模型进行迭代计算，每次计算后，利用计算结果再次修正材料参数，直至两次计算获得应力场差值小于之前设定的容差，则本层填筑计算完毕。

(4)填筑第二层分层模拟，通过上级计算获得本级的第一层应力场及边界条件，复活第二层填筑体单元，将第二层单元视为第一层，重复(2)、(3)步，进行新填筑部分和已有填筑本体共同作用的非线性有限元计算。

(5)重复(4)步操作，进行分层填筑施工过程模拟，直至分层填筑结束。

3 数值分析算例

3.1 数值分析模型

以填高60 m路基为例建立有限元计算模型，路基面宽20 m，边坡坡率1∶1.5。考虑计算域效应，计算范围取左右坡脚外125 m，竖直方向取填方基底向下50 m。三维有限元计算模型采用空间八节点等参单元进行离散。

采用邓肯-张(E-B)模型，并按本文方法进行修正。基底以下按基岩考虑，弹性模量取1.2 GPa，粗粒填料密度为 2 000 kg/m3，模型材料计算参数如表1所示。

表1 模型材料计算参数表

计算过程中只考虑填筑体自重。路基按照分层填筑的施工工序共分为5步，每层填筑高度相同。

3.2 计算结果分析

(1)填筑体位移量值整体随施工分级加载过程而增加，路基填筑完毕后的位移等值线如图3所示。从图3可以看出，完建工况下，位移、竖向位移以路基面中心整体呈对称分布，其中横向最大位移出现在左右侧路基边坡中部高程处，最大位移量值为28 cm，方向向外。竖向最大沉降值52 cm，约占填高的0.87%，出现在填筑体中心墙距基底面2/3填高处。路基面最大沉降46 cm，基底最大沉降10 cm。相比于一次性加载，本体最大沉降出现在路基面，采用分层模拟的计算办法，更接近于实际情况。

(2)填筑体大小主应力整体随施工分级加载过程而增加，填筑体大主应力分布接近自重应力。填筑体大、小主应力均呈对称分布，完建工况下，大、小主应力均随高程降低逐渐变大。大主应力最大值为 0.92 MPa，出现在路基填筑体基底中心处，小主应力最大值为0.67 MPa，同样出现在路基填筑体基底中心处，填方中未出现拉应力。小主应力等值线如图4所示。

图3 位移等值线图(cm)

图4 小主应力等值线图(MPa)

(3)路基填筑完毕后，应力水平等值线如图5所示。从图5中可以看出，填筑体内部应力水平整体呈对称分布，应力水平整体较低，集中在0.1～0.5之间。最大值分别出现在路堤边坡中部高程的表层部分和两侧填方的坡脚部分。其中，边坡浅表层应力水平达0.5，坡脚部位应力水平达0.6。可见，填筑过程向两侧的水平变形和刚性基底对坡脚的约束，使上述两个部位的应力水平有所增加，但在此坡率下，路基整体处于稳定状态。

图5 应力水平等值线图(MPa)

4 结论

(1)在非线性计算条件下，通过分步、分层加载的方式模拟路基分层填筑，较一次加载更加符合实际情况。

(2)本文提出的采用非线性有限元模拟高填方路基分层填筑施工过程的方法可有效模拟路基的分层填筑过程，充分考虑填筑过程中粗粒填料的应力路径变化，计算结果符合逐层加载的计算规律。

(3)计算所得填筑体完建工况下的应力状态和应力水平分布，可为后期进行稳定性分析、工后沉降计算等提供更加准确的基础数据。