非均质黏土隧道开挖环向稳定性的刚体平动有限元上限分析

刘国良， 黄飞衡， 罗静静， 杨 峰

(1. 中铁隧道勘察设计研究院有限公司， 广东 广州 511458； 2. 中南大学土木工程学院， 湖南 长沙 410075)

0 引言

黏土地层普遍存在非均质、各向异性等特点[1-3]。其中，非均质性体现在土体抗剪强度沿深度方向线性增加，这将影响隧道开挖过程中的地层稳定性。通常环向角度下的隧道开挖面稳定性，适用于分析盾构法隧道开挖工作段较长以及矿山法隧道刚性支护未施作的情况。在非均质黏土地层条件下，关于隧道环向开挖面稳定性的研究已有不少文献报道。黄茂松等[4]应用多块体极限分析上限法揭示了非均质黏土地层圆形隧道开挖面稳定性最优上限解及其规律。Wilson等[5]采用刚性块体上限法和有限元极限分析法，揭示了抗剪强度随深度线性增加时黏土地层方形隧道开挖面稳定性。Ukritchon等[6]以三维有限元极限分析方法研究了抗剪强度沿深度线性增加条件下的隧道稳定性，总结了埋深比、土体非均质强度、土体重度等因素的影响规律。Chen等[7]利用离散元技术改进的三维破坏模式，在极限分析理论框架下结合土拱效应评估了方形隧道开挖面稳定性。Zhang等[8]运用极限分析运动学方法，构建出马蹄形隧道地层连续速度场，探讨了不排水条件下黏土地层隧道开挖面失稳破坏机制。

对于抗剪强度沿深度线性增加的非均质纯黏土地层隧道纵向开挖面稳定性问题，文献[9]应用刚体平动运动单元上限有限元进行了多因素综合影响分析，揭示了滑移线网破坏模式的形态特征。对于非均质黏土地层隧道而言，沿着横断面的环向开挖面稳定性同样值得关注，特别是不同横断面形状的影响规律。此外，地层潜在破坏特征和扩展范围的定量分析，对于稳定性评价及工程措施的制定具有重要的理论意义，现阶段这方面的探讨仍然少见。刚体平动运动单元上限有限元采用非线性数学优化方法，具有搜索获取精细滑移线网破坏模式的特色功能，可作为破坏区域定量分析的有力手段[9-10]。

本文沿用刚体平动运动单元上限有限元方法[9-10]，对不排水抗剪强度沿深度线性增加的非均质黏土地层隧道环向开挖面极限状态下的稳定性进行研究，重点探讨圆形、马蹄形及椭圆形等4种代表性隧道断面形状下，不同的隧道埋深比、土体非均质强度和土体重度等参数与临界荷载比上限解及失稳破坏模式的关联规律，揭示滑移线网破坏模式的特征规律，给出破坏范围的定量几何尺度，以期为地层加固等工程措施提供参考与借鉴。

1 问题描述与计算模型

1.1 非均质黏土地层隧道环向开挖面稳定性模型

将非均质黏土地层隧道环向开挖面稳定性问题简化为沿横断面剖开的二维平面应变模型，即沿着隧道纵向力学状态一致，如图1和图2所示。如此简化与三维开挖面空间问题有所差异，不过用来分析沿横断面的环向开挖面潜在破坏模式的形态和扩展范围是合适的，从稳定性评价结论方面也是保守的。

(a) 圆形隧道

(b) 椭圆形隧道 (c) 椭圆形隧道 (d) 马蹄形隧道

图2 非均质黏土地层圆形隧道环向开挖面稳定性分析模型

工程现场常见圆形和马蹄形隧道断面形式，而矢跨比小的瘦高型隧道(如单线铁路隧道、连拱隧道中导洞等)以及矢跨比大的大跨型隧道(如3车道、4车道公路隧道等)，可近似为不同跨高比的椭圆断面，以此简化研究不同断面形式隧道环向开挖面稳定性和破坏形态的演化规律，可为隧道断面设计、变截面隧道施工等问题提供借鉴。

如图1所示，考虑4种断面形状的隧道环向开挖面： ①椭圆形断面(B/D=0.5)； ②圆形断面(B/D=1)； ③椭圆形断面(B/D=1.5)； ④马蹄形断面(上部半圆，下部矩形)。

由于隧道形状和计算模型受力均沿竖向中心线对称，故选取右侧一半开展计算分析(以圆形为例，如图2所示)。

本文考虑的非均质黏土地层为不排水抗剪强度变化规律仅沿深度线性增长的情况，参考文献[9]定义如下：

cu(y)=cu0+ρ(C+D/2-y)。

(1)

式中：cu(y)为竖向坐标y处对应的土体不排水抗剪强度；ρ为抗剪强度随深度变化率。

1.2 变量分析与参数确定

如图2所示，以圆形断面为例，隧道环向开挖面稳定性分析涉及7个计算参数：{C,D,σs,σt,cu0,γ,ρ}。将这些参数无量纲化[9]： 埋深比C/D，临界荷载比(σs-σt)/cu0，无量纲土体重度γD/cu0，无量纲土体非均质强度ρD/cu0。其中，临界荷载比(σs-σt)/cu0为待求解量，可视为稳定性评价指标。为方便计算，将(σs-σt)/cu0中的σt置为0；对于σt不为0的情况，将σs或σt具体值代入(σs-σt)/cu0(即σt=0时的σs/cu0值)，即可获得对应的σt或σs。计算参数取值见表1。

表1 隧道环向开挖面稳定性计算参数

1.3 运动单元上限有限元模型

以圆形隧道为例，非均质黏土地层隧道环向开挖面稳定性刚体平动运动单元上限有限元模型如图2所示。其中，地层不排水抗剪强度沿深度线性增长。上限有限元模型初始三角形网格类似文献[10]，本文未显示。利用初始网格可获取σs/cu0初始解，再进行一系列网格更新操作得到最终较优解答。求解流程与文献[9]一致，即依据上限定理将问题转换为非线性规划，在一系列约束条件下搜索σs/cu0最优(最小)上限解。

上限有限元模型非线性规划目标函数表达式为：

(2)

2 临界荷载比σs/cu0上限解

2.1 临界荷载比σs/cu0上限解结果

采用自编刚体平动运动单元上限有限元程序[9]，对表1无量纲参数的多种组合工况进行系列计算，获得4种断面形状的隧道于不同参数影响下的环向开挖面稳定性临界荷载比σs/cu0上限解，如表2所示。

2.2 临界荷载比σs/cu0上限解对比

为检验程序计算结果的可靠性，选取非均质地层圆形隧道环向开挖面稳定性临界荷载比σs/cu0上限解(ρD/cu0=0.5,γD/cu0=0～2)，与黄茂松等[4]、Wilson等[11]求解对比如图3所示。

极限分析上限解精度与破坏模式合理性直接关联[4]，UBRB(Wilson等[11])和多块体上限法(黄茂松等[4])需预先构造破坏机制进行上限分析，通常应用于求解圆形隧道时临界荷载比σs/cu0上限解的数值稍大。如图3所示，本文σs/cu0上限解曲线处于FEM(Wilson等[11])方法上、下限解之间，与LBFEM方法下限解平均相对差值仅2%左右，说明程序可靠性良好。

2.3 临界荷载比σs/cu0上限解的分析与讨论

以圆形隧道为例，选取典型数据绘制σs/cu0上限解曲线如图4所示，以此分析无量纲土体重度γD/cu0、无量纲土体非均质强度ρD/cu0以及埋深比C/D对于隧道环向开挖面稳定性的影响规律。

由图4(a)可知，C/D=1时σs/cu0上限解与γD/cu0呈显著的线性负相关，γD/cu0增加，地层稳定性随之下降；ρD/cu0=0.05～1.0对应曲线由下至上依次等间距平行排列，说明σs/cu0上限解与ρD/cu0呈线性正相关。图4(b)和4(c)分别对应于C/D=3、5的情况，此时σs/cu0上限解曲线规律与图4(a)相同，而对应的σs/cu0上限解数值增加。

表2 隧道环向开挖面稳定性临界荷载比σs/cu0上限解

图3 圆形隧道环向开挖面稳定性临界荷载比σs/cu0上限解对比 (ρD/cu0=0.5)

由图4(d)可知，ρD/cu0=0.05时，σs/cu0上限解与C/D呈近似线性的增加或减小规律，其中γD/cu0较大时近似线性负相关，而γD/cu0较小时近似线性正相关。图4(e)和4(f)分别对应于ρD/cu0=0.5、1.0的情况，此时σs/cu0上限解与C/D近似呈线性正相关规律。

由图4(g)可知，γD/cu0=0时，随着ρD/cu0增大，σs/cu0上限解表现为线性增加[12]。比较图4(g)—(i)发现，无论γD/cu0取何值，C/D较大时σs/cu0上限解增长速度更加明显。图4(h)—(i)各曲线近似交于一点，说明存在一种特定状态，σs/cu0上限解随C/D的发展趋势在此刻发生转折。

(a) C/D=1 (b) C/D=3 (c) C/D=5

(d)ρD/cu0=0.05 (e)ρD/cu0=0.50 (f)ρD/cu0=1.00

(g)γD/cu0=0 (h)γD/cu0=2 (i)γD/cu0=3

其他条件相同时，对比表2中不同断面形状隧道的临界荷载比σs/cu0上限解，发现隧道环向开挖面稳定性规律如下： 椭圆形隧道和圆形隧道临界荷载比σs/cu0上限解由大到小依次为①椭圆形隧道(B/D=0.5)、②圆形隧道、③椭圆形隧道(B/D=1.5)； 其中，①椭圆形隧道较③椭圆形隧道临界荷载比σs/cu0上限解高出15%～40%，且随着埋深比C/D增加，差值进一步增大，反映出隧道跨度是影响稳定性的关键因素之一。此外，与圆形隧道轮廓接近的④马蹄形隧道，其临界荷载比σs/cu0上限解普遍稍小一些。

3 隧道环向开挖面破坏模式

3.1 破坏形态对比

图5示出C/D=4、ρD/cu0=1、γD/cu0=0时2种极限分析上限法得到的圆形隧道环向开挖面潜在破坏模式。其中，图5(a)示出UBFEM(Wilson等[11])方法对应的破坏形态，以速度矢量形式展示。图5(b)示出本文刚体平动运动单元上限有限元计算得到的滑移线网破坏模式。需要说明的是，这里的滑移线即为模型里面活动的或者有效的速度间断线，与主应力的迹线关联。而模型中的大量滑移线通常表现为网状的形式，可等效为破坏区域中的塑性区。

(a) UBFEM法 (b) 本文方法 (c) 速度矢量图

由图5可知，2种上限法展示的隧道环向开挖面破坏形态相似，表现为： 图5(a)洞顶正上方速度矢量一致的规则区、拱顶左上方速度矢量过渡区，分别对应于图5(b)洞顶正上方整体下沉区及右拱腰上方网状滑移区。此外，图5(c)为速度矢量闭合图，以模型中最大数值V为参考，对速度矢量值进行了归一化处理。速度矢量闭合图中每条线段对应于图5(b)中特定的有效间断线，而图中右上方交点与下方任意交叉点的连线，反映了图5(b)中特定刚性单元的速度矢量。

由图5(b)可知，本文圆形隧道环向开挖面的失稳破坏形态，与太沙基隧道土压力理论假定的塑性剪切带和扇形塑性变形区的破坏形态基本吻合(此时内摩擦角为0)。图5(b)展示的潜在失稳破坏模式，可提供破坏区域界限及面积、滑移面角度等信息，为破坏机制的量化研究提供数据。

3.2 破坏模式演变规律分析

选取圆形隧道，绘制典型参数对应的环向开挖面潜在破坏模式如图6所示，以分析埋深比C/D、无量纲土体重度γD/cu0和无量纲土体非均质强度ρD/cu0等参数对破坏模式演变规律的影响。

图6(a)示出C/D=2、γD/cu0=3、ρD/cu0=0.5时的滑移线网破坏模式，其形态同图5(b)一致。即上方为整体下沉规则区，外侧滑动面竖直；右拱腰附近为网状滑移区，其2簇滑动面交线近似垂直，而外侧滑动面构成平滑曲线过渡。将相同条件下γD/cu0=0～3对应的主要滑动面叠加绘制如图6(b)所示，发现γD/cu0增大时，破坏区域面积有所增加，且稍有下移趋势。

图6(c)示出C/D=2、γD/cu0=1、ρD/cu0=0.5时的滑移线网破坏模式。将相同条件下ρD/cu0=0.05～1.0对应的主要滑动面叠加绘制如图6(d)所示，发现随着ρD/cu0增大，破坏区域面积反而明显减小，破坏区域有向上和向中心线靠拢的趋势，这符合土体强度增加、破坏范围随之减小的规律。

图6(e)示出C/D=3、γD/cu0=1、ρD/cu0=0.5时的滑移线网破坏模式。将相同条件下C/D=1～3对应的主要滑动面叠加绘制如图6(f)所示，发现随着C/D增大，破坏区域面积显著增加，破坏区域有向下和向外扩张的趋势，说明对于非均质纯黏土地层(内摩擦角为0)，隧道埋深是影响破坏范围的关键因素之一。

选取C/D=2、ρD/cu0=0.5、γD/cu0=1时，绘制不同断面形状对应的隧道环向开挖面潜在破坏模式如图7所示。由图7(a)可知，直墙马蹄形隧道破坏范围易于扩展至墙角处； 而图7(b)显示的跨度较大的椭圆形隧道(B/D=1.5)，其潜在破坏范围反而局限于隧道拱部区域。在相同参数条件下，将4种不同断面形状隧道对应的主要滑动面叠加绘制如图7(c)所示，可以看出水平方向破坏区域范围由大到小排序为： ④直墙马蹄形隧道>①椭圆形隧道(B/D=0.5)>②圆形隧道>③椭圆形隧道(B/D=1.5)。

实际上，破坏区域的范围仅反映了非均质黏土地层隧道环向开挖面破坏特征的一个方面，其他方面如破坏状态下的速度场也是关键要素。图7(d)和7(e)分别示出相同参数条件下圆形隧道(B/D=1.0)和椭圆形隧道(B/D=1.5)滑移网破坏模式对应的速度矢量闭合图。图中均以最大矢量值进行归一化处理，可以看出圆形隧道对应的速度矢量值普遍大一些，特别是水平方向的速度，由此带来临界状态下耗散能的增加，反衬出跨度较大的椭圆形隧道(B/D=1.5)尽管破坏范围小，但速度矢量值稍小，系统耗散能降低，反而不利于地层稳定。

(a) (b) (c) (d) (e) (f)

(a) (b) (c) (d) (e)

3.3 破坏区域面积定量分析

破坏范围的定量分析，结合对应速度场的特征规律，可为地层加固范围的确定提供数据支撑。为进一步说明隧道失稳极限状态下破坏区域面积随综合因素的演变规律，定义无量纲参数Sf/Sc为相对破坏范围，Sf为破坏区域面积(筛选模型速度不为0的刚体单元，累计单元面积之和)，Sc为隧道轮廓面积，Sf和Sc均考虑右侧的半模型。以圆形隧道为例，不同参数条件下的相对破坏范围Sf/Sc计算值见表3，选取C/D=2、4的数据绘制曲线如图8所示。

由图8和表3可以看出，圆形隧道相对破坏范围Sf/Sc随着C/D和γD/cu0的增加而增大，并随着ρD/cu0的增加而减小，这直接以定量数据反映了图6所示的演变规律。

此外，埋深比C/D变化时，Sf/Sc随之变化非常明显。如C/D=5较之C/D=1的破坏范围普遍增长了10倍以上，说明隧道埋深较大时，受地表超载作用，破坏面水平向延伸范围较大； 同时，对于不排水抗剪强度沿深度线性增长的非均质黏土地层(内摩擦角为0)，潜在破坏面易于扩展至地表，于是Sf/Sc呈现出与(C/D)2近似线性相关的特征。

Sf/Sc随γD/cu0增加而线性增长，说明自重使得隧道附近土体更易于产生剪切滑移破坏；不过该增长比较缓慢，尤其是ρD/cu0较大时。由图8可知，Sf/Sc随着ρD/cu0增加而非线性减小，即土体抗剪能力提高，破坏范围受到极大的限制从而缩小。

4 σs/cu0上限解拟合公式

如图4所示，环向开挖面稳定性临界荷载比σs/cu0上限解与γD/cu0、ρD/cu0呈线性关系，且线性斜率值和常数项仅与埋深关联，故C/D固定时σs/cu0上限解可认为是无重度均质土上限解在γD/cu0、ρD/cu02变量作用下的线性变化。于是，σs/cu0上限解可写为：

σs/cu0=Bγ·γD/cu0+Bρ·ρD/cu0+B0。

(3)

式中：Bγ为σs/cu0上限解与γD/cu0线性关系斜率；Bρ为σs/cu0上限解与ρD/cu0线性关系斜率；B0为无重度均质土σs/cu0上限解。

表3 圆形隧道相对破坏范围Sf/Sc计算值

以圆形隧道为例，对不同埋深对应的σs/cu0上限解公式拟合如下：

(4)

可以看出Bγ、Bρ、B0随C/D变化规律性较强，数据拟合可得：

Bγ=-1.034·C/D-0.141；

(5)

Bρ=4.526·C/D-3.154；

(6)

B0=-0.096·(C/D)2+1.229·C/D+1.372。

(7)

故圆形隧道σs/cu0上限解最终拟合公式为：

σs/cu0=(4.526·C/D-3.154)·ρD/cu0+(-1.034·C/D-0.141)·γD/cu0-0.096·(C/D)2+1.229·C/D+1.372。

(8)

同样地，拟合其余3种断面隧道σs/cu0上限解公式见表4。其确定系数R2均大于0.99，表明拟合精度较好。上限解拟合公式只需少量参数(C,D,γ,ρ,cu0)即可初步评估抗剪强度随深度线性增长的纯黏土地层隧道环向开挖面稳定性，方便工程应用。

(a) C/D=2

(b) C/D=4

5 试验验证

为进一步检验刚体平动运动单元上限有限元方法的适用性，选取基于离心机模型试验[13]的均质软黏土地层单孔隧道临界荷载比σs/cu0结果进行对比验证。其试验采用直径D=100 mm的圆形隧道模型，土体重度γ=18.1 kN/m3，均质地层中抗剪强度随深度变化率ρ=0。试验过程中控制土体孔隙水压力变化不超过1%以保证不排水状态。试验考虑不同埋深C、不同抗剪强度cu0条件下隧道环向开挖面失稳破坏，临界荷载比σs/cu0试验结果与本文运动单元上限解及拟合公式上限解对比如表5所示。

表4 4种断面隧道环向开挖面稳定性临界荷载比σs/cu0上限解拟合公式

表5 临界荷载比σs/cu0试验结果与本文上限解对比

由表5临界荷载比σs/cu0数据对比可知，尽管数值上存在不同程度的差异，刚体平动运动单元上限解与离心模型试验结果总体上吻合良好，从而验证了本文计算方法的可靠性和计算结果的精度； 同时，拟合公式对应的σs/cu0上限解也具有较好的一致性。

6 结论与讨论

1)在抗剪强度沿深度线性增加的非均质黏土地层条件下，不同断面形状隧道环向开挖面稳定性差异显著，椭圆形隧道(B/D=0.5)较之椭圆形隧道(B/D=1.5)临界荷载比σs/cu0上限解高出15%～40%，且埋深增加将增大这一差值。σs/cu0上限解与ρD/cu0正相关而与γD/cu0负相关，ρD/cu0和γD/cu0均影响C/D对上限解的作用趋势。

2)由隧道环向开挖面潜在破坏模式形态演变及范围定量分析发现，破坏范围主要取决于断面形状和埋深大小。较之圆形隧道，直墙马蹄形隧道失稳破坏区域易扩展至边墙底部。此外，埋深成为影响隧道失稳破坏范围的关键因素。

3)基于刚体平动运动单元上限有限元数据，获得了非均质黏土地层不同形状隧道的临界荷载比σs/cu0上限解拟合公式。对比已有离心机模型试验结果，验证了运动单元上限解及拟合公式的可靠性，可为工程应用提供参考。

4)考虑到现阶段分析模型与工程实际比照存在不少简化，下一步的研究应考虑荷载的非均匀分布以及地层成层条件等因素的影响。