基于CDEM的砂卵石地层盾构开挖面稳定性分析

薛亚东 ，张 森 ，李 兴 ，葛嘉诚

（1.同济大学岩土及地下工程教育部重点实验室，上海 200092；2.同济大学土木工程学院，上海 200092）

随着盾构法在我国隧道工程中的广泛应用，以成都为代表的砂卵石地层对盾构施工技术的适用性提出了挑战.在成都地铁1 号线盾构施工过程中，多个区间出现了明显的滞后沉降现象，造成了巨大的经济损失和不良的社会影响.砂卵石地层与均质地层相比具有结构变异性大、力学性质复杂等特点[1]，因此深入研究砂卵石的细观力学特性，探明砂卵石地层中盾构开挖面的失稳机理，对于我国城市地铁隧道建设具有重要意义.

针对砂土、砂卵石土的力学特性及其地层条件下的盾构施工等问题，国内外学者开展了广泛的研究.模型试验方面：范祚文等系统开展了砂卵石地层中盾构开挖对地层影响的试验工作，分析得出了开挖面支护压力分布与地表沉降规律的内在联系[2]；Chen Renpeng 等对砂土中盾构开挖面失稳机理进行了研究，发现开挖面失稳分为局部失稳和整体失稳，并通过监测应力集中比得到了土拱的高度[3]；胡雄玉等通过缩尺试验和PFC3D (particle flow code in 3D)数值模拟，对砂卵石地层和砂土地层中盾构开挖引起的地层变形做了对比研究，发现砂卵石地层地表沉降曲面自上而下呈现逐渐收缩的“V形河谷”状，而砂土地层地表沉降曲面呈“圆形漏斗”状[4].数值模拟方面：考虑砂卵石地层的结构特征与变形不连续性，Wang Yunnian 采用CDEM （continuum-based discrete element method）离散元方法，建立砂卵石混合土体数值模型并进行单轴加载试验，数值模拟的结果与原位测试在非线性变形和破坏模式方面具有很好的一致性[5]，为砂卵石地层的数值模拟提供了一种新手段；白永学等基于现场调查，结合FLAC3D （fast lagrangian analysis of continua 3D）数值模拟，发现盾构开挖面支护压力不足导致地层局部区域失稳滑动，是形成地层空洞的主要原因[6]；江英超等采用PFC 颗粒流离散元对砂卵石地层中滞后沉降的形成机理开展了细观层面的研究，得出滞后沉降的形成过程主要包括：地层扰动形成松散区域——松散区域贯通形成空洞——空洞失稳——向地表扩展——地表塌陷[7].

本文以成都地铁隧道盾构施工为工程背景，采用CDEM 方法开展砂卵石土三轴压缩试验的数值模拟计算，定量研究不同围压、含石率对其力学特性的影响.基于“颗粒流动”、“土拱效应”、“超挖出土”等盾构施工特点，建立模拟砂卵石地层盾构开挖面超挖出土失稳的二维动态离散元模型，研究开挖面失稳的渐进演化机制.

1 CDEM 原理及方法

CDEM 是由中国科学院力学研究所开发的基于连续介质力学的可变形块体离散元方法，将有限元与离散元进行耦合计算，通过块体边界的断裂实现材料从连续到非连续的渐进破坏，因此块体及界面是CDEM 的两个基本要素[8].块体由若干有限元单元组成，可被设置为弹塑性体，在其内部根据力边界条件基于有限元法求解单元节点位移，从而表征材料的连续变形特征；接触界面采用基于结构层缩进技术的接触模型，引入可断裂的法向和切向弹簧传递接触力，界面受压时接触，受拉时分离，界面切向依据摩尔-库伦强度准则判断块体单元状态；界面破坏前与有限元方法求解弹塑性问题一致，当界面满足相应的破裂准则时，接触界面发生滑移，而块体单元仍只发生变形，从而实现破裂、大位移或大转动的模拟.

2 砂卵石三轴压缩数值模拟

2.1 砂卵石的物理力学特性

砂卵石土是一种典型的土石混合体，表现出显著的二元结构性.已有的研究表明：薄弱接触界面增加了材料的不连续性[9]；含石率是影响混合体力学特性的重要参数[10]；粗颗粒间依靠点对点接触、摩擦、咬合等提高了地层的抗剪强度，在非扰动情况下具有较好的自稳性，而扰动状态下接触面易失效滑移而导致失稳[11-12].基于砂卵石土的物理力学特性，通常将相关力学问题划分为两类：以应力主导的小变形问题（如三轴压缩试验），宜采用连续或准连续介质方法模拟；以颗粒流动主导的大变形问题（如盾构开挖问题），宜采用非连续介质方法模拟.CDEM不仅适用于上述两类问题的精确模拟，还具有较高的计算效率.

2.2 试验模型及参数

砂卵石土在三轴压缩条件下属小变形范畴，故采用弹塑性本构的连续介质模型描述；卵石在三轴压缩条件下一般不发生屈服，故采用线弹性本构的连续介质模型描述；采用脆断弹簧模型描述土石接触界面，形成准连续介质模型.

基于小变形连续介质假设，在满足应变梯度精度要求的条件下，砂土单元的尺寸一般是任意的.成都地铁现场勘查资料显示，“卵石粒径以30～70 mm为主，局部为80～120 mm”.结合《土工试验规程》[13]，确定本试验模型为直径（B）300 mm、高（H）600 mm的圆柱体，采用四面体单元进行自由网格划分.划分控制尺寸与卵石粒径相匹配，分别划分0.02、0.03、0.04 m 3 种尺寸的网格，经单元组集后卵石块体平均粒径分别为0.04、0.06、0.08 m.对基本网格随机抛撒指定大小、密度的虚拟球体域，当任意四面体单元节点处于随机球体域内，则该单元被指定为卵石属性，其余单元则被指定为砂土属性.将土石单元接触界面离散化，建立界面接触模型.试验采用位移控制式加载，在试样上表面设置刚性板并设定合适的加载速度，下表面设置固定位移边界.圆周方向采用径向压力边界模拟围压，分别设置100、200、300、400 kPa 4 种围压.根据现场地勘调查报告并参考《工程地质手册（第4 版）》[14]，设置单元参数见表1，数值计算模型如图1所示.

表1 单元力学参数Tab.1 Mechanical parameters of elements

图1 数值计算模型Fig.1 Numerical model

2.3 试验结果及影响因素分析

2.3.1 典型应力-应变曲线

以0.03 m 单元尺寸，40%体积含石率（WV），200 kPa围压的试验结果为例，绘制偏应力关于轴应变关系曲线，如图2所示.其受力变形过程可分为3 个阶段：（1）线弹性阶段，偏应力与轴应变呈线性关系；（2）弹塑性阶段，曲线出现拐点，偏应力随轴应变增长的速率大幅下降，但仍呈上升趋势；（3）理想塑性阶段，偏应力基本不随轴应变增加而变化.该结果与滕丽等得到的成都地区砂卵石地层真三轴试验典型应力应变曲线基本一致[15].

为研究卵石对砂土的骨架增强作用以及土石接触界面对砂卵石宏观强度的弱化作用，在相同围压和卵石粒径条件下分别设置3 组试样进行试验：均质砂土、砂卵石土（考虑界面接触）、砂卵石土（无界面接触），试验结果见表2.砂卵石土（考虑界面接触）相较于均质砂土峰值应力和初始切线模量显著提高，弹性极限应力稍有提高，表明卵石对砂土具有一定的骨架增强作用；砂卵石土（考虑界面接触）较砂卵石土（无界面接触）峰值应力显著降低，是由于土石界面发生滑移导致两种材料无法在界面上充分传递荷载；砂卵石土（考虑界面接触）最终进入理想塑性阶段，而砂卵石土（无界面接触）持续处于弹塑性阶段，但弹性阶段两者应力应变路径几乎相同.

图2 典型偏应力-轴应变关系曲线Fig.2 Typical deviation of stress-strain curve

表2 不同模型试验结果Tab.2 Results of different models

2.3.2 围压、含石率的影响

取0.03 m 单元尺寸，40%体积含石率，分别在不同围压（σ3）下对试样进行固结加载，得到各围压下偏应力-轴应变关系曲线，如图3所示.随着围压增加，试样峰值强度增加，线弹性阶段延长，但围压并不影响试样进入理想塑性阶段的时机.

在0.03 m 单元尺寸条件下分别设置含石率为40%～60%、20%～40%、0～20%的试样，在不同围压下进行试验，由试验结果（如图4所示）可知：在一定范围内随着含石率增加，砂卵石的初始切线模量和峰值应力增大，但当含石率增加到一定程度后该作用减弱；尽管更多的卵石单元增加了材料的宏观强度，然而更多的接触界面使砂土单元呈现“支离破碎”的形态，难以有效传递应力；接触界面上能够承受的最大剪应力遵循摩尔-库伦准则，当实际剪应力超过最大剪应力后界面破坏，因此过多的接触界面会削弱砂卵石整体承载能力.

图3 不同围压下偏应力-轴应变关系曲线Fig.3 Deviation curves of stress versus strain under different confining pressures

图4 不同含石率条件下偏应力-轴应变关系Fig.4 Deviation curves of stress versus strain for different cobble content rates

3 砂卵石地层盾构开挖面模型

3.1 平面三角块体离散元摩擦角标定

由于平面三角形单元与砂土颗粒几何相似，并且在细观层面上能够较好地表现转动状态，故选择其作为基本离散单元.当考虑计算经济性导致模型单元尺寸远大于实际颗粒尺寸时，尺寸效应会对离散系统的流动性与结构性带来较大影响.综合考虑后拟取0.08 m 为基本单元尺寸，该尺寸与真实卵石尺寸接近，但与真实砂土尺寸差异悬殊，因此需考虑数值计算中砂土的尺寸效应.

颗粒体的内摩擦角反映其摩擦特性，一般认为它由两个部分构成：接触摩擦作用和角粒效应.当离散单元尺寸较小时角粒效应较弱，可认为宏观内摩擦角与接触界面摩擦角近似相等；但当离散单元尺寸较大时，角粒效应对内摩擦角的贡献不可忽略，因此需通过一定的相似条件进行等效内摩擦角的标定.基于实际的料斗与散粒体物料的相互作用，设计Hopper Flow 料斗试验来标定大尺寸单元的“等效内摩擦角”.设置一个倒梯形料斗，料斗侧壁保持固定，在料斗底部固定刚性挡板.使散体物料重力固结，待系统稳定后，以一定速度释放底部挡板，散体物料由下往上开始发生塌落流动，监测对称轴上5 个测点的位移、水平应力的动态变化情况.各监测点均处于对称轴上，所测水平应力即主固结应力，反映了不同高度处介质成拱效应的程度，可作为标定界面摩擦角的标准.假设真实砂土的离散单元尺寸为0.002 m，已测得成都地区砂卵石地层中砂土的内摩擦角为35°，将此作为基准单元.以0.08 m 尺寸的砂土单元作为目标单元.通过引入尺寸为0.012 m的中间单元减少计算量，经两级Hopper Flow 标定试验得到目标单元的等效内摩擦角为11°，用于盾构开挖面工程尺度的模拟计算.

3.2 盾构开挖面模型

建立图5所示的盾构开挖面二维超挖失稳模型，包含有限范围的地层、已支护的部分混凝土衬砌、开挖面挡板、底部的螺旋出土器.其中螺旋出土器的出土构造作了相应简化：在该处设置可动出土挡板，通过其水平位移和土体单元置零模拟出土过程.按照上述模型，选取具有一定颗粒级配的组合卵石模型（颗粒半径r= 0.046～0.080 m），建立浅埋盾构开挖面的二维网格模型.为研究不同体积含石率下地层渐进破坏过程的异同，建立0、20%、40%、60% 4 种含石率的模型.模型基本参数、材料基本参数及界面接触单元基本参数分别见表3、4、5（表5中砂土-砂土的摩擦角基于3.1 节的标定结果）.

图5 盾构开挖面超挖失稳模型Fig.5 Failure model of tunnel face

表3 模型基本参数Tab.3 Basic parameters of the model

表4 材料基本参数Tab.4 Basic parameters of the materials

表5 界面接触单元基本参数Tab.5 Basic parameters of the jointed elements

3.3 模拟过程

模拟过程主要包括以下阶段：（1）对模型外部边界、衬砌、各挡板分别施加位移约束；（2）在重力作用下迭代计算至静力稳定；（3）出土挡板向左侧移动1.2 m，开挖面土体塌落出土；（4）地层迭代计算至稳定；（5）通过特定区域（出土口至左侧挡板之间0.5 m × 1.0 m 的矩形区域）单元置零挖土；（6）地层重新迭代计算至稳定；（7）重复（5）、（6），共挖土25 次，最终地层迭代计算至稳定，模拟过程结束.

4 盾构开挖面超挖失稳分析

4.1 开挖面超挖土体渐进破坏

盾构开挖面超挖引起的破坏是渐进的，主要是由于超挖量较小时砂卵石土体内部存在显著的土拱效应，能够在一定范围内通过颗粒接触与位移调整形成承载土拱来保证上部地层与开挖面的稳定.随着超挖量增大，开挖面附近土体减少导致地层内局部区域产生临空面，砂卵石介质随之产生滑移、滚动，空洞区开始形成并逐渐向上方与前方发展，最终可能导致开挖面塌落失稳.研究不同含石率情况下土体超挖渐进破坏过程，并重点分析对开挖面稳定有利的“土拱效应”的发展和对开挖面稳定不利的“空洞区”的演变规律.

4.2 土拱效应发展规律

根据土拱效应理论，当土体形成较稳定拱应力条件时，主应力方向会发生偏转，拱轴线即为主应力迹线.因此，最小主应力云图能够直观地显示土拱效应的发展情况.以40%体积含石率模型为例，比较3 个开挖阶段最小主应力云图的变化（如图6所示），其清晰地显示出地层内部土拱效应随着超挖出土的发展进程：承载土拱形成于出土口底部——向上发展巩固——上部拱角发展至开挖面顶端时进入极限状态——向地表发展并趋于弱化——达到地表时破坏消散.

随着含石率降低土体成拱能力弱化，承载土拱出现时机推迟，失效时机提前，拱顶的极限高度降低；高含石率地层中土拱效应趋于垂直向上发展，低含石率地层中趋于斜向上发展.

4.3 空洞区演变规律

传统的颗粒离散元采用接触力线反映颗粒堆体内部的渐进破坏过程，但其存在一定的主观性，CDEM方法中界面弹簧的失效状态即反映土体内部滑移错动及空洞情况.在松散介质地层扰动问题与稳定性研究中，基于界面接触状态分析的CDEM 方法能够从物理角度本质地解释地层内部的变形状态及松散或空洞区的发展规律.

图6 不同开挖阶段最小主应力云图Fig.6 Minimum principal stress contours in different excavation stages

以40%体积含石率模型为例，各开挖阶段单元接触状态云图如图7所示.超挖出土在细观层面上会带来两类不同的扰动影响：其一是界面受剪破坏，表示界面发生接触滑移，随着开挖出土其产生频率逐渐增加，总体分布范围较广，以红色表示；其二是界面受拉破坏，表示界面发生接触分离，随着开挖出土逐步从出土口斜向上渐进扩展，集中分布在开挖面附近，以黄色表示.这两种扰动对地层稳定性的影响有明显区别：界面受剪破坏仅表示原有接触关系发生更新，仅可能引起地表沉降；而界面受拉破坏是一种不稳定状态，在偶然因素的作用下易引起地层失稳，蓝色表示界面接触完好.

图7 不同开挖阶段单元接触状态云图Fig.7 Contact condition contours in different excavation stages

根据模拟结果总结出空洞区的一般发展规律为：原始低位土拱破坏——形成分散的空洞区域——土拱逐渐上移——空洞区域逐渐向上扩展（一般处于新土拱的下部）——土拱发展到开挖面顶端并进入极限状态——空洞区继续扩展直到穿透极限土拱——地层塌落，地表出现“水滴状”空洞区.比较不同含石率情况的模拟结果后发现：高含石率地层空洞区较窄且向地表垂直方向发展，低含石率地层空洞区稍宽且斜向上发展，这与胡雄玉等[4]结论一致.

4.4 地层变形分析

以 40% 体积含石率模型为例，在开挖面至开挖面前方 9 m 范围内以 0.25 m 间隔均布 37 个监测点，汇总各监测点监测到的地表沉降随超挖量（G）的时程变化，由于对数横坐标对于曲线的拐点具有更好的显示作用，将其时程变化在对数横坐标下绘出，如图8所示（x为监测点离开挖面距离）.由此，地表各点沉降速率大致可分为 3 个阶段：线性沉降阶段、缓慢沉降阶段、加速沉降阶段.图8中黑色纵线为 3 个沉降阶段的分界线.各阶段分界点与土拱发展有一定的内在联系：随着超挖出土承载土拱上升，覆盖面积扩大，减弱了地表沉降对于超挖的敏感性，进入缓慢沉降阶段；随着超挖继续，空洞区逐渐击穿极限土拱，土拱效应破坏，进入加速沉降阶段，即代表地层失稳.

进一步分析可知，加速沉降阶段中部分测点的沉降速率明显高于其它测点，这说明各区域对于超挖失稳的敏感性有一定的差异，可分为敏感失稳区与一般失稳区，这与实际工程中不同区域受到盾构开挖影响不同的事实相符.图8（a）为敏感失稳区，其测点均位于开挖面前方 3.25 m 范围以内，曲线在纵线处出现了显著的拐点，代表该测点下方土体受到土拱效应的影响；图8（b）为一般失稳区，其测点均位于开挖面前方 3.5 m 以外，总体基本呈线性沉降.随着含石率降低，缓慢沉降阶段的开始与结束分别延后和前，即土拱效应维持开挖面稳定的效果减弱；敏感失稳区向远离开挖面方向偏移；纯砂模型整体接近线性沉降，是因为砂土地层中的土拱效应远弱于砂卵石地层.

图8 地表沉降-超挖量时程曲线Fig.8 Curves of ground settlements versus over-excavation volumes

4.5 超挖量控制标准

开挖面土体失稳是由于空洞区穿越极限承载土拱造成的，据此从3 种角度提出3 类超挖量控制标准（严格程度递增）：（1）从稳定开挖面土体的角度：超挖量不破坏极限承载土拱；（2）从施工后治理沉陷的角度：超挖量不使空洞区发展高度超过开挖面顶端；（3）从控制地表沉降的角度：超挖量不使地表沉降超过环境变形要求.

图9汇总了60%、40%、20% 3 种含石率开挖面模型的超挖次数与超挖量的关系（砂土由于极限承载土拱作用不显著，不将其列入）.分别以极限承载土拱破坏、空洞区发展至开挖面顶端，开挖面顶端承载土拱形成作为关键超挖节点的判断准则，各关键超挖节点对应的超挖量分别为第1 超挖控制量、第2 超挖控制量和第3 超挖控制量.根据模拟结果，高含石率地层（体积含石率大于40%）第1 超挖控制量为5.5 m3，第2 超挖控制量为4.5 m3；低含石率地层（体积含石率小于20%）由于超挖出土现象不易发生，建议采用开挖面极限支护压力模型确定实际支护压力控制标准.

图9 关键超挖节点汇总Fig.9 Significant over-excavation stages

5 结 论

（1）首次将砂卵石土置于三维细观尺度，采用考虑土石界面的准连续介质模型，利用CDEM 模拟了砂卵石土的三轴压缩试验，研究其细观渐进破坏规律.

（2）基于土拱效应开展Hopper Flow 标定试验，获得了用于工程尺度计算的砂土模型单元的等效内摩擦角.基于“颗粒流动”等盾构施工特点，建立模拟砂卵石地层盾构开挖面超挖出土失稳的二维动态离散元模型，得到了土拱效应、空洞区发展等开挖面失稳的渐进演化机制.

（3）地表沉降随超挖量的时程曲线可分为3 个阶段.开挖面土体的最终失稳破坏是由于空洞区范围穿越极限承载土拱造成的，据此提出3 类超挖量控制标准.

（4）本文未考虑水的影响，后续研究可考虑地下水及其渗流对于砂卵石土的力学特性与砂卵石地层盾构施工的影响.