基于PFC2D的高填黄土减载明洞土体固结蠕变分析

后浩斌，李盛，2*，尤著刚，何永泽

(1.兰州交通大学土木工程学院，兰州 730070；2.兰州交通大学道桥工程灾害防治技术国家地方联合工程实验室，兰州 730070)

高填黄土明洞因能较好地适应黄土高原地区山高谷深、沟壑纵横的地形限制，被广泛应用于西北地区工程建设中。然而，由于其采用的填料多为黄土，而黄土具有明显的固结蠕变特性，这将导致土体内部的应力重分布，并对减载材料的减载效果产生影响。同时，若减载材料自身也具有固结蠕变特性，其减载效果的长期有效性也将无法保证。因此，探究黄土固结蠕变对高填减载明洞受力的影响十分必要。

有关固结蠕变的研究多集中于探究其发生的微观机理及对黄土物理指标变化的敏感性方面。朱才辉等[1]通过对压实Q3马兰黄土进行不同压实度和不同含水率的蠕变试验，证明了黄土固结蠕变特性对含水率的敏感度较大，对压实度次之。葛苗苗等[2-3]和张豫川等[4]分别对压实黄土进行固结蠕变试验，验证了压实黄土的长期变形特性，分析了黄土固结蠕变对含水率、压实度、固结应力等单一影响因素变化的响应。周远强等[5]对黄土固结蠕变的快慢与填土速率、填料压实度、强夯深度的相关性进行研究，发现黄土的固结蠕变随填土速率的加快呈线性增加，随填料压实度、强夯深度的增加呈线性衰减。Kang等[6]和Wang等[7]分别开展数值模拟和固结蠕变试验，研究了不同土体在固结蠕变过程中的宏微观响应。张凯旋[8]从细观组构演化的角度对黄土的固结蠕变特性及机理进行了深入分析。与此同时，众多学者在研究高填方结构时也开始关注土体固结蠕变对结构所受土压力的影响。Mcaffee等[9]采用现场试验、数值模型等方式研究了管涵填方试验中，路堤内部土压力及沉降随土体固结蠕变的变化规律。文献[10-12]通过现场试验研究了高填黄土减载明洞土压力及结构内力的分布规律，并推导了高填明洞土压力计算公式；姚玉相等[13]利用数值模拟确定了明洞采用聚苯乙烯泡沫板(expanded polystyrene boar，EPS)板减载时的最优参数取值。赵亮亮等[14]分析了槽宽比、边坡坡角等的改变对高填减载明洞土压力变化的影响。

然而，针对黄土固结蠕变特性对高填减载明洞长期土压力分布及对减载措施实际效果的影响鲜有报道。为此，以离散元软件PFC2D建立虚土减载的高填黄土明洞数值模型，探究了黄土固结蠕变对虚土实际减载效果的影响；分析了不同填方深度处土压力及竖向位移随黄土固结蠕变的发展规律；并通过土颗粒间接触力链的演化规律，从细观层面对黄土固结蠕变过程中土颗粒间的接触变化及虚土减载的荷载传递机理进行分析；最后，通过参数化分析进一步探究了土体固结蠕变条件下，虚土厚度、宽度和矢跨比的改变对明洞受力的影响，得到虚土减载时的最佳参数，以期为高填黄土减载明洞在设计、施工中考虑土体固结蠕变的影响以及采取合理的设计参数以优化明洞设计提供参考。

1 基于离散元法的数值模拟

1.1 模型的建立

所建立的高填黄土减载明洞模型如图1所示。模型总高H=41 m，底宽B=13.8 m。明洞高h=11 m，跨径b=12.8 m，边坡坡脚70°，虚土位于明洞顶，厚度T=5 m，宽度W=12.8 m，沿明洞中轴线对称设置；明洞顶以上覆土高度H1=30 m。模型边坡、明洞、地基均采用wall单元模拟，填土和虚土为黄土，用ball单元模拟，填土压实度设置为96%，虚土压实度为85%。为了使数值模拟的结果符合实际情况，在建立模型时，采用改进的分层压实法(improved multi-layer compaction method，IMCM)[15]将填土分7层填筑，以避免生成初始模型时产生过大的水平应力。为保证顶面土体被压实，留白区2 m。

图1 层次模型

为了监测伴随土体固结蠕变，明洞洞顶、明洞不同填方深度处的竖向土压力和竖向位移以及填土表面沉降的变化情况，分别在A-A、B-B、C-C、D-D截面布置了一定数量的编号测量圆。其中，A-A截面处1～9号测量圆监测明洞顶竖向土压力随土体固结蠕变的变化及洞顶平均竖向土压力大小；C-C截面处5、10、11号测量圆及D-D截面处12、13、14号测量圆分别监测不同填方深度处的土压力及竖向位移随土体固结蠕变的变化；B-B截面处水平布置的一排测量圆监测填土表面沉降随固结蠕变的变化。测量圆具体布置如图2所示。

红色矩形区域为设置虚土区域

1.2 本构模型的选择

选用PFC2D内置的Burgers本构模型来描述荷载作用下填土和虚土的固结蠕变行为，其本构方程为

(1)

式(1)中：EM为Maxwell体弹性模量，MPa；ηM为Maxwell体黏性系数，MPa·h；EK为Kelvin体弹性模量，MPa；ηK为Kelvin体黏性系数，MPa·h；EM、ηM、EK、ηK均为本构方程中表征填土和虚土固结蠕变行为的宏观参数；t为固结蠕变的时间，h；σ0为初始应力，kPa；ε(t)为初始应力σ0在t时间内产生的应变。

1.3 材料参数的获取

使用的相关概念说明如下：①“土体”指代“填土及虚土”；②以第7层填土填方完成时刻(S0阶段)作为土体固结蠕变前后的分界;③填方完成时刻以前为固结蠕变前阶段，未考虑土体固结蠕变；④填方完成时刻之后为固结蠕变阶段。

与有限元软件能够直接识别控制材料力学特性的宏观参数不同，PFC2D程序是通过识别与宏观参数对应的细观参数进行模拟分析的。因此，研究过程中需依次获取三项材料参数：①土体固结蠕变前的细观参数；②填方完成时刻之后，土体固结蠕变时Burgers本构模型中的宏观参数；③与Burgers模型中的宏观参数相对应的表征土体固结蠕变的细观参数。具体过程为：土体固结蠕变前的细观参数通过建立数值双轴压缩试验模型，利用参数标定的方法获取的。即通过不断调整控制颗粒整体宏观力学特征的细观参数值，直到颗粒集合表现出与实验室测得的填土相同的宏观性质时即可。最后，利用PFC2D程序内置的FISH语言将细观参数赋予图1模型并运行程序，可确定各层填土层及虚土层的初始应力σ0。标定土体细观参数所需的宏观参数(土体密度、压实度、弹性模量等)参考文献[4]，细观参数标定结果如表1所示。

表1 土体(填土和虚土)细观参数

考虑到实际明洞地基一般由钢筋混凝土浇筑而成，变形相对于土体要小得多，可将其视为刚性。根据文献[13]的经验，采用wall单元模拟明洞地基时，切向与法向刚度可取为1.0×1010N/s；而边坡一般经过加固处理，墙体刚度取为土颗粒刚度的10倍。故除土体外，未考虑其他材料的固结蠕变特性。

在确定Burgers本构模型的宏观参数时，首先对文献[4]中含水率ω=15.7%，压实度分别为k1=96%、k2=85% 的黄土在不同加载等级下进行一维固结压缩试验记录的应变—时间曲线进行拟合，得到一组宏观参数(EM、ηM、EK、ηK)作为基础；然后对初始宏观参数进行数值内插，并根据前述确定的各层填土和虚土层的初始应力σ0得到对应的初始宏观参数。

考虑土体固结蠕变的细观参数的标定，以填土为例详述如下：①将前述确定的各层填土的初始应力σ0及对应的宏观参数(EM、ηM、EK、ηK)代入Burgers本构方程[式(1)]，确定各层填土固结蠕变达到稳定时的理论应变值ε0；②建立蠕变参数标定数值模型标定填土固结蠕变时的细观参数，即：根据经验先假定一组细观参数作为参照赋予填土并运行PFC2D程序；运用“试错法”不断调整细观参数值，直到样本在不同初始应力σ0下的时间—应变曲线在经历初期的瞬时弹性变形阶段后平缓延伸，曲线的峰值应变ε1与理论值ε0基本相符，应变速率趋近于0，此时认为填土固结蠕变达到基本稳定。各层填土及虚土固结蠕变的细观参数如表2所示。部分土层蠕变细观参数的标定曲线如图3所示。

表2 土体(填土和虚土)固结蠕变的细观参数

图3 填土部分标定曲线

2 高填减载明洞土体固结蠕变结果分析

根据图2布置的测量圆，对土体处于不同固结蠕变阶段[Stage 0(S0)～Stage 10(S10)]时，明洞顶的竖向土压力、不同填方深度处的竖向土压力、竖向位移及填土表面沉降进行监测，结果如图4～图9所示。

图4 洞顶竖向土压力

2.1 洞顶竖向土压力分析

由图4可知，填方完成时刻，明洞洞顶(5号测量圆处)竖向土压力为390.4 kPa，相比竖向静止土压力(γH= 544 kPa，其中，γ为填土容重，H为填土表面至明洞顶的距离)减载效果明显，减载率为28.2%；虚土宽度范围内(1～9号测量圆)的平均竖向土压力为386.9 kPa，平均减载率为29.0%。这是由于明洞顶设置了压实度相对较小的虚土，其颗粒结构较为松散，压缩性较两侧填土更大，在上层土体荷载作用下，较容易产生竖向变形，从而使内外土柱之间产生沉降差，激发了土拱效应，致使明洞洞顶上方土体荷载转移向明洞两侧。而伴随着土体固结蠕变的发展，虚土颗粒间有效应力增加，宏观表现为虚土被逐渐压实，压缩性减小，内外土柱沉降差降低，土拱效应被削弱。此时，洞顶竖向土压力又回调增大至约530 kPa，相比竖向静止土压力(γH=544 kPa)仅相差2.57%，说明固结蠕变的发展明显削弱了虚土的减载效果。

2.2 不同填方深度处竖向土压力分析

图5为C-C截面和D-D截面不同填方深度处土压力随土体固结蠕变的变化，S0阶段代表填方完成时刻。由图5(b)可知，各层土体的初始土压力(S0阶段)均小于竖向静止土压力γH，也小于图6(a)中对应位置的初始土压力，且随着土体的固结蠕变呈现出“先上升(S0～S2阶段)、后回调(S2～S3阶段)、最终稳定(S3～S10阶段)”的特征。原因如下。

图5 截面竖向土压力

(1)填方完成时刻，C-C截面的虚土具有减载作用。但由于土体固结蠕变速率在填方完成后的初期(S0～S2阶段)发展较快，土体内部颗粒间有效压力迅速增大，虚土密实度增加，减载效果被削弱，致使明洞洞顶土压力显著上升。同时对比图1分析可知，土层位置越低，上覆填方高度越大，所受上覆土体竖向荷载更大，固结蠕变也更加显著，因而位置较低的土层较位置较高处的土层竖向土压力增加明显。

(2)土压力达到峰值后回调(S2～S3阶段)是由于虚土的残余减载效果被增大的土压力进一步激发导致的，但由于此时虚土的密实度增加，其残余减载效果有限，故而土压力的回调程度较小，且由于土层位置越高，其上覆填方高度越小，所受上覆土体竖向荷载相对较小，故回调趋势的变化也随土层位置由低处向高处逐渐递减；S3～S10阶段，土压力变化基本稳定，说明土体的固结蠕变在S3阶段以后趋于稳定。图5(b)表明，由于D-D截面处于虚土外侧，填土土压力未受到虚土残余减载效果的影响，且由于压实度较大，固结蠕变阶段(S3～S10阶段)各层土压力逐渐增加，至S3阶段以后也达到稳定。

2.3 不同填方深度处竖向位移分析

图6为C-C截面、D-D截面不同填方深度处竖向位移随土体固结蠕变的变化规律。对比图5可知，各层土体竖向位移的变化规律与土压力的变化相对应，C-C截面竖向位移随土体的固结蠕变同样呈现出“先上升(S0-S2阶段)、后回调(S2～S3阶段)、最终稳定(S3～S10阶段)”的特征；D-D截面竖向位移则持续增加，且各层土体竖向位移在S3阶段以后基本稳定。同时可以看出，在固结蠕变初期，由于虚土压实度较低，压缩变形明显，故C-C截面各层的竖向位移增量相较D-D截面对应各层较大。而随着土压力的回调[图5(a)中S2～S3阶段]，土颗粒之间咬合嵌套作用减弱，颗粒体系无法保持原有状态下的平衡而自发的调整相对位置，产生相互滑移、错动。故而C-C截面处土体各层竖向位移也在S2～S3阶段出现小幅回调。

2.4 填土表面沉降分析

填土表面沉降是表层填土竖向位移和低处土层竖向位移共同累积的结果。图7为图2中B-B截面布设的位移测量圆监测的土体(填土和虚土)固结蠕变前后填土表面沉降变化的对比情况。可以看出，填方完成时刻，填土表面沉降量较小且分布均匀，平均沉降量仅为35.4 mm；当固结蠕变基本稳定时，填土表面沉降量相较于填方完成时显著增加，且由于虚土的存在，洞顶宽度范围正上方的填土表面累积沉降量最大(830 mm)，最小沉降量(600 mm)则分布在两侧，总体呈现出“中间大，两头小的特征”。由固结蠕变引起的沉降量分别占累积沉降量的95.7%和94.1%。填土表面沉降量随固结蠕变的持续增大将对土体内部土压力环境的稳定性产生扰动，并对明洞衬砌结构产生危害。

图7 填土表面沉降

2.5 接触力链分析

在PFC2D中，接触力链(contact force chain)的分布实质上反映了土体中荷载的传递，力链的细密程度代表土颗粒间接触力的大小，线条组成的网络结构代表接触力的分布集度和荷载传递方向，网络结构越密集说明颗粒接触越紧密。

图8(a)～图8(c)为颗粒间接触力链伴随土体固结蠕变的细观演化过程。分析图8(a)可知：填方完成时刻(S0阶段)，由于虚土的压实度较小，结构松散，土颗粒相互之间接触较弱，故而力链网络稀疏且与周围填土的力链网络界限分明；其与两侧填土的刚度差异激发了土拱效应，从而使虚土上方的土体荷载逐渐向两侧转移并向下传递，降低了洞顶承受的土压力。虚土上方土体颗粒间的接触力链呈“拱形”，其在虚土两侧分布密集、稠密而在虚土宽度范围内分布稀疏也从细观层面与这种荷载传递规律相对应。

图8 颗粒间接触力链分布

图8(b)表明：随着土体固结蠕变的发展(S0～S3阶段)，虚土与两侧填土颗粒间的接触显著增强，力链线条由图8(a)中的“细而稀疏”变为“粗且密集”，并向两侧延伸，分布于虚土上方的“土拱”形态逐渐消失，说明土体的固结蠕变使虚土与两侧填土间的刚度差异降低，虚土颗粒被压密。图8(c)为土体固结蠕变趋于稳定时(S3阶段之后)的接触力链图。可以看出，此时虚土上方土体颗粒间力链的“土拱”形态已经消失，各层土体接触力链分布更加均匀。虚土范围内颗粒间接触力链更加紧密，颗粒间的接触力增大。说明土体的固结蠕变使虚土的减载效果大幅降低，由虚土激发的土拱效应被严重削弱。这与图4分析的减载明洞洞顶宏观土压力的变化规律相一致。

3 高填减载明洞土体固结蠕变影响因素分析

为方便设计者采取合理的设计参数优化明洞设计，以图1建立的数值模型(T=5.0 m；W=12.8 m；f/b=0.5，其中f为矢高)为参考，分别增设虚土厚度T= 3 m(1/10H1)、4 m(1/7.5H1)、6 m(1/5H1)；宽度W= 6.4 m(0.5b)、19.2 m(1.5b)；矢跨比分别为0.4、0.6、0.7，从洞顶平均竖向土压力(图2中A-A截面1～9号测量圆)、填土表面沉降的角度进一步探讨虚土的厚度T、长度W和明洞结构尺寸等参数变化对土体固结蠕变的影响。矢跨比f/b定义为明洞拱形衬砌的矢高f与跨径b之比。

3.1 虚土厚度的影响

图9为不同虚土厚度T减载时，洞顶平均竖向土压力和填土表面沉降的分布。分析图9(a)可知，填方完成时刻，洞顶平均竖向土压力随T的增大而减小，且差值显著。这是因为T越大，土颗粒间的空隙越多，压缩储备空间增大，激发的土拱效应越强，减载效果更好。随着土体固结蠕变趋于稳定，虚土被逐渐压密，土拱效应减弱，不同虚土厚度T减载下的竖向土压力均回调增大，洞顶平均竖向土压力也接近γH(544 kPa)，但相对大小仍与填方完成时刻一致，说明T越大，虚土经土体固结蠕变后的残余减载效果越大。

由图9(b)可知，填方完成时刻不同虚土厚度T减载时，填土表面沉降均较小且分布均匀，明洞顶正上方最大沉降量在39.43～41.84 mm。土体的固结蠕变使填土表面沉降大幅增加，且T越大，虚土压缩变形越大，各层回填土的沉降也随之增大，最终填土表面的沉降也随T的增大而增加。综上分析可知，高填黄土减载明洞考虑固结蠕变时，虚土厚度T越大，残余减载效果越佳。但在实际工程中，应综合考虑明洞顶以上填土高度、工程实际要求等选择适宜的虚土厚度。可得虚土厚度T=6.0 m(1/5H2)时最优。

图9 虚土厚度对下洞顶平均竖向土压力及填土表面沉降的影响

3.2 虚土宽度的影响

图10为不同虚土宽度W减载时，洞顶平均竖向土压力和填土表面沉降的分布。分析图10(a)可知，填方完成时刻，当W=6.4 m(0.5b)小于明洞宽度(b=12.8 m)时，通过虚土的减载作用向两侧转移的上覆土体荷载仍有部分作用于洞顶，相比W=12.8 m(1.0b)和W=19.2 m(1.5b)时上覆土体荷载完全转移向明洞宽度外侧，洞顶平均竖向土压力较大。此外，W=19.2 m(1.5b)时，由其激发的土拱被分布在更大的范围，减载范围扩大，反而降低了虚土的减载效率，故而洞顶平均竖向土压力随W的增大先减小，后小幅增大。随着土体固结蠕变趋于稳定，洞顶竖向土压力又回调增大，W=0.5b、1.0b、1.5b时对应洞顶平均竖向土压力较填方完成时刻均显著增大，分别为537.8、525.3、526.1 kPa，此时由于虚土的减载效果大幅减弱，土压力的差异主要受土体自重荷载的影响，即在洞宽范围内，当W=6.4 m(0.5b)时，平均竖向土压力是虚土和部分填土自重共同作用的结果，较W=12.8 m(1.0b)和W=19.2 m(1.5b)时荷载全部源于虚土更大。综上所述可知，填方完成时刻，虚土宽度W=12.8 m(1.0b)时减载效果较好；土体固结蠕变趋于稳定时，W=12.8 m(1.0b)的洞顶平均竖向土压力亦相对较小。

图10 虚土宽度对洞顶平均竖向土压力及填土表面沉降的影响

由图10(b)可知，虚土宽度W由6.4 m(0.5b)增大至19.2 m(1.5b)时，伴随着土体的固结蠕变，填土表面沉降较填方完成时刻均显著增大，但数值接近，明洞顶正上方对应最大沉降量分别为844.5 mm、833.9 mm、836.6 mm，平均沉降838.3 mm，最大沉降差7.9 mm，仅占平均沉降量的0.94%，W的增大并未提高其竖向可压缩性，而是使填土表面发生显著沉降的范围进一步扩大。

3.3 矢跨比的影响

图11为矢跨比f/b分别为0.4、0.5、0.6、0.7时洞顶平均竖向土压力和填土表面沉降的分布。分析图11(a)可知，填方完成时刻，不同矢跨比f/b下洞顶平均竖向土压力均在386.7 kPa。这表明虚土紧邻明洞顶布置时，其减载作用消解了f/b变化对洞顶竖向土压力的影响。随着土体固结蠕变，土拱效应被削弱，不同矢跨比f/b对应的洞顶平均竖向土压力相较填方完成时刻均回调增大，总体上随f/b的增大逐渐减小，这主要是因为当明洞顶以上填土高度H1不变时，f/b的增大导致上覆回填土高度降低，覆土荷载减小。由图11(b)可知，无论填方完成时刻或土体固结蠕变趋于稳定，不同矢跨比f/b时填土表面沉降曲线分布均具有较高的重合。综上所述，无论是否考虑回填土的固结蠕变作用，采用虚土减载的高填黄土明洞矢跨比f/b的变化对填土表面沉降几乎没有影响。

图11 矢跨比对洞顶竖向土压力及填土表面沉降的影响

4 结论

为探究黄土的固结蠕变对高填减载明洞长期土压力分布及虚土减载实际效果的影响，本文利用PFC2D建立数值模型，研究了伴随土体的固结蠕变，虚土的减载效果、不同填方深度处土压力及竖向位移的发展规律；通过土颗粒接触力链的演化规律，从细观层面对黄土固结蠕变过程中土颗粒间的接触变化、荷载传递机理进行分析；并从洞顶平均竖向土压力、填土表面沉降的角度探究了虚土厚度T、宽度W和明洞矢跨比f/b等的变化对考虑土体固结蠕变条件下明洞受力的影响。得出如下结论。

(1)高填黄土明洞采用虚土减载能够降低明洞洞顶竖向土压力。但当土体固结蠕变达到基本稳定时，虚土的减载能力会明显降低，洞顶土压力值将回调至接近静土压力γH。故在进行高填黄土减载明洞的承载能力设计时，应充分考虑由黄土固结蠕变导致的减载明洞洞顶土压力回调及其对明洞衬砌结构的影响。

(2)土颗粒间接触力链的演化规律从细观层面表明了虚土的减载效果随土体的固结蠕变被削弱，其上方的“土拱”逐渐消失是因为：土体的固结蠕变使虚土颗粒间的接触力链由疏变密，由细变粗，并与两侧填土颗粒间的接触增强，刚度差异降低，土拱效应明显减弱。当土体固结蠕变达到基本稳定时，土颗粒间的应力分布趋于均匀。土体的固结蠕变将使填土表面沉降持续增加，沉降增量将达累积总沉降的90%以上。

(3)明洞虚土厚度T越大，减载效果越好；土体固结蠕变后，其残余减载效果与T的大小呈正相关。填土表面沉降亦随T的增大而增加，最佳厚度为T=6.0 m(即1/5H1)；填方完成时刻，洞顶竖向平均土压力随虚土宽度W的增大先减小后增大，填土表面沉降均较小；土体固结蠕变后，不同W减载下的平均竖向土压力差值显著减小，填土表面沉降的范围随W的增大而增大，最佳宽度为W=12.8 m(即1.0b)；虚土的减载作用会抵消f/b变化对洞顶竖向土压力的影响；土体固结蠕变后，洞顶平均竖向土压力随f/b的增大而减小，f/b变化对填土表面沉降几乎没有影响。