许梦飞，姜谙男，史洪涛，李德生，万友生，程利民

(1.大连海事大学 道路与桥梁工程研究所，辽宁 大连 116026；2.中铁建大桥工程局集团第一工程有限公司，辽宁 大连 116021;3.南昌轨道交通集团有限公司，江西 南昌 330013)

1 概述

随着城市地下空间开发规模的逐渐增大，地铁隧道施工面临着更为复杂的周边环境。开挖造成的土体扰动必然会对既有建构筑物的安全稳定性产生威胁。如何对开挖扰动区进行预测，优化施工参数，减小对周边环境的扰动程度，是当前研究的热点问题。

在理论公式方面，吴昌胜[1]等基于Mindlin解，推导了大直径泥水盾构掘进引起的地表变形公式，分析了施工参数对地表沉降的影响规律；苟长飞[2]等针对盾构隧道建设过程中的同步注浆问题，推导了浆液充填压力沿盾尾间隙环向分布模型;魏刚[3]等基于随机介质理论推导了双圆盾构隧道施工引起的地表沉降公式；在盾构施工造成的对周边环境影响方面，童建军[4]等运用力学方法分析了下穿机场盾构施工过程中对上覆跑道变形的影响特征；史江伟[5]等、王正兴[6]等采用室内试验和数值模拟的方法研究了盾构隧道开挖对既有管线的影响规律。戴轩[7]等通过现场实测和有限元数值模拟的方法，分析了盾构施工过程中临近建筑物的变形规律。

综上，相关研究中，较少考虑到隧道开挖引起的应力-损伤-渗流耦合问题[8-9]。对此，本文首先建立了基于M-C准则的岩体弹塑性耦合损伤模型，同时引入损伤-渗透系数演化方程；其次，利用完全隐式的向后欧拉算法解决了积分过程中的飘零问题；通过在主应力空间对问题进行分区域讨论，克服了M-C准则的数值奇异点问题；最后，将该模型应用于南昌地铁三号线工程，分析了盾构施工对临近暗涵变形特征的影响规律，对比了施工过程中围岩位移和损伤值的变化趋势，给出了合理的施工建议。本文结论对类似工程的安全实施提供了一定理论依据。

2 岩石弹塑性损伤-渗流耦合损伤模型

2.1 基于M-C准则的岩石弹塑性损伤模型

考虑损伤的M-C准则表达式为：

(1-ω)ccosφ

(1)

式中：I1为应力张量第一不变量；J2为应力偏量第二不变量；θ为罗德角；φ为内摩擦角；c为黏聚力。与黏聚力c相比，损伤对内摩擦角的影响较小，可以忽略，由式(1)可以看出黏聚力随着损伤的累积逐渐减小。

相关研究表明，岩石内部损伤的演化通常伴随着塑性应变的累积，而幂指数形式又比较符合岩土材料的损伤演化规律，因此本文给出了损伤演化方程如下：

(2)

(3)

式中：εp1、εp2、εp3为3个主应力；κ为常数，可以通过室内试验或反分析求得；损伤变量演化规律如图1所示。

图1 不同κ值下损伤变量演化规律

此外，由Lemaitre有效应力和应变等效假设理论可知损伤材料的应力-应变关系为：

σ=(1-ω)Dε

(4)

式中：D为材料的弹性刚度矩阵。由式(4)可看出，当材料内部出现损伤时，其弹性模量也会发生弱化，损伤材料的应力应变曲线如图2所示。从图中可以看出，当岩石材料应变ε达到一定阈值ε0时，其内部产生损伤，致使其弹性模量由初始的E0变为E0(1-ω)。

图2 损伤作用下岩石应力-应变曲线

2.2 岩石应力-渗流耦合方程

岩石内部存在的孔隙、微缺陷等为地下水的赋存和发育提供了条件，由有效应力原理可知，渗流场产生的孔隙水压值p对岩石应力场造成影响：

(5)

岩石内部应力场改变的同时，其骨架结构也会发生变化，从而导致渗透系数的变化，根据Kozeny-Carman公式，非破裂区岩石渗透系数可以表示为体积应变的函数：

(6)

式中：K0为渗透率；n0为初始孔隙度；εv为体积应变。

对于破裂损伤区的岩体，其渗透系数会发生数量级上的突增，其可以表示为损伤变量的函数[10]：

K=K0·10ξ(A′e-D/α+B')

(7)

式中：ξ为跳突系数，用于描述损伤岩体渗透系数数量级上的突增；A'=1/(e-1/α-1)；B′=-1/(e-1/α-1)；α为经验参数。

3 完全隐式的向后欧拉算法

由于Mohr-Coulomb准则在屈服面上存在棱线、尖点等数值不连续特征，造成了应力积分的困难。本文为保证计算精度，摒弃了角点光滑法等近似解法，从主应力空间出发，建立M-C准则的完全隐式向后欧拉算法，避免了数值“奇异点”的问题,见图3。

图3 M-C准则在π平面上的图形

为了减少应力维数，简化计算过程，在主应力空间中对问题进行讨论。由于空间应力的对称性，只在σ1>σ2>σ3区域分析问题即可。

主应力空间中M-C准则的表达式为：

f1=σ1-σ3+(σ1+σ3)sinφ-2(1-ω)ccosφ

f2=σ2-σ3+(σ2+σ3)sinφ-2(1-ω)ccosφ

f3=σ2-σ1+(σ2+σ1)sinφ-2(1-ω)ccosφ

f4=σ3-σ1+(σ3+σ1)sinφ-2(1-ω)ccosφ

f5=σ3-σ2+(σ3+σ2)sinφ-2(1-ω)ccosφ

f6=σ1-σ2+(σ1+σ2)sinφ-2(1-ω)ccosφ

(8)

用膨胀角ψ代替式(8)中的内摩擦角φ可得到与屈服函数具有相同形式的塑性势函数gi。

整个应力计算过程分为弹性预测、塑性修正和损伤修正3个部分：

a.弹性预测：已知tn时刻的计算应力σn和当前迭代步的应变增量Δε，计算tn+1时刻的应力值为：

(9)

(10)

b.塑性修正阶段，映射应力可能存在的区域为应力平面f，棱线l1或l2和尖点P处。为判断应力所在区域，利用边界面法对应力区域进行划分，边界面方程可以定义为[15]：

(11)

式中：rp为塑性修正方向，其表达式为：

(12)

(13)

式中：p为向量长度，计算时取1。

σa的表达式为：

(14)

同理可以建立边界面：

(15)

建立边界面后，将试算应力带入，当pΙ-Ⅱ≥0且pΙ-Ⅲ≤0时，应力返回至屈服面；当pΙ-Ⅱ<0且pΙ-Ⅲ<0，应力返回至棱线l1;当pΙ-Ⅱ>0且pΙ-Ⅲ>0时，应力返回值棱线l2；否则应力返回至尖点。

判定返回区域后，根据塑性增量理论，分别建立N-R方程对更新应力进行求解：

(16)

式中：Δλ为塑性因子。

c.损伤修正。

利用求得的塑性应力求解塑性应变：

(17)

将求得的塑性应变代入式(3)和式(2)中计算损伤变量Dn+1。在新求得的损伤变量基础上，对应力进行再次修正：

(18)

基于ABAQUS软件中的UMAT子程序接口，利用FORTRAN语言实现编写了上述应力计算流程；每个应力步计算完成后，输出损伤值和体积应变，利用USDFLD子程序接口更新当前应力状态下积分点上的渗透系数，并以场变量的形式进行输出，最终实现流固耦合的完全计算。计算流程如图4所示。

图4 流固耦合计算流程图

4 工程应用

4.1 工程概况

南昌地铁3号线起点为莲塘站，终点为京东大道站。为联系莲塘、青云谱区、绳金塔、老城核心区、青山路、二七北路、青山湖风景区、火炬大街、高新区等重要区域的联系通道。线路全长28.5 km，全线设车站22座，均为地下站，平均站间距1.36 km。

南昌市气候为亚热带季风气候，气候比较温暖，降雨量多。年降雨量具备分派不均的特点，每一年4～6月份(梅雨季)为丰水期，降雨量较大，降水量占整年总量的51%；10月至翌年2月为少雨季节，为枯水期，降水量仅占总量的12.8%。因此，在市内进行地铁施工时，应做好雨季施工准备，做好防水措施。

本文所研究区域段CK34+660-885下穿玉带河暗涵，垂直净距10.5～12.3 m，该暗涵为混凝土框架结构。该段地质条件相对复杂，地层以上至下依次为杂填土、素填土、粉质黏土层、粗砂、砾砂、强风化泥质粉砂岩、中风化泥质粉砂岩、弱风化泥质粉砂岩。地层条件复杂，盾构隧道开挖对暗涵影响较大，因此需要在施工过程中对暗涵进行处理和隧道施工方法进行优化，见图5和图6。

图5 暗涵与盾构隧道相对位置图

图6 施工现场

4.2 有限元计算模型

根据工程实际情况建立三维有限元计算模型如图7所示。模型分为51 613个单元，73 706个节点。土体左右两侧施加水平方向约束，底部施加固定约束。

图7 有限元计算模型

按照地勘资料，地层由上而下分别为素填土1、素填土2、粉质黏土、粗砂、强风化泥质粉砂岩(占比较少暂不予考虑)、中风化泥质粉砂岩和弱风化泥质粉砂岩。除弱风化岩层外，其他岩土体采用理想弹塑性Morh-Coulomb准则，其参数选取如表1所示，弱风化岩层损伤参数取κ=15。在土体内部施加随深度呈梯形分布的孔隙水压，并设地表处孔隙水压为零。

盾构支护形式主要为管片和注浆层，按照施工设计，设置支护压力为0.23 MPa，注浆压力为0.3MPa，进尺为1.2 m。在有限元计算时，将盾构工法简化为以下步骤：移除(开挖)土体，激活当前环的盾构机模型，掌子面施加支护压力；进行下一步开挖时，移除开挖土体，激活当前环盾构机模型，移除上一环盾构机模型，同时在上一环添加注浆层和管片，并在对应位置处施加注浆压力，将上一环的支护压力撤去；如此反复直至开挖完毕。

表1 模型计算参数Table 1 Model calculation parameters材料名称重度/kN·m-3 弹性模量/MPa泊松比内摩擦角黏聚力/kPa渗透系数×10-9/m·s-1 杂填土117 3.5 0.2510527杂填土2173.50.2510527粉质黏土19.660.316204.7中沙19.5100.3280.001900中风化岩层23.47000.26321 15018弱风化岩层25.21 7000.3342 20014加固区(C20)2425 5000.25— — — 管片(C50)2534 5000.25— — — 注浆层(C20)2425 5000.25— — —

4.3 计算结果分析

4.3.1地表与暗涵变形模式分析

分别对未进行注浆加固，加固长度为2 m和加固长度为3 m共3种工况下的地表沉降值进行提取，结果如图8所示。

图8 不同加固条件下地表沉降曲线

由图8可以看出，未进行注浆加固时，地表最大沉降量达到14.08 mm左右；在衬砌环顶部中心线左右60°范围进行注浆加固，加固长度为2 m时，地表最大沉降量仅为9.58 mm左右；在衬砌环顶部中心线左右60°范围进行注浆加固，加固长度为3 m，地表沉降减少为8.08 m左右。随着距离隧道中心线距离越近，地表的沉降值越来越大，是因为距离中心线越近受隧道开挖的影响越大。

此外，随着加固长度的增加，地表沉降量的差值减少，不均匀沉降的影响减少。因此，当盾构机穿越暗涵时，对土体进行注浆加固处理，能够有效地增强土体力学参数，降低地表沉降，减小开挖扰动程度，保障施工的安全进行，增加隧道开挖的稳定性。由图8我们还可以看出，测线上距离侧线中点左右一定范围的测点地表沉降量相差不大，这是由于该范围内测点位于暗涵上方，受暗涵影响，因此沉降比较均匀，各测点沉降差值较小。

对未加固条件下，不同开挖距离下地表测线沉降值进行计算得到结果如图9所示。

图9 不同开挖距离对地表沉降的影响

由图9可以看出，随着隧道开挖距离中心测线越来越近，测线上地表沉降量逐步变大，并且增加趋势越来越明显。当隧道开挖至10 m左右时，测线位置处沉降量整体较为平缓，最大沉降量为3.9mm左右；当隧道开挖至20 m左右时，测线上各测点的沉降量变化较为明显，距离测线中点越近变化越明显，测线上最大沉降量为6.4 mm左右；当隧道开挖至30 mm左右时，此时隧道开挖至测线下方附近，测线上各个测点的变化越来越明显，测线上最大沉降量为9.0 mm左右。

为研究盾构隧道开挖对上部暗涵的影响，在该暗渠模型上选取6个测点，对不同加固条件下暗渠的变形情况进行计算，结果如图10所示。

图10 注浆加固对暗涵变形的影响

由图10可知，未进行注浆加固时，暗涵最大沉降值为14.08 mm，发生在测点2处；进行注浆加固后，加固长度为2 m时，暗涵最大沉降值减小为9.61 mm，加固长度为3 m时，暗涵的最大沉降量为7.98 m，同样发生在测点2处。整体注浆加固对开挖过程中，暗涵的稳定性起到了重要作用，且加固区的长度越大，加固对暗涵的稳定作用更加明显。

另外从图5还可以看出，位于暗涵外部测点最终沉降值相对于暗涵中部测点沉降值较小，暗涵底部上测点(测点5)沉降值明显小于周边测点沉降值。4.3.2不同开挖步下孔隙水压变化情况

为了研究隧道开挖过程中地层中渗流的变化情况，对未进行注浆加固条件下，盾构隧道开挖过程中孔隙水压分布规律进行分析如图11所示。

(a) 开挖0 m

由图11可以看出，在隧道开挖初期，地层中的孔隙水压力变化非常明显，渗流场成“漏斗状”分布，这是由于隧道开挖初期，土体突然受到扰动，进而导致地层中的渗流过程发生了改变。随着开挖的进行，由于管片和注浆层的止水作用，孔隙水压分布规律逐渐趋于为开挖时的分布规律，各个地层中的渗流过程基本趋于稳定，渗流场逐渐变化为原来地层的分布状态。

由此可以看出，隧道开挖初期对地层中孔隙水压的影响最大，对渗流过程的改变最大，隧道开挖最不稳定。因此，在隧道开挖初期应密切关注洞内涌水情况，避免发生灾害。

4.3.3开挖过程中围岩损伤值变化规律

对未进行加固，不同孔隙水压下，洞周关键点的位移与损伤值变化情况进行计算，结果如图12所示。

(a) 拱顶点

由图12可以看出，随着孔隙水压的增大，洞周监测点的位移与损伤值均随之增大，对于拱顶点，不同孔隙水压下的最大沉降值分别为17.3、18.3和20.4 mm;对于拱底点，最大隆起值分别为16.65、17.2和19 mm；对于拱腰点，最大收敛值分别为3.03，14.5和17.45 mm。相对而言，拱腰收敛值对孔隙水压变化较为敏感。损伤值与位移变化规律趋于一致。由计算结果可知，当雨季来临时，应加强监测手段，防止出现工程事故。计算结果也表明了本文模型具有一定的工程适用性。

4.3.4耦合模型验证分析

为进一步验证本文所建耦合模型的准确性和工程适用性，分别对不考虑渗流、不完全耦合(渗透系数不发生变化)和完全耦合[渗透系数按式(6)和式(7)进行演化] 3种条件下地表沉降值以与拱顶沉降值的计算结果进行对比分析，结果如图13和图14所示(根据现场实际情况， 拱顶沉降值取孔隙水压为0.2 MPa时的计算结果， 为便于与实测数据对比， 只选取了安装监测点后， 即开挖距离30 m以后的数值)。

图13 不同条件下地表沉降计算值

图14 不同条件下拱顶沉降计算值

由图13、图14可以看出，当不考虑渗流作用时，地表沉降值和拱顶沉降值要远小于考虑渗流条件下的数值。由于本文所建立的完全耦合模型考虑了渗透系数的演化，随着开挖的进行，渗透系数会不断增大，渗流作用愈发明显，造成位移的变化量也较不完全耦合模型的计算结果大。将计算结果与现场监测值进行对比，可以看出由完全耦合模型计算而得的结果与监测值更接近，进一步说明了本文所建模型的正确性。

5 结论

本文建立了基于M-C准则的弹塑性损伤-渗流耦合模型，通过UMAT子程序接口编写了模型的求解程序，最后利用该模型对南昌下穿暗涵盾构隧道开挖稳定性进行了数值计算，分析了注浆加固与未加固条件下地表沉降与暗涵沉降规律，开挖过程中的孔隙水压变化情况与围岩损伤变化规律，提出了合理的施工建议，验证了该模型的正确性和工程适用性。