卵石倾角对砂卵石地层隧道开挖影响的细观力学研究

2022-10-11刘泽昊齐吉琳杜修力

工程力学 2022年10期

张佩，刘泽昊，齐吉琳，杜修力

(1. 北京建筑大学土木与交通工程学院，北京 102616；2. 香港理工大学土木与环境工程系，香港 999077 ；3. 北京工业大学城市与工程安全减灾教育部重点实验室，北京 100124)

砂卵石地层是由具有一定尺寸、强度较高的岩块和强度较低的土体构成的不均匀松散岩土介质，为典型的土石混合体地层[1-3]。由于组成成分的物理力学性质差异很大，砂卵石地层表现出显著的非均质性，其材料力学行为和场地力学行为将显著地受到内部细观结构的影响。

对于单独的卵石颗粒，在三维分析时通常呈椭球形或球形，在二维分析时呈椭圆形或亚圆形。在实际地层中，由于河流的冲洪积作用及沉积作用，椭球形或亚圆形的卵石将在地层中水平分布或与水平面呈一定夹角，呈现出定向性排列的现象。卵石在地层中的沉积方向通常采用卵石长轴倾角来描述，其是指卵石长轴与水平面之间的夹角。王振峰[4]对西安地铁2 号线沿线卵石进行了统计，结果表明：卵石含量随长轴倾角增加呈近似指数下降的趋势。基于数值图像技术，马辉等[5]对成都地铁1 号线沿线卵石的统计也给出了相同的结论，并得到卵石长轴倾角平均值为32°。牟迪[6]对成都地铁4 号线沿线卵石进行了统计，结果表明卵石长轴倾角集中在30°～60°区间。通过对三峡库区的土石混合体进行统计，张抒和唐辉明[7]指出，块石含量与块石长轴倾角间呈现出均匀分布、正态分布、伽马分布三种形式。XU 等[8]研究表明：块石数目随倾角的增加呈现出“M”形的变化规律，长轴倾角集中在10°～40°与130°～170°两个区间。

目前，各国学者已在卵石倾角对材料力学行为的影响方面展开深入研究[9-11]，而关于卵石倾角对场地力学行为方面的研究鲜有报道。层状岩体隧道施工力学方面的研究表明，当岩层倾角发生变化时，围岩变形的分布特征和破坏模式将发生显著变化[12-15]。因此，基于在层状岩体施工力学方面的认识，考虑砂卵石地层材料的细观结构特征，开展不同卵石倾角下的隧道开挖计算，探讨卵石倾角对围岩应力、围岩变形的影响规律，具有重要的工程意义。

当采用数值方法开展隧道开挖模拟时，目前考虑材料细观结构特征的计算模型主要分为两类：基于非连续介质力学的离散元模型[16-18]和基于连续介质力学的细观有限元模型[1,2,19]。离散元模型将砂卵石地层材料视为圆形或球形颗粒集合体，并结合单元体试验结果标定土体颗粒接触参数。该类模型可以再现围岩灾害发生过程，但其计算代价昂贵，计算规模有限。细观有限元模型将砂卵石地层介质视为由土体基质和块石等形成的多相材料，分别赋予各自的材料属性，并将块石的细观结构特征，如块石形状、尺寸、倾角以及空间分布形式等都进行显式地表征[1-2]。该类模型能准确地反应围岩变形发展和应力传递路径的细观机制，特别是大块石对变形模式的影响。

本文从细观角度出发，考虑地层材料的非均质性，建立砂卵石地层细观有限元分析模型，开展不同卵石倾角下的隧道开挖模拟。通过对比围岩应力计算值与理论值，对细观分析模型的适用性进行验证。结合PECK 经验公式，验证高斯分布函数在描述砂卵石地层和地中沉降曲线的适用性；在此基础上，分析了沉降最大值和沉降槽宽度系数随卵石倾角的变化规律；提出地层沉降最大值、沉降槽宽度系数沿深度方向变化规律的数学描述公式，并对隧道断面变形模式进行了分析。

1 PECK 经验公式

基于隧道工程实践经验和大量监测数据，PECK[20]于1969 年提出地表沉降曲线可以采用高斯分布函数进行描述。地表沉降S可表述为：

式中：Smax为地表沉降最大值；i为地表沉降槽宽度系数，其物理含义为沉降曲线反弯点至隧道中心线的水平距离；x为地表任一位置至隧道中心线的水平距离。

PECK 经验公式现已成为评价隧道开挖时地表变形规律的基本方法。通过现场监测和离心机模型试验，MAIR 等[21]发现地中沉降曲线也可采用高斯分布函数进行描述，并将式(1)拓展至如下形式：

式中：Smax(z)为z深度处的地中沉降最大值；i(z)为z深度处的沉降槽宽度系数；z为地表到地中某位置处的竖向距离，如图1 所示。从式(2)可以看出，当采用高斯函数描述地中沉降曲线时，存在地中沉降最大值和地中沉降槽宽度系数两个关键参数。

图1 地表与地层沉降槽(横断面)Fig. 1 Ground settlement trough in transverse section

已有研究[22]表明：在砂卵石地层中，PECK公式的适用性与地层中的块石含量密切相关，界限含石量通常在50%左右。当地层中的块石含量小于50%时，地表和地层沉降曲线仍可采用高斯函数进行描述；当含石量大于50%时，地表沉降槽形态呈“V”型分布，采用高斯曲线描述将存在一定的误差。

2 砂卵石地层隧道开挖细观分析模型

2.1 二维细观有限元模型

考虑到计算量的限制，采用砂卵石地层二维细观有限元模型开展隧道开挖模拟。图2 所示为建立的砂卵石地层隧道开挖细观分析模型，其计算区域长度为40 m、高度为30 m。隧道开挖断面为圆形，直径D为6.0 m、洞顶埋深H为12.0 m。隧道开挖方式采用全断面一步开挖，不设衬砌支护，从而可以最大程度地体现出开挖作用对围岩的影响。

图2 砂卵石地层隧道开挖细观分析模型 /mFig. 2 Mesoscale model of sandy cobble stratum

在全部计算区域内，均考虑地层材料的非均质特性，将材料视为由土体与块石组成的两相介质，暂不考虑土体与卵石间界面过渡区的影响。为方便模型的建立，卵石形态均采用椭圆形描述，并采用“取与放”的方法[23-25]将卵石均匀地放置于计算区域中。椭圆形卵石具有长轴长度、扁平度和长轴倾角等细观结构因素。为避免细观因素间的耦合影响，在分析长轴倾角对隧道开挖围岩反应的影响时，将卵石扁平度及卵石长轴长度设为定值。监测路径Path-1 沿竖直方向分布，位于模型的纵向对称线上，如图2 所示。

2.2 工况设计

结合实际工程，研究者对卵石的细观结构特征进行了统计分析。通过对西安地铁2 号线沿线卵石进行统计，王振峰[4]得到卵石扁平度的平均值为1.75，分布范围为1.0～2.50。杨期祥[26]对成都光华公园站的卵石特征进行了分析，结果表明：卵石扁平度的平均值为2.35，分布范围集中在1.30～2.60。GAO 等[10]对成都地铁1 号线沿线卵石进行了统计，结果表明：卵石扁平度的平均值为1.50。鉴于此，本文选取卵石的扁平度为1.50。

在卵石尺寸方面，统计结果表明，对于北京市砂卵石地层，颗粒粒径由北向南逐渐增大。其中，北京地铁7 号线达官营站至湾达站区间，颗粒尺寸分布范围是20 mm～500 mm，超过200 mm的漂石含量约为15%～45%[27]；地铁9 号线局部区域，颗粒尺寸大于200 mm 的漂石含量大于20%[26]；地铁10 号线西南段，颗粒尺寸集中在20 mm～300 mm，最大粒径为1630 mm，其中100 mm～200 mm 的颗粒含量约占40%[28]。成都市也属于典型的砂卵石地层。其中，成都地铁1 号、2 号线工程，颗粒尺寸集中在20 mm～80 mm，局部区域有超过300 mm 的漂石[7]；地铁四号线西延线工程，颗粒尺寸分布范围为20 mm～660 mm，超过200 mm的漂石含量为4%～20%[6]。对于兰州市砂卵石地层，大粒径颗粒尺寸集中在150 mm～400 mm，最大粒径为720 mm[28]。

通过对Franciscan 的土石混合体进行研究，LINDQUIST[29]和MEDLEY[30]发现，对于隧道工程，当地层中块石尺寸小于0.05 倍隧道直径时，可将地层材料视为均匀介质；而当块石尺寸大于0.05 倍隧道直径时，需考虑块石细观结构特征的影响，视为非均质材料。也就是说，对于直径为6.0 m 的隧道，当地层中块石粒径大约300 mm时，就需考虑其在地层中的细观分布特征。鉴于此，本文选取的卵石含量cr为45%，卵石长轴尺寸为300 mm，则卵石短轴的尺寸为200 mm。网格尺寸取为0.10 m，该模型单元数为611 340，节点数为306 371。

对于所有的细观分析模型，卵石扁平度、卵石含量和长轴长度均保持不变，仅改变长轴倾角θ。在同一个模型中，所有的卵石均采用同一长轴倾角。对于隧道、地下管廊等长线型工程，其纵向跨度较大，穿越地层中的卵石倾角将随冲积方向和沉积作用程度的不同而发生变化[9]。因此，本文选取的长轴倾角θ 变化范围为[0°, 90°]，以10°为单位变化。

2.3 细观组分材料参数

由于块石拉/压强度较大，在隧道开挖荷载作用下通常不发生破坏，可视为弹性材料。土体基质采用经典的Mohr-Coulomb 弹塑性本构模型来描述其力学性质和变形行为。砂卵石地层细观组分材料参数如表1 所示。细观组分相关参数是由杜修力等[31]结合中型三轴试验标定得出的，详细的标定和验证过程见文献[31]。

表1 砂卵石地层细观组分材料参数[31]Table 1 Material parameters for meso-constituents

2.4 边界条件

砂卵石地层隧道开挖细观分析模型的边界条件为上表面自由，底部设置水平与竖向位移约束，侧面设置法向水平位移约束。

2.5 模拟过程

模拟过程共包含2 个分析步：地应力平衡分析步、隧道开挖分析步。地应力平衡分析步中，仅考虑土体自重应力的作用，通过给定有限元模型地表和底部的竖向应力值及其对应的坐标值，形成初始应力场。计算模型中其他部位的竖向应力值则根据地表和底部的差值得到，水平应力通过竖向应力乘以侧压力系数获取，侧压力系数K0=0.5。隧道开挖分析步中，采用“生死单元”的方法移除相应于隧道断面位置处的单元，来实现隧道的开挖模拟。

3 围岩应力结果与分析

3.1 计算值与设定值的对比

地层初始应力场是隧道开挖计算的初始条件。对于初始应力，将细观分析模型计算值与理论值进行对比，对细观分析模型的可靠性进行判断。

在本文研究中，初始应力理论值的计算公式为：

图3 为不同卵石倾角下，监测路径Path-1 竖向应力和水平应力计算值与理论值的对比。从图可见，不同卵石倾角下围岩竖向应力与理论值均拟合较好，其值随地层深度线性增加。对于围岩水平应力，模拟值与理论值也吻合较好，且其值等于竖向应力的0.50 倍，与设定值保持一致。对比结果表明，细观有限元模型能准确地描述地层中的围岩应力分布。

图3 初始状态下Path-1 应力计算值和理论值的对比Fig. 3 Stress comparison of simulation results and theoretical values in initial state for monitoring path

3.2 竖向应力云图

隧道开挖前后，围岩竖向应力云图如图4 所示。从图4 中可以看出，在隧道开挖前，竖向应力沿地层深度不断增大。当隧道开挖后，竖向应力在断面周围发生显著变化。

将不同倾角下的围岩竖向应力云图进行对比，可以发现，当卵石倾角从0°增加到90°时，拱腰处的竖向应力逐渐增大，且竖向应力增大的区域是不断扩大的，如图4 中拱腰处蓝色区域所示。当卵石倾角为0°和90°时，竖向应力分布表现出对称性，当卵石倾角为10°～80°时，竖向应力呈现出不同程度的非对称性，且当卵石倾角为40°～60°时，非对称性最为明显。竖向应力云图沿顺时针方向发生偏转，向模型对称轴右侧区域倾斜。

由以上分析可知，当地层中卵石倾角发生变化时，围岩竖向应力大小和分布范围均发生显著变化。对于由卵石和土体基质组成的砂卵石地层，荷载由卵石骨架和土体基质共同承担。当卵石倾角由0°逐渐变化至90°时，其在水平方向上荷载传递面积逐渐减小，而在竖直方向上荷载传递面积逐渐增大。因而，在相同的开挖荷载作用下，拱腰处的应力将逐渐增大，围岩应力分布会呈现出显著差别。

3.3 塑性区

塑性区为围岩土体应力达到屈服状态的区域，可作为分析围岩扰动区的参考。隧道开挖后，围岩塑性区分布如图5 所示。从图5 中可以看出，与均质地层类似，在非均质砂卵石地层中，围岩塑性区也主要分布于开挖断面周围。隧道开挖后，围岩塑性区呈蝴蝶状，分布于隧道拱肩、拱腰及拱脚位置。

图5 不同卵石倾角下的围岩塑性区Fig. 5 Plastic zones under different rock orientation angles

将不同倾角下的围岩塑性区进行对比，可以发现，当卵石倾角为0°和90°时，围岩塑性区呈对称分布，对称轴为过隧道中心的竖直线。当卵石倾角为0°时，拱肩、拱腰、拱脚处的分布范围分别为3.75R、1.75R和1.75R；当卵石倾角为90°时，相应的分布范围分别是4.0R、1.75R和1.75R。当卵石倾角在10°～80°时，围岩塑性区呈非对称性分布，且当卵石倾角为40°时，非对称特征最为显著。围岩塑性区沿逆时针方向发生偏转，向模型对称轴左侧倾斜，左拱肩处的塑性区面积大于右拱肩，而左拱脚处的塑性区面积小于右拱脚，拱腰处的塑性区面积未发生显著变化。当卵石倾角为40°时，左拱肩处塑性区面积约为4.0R×3.0R，右拱肩处为4.0R×2.0R，左拱脚处约为π ×(0.75R)2，右拱脚处约为π ×(1.25R)2。

4 围岩变形结果与分析

4.1 变形云图

地层变形是地下工程施工模拟的关键问题，也是对工程环境影响评价的基础。隧道开挖前，通常认为地层中围岩变形为零。隧道开挖后，地层中的围岩变形发生显著变化。

4.1.1 竖向位移

图6 为不同卵石倾角时的围岩竖向位移等值线图。从图6 中可以看出，隧道开挖后，洞顶围岩表现为沉降变形，洞底围岩表现为隆起变形。将不同倾角下的围岩竖向位移进行对比，可以看出围岩竖向位移的分布形态和大小因卵石倾角的不同而存在差异。当卵石倾角为0°和90°时，围岩竖向变形呈对称分布，对称轴为过隧道中心的竖直线。当卵石倾角在10°～80°时，围岩竖向变形呈现出非对称性分布的现象。在本文给定条件下，当卵石倾角在0°～90°变化时，洞顶沉降值在156.87 mm～192.94 mm 范围内变化，洞底回弹值在84.76 mm～105.65 mm 范围内变化。

图6 不同卵石倾角下的竖向位移等值线图Fig. 6 Vertical displacement contour plots under different rock orientation angles

4.1.2 水平位移

图7 为不同卵石倾角时的围岩水平位移等值线图。从图7 中可以看出，隧道开挖后，卵石倾角的变化对水平位移的分布形态和大小均产生影响。在本文给定条件下，当卵石倾角在0°～90°变化时，左拱腰处水平位移值在67.74 mm～148.83 mm范围内变化，右拱腰处水平位移值在67.08 mm～146.88 mm 范围内变化。

图7 不同卵石倾角下的水平位移等值线图Fig. 7 Horizontal displacement contour plots under different rock orientation angles

4.2 地表沉降槽

图8 为不同卵石倾角下的地表沉降曲线。从图8 中可以看出，当地层中的卵石倾角发生变化时，地表沉降值呈现出显著差异，沉降最大值和沉降槽宽度均有所不同。同时，沉降最大值的位置也会偏离模型竖向对称轴。将不同倾角下的地表沉降曲线采用PECK 高斯公式进行拟合，如图8所示。从图8 中可以看出，不同卵石倾角的地表沉降值和拟合曲线吻合较好，相关系数均达到0.99。分析结果表明：当块石含量为45%时，PECK 公式可以较好地描述砂卵石地层不同卵石倾角下的隧道开挖地表沉降特征。

图8 不同卵石倾角下的地表沉降曲线Fig. 8 Ground surface settlements under different rock orientation angles

提取地表沉降曲线上的最大值，并根据拟合曲线得到沉降槽宽度系数，分别绘制沉降最大值、沉降槽宽度系数随卵石倾角的变化图，如图9(a)和图9(b)所示。将沉降最大值位置偏离隧道中心的水平距离定为δ，绘制δ 随倾角的变化图，如图9(c)所示。

图9 关键参数随卵石倾角的变化Fig. 9 Variation of key parameters with rock orientation angle

从图9(a)可以看出，当卵石倾角从0°变化到90°时，地表沉降最大值Smax(0)整体上呈现出线性增加的趋势。当卵石含量为45%时，在本文给定条件下，地表沉降最大值由88.2 mm 增加至129.2 mm。该变化规律相当于卵石倾角每增大10°，沉降最大值增加4.56 mm，增幅比例为5.17%。图9(b)为地表沉降槽宽度系数i(0)随倾角的变化规律，从图中可以看出，当卵石倾角从0°增加到90°时，地表沉降槽宽度系数呈现出先增大后减小的规律，呈下开口的抛物线型变化形式。在卵石倾角为30°时，其值达到最大值；在卵石倾角为90°时，其值达到最小值。当卵石含量为45%时，地表沉降槽宽度系数在4.71 m～6.24 m 范围内变化。图9(c)为不同卵石倾角下，地表沉降最大值偏移隧道中心的水平距离δ。从图中可以看出，当卵石倾角从0°增加到90°时，δ 的变化形式呈现出先减小后增大的变化规律，呈上开口的抛物线形式。当卵石含量为45%时，偏移距离δ 在-0.45 m～0.07 m 范围内变化。相对隧道直径，偏移量δ 均小于7.5%D，可忽略不计，仍可认为沉降最大值位于模型纵向对称线上。

对于平面应变条件，在重力荷载作用下，地层初始应力状态以竖向应力为大主应力，水平应力为小主应力。隧道开挖后，不同位置处的围岩应力发展路径不同，因应力重分布而产生的承压、卸载状态也不同。其中，在开挖断面周围，拱顶和洞底处的围岩土体主应力方向发生旋转，以水平应力为大主应力，竖向应力为小主应力。研究表明[32]：对于此类应力状态，当颗粒倾角或沉积方向从0°增大至90°时，材料内摩擦角呈现出单调递减的规律。因此，地表沉降值将呈现出逐渐增加的现象。

4.3 地中沉降槽

地中变形决定了隧道开挖对地下结构的影响程度，是围岩变形分析的关键问题之一。图10 所示卵石含量为45%、卵石倾角为0°时的地中沉降曲线，地层深度分别为0.0 m、3.0 m、6.0 m、9.0 m 和12.0 m。由图10 可见，随着深度变化，地层沉降最大值和沉降槽宽度均发生变化。采用Peck 高斯公式对地表沉降曲线进行拟合，可以看出，不同深度处的地层沉降曲线和拟合曲线均吻合较好，相关系数均达在0.98 以上。分析结果表明：当块石含量为45%时，高斯公式可以很好地描述砂卵石地层不同深度处的地层沉降特征。

图10 卵石倾角为0°时的地层沉降曲线Fig. 10 Subsurface settlement curves at rock orientation angle of 0°

从式(2)可以看出，Smax(z)和i(z)是地中沉降变形曲线的两个关键参数，理清Smax(z)和i(z)随深度z以及卵石倾角θ 的变化规律，是确定整个地层沉降变形的基础。

4.3.1Smax(z)的变化规律及数学描述方法

提取路径Path-1 上的竖向位移值Smax(z)，并将其整理在(z/z0)-(Smax(z)/Smax(z0))的坐标系中。Smax(z0)为深度z0处的竖向沉降最大值，即洞顶竖向位移。图11 为不同卵石倾角时，Smax(z)/Smax(z0)随深度的变化规律。由图11 可见，其变化具有以下特征：

图11 不同卵石倾角Smax(z)模拟值与式(4)预测值的对比Fig. 11 Comparison between simulated Smax(z) and predicted curves of Eq.(4) under different rock orientation angles

1) 在地表处(z=0)，Smax(z)/Smax(z0)等于Smax(0)/Smax(z0)；

2) 在洞顶处(z=z0)，Smax(z)/Smax(z0)=1.0；

3) 当地层深度z从0 增加至z0时，Smax(z)/Smax(z0) 呈现出非线性增加的趋势；

4) 在不同卵石倾角下，Smax(z)/Smax(z0)与深度z均表现出以上的变化特征。

基于以上认识，本文提出以下幂函数来描述Smax(z)/Smax(z0)与z/z0间的变化规律：

式中，n为地中沉降最大值影响参数。

采用式(4)对Smax(z)进行拟合，如图11 所示。从图11 中可以看出，当不同卵石倾角下，细观模型模拟值与式(4)计算值吻合较好，相关系数均大于0.98。当卵石倾角从0°增加到90°时，n值呈现出先减小后增大的趋势，当卵石倾角为30°时取得最小值，如图12 所示。当卵石含量为45%时，在本文给定条件下，n值在0.59～0.81 范围内变化。

图12 不同卵石倾角下的沉降最大值影响参数nFig. 12 n value under different rock orientation angles

4.3.2i(z)的变化规律及数学描述方法

根据地中沉降槽拟合曲线，得到不同深度处沉降槽宽度系数i(z)，并将其整理在(z/(z0+R))-(i(z)/i(0))的坐标系中。图12 为不同卵石倾角时，i(z)/i(0)随深度的变化规律。由图可见，其变化具有以下特征：

1) 在地表处(z=0)，i(z)/i(0)=1；

2) 在洞顶处(z=z0),i(z)/i(0)=i(z0)/i(0)；

3) 当地层深度z从0 增加至z0时，i(z)/i(0)呈现出非线性减小的趋势；

4) 在不同卵石倾角下，i(z)/i(0)与深度z均表现出以上的变化特征。

基于以上认识，本文提出以下幂函数来描述i(z)/i(0)与z/(z0+R)间的变化规律：

式中，f为沉降槽宽度系数影响参数。

采用式(5)对i(z)进行拟合，如图13 所示。由图13 可见，不同卵石倾角下，细观模型模拟值与式(4)计算值吻合较好，当卵石倾角为0°～70°时，相关系数均大于0.90。当卵石倾角从0°增加到90°时，f值在整体上呈现出逐渐增大的趋势，如图14 所示。当卵石含量为45%时，在本文给定条件下，f值在2.88～4.41 范围内变化。

图13 不同倾角下i(z)模拟值与式(5)预测值的对比Fig. 13 Comparison between simulated data of i(z) and predicted curves of Eq.(5) under different rock orientation angles

图14 不同倾角下的沉降槽宽度影响参数fFig. 14 f value under different rock orientation angles

4.4 隧道断面变形模式

图15 为不同卵石倾角时的隧道开挖断面变形图，变形系数统一取为2.0。图15 中圆形代表隧道开挖前断面形状，不规则线框为隧道开挖后断面形状。由图15 可见，隧道开挖后，在无支护状态下，断面向隧道中心收敛。不同于均质地层，由于卵石的存在，变形模式呈现非光滑的锯齿状。同时，从图15 中可以看出卵石倾角的变化对开挖断面变形模式具有一定影响。在卵石倾角为0°和90°时，变形形态呈对称分布，对称轴为通过隧道中心的竖直线。当卵石倾角在10°～80°间变化时，变形模式沿着顺时针方向发生偏转，形态逐渐由水平扁平状向竖直扁平状发展。