考虑黏土土性参数与支护压力变异性的盾构掌子面稳定性分析

杨文钰，郑俊杰，章荣军，乔雅晴

(1.华中科技大学 土木与水利工程学院，武汉 430074；2.武汉大学 土木建筑工程学院，武汉 430072)

随着城市的发展，盾构隧道成为充分利用地下空间最为有效的手段之一。掌子面稳定性是盾构隧道的经典问题，合适的掌子面支护压力不仅能够保证盾构施工的安全，同时也保障了盾构隧道周围建构筑物的安全。

对于盾构隧道开挖面失稳问题，大量的学者运用理论分析法[1-3]、数值模拟法[4-6]、模型试验法[7-9]已经得到了丰富的研究成果。然而，土体的天然变异性与试验或场地的限制导致了人们对岩土体的认知缺乏，这决定了岩土体参数的不确定性[10]。以上这些研究都把土体视作均质且各向同性的材料，在预测极限支护压力时可能会存在偏差。近年来，有学者将土性参数的不确定性考虑在盾构隧道掌子面稳定性的研究中，Mollon等[11-14]基于随机响应面等方法，研究了抗剪强度参数的变异性对掌子面稳定性的影响，揭示了土性参数变异性对极限支护压力的影响的基本规律。以上研究将土性参数视作随机变量，为进一步考虑土性参数的空间变异性，Mollon等[15]首先基于极限分析理论，得到了考虑抗剪强度参数空间变异性二维掌子面失效机制，Cheng等[16-19]考虑了砂土、黏性土以及成层土抗剪强度参数的空间变异性，研究表明，抗剪强度空间变异性对掌子面稳定性有重要影响，忽略这一特性可能高估掌子面的稳定性；极限支护压力与抗剪强度参数的变异系数、自相关距离密切相关。以上研究较全面地揭示了土性参数的不确定性对盾构隧道掌子面稳定性的影响规律，然而，鲜少有研究考虑掘进参数的不确定性。

在盾构隧道施工过程中，掌子面土体经过刀盘的切削作用进入土舱内部，掌子面处土体的稳定靠土舱压力与掌子面处水土压力相平衡来保持，盾构机推进油缸的推力、推进速度、螺旋输送机的出土量[20]与刀盘开口率[21]等盾构机的参数都对土舱压力有一定的影响。通常情况下，为保持盾构掌子面的稳定，先根据地质条件设定土压力值，在盾构掘进过程中再根据土压传感器的变化人为做出调整[22-23]。但目前土舱压力值的计算尚无固定方法[24]，一般参考工程地质、盾构机型等因素，由施工经验确定。工程地质条件的复杂性决定了土舱压力的值并非定值[25]；另外，盾构机的人为控制可能存在操作不当、违规操作等问题，也会导致盾构施工中土舱压力具有一定的变异性。因此，忽略支护压力(本文中提及的支护压力仅限于盾构中的土舱压力)变异性可能无法反映工程实际的复杂性，相应地可能会高估掌子面稳定性。基于此，笔者采用随机场理论与数值模拟分析相结合的方法，研究了黏性土内摩擦角与黏聚力空间变异性共同作用下对盾构隧道掌子面稳定性的影响，分析了黏聚力与内摩擦角变异性共同作用下对掌子面失稳模式和极限支护压力的变化规律，并结合概率分析法探讨了同时考虑土性参数空间变异性与支护压力变异性的可靠度分析方法。

1 数值模型的建立

1.1 问题描述

建立数值模型的核心问题是盾构隧道在黏性土掘进过程中的掌子面稳定性。盾构隧道掌子面稳定性的数值模型参考了Mollon等[12]的研究，数值模型如图1所示。采用的数值模型隧道直径(衬砌管片外径)D为10 m，埋深C(隧道拱顶距地表距离)为10 m，为简化分析，不考虑地下水的影响。整个模型的计算域为50 m×40 m×26 m(长×宽×高)，经验证，计算域的大小能使计算结果不受边界的影响[12]。数值模型共有52 240个单元。由于采用应力控制法分析盾构隧道掌子面稳定性时，掌子面上会出现较高的应力梯度，所以对掌子面以及掌子面后方土体单元进行了网格的加密，掌子面上共分为了198个单元。模型的底部为固定边界，四周为法向位移约束边界，顶部为自由边界。为了简化分析，数值模型仅模拟了黏性土与衬砌管片两种材料，本构模型采用Mohr-Coulomb模型，衬砌管片采用线弹性模型，衬砌管片的厚度为0.4 m。黏性土与衬砌管片材料参数如表1所示。

图1 数值模型示意图Fig.1 Schematic diagram for the numerical modelling of the problem in this

表1 黏土与衬砌材料参数Table 1 Properties for the clay and the lining segment

1.2 改进二分法

由于掌子面上的位移事先难以确定，所以采用较为常用的应力控制法进行计算。考虑到不确定性计算中的Monte-Carlo策略需要多次计算，选择了计算效率较高的改进二分法。

改进二分法基于简单二分法，简单二分法的流程可简单分为：1)确定计算上下限；2)将上下限的中间值代入计算，并根据计算结果更替上下限；3)不断重复前两个步骤直到计算结果达到特定的精度值。简单二分法的本质是进行几个稳定与不稳定的循环计算。然而，塑性流动往往出现在大量的计算步之后，尤其是对结果的精度要求苛刻时，将会耗费大量的时间。针对这一不足，Mollon等[12]提出了改进二分法，改进二分法的步骤如下：首先将土体的黏聚力赋予一个很大的值，这让土体变成弹性材料；接下来，手动设置内部应力为初始值的两倍，统计系统重新回到平衡状态所需要的计算步数N，根据Mollon等[12]的研究，N的值接近于3 000步；N的值确定后，将初始黏聚力重置为真实值，若循环N步后系统仍不平衡，则可认为该工况会进入塑性流动状态。由此可见，相比简单二分法，改进二分法大大的缩短了塑性流动状态工况的计算时间。所以，笔者采用改进二分法分析盾构隧道掌子面稳定性。

1.3 数值模拟的验证

在二分法中，选取合适的上下限对计算准确的极限支护压力非常重要。选择60 kPa与20 kPa作为改进二分法的上下限。图2展示了当计算步数N为2 700时，掌子面中心点A(A点的几何位置见图1)位移与速度随计算步数的变化。支护压力为60 kPa时，在2 000步左右时，A点的水平位移已经趋于定值，A点的水平速度降为0；而当支护压力为20 kPa时，计算2 700步后，A点的水平位移与水平速度仍在增长。说明对于表1参数的黏性土而言，分析掌子面稳定性的改进二分法中N取2 700步，上下限取60 kPa与20 kPa是合理的。

图2 掌子面A点水平位移与水平速度随计算步数的变化Fig.2 Variation of the horizontal displacement and the horizontal velocity of point A on the tunnel face against the calculation

为验证数值模型的合理性，首先运用了逐级减小支护压力的应力控制法。支护压力从100 kPa开始逐步减小，为使系统达到塑性流动的状态，在每一支护压力所对应的工况下计算10 000步。图3是掌子面中心点水平位移-支护压力变化曲线，如图所示，曲线的拐点在34 kPa处，即极限支护压力为34 kPa。

图3 掌子面A点水平位移支护压力变化曲线Fig.3 Variation of the horizontal displacement and of point A on the tunnel face against support

运用改进二分法进行计算，设置精度为1 kPa，计算结果为34.53 kPa，此时相应的速度场如图4所示。这个结果不仅与应力控制法的计算结果相同，也与理论解[26]与数值计算结果[12]相符，这证明了模型与改进二分法的正确性。

图4 支护压力为34.53 kPa时的速度场Fig.4 Velocity field when the support

2 不确定性分析

2.1 基于随机场理论的盾构隧道掌子面稳定性分析

运用随机场理论对盾构掌子面的稳定性进行不确定性分析。采用随机场分布光滑度与连续性较好的高斯型自相关函数[27]，能够有效的描述土性参数的空间自相关性。高斯型自相关函数主要的形式如式(1)所示。

(1)

式中：τx、τy、τz分别为空间两点在x、y、z方向上的相对距离；δx、δy、δz分别为x、y、z方向的波动范围。在本文中，各向异性随机场的水平向波动范围δx、δy为20 m，竖向波动范围δz为2 m。中心极限定理表明[28]，受大量不确定性因素影响的因变量通常近似服从正态分布或对数正态分布，同时对数正态分布严格非负[29-30]，这与岩土体参数的概率分布相符合，因此，采用对数正态分布描述黏聚力与内摩擦角的不确定性。

K-L级数分解法具有运算效率较高、生成随机场精度较好、对于各向异性随机场的生成有较强的适应性等优点，被广泛应用。K-L级数分解法将土体参数随机场H(A,θ)(其中连续坐标A∈Ω⊆Rn，目标空间中坐标θ∈Θ)的离散转化为求解Fredholm积分方程的特征值问题，这一特征值问题如式(2)。

(2)

式中：A1、A2为目标离散空间Ω中的任意两点坐标，ρ(A1,A2)为这任意两点处随机场特征值之间的相关函数值，λi和fi分别为与相关函数对应的特征值和特征函数。关于该特征值的求解具体可参考文献[31]。因此，随机场相关函数的特征值计算可以简化为相关函数特征值的乘积，随机场即可离散为

(x,y,z∈Ω)

(3)

式中：ξi(θ)为独立标准正态随机向量，θ∈Θ为外部空间坐标；(x,y,z)为随机场区域中的任意坐标点，与A1、A2对应；μ、σ分别为随机场的均值与标准差。为实现高效计算，通常在保证精度的前提下，截取式(3)的前n项来提高计算效率，n的选取参考文献[31]中比率因子ε的大小进行调整，其定义为

(4)

通常ε大于0.95即可认为随机场精度已经满足计算要求。截取后的前n项为

(x,y,z∈Ω)

(5)

基于此，采用K-L级数分解法建立黏聚力与内摩擦角的各向异性随机场。

采用数值模型及改进二分法，模型中仅考虑黏性土内摩擦角与黏聚力空间变异性对盾构掌子面稳定性的影响，其他参数均为常量。保持内摩擦角与黏聚力的均值与确定性计算的参数一致，重点研究内摩擦角与黏聚力共同变异的情况下对掌子面失稳模式、极限支护压力的影响。分别选取黏聚力的变异系数COV(c)为0.1、0.2、0.3、0.4，内摩擦角的变异系数COV(φ)为0.05、0.01、0.15、0.2，组合得到16组工况建立随机场进行计算分析。为方便后续分析，工况名命名如表2所示。

表2 工况名称表Table 2 The names of numerical modelling cases

每组工况分别进行500次随机计算。由于研究结果表明[18]，抗剪强度参数的负相关关系相对其变异系数对极限支护压力的影响较小，因此，不考虑抗剪强度参数的负相关。

2.2 掌子面失稳模式

以黏聚力变异系数COV(c)为0.4，内摩擦角变异系数COV(φ)为0.2这组工况为例，对掌子面失稳的模式进行分析。

选取3组典型工况进行分析。图5为3组典型工况黏聚力、内摩擦角的随机场分布。图6是3组典型工况的最大剪应变增量云图。可以发现，考虑土性参数的空间变异性时，虽然土性参数随机场的均值与变异系数相同，不确定性分析计算得到的开挖面失稳模式也不尽相同，失稳区域的大小和形式也都有一定的区别。很显然，掌子面失稳模式与掌子面局部区域土性参数随机场的分布相关。几组工况中，工况a、工况c发生整体破坏，工况b发生局部破坏。工况b的掌子面前方内摩擦角随机场呈现明显的分层，隧道轴线上方内摩擦角大、下方内摩擦角小，土体沿着强度软弱处开始破坏，逐渐扩展、传递，导致隧道轴线下方出现局部破坏。

图5 典型工况的黏聚力与内摩擦角随机场分布Fig.5 Thecohesion and the friction random field

图6 典型工况的最大剪应变增量云图Fig.6 Maximum shear strain increment contour

为定量刻画掌子面局部区域土性参数随机场分布与掌子面失稳模式间的关系，引入归一化极限支护压力σc_nor,掌子面局部归一化均值cave_nor、φave_nor、标准差cdev_nor、φdev_nor与变异系数COV(c1D)、COV(φ1D)作为量化指标进行探究，这些量化指标按式(6)～式(9)计算。

(6)

(7)

(8)

(9)

式中：cave_1D、φave_1D为c、φ掌子面前方一倍直径内的均值，cdev_1D、φdev_1D为c、φ掌子面前方一倍直径内的标准差。cave、φave分别为7 kPa、17°。σc_det为极限支护压力的确定性计算结果。

表3列出了3组典型工况的归一化统计参数。不难发现，归一化均值指标决定了σc_nor的大小，cave_nor与φave_nor越大，σc_nor越小。而σc_nor与归一化标准差cdev_nor和φdev_nor、归一化变异系数COV(c1D)和COV(φ1D)无明显关系。为进一步说明以上结论，分析500组工况的归一化极限支护压力与局部归一化黏聚力与内摩擦角随机场的归一化指标间的关系，绘制图7。对比图7(a)～(e)可以发现，极限支护压力归一化指标与σc_nor的相关系数均为负数，σc_nor与归一化均值指标呈现较明显的负相关关系，而与归一化标准差和归一化变异系数相关系数很小。

表3 典型工况的掌子面局部区域归一化参数Table 3 Local normalized parameters for tunnel face of typical cases

图7 归一化极限支护压力与归一化掌子面局部区域参数间的关系散点图Fig.7 Relationship between normalized support pressure and normalized local parameters for tunnel

结合图5～图7与表3分析归一化统计参数对失稳模式的影响可以发现，工况b中φdev_nor与COV(φ1D)较工况a、工况c大，掌子面发生局部破坏，说明φdev_nor，COV(φ1D)的值能够一定程度上解释该工况发生局部破坏的原因。这是因为，当掌子面局部随机场分布较为分散时，可能会出现如图5(b)中的内摩擦角局部随机场的上下分层情况，这种情况很有可能导致局部破坏。

2.3 极限支护压力

为进一步分析土性参数的空间变异性对极限支护压力σc的作用，对比了16组黏聚力变异系数COV(c)与内摩擦角变异系数COV(φ)不相同的工况。

图8绘制了工况1～工况4以及工况13中σc的概率分布直方图。所有工况中，σc的中位值与均值均大于σc_det，很显然土性参数的空间变异性对σc的影响不容忽视。将工况1视作基准工况进行分析。工况4与工况13的COV(φ)与COV(c)分别在基准工况的基础上增大了3倍。相应地，相比工况1，工况4的σc中位值增大了6.27%，均值增大了5.51%，标准差与变异系数均增大了2倍左右，波动范围从29.33～39.19 kPa增大到22.68～50.35 kPa；工况13的中位值则增大了2.26%，均值增大了1.29%，标准差与变异系数均增大了1倍左右，波动范围从29.33～39.19 kPa增大到21.02～45.87 kPa。很显然，随着c、φ的变异性增大，σc的统计参数均相应地增大，但φ变异性增大对σc的统计参数造成的影响更甚。分析所有工况σc的统计特征值可以发现，σc的统计特征值与COV(c)、COV(φ)均成正比关系，c、φ的变异性共同作用时，对σc的影响较某一土性参数单独作用时更大，会使掌子面的稳定性更趋于不安全，所以，对于内摩擦角与黏聚力同时具有变异性的场地地质条件，应当给予重视。

图8 极限支护压力频率分布直方图Fig.8 Relationship between normalized support pressure and normalized local parameters for tunnel

为进一步量化土性参数空间变异性对σc的影响，将现有的16组工况中的数据进行拟合如图8，拟合结果表明，σc服从正态分布，未在图8中出现的工况的拟合参数见表4。

表4 各工况极限支护压力正态分布拟合参数Table 4 Fitting parameters for normal distribution of the critical support pressure for every cases

3 考虑支护压力变异性的盾构掌子面可靠度分析

由于忽略支护压力的变异性可能会高估掌子面的稳定性，采用蒙特卡洛策略将支护压力的变异性考虑进盾构隧道掌子面可靠度分析中，同时为将土性参数的变异性考虑进来，将实际的掌子面极限支护压力视作服从表4中的拟合分布结果的随机变量。对于支护压力的概率分布，由于现有的研究中对支护压力的统计规律并无统一定论，采用某一工程实测数据[32]的统计结果作为依据进行分析。该工程中土舱压力的变异系数为0.114，并且符合正态分布，故支护压力服从正态分布，且变异系数的上限设置为0.1。为使结果有效，蒙特卡洛模拟进行100万次，失效概率定义为

(10)

式中：N为蒙特卡洛模拟的总计算次数，即1×106次，Nf指的是支护压力小于极限支护压力的次数。

图9表示支护压力服从正态分布、变异系数为0.1时，支护压力均值的变化对失效概率的影响。随着支护压力均值的增大，失效概率逐渐减小为0%，支护压力均值小于24 kPa或大于60 kPa，此时土性参数的变异性对失效概率的影响不大。观察工况1～工况4或工况4、工况8、工况12、工况16的失效概率随支护压力均值的变化，当COV(c)与COV(φ)在增大时，失效概率随支护压力均值的增大而减小的速率在慢慢变缓，但减缓的路径各不相同。这种现象在工况1～工况4中观察得更明显，在工况4、工况8、工况12、工况16中，随着黏聚力的变异系数增大，失效概率-支护压力均值的变化曲线虽然在减缓，但是减缓并不明显。这说明在考虑支护压力变异性时，内摩擦角的变异性对掌子面稳定性失效概率的影响更甚。通过以上的分析可知，支护压力的均值对失效概率的影响较大，所以在实际工程中，需要严格把控支护压力的均值来控制掌子面的稳定性。为选取合适的支护压力均值，引入支护压力均值特征值的概念，借鉴概率统计的方法，将失效概率与支护压力均值特征值结合起来，在此规定失效概率为5%时为支护压力均值特征值σk。在实际应用中，可在盾构的不同时段通过传感器监测土舱压力，将数据进行正态分布的拟合，如果该分布中的支护压力的均值小于支护压力均值特征值σk，那么应当引起重视。图9中虚线与曲线的交点代表的就是各工况支护压力特征值σk。

图9 支护压力均值对失效概率的影响Fig.9 The impact of critical support pressure mean value

实际工程中，支护压力受工程地质条件与人为因素的限制而非定值，具有一定的变异性，接下来分析支护压力变异系数对掌子面稳定性的影响。图10表示支护压力变异系数对支护压力均值特征值σk的影响。总体而言，σk随着支护压力变异性增强而增大，近似线性关系。对于工况1～工况4，σk都随支护压力的变异系数的增大而增大，支护压力变异系数从0.01增大到0.1，4组工况的σk分别增大了12.88%、9.57%、8.06%、6.49%，增幅随内摩擦角变异系数增大而逐步减小。这说明，随着内摩擦角的变异性增强，支护压力的变异性对σk的影响会被减弱。对于工况4、工况8、工况12、工况16，σk都随支护压力的变异系数的增大而增大，但增幅并不随工况的黏聚力变异系数变化而变化，支护压力的变异系数从0.01增大到0.1，4组工况的σk分别增大了6.49%、6.15%、5.98%、6.00%。这说明，支护压力的变异性对σk的影响主要受到内摩擦角变异性的限制，而受黏聚力的变异性影响较小。

图10 支护压力变异系数对支护压力特征值的影响Fig.10 The impact of critical support pressure coefficient of variation on thecritical support pressure characteristic

综上所述，考虑支护压力变异性的盾构隧道掌子面稳定性分析中，确定合理的支护压力均值特征值是保证掌子面稳定性的关键。σk受土性参数与支护压力变异系数的影响，将σk的值与确定性计算的结果、土性参数、支护压力的变异系数结合起来，能够得到式(11)～式(13)。

σk=[k·COV(σ)+b]·σc_det

(11)

k=a1+a2·COV(c)+a3·COV(φ)·ln(COV(φ))

(12)

b=a4+a5·COV(c)+a6·COV(φ)·ln(COV(φ))

(13)

图10表明σk与支护压力的变异系数呈线性关系，其截距k、斜率b与土性参数变异系数COV(c)、COV(φ)相关，可以由式(12)与式(13)表示，其中常量a1=1.828、a2=-0.638、a3=2.234、a4=0.844、a5=0.297、a6=-1.251。图11展示了通过式(11)预测的支护压力均值特征值σk与通过数值模拟与蒙特卡洛模拟得到的σk的对比，能够发现大部分数据落在了1∶1线附近，拟合精度为0.89，表示拟合公式的预测效果不错。

图11 σk的数值计算值与式(11)的预测值对比Fig.11 Comparison between the values of σk obtained from numerical modelling and Eq.

4 结论

运用随机场理论与数值模拟结合的盾构隧道开挖面稳定性分析方法，同时考虑了黏性土黏聚力与内摩擦角的空间变异性，进一步考虑了服从正态分布的支护压力的变异性，研究了土性参数与支护压力变异性对盾构隧道掌子面稳定性的影响，得到以下结论：

1)同时考虑黏性土黏聚力与内摩擦角变异性时，盾构隧道掌子面的失稳模式与掌子面局部区域的土性参数随机场分布相关，土性参数在掌子面前方一倍直径处的局部随机场均值越小，掌子面越不稳定，而土性参数在掌子面前方一倍直径处的局部随机场的标准差或变异系数越大，发生局部破坏的可能性越大。

2)黏性土黏聚力与内摩擦角的空间变异性对盾构隧道掌子面稳定性有重要的影响，其中内摩擦角对掌子面稳定性的影响更甚。黏聚力与内摩擦角对极限支护压力有正向的作用，随着黏聚力与内摩擦角的变异性增大，对掌子面稳定性的影响逐渐增强，掌子面越不稳定。

3)支护压力的不确定性对盾构隧道掌子面稳定性也有一定的影响，支护压力均值越大，变异系数越小，掌子面越稳定。然而，支护压力的变异系数对掌子面稳定性的影响受到内摩擦角变异性的限制，内摩擦角的强变异性会减弱支护压力变异性对掌子面稳定性的影响。

4)提出了掌子面支护压力均值特征值的概念，结合掌子面的失效概率、极限支护压力的确定性结果、土性参数的变异系数以及支护压力的变异系数，对掌子面支护压力均值特征值给出了初步的确定方法。