岩石弹性模量应力相关性(stress-dependent elastic moduli, SDEM)是岩石弹性模量随应力变化的性质^{[1 − 4]}。KULHAWY^[5] 对岩石弹性模量与围压的关系开展研究，提出了幂函数模型，即$E = {E_0}{\sigma _3^m}$，式中：${E_0}$为单轴条件下岩石线弹性阶段的切线模量；m介于0~1之间。BROWN 等^[1]利用砂岩三轴压缩试验测得0 MPa~80 MPa范围内不同围压下试件的应力-应变曲线，分析发现泊松比约为0.2且几乎不受围压影响；以峰值强度50%处的切线模量为岩石弹性模量E，E随围压增加(减小)而增加(减小)；建立了$E$与${\sigma _3}$修正的幂函数关系模型，即$E = {(A{\text{ + }}B{\sigma _3})^C}$。在幂函数模型中，$E$可以随${\sigma _3}$无限增加的特征与实际不符。为此，SANTARELLI和BROWN^[6]提出了一种$E$与${\sigma _3}$的负指数关系模型，即$E{\text{ = }}{E_\infty } - ({E_\infty } - {E_0}) {{\mathrm{e}}^{ - A{\sigma _3}}}$，式中，${E_\infty }$、${E_0}$分别为当${\sigma _3} \to \infty$、${\sigma _3} = 0$时的$E$。显然，当${\sigma _3} \to \infty$时，${E_\infty }$是一个常数，它消除了$E$随${\sigma _3}$无限增加的情况，是一个较为合理的数学模型。另外，还有不少学者建立并使用了其他形式的数学模型，例如：$E{\text{ = }}A{\sigma _3} + B$^{[7 − 8]}、$E{\text{ = }}A\sigma _3^2 + B{\sigma _3} + C$^{[9 − 14]}、$E{\text{ = }}A\ln {\sigma _3} + B$^[15]、$E = {A / {(1 + B{{\mathrm{e}}^{ - C{\sigma _3}}})}}$^[16] 等。

众所周知，隧洞开挖后，径向卸荷，在径向上(即垂直于洞壁→围岩内部的方向)${\sigma _3}$逐渐增加。根据岩石的SDEM，随着${\sigma _3}$径向递增，围岩弹性模量也将径向递增。岩石SDEM影响的显著性分2种情况：① 当隧洞深埋时，开挖卸荷作用较强，围岩弹性模量受围压影响明显。因为，此时围压的径向变化量较大，导致围岩弹性模量的径向变化量也较大，弹性模量不再为常数，围岩为径向非均质体。② 当隧洞浅埋时，开挖卸荷作用较小，围岩弹性模量虽然受围压影响，但影响较小。此时，弹性模量可以近似为常数，围岩可近似为均质体。可见，在计算深埋隧洞围岩应力时，隧洞围岩不再是均质体，而是径向非均质体，岩石的SDEM是不可忽略的因素。

根据隧洞断面和地应力条件，考虑SDEM的深埋隧洞围岩应力计算问题，可分为静水压力下圆形隧洞、非静水压力下非圆形隧洞2种情况。其中，前者最为基础，可简化为“考虑SDEM的厚壁筒”，存在理论解。对此，已有研究大致可分2类：弹性应力解析解和弹塑性应力解析解。① 对于弹塑性应力解析解，一些学者将围岩视为弹塑性体，通过建立弹性模量与极半径的函数关系并结合统一强度理论^[17]，或者考虑围压依赖的应变软化行为和Hoek–Brown强度准则^[18]，再或者通过建立弹性模量与围压的函数关系并结合摩尔-库仑强度准则^[19]，分别推导了圆形隧道围岩弹塑性应力半解析解。这些弹塑性应力的半解析解所用强度准则多与剪应力强度有关，而深部围岩破坏(板裂、分区破裂、岩爆等)多为拉张性质的破坏，难以用来解释。② 对于弹性应力解析解，BROWN等^[1]将岩石视为线弹性体，根据应变相容方程、平衡微分方程、物理方程推导了考虑SDEM的厚壁筒切向、径向应力的关系式，以及径向应力关于极半径的积分方程，并利用Guass-Legendre数值积分方法求出了厚壁筒弹性应力半解析解。为便于对比，将该方法称为“积分方程法”。LIONÇO 等^[20]利用积分方程法研究了不同埋深圆形隧洞的围岩应力分布特征，发现当埋深较小时，围岩最大应力在洞壁，而当埋深超过到一定数值，则最大拉应力将出现在洞壁附近的围岩内部。然而，① 积分方程法所得半解析解形式相对复杂(见附录A)，应用较为不便，而且缺乏验证。而根据几何方程、平衡方程和物理方程，结合微分方程数值求解方法或能推出别的半解析解，这种方法也值得尝试。② 以往弹性应力半解析解，多用于围岩主应力分析，与格里菲斯强度理论(Griffith's strength theory)结合的围岩拉应力分布及其影响因素(与深部围岩板裂联系密切)仍需继续研究。

因此，本文根据几何方程、平衡方程和物理方程，联合微分方程数值求解方法，提出了一种考虑SDEM的厚壁筒应力半解析解，并结合格里菲斯强度理论，对厚壁筒拉应力及其影响因素进行了分析。结果有望对认识深部围岩破坏机理提供参考。

1   考虑SDEM的厚壁筒应力半解析解

问题描述：在处于静水压力下的无限大弹性岩体中，开挖圆形隧洞，水平地应力和垂直地应力均为p，隧洞半径为r₀，圆筒内壁受有均布支护压力p₀，考虑岩石SDEM的影响，求围岩应力分布。

问题分析：静水压力下圆形隧洞的受力情况，可简化为无限大厚壁筒平面应变模型。根据对称性，其围岩剪应力为0，弹性模量E、径向应力σ_r和切向应力σ_θ只随极半径r变化。上述问题，可等价为“考虑SDEM的已知边界条件下无限大厚壁筒应力分布问题”。边界条件为：

岩石弹性模量与σ₃的关系^{[1, 19]}：

式中：a、b和m为模型系数，可由常规三轴试验确定；${\sigma _3}$为正数，并且只能为压应力。

1.1   微分方程法

根据厚壁筒的几何微分方程、平衡微分方程和物理方程推导出了考虑SDEM的厚壁筒应力半解析解的微分方程法，下文简称为微分方程法(differential equation methods, DEM)。

求解思路：首先，列出厚壁筒的基本微分方程，包括几何方程、静力平衡微分方程和平面应变模型的物理方程；其次，通过这些方程联立、变换，消去位移量，建立σ_r和σ_θ的关系式，以及σ_r关于r的微分方程；最后，通过求解微分方程，得到σ_r及σ_θ的解析解或半解析解。

1.1.1   基本微分方程

考虑SDEM的厚壁筒依然满足弹性力学中各向同性、连续、小变形的假设，因此它的几何方程、平衡方程与均质厚壁筒的相同，分别为：

受到SDEM的影响，厚壁筒已不再是均质体，而是弹性模量随应力变化的非均质体，但仍然是线弹性体，服从胡克定律。由于在静水压力条件下σ₃ = σ_r，在式(3)中可将σ₃替换为σ_r；在均质厚壁筒物理方程中，将E替换成${{E}}({\sigma _r})$得到考虑SDEM厚壁筒的物理方程，即：

或者：

1.1.2   微分方程变换

根据基本微分方程，通过方程联立、化简，求解径向位移，进而建立${\sigma _r}$与${\sigma _\theta }$的关系式。

1) 径向位移求解

首先，联立式(4)、式(6)，得出位移-应力关系：

其次，将式(8)的两式相减得到：

再次，将式(9)代入(5)，消去${\sigma _\theta } - {\sigma _r}$，得到：

随后，为了便于开展积分运算，将式(10)变形，将含有${u_r}$、${\sigma _r}$的项分别移到等号两边，得到：

然后，将式(11)等号两端同时对$r$求积分：

由于：

则：

式中，${C_0}$为积分常数。将式(14)代入式(12)得到：

由于式(15)的等号右侧有积分项，等号左侧的积分常数${C_0}$可以消去，整理后得到：

令${{f}}({\sigma _r}) = \int {{1 \mathord{\left/ {\vphantom {1 {{{E}}({\sigma _r})}}} \right. } {{{E}}({\sigma _r})}}{\text{d}}{\sigma _r}}$，则式(16)可变换为：

2) σ_r和σ_θ的关系式建立

将式(17)代入式(4)可得：

将式(18)代入式(7)，消去${\varepsilon _r}$和${\varepsilon _\theta }$得到：

由于：

将式(20)代入式(19)，化简得σ_r、σ_θ表达式：

式中：

由式(21)~式(23)联合推知σ_r和σ_θ的关系式：

虽然本文推导过程未用到相容方程，但推导结果式(24)与文献[1]中根据应变相容方程、平衡方程和物理方程推导出的结果是一样的，见附录A式(A7)。

将式(23)代入式(21)得到：

将式(25)中${{{\text{d}}{\sigma _r}} \mathord{\left/ {\vphantom {{{\text{d}}{\sigma _r}} {{\text{d}}r}}} \right. } {{\text{d}}r}}$移到等号左边，变形化简得到${\sigma _r}$关于r的微分方程，即：

理论上，解微分方程式(26)即可得${\sigma _r}$关于r的函数关系式，再将${\sigma _r}$代入式(22)，即可得到${\sigma _\theta }$。而式(26)求解难度主要取决于${{E}}({\sigma _r})$表达式的复杂程度：当表达式简单时，存在理论解；当它复杂时，理论求解困难，可用数值计算方法求数值解。

1.1.3   微分方程求解

如果厚壁筒是均质的，则它的弹性模量是一个常数，即${{E}}({\sigma _r}) = E$。那么根据式(23)计算可得${{f}}({\sigma _r}) = {{\left( {A + {\sigma _r}} \right)} \mathord{\left/ {\vphantom {{\left( {A + {\sigma _r}} \right)} E}} \right. } E}$(其中，$A$为不定积分的常数项)，微分方程式(26)存在解析解。该条件下，式(26)可化简为：

结合应力边界条件式(1)、式(2)，令${p_0} = 0$，根据高等数学，求解一阶非齐次微分方程式(27)得：

将式(28)、式(29)代入式(21)得到：

式(29)、式(30)即为无支护(${p_0} = 0$)条件下无限大均质厚壁筒应力分布函数。它们与弹性力学的解析解是一样的。这说明前述公式的推导过程是正确的。

如果考虑SDEM的影响，并且弹性模量与应力的关系服从式(3)，那么通过计算可以得到${{f}}({\sigma _r}) = {{\ln (a{{\mathrm{e}}^{m{\sigma _r}}} - b)} \mathord{\left/ {\vphantom {{\ln (a{{\mathrm{e}}^{m{\sigma _r}}} - b)} {am}}} \right. } {am}} + C$(其中，$C$为不定积分的常数项)。该条件下，式(26)可化简为：

由式(31)可知，该微分方程较为复杂，求解析解较为困难，但可进行数值求解。由于在微分方程式(31)中含有待定系数$C$，直接数值解会非常麻烦，需要对它提前确定。为确定$C$，将式(31)等号右侧的$r$移到等号左侧，变形为以下形式：

在式(32)的等号两端同时对$r$求极限，即：

结合边界条件式(2)${\left. {{\sigma _r}} \right|_{r \to \infty }} = p$对式(33)化简得：

由于式(34)等号的右边为常数，设它为$B$，则：

下面利用反证法来证明$B$只能等于0。假设$B \ne 0$，将B移到式(35)的等号左侧得到：

由式(36)可知，${\sigma '_r}$是${B \mathord{\left/ {\vphantom {B r}} \right. } r}$的等价无穷小。根据等价无穷小的充要条件，${\sigma '_r}$可表示为：

式中，$o({B \mathord{\left/ {\vphantom {B r}} \right. } r})$为${B \mathord{\left/ {\vphantom {B r}} \right. } r}$的高阶小量。根据高阶小量积分后仍然为高阶小量的性质(在弹性力学中，根据小变形假设构建几何方程也是利用了该性质，是忽略高阶小量的结果)，忽略高阶小量$o({B \mathord{\left/ {\vphantom {B r}} \right. } r})$，并对式(37)积分得到：

再由式(38)计算得到：

式中，D为积分常数。根据式(39)计算${\sigma _r}$的极限：

式(40)所示结果与式(2)${\left. {{\sigma _r}} \right|_{r \to \infty }} = p$的边界条件是矛盾的。因此，$B \ne 0$的假设不成立，B只能为0。根据$B = 0$，结合式(34)、式(35)解出C：

为便于对微分方程式(31)进行数值求解，令：

则式(31)可改写成：

利用四阶Runge-Kutta微分方程数值求解法，直接套用其公式，对于微分方程式(43)进行数值求解，求解公式如下：

式中：$\Delta l$为对$r$进行离散的步距；$n + 1$为离散点的数量；${r_n}$、${\left[ {{\sigma _r}} \right]_n}$分别为第$n$个点对应的$r$、${\sigma _r}$，其中，${r_n} = {r_0} + \Delta l(n - 1)$，根据边界条件式(1)，当$n = 1$时，${r_1} = {r_0}$，${\left[ {{\sigma _r}} \right]_1} = {p_0}$。利用式(44)可解出各点的径向应力${\sigma _r}$的数值解，解的集合表示为：

将式(45)代入式(24)，求出${\sigma _\theta }$得：

该“半解析解”与“解析解”不同。解析解可以确定${\sigma _r}$、${\sigma _\theta }$与$r$的函数表达式，但是半解析解却不能。半解析解的形式是点的集合，形式如式(45)，它解决的是当${\sigma _r}$与$r$的函数关系难以确定时，用数值计算的方法得到近似理论解的问题，而且这种解的精度又高于ANSYS等数值模拟软件的计算结果。一般地，对于半解析解，它只能根据给定的$r$，计算出对应的${\sigma _r}$、${\sigma _\theta }$数值，无法显示出${\sigma _r}$、${\sigma _\theta }$与$r$的数学函数表达式及各物理参数与应力的函数关系式。因此，即使得到半解析解，各因素与拉应力的函数关系也是难以建立的。这是许多半解析解和数值解均会遇到的共性问题。

1.2   微分方程法的验证

以上微分方程法所求数值解是否正确，还需要进行验证。一般地，就同一个问题，采用两种方法计算，如结果一致，则能够说明两种方法都是可靠的；否则，则说明至少有一种方法是存在问题的。

1.2.1   与积分方程法计算结果的对比

微分方程法的验证主要基于“积分方程法”(推导过程见附录A)。积分方程法来源于文献[1]，是另外一种通过用求解积分方程来获得考虑SDEM的厚壁筒应力半解析解的方法。对1.1节所提问题，分别采用微分方程法、积分方程法求解，通过对比计算结果来对微分方程法进行验证。

两种方法所得应力结果都是数值解，在计算过程中，对涉及到的厚壁筒物理参数进行量化是必要的。物理参数主要包括：r₀、p、p₀、μ、a、b、m 共7个，它们的取值如表1所示。其中，式(3)的模型系数$a = 59.96$ GPa、$b = 50.31$ GPa、$m = 0.044$来源于文献[19]。利用两种方法做出${\sigma _r}$、${\sigma _\theta }$随$r$的变化曲线，如图1所示，由图可以看出，两种计算方法所得曲线重合。这说明，两种方法均可以用于求解考虑SDEM的厚壁筒的应力分布，求解结果都是较为准确的。

表  1  考虑SDEM的厚壁筒物理参数取值表

Table  1.  physical parameters of thick-walled cylinder considering SDEM

内径 r0/m	外壁压 应力p/MPa	内壁压 应力p0/MPa	泊松比 μ	式(3)参数
				a/GPa	b/GPa	m
4	54	0	0.25	59.96	50.31	0.044

下载: 导出CSV 

| 显示表格

图  1  积分方程法与微分方程法计算结果对比

Figure  1.  Comparison between the integral equation method and the differential equation method

下载: 全尺寸图片 幻灯片

1.2.2   与积分方程法的区别

微分方程法和积分方程法都是以厚壁筒外径是无限大的假设下推导的，符合地下岩体中开挖圆形巷道的条件，计算精度均较高，可相互验证。它们也存在以下区别。

1)求解思路不同

微分方程法根据几何方程、平衡方程和物理方程(未用到相容方程)，构建${\sigma _r}$关于$r$的微分方程，并结合微分方程数值求解方法(比如龙格-库塔法)求解该微分方程，来获得${\sigma _r}$关于$r$的数值解(见式(45))；积分方程法根据相容方程、平衡方程和物理方程(未用到几何方程)，构建${\sigma _r}$关于$r$积分方程，并结合数值积分方法(比如高斯-勒让德法)求解该积分方程，来获得${\sigma _r}$关于$r$的数值解。

2)求解精度的影响因素不同

微分方程法的求解精度，受微分方程数值求解方法的影响。一般地，龙格-库塔法的求解精度高比欧拉法的要高，龙格-库塔法的阶数越高，求解精度就越高，四阶龙格-库塔法可满足要求。积分方程法的求解精度，受积分方程数值求解方法的影响。一般地，高斯-勒让德法的求解精度高比梯形积分法的要高，高斯-勒让德法的节点数越多，求解精度就越高，六节点高斯-勒让德法可满足要求。

3) ${\sigma _r}$-r曲线的绘制方法不同

微分方程法在绘制${\sigma _r}$-r曲线时，通过对一定范围内的极半径r离散化，将各点的r坐标代入数值求解式(44)求解对应${\sigma _r}$，此时r为自变量，${\sigma _r}$为因变量，与常规的思路是一致的。而积分方程法在绘制${\sigma _r}$-r曲线时，通过对一定范围内的${\sigma _r}$离散化，将各点的${\sigma _r}$坐标代入数值求解公式(附录A式(A16))求解对应的r，此时r为因变量，${\sigma _r}$为自变量，与常规的思路相反。相对而言，前者使用更加方便。

4)微分方程法的数值求解过程相对更简单

与积分方程法相比，微分方程法的计算过程相对更加简便。因为积分方程法的数值求解，需要提前确定高斯节点和积分系数，这个过程相对复杂，需要用到更专业的数学知识，而龙格-库塔微分函数数值解法不需要，直接套用公式求解即可。

5)微分方程法的作用在于验证

在文献[1]的积分方程法的推导过程中，几处关键环节的推导不够详细(可参考附录A)，而且如何根据附录A中式(A16)得到切向和径向应力也未做出详细的说明。这容易使读者产生困惑，对积分方程法的可靠性产生质疑。而微分方程法，从推导过程到求解方法都与积分方程法存在较大的区别；它的出现不仅验证了积分方程法可靠，也证实了考虑SDEM的厚壁筒应力半解析解的可靠性，可以为其后期的应用奠定更加坚实的基础。

1.2.3   微分方程法的适用条件和局限性

微分方程法可用于求解考虑SDEM的厚壁筒应力，但并不是所有的考虑SDEM的厚壁筒的应力都可以用它来求解。一般地，厚壁筒应同时满足以下几个方面的条件。

1) 在物理性质方面

① 厚壁筒应同时具备连续、各向同性、小变形、线弹性的性质。因为这些条件是推导过程所用到的基本微分方程的内在要求。例如，几何方程式(4)就是以连续和小变形为前提推导出来的；物理方程式(6)是以各向同性和线弹性为前提得到的，仍然要求厚壁筒任意微元应力与应变的关系服从胡克定律。② 厚壁筒的弹性模量应该具有应力相关性。而且弹性模量与最大、中间主应力无关，只与最小主应力相关，是关于最小主应力的一元函数，否则式(23)中的积分运算将会变得相对复杂。另外，在SDEM的影响下，厚壁筒的弹性模量还会呈现出径向的递变性。③厚壁筒的泊松比必须为常数，处处相等，如果泊松比也是空间变化的，则该方法不适用。

2) 在几何尺寸方面

① 厚壁筒的外壁半径要足够大。因为在求解积分常数项C(见式(41))时，要进行$r \to \infty$过程的极限运算。如果厚壁筒的内径和外径相差不大，利用该方法计算，将会出现一定的误差。② 厚壁筒在轴向上要足够长。因为只有足够长，垂直于轴向的剖面应力分布问题才能转化为平面应变问题求解。不难看出，式(6)或者式(7)就是基于平面应变问题的物理方程得到的。如果厚壁筒在轴向上较短，则可简化为平面应力问题，其求解过程是类似的，只是需要将物理方程式(6)或式(7)替换为平面应力问题的物理方程。

3) 在受力条件与应力状态方面

① 厚壁筒外径和内径所受力必须为均布的压力，且厚壁筒整体处于静力平衡状态。因为平衡方程式(5)是以厚壁筒处于静水压力下的平衡状态为前提得到的^[21]。之所以要求所受力必须为压力，是因为：对于推导过程所用到的式(3)，它是弹性模量随围压增加而增加的定量表征，其内在的要求是${\sigma _3}$要为正数，并且只能为压应力。如果${\sigma _3}$为负数，将呈现为弹性模量随着${\sigma _3}$的增加而减小，这是不符合要求的。因此，推导过程默认压应力的符号为正、拉应力为负(图1即是按此符号规定下的结果)。② 但是在弹性力学、数值模拟软件中，正应力如果数值为正，则为拉应力，否则为压应力。如果按此符号规定，微分方程法依然是适用的。只是不能将为负数的p和p₀直接代入前述公式求解，需要先对它们取绝对值，再代入计算，最后还需要在计算出来的径向和切向应力前添加负号。③ 厚壁筒的任意微元的最小主应力不能出现拉应力。因为式(3)中，代入的最小主应力必须为压应力。如果计算结果中出现了最小主应力为拉应力的情况，则不适用。显然，对于静水压力下的圆形隧洞，是不存在这个问题的。但是，如果在非静水压力和非圆形隧洞的数值模拟研究中，如果用到式(3)，这一点是需要考虑的。

微分方程法与积分方程法存在相同的局限性，局限性主要体现在适用于静水压力下连续-各向同性-小变形-线弹性-弹性模量只与最小主应力相关-无限大厚壁筒的应力求解，难以用于复杂断面形式和复杂应力环境下隧洞围岩应力分析。

1.3   考虑SDEM的厚壁筒应力分布特点

为了形象说明考虑SDEM的厚壁筒应力分布特征，对厚壁筒的参数赋值，大部分参数与表1一致。但是为了与弹性力学对正应力符号的规定保持一致，规定：如果正应力为正，则该应力为拉应力，否则为压应力。因此，边界应力条件需要改为$p = - 54$ MPa。绘制考虑SDEM的厚壁筒应力分布曲线如图2所示。图中，${\sigma _\theta }$、${\sigma _r}$的数值均为负数，负号只表示应力为压应力。通过分析图中应力大小的变化趋势可以得到：① 对于考虑SDEM的厚壁筒，其${\sigma _r}$的大小随$r$的增加而单调递增；② ${\sigma _\theta }$的大小随$r$的增加呈现出了先增加后减小的趋势，存在${\sigma _\theta }$峰值；③ 以峰值为界，厚壁筒分为应力上升区(位于${\sigma _\theta }$峰值的左侧)和应力下降区(位于${\sigma _\theta }$峰值的右侧)；④ ${\sigma _\theta }$始终大于${\sigma _r}$，分别为最大、最小主应力${\sigma _1}$、${\sigma _3}$；⑤ 考虑SDEM的厚壁筒与均质厚壁筒应力分布存在显著区别，主要集中在内径附近的应力上升区。${\sigma _\theta }$极值及其对应$r$的理论表达式也许是重要的，但是由于只有数值解，${\sigma _r}$、${\sigma _\theta }$与$r$的函数关系，还难以给出。

图  2  考虑SDEM的厚壁筒应力分布曲线

Figure  2.  Radial distribution curves of stresses in the thick-walled cylinders with SDEM

下载: 全尺寸图片 幻灯片

计算出${\sigma _r}$后，将其代入式(3)，替换掉${\sigma _3}$，即可以得到厚壁筒弹性模量分布函数$E(r)$：

式中，${\sigma _r} = {\sigma _r}(r)$。

根据式(47)可知，厚壁筒弹性模量的分布函数具有以下特点：① 弹性模量$E$是关于径向应力${\sigma _r}$的函数，而径向应力${\sigma _r}$又是$r$的函数，$E$是关于$r$的复合函数；② 在内壁${\sigma _r} = 0$，${{E}}(0){\text{ = }}a - b$；当$r \to \infty$， ${\sigma _r} \to p$，厚壁筒弹性模量无限趋向于一个常数$a-b\mathrm{e}^{-mp}$，即E-r曲线存在一条水平渐近线。图3为考虑SDEM的厚壁筒弹模径向分布曲线。由图可以看出，考虑了SDEM的厚壁筒弹性模量不再是一个常数，而是随位置变化的函数；$E$随$r$的增加而单调递增，增长速率逐渐递减，$E$随$r$的增加呈现出不断趋近于一个常数的特征。

图  3  考虑SDEM的厚壁筒弹性模量径向分布曲线

Figure  3.  Radial distribution curve of elastic moduli in the thick-walled cylinder with SDEM

下载: 全尺寸图片 幻灯片

以上分析说明：SDEM的作用可导致厚壁筒弹模径向递变，并对厚壁筒应力分布产生影响，影响范围主要集中在筒壁附近。因此，在研究对象方面，考虑SDEM的围岩，既不同于常规的均质围岩，又与常规的非均质围岩(弹性模量随机变化而非定向递变)^[22]有所区别，是弹性模量空间递变的围岩。

2   考虑SDEM的厚壁筒应力半解析解的应用

围岩板裂是深部工程岩体典型的破坏形式^{[23 − 24]}，表现为完整脆性岩体被平行的裂隙切割，形成多层近似平行于开挖面的岩板^{[25 − 26]}，严重威胁着深部工程施工安全^{[27 − 28]}。例如，我国西藏自治区某大断面隧洞拱顶围岩板裂严重，钻孔凿岩台车、身着防弹衣在掌子面前装填炸药的工人，时刻承受来自拱顶岩板掉落的风险，给隧洞安全高效施工带来了巨大的挑战^[29]。这种围岩破坏形式难以用传统的“三圈”(破碎圈、破裂圈和弹性圈)理论进行解释，其诱发机理至今仍然未得到很好的解释。

如果考虑到岩石SDEM的影响，隧洞开挖后，径向卸荷，${\sigma _3}$在径向(洞壁→围岩内部方向)不断增加，导致围岩弹性模量呈现出径向递变性，进而影响围岩应力分布。SDEM可能是造成围岩板裂的重要诱发因素。又因为深部围岩板裂属于张拉破坏，其板裂特征与拉应力分布有关。因此，需要重点关注围岩拉应力，而前述微分方程法求出的是最大、最小主应力，还需要进一步结合格里菲斯强度准则，求出围岩拉应力。研究围岩拉应力的影响因素及其影响效果，对认识深部围岩板裂机理具有启示意义。

2.1   隧洞围岩拉应力分布

格里菲斯强度理论认为，玻璃、岩石等脆性固体材料内部通常含有大量微裂纹。在压应力状态下，裂纹周边也可能出现较高的拉应力，同样可导致裂纹的不稳定扩展。岩石的裂纹扩展应满足格里菲斯强度准则。根据格里菲斯强度准则，围压拉应力计算公式如下：

式中：${\sigma _1}$、${\sigma _3}$分别为最大、最小主应力(正号表示拉应力，反之为压应力；对于静水压力下的圆形隧洞，${\sigma _1}$、${\sigma _3}$均为压应力，拉应力应由式(48)的第一个公式计算)；${\sigma _{\mathrm{t}}}$为岩石的拉应力。如果${\sigma _1}$、${\sigma _3}$满足${\sigma _1} + 3{\sigma _3} {\leqslant} 0$的条件，则拉应力为$- {{{({\sigma _1} - {\sigma _3})}^2}}/ {8({\sigma _1} + {\sigma _3})}$，否则${\sigma _{\text{t}}}{\text{ = }}{\sigma _3}$。该公式反映的是岩石内部孔隙周边的最大拉应力。在此规定下，${\sigma _3}$即使为正(拉应力)也不能作为围岩的拉应力，除非满足${\sigma _1} + 3{\sigma _3} > 0$的条件。根据${\sigma _{\text{t}}}$集中区域，可以判断最容易发生张拉破坏的位置。

为便于说明围岩拉应力分布特征，做出围岩拉应力分布云图(计算所用模型参数与1.3节相同)，如图4(a)所示。由图可知，最大拉应力未在洞壁，而是位于到洞壁一定距离的围岩内部。为了便于定量分析拉应力分布的影响因素，引入两个指标，即围岩拉应力峰值${\sigma _{{\text{tmax}}}}$、${\sigma _{{\text{tmax}}}}$所在位置到洞壁的距离$d$。前者决定了围岩拉应力是否大于围岩抗拉强度，决定了是否会发生板裂；后者决定了围岩发生板裂的起裂位置和所产生岩板的厚度。为了更形象表示这两个指标，根据图4(a)，做出围岩拉应力径向分布曲线如图4(b)所示，并在图中对指标做了标注。

图  4  围岩拉应力分布

Figure  4.  Distribution of tensile stress in surrounding rocks

下载: 全尺寸图片 幻灯片

2.2   隧洞围岩拉应力分布的影响因素

虽然半解析解无法展现出${\sigma _r}$、${\sigma _\theta }$与$r$及各个参数之间的函数关系式。但通过表1不难看出，即使是半解析解，围岩应力依然受到r₀、p、p₀、μ、a、b、m 7个参数的影响，围岩拉应力也是如此。其中，${r_0}$对应隧洞半径；${p_0}$对应于隧洞支护应力；$p$对应于地应力，它与埋深h的关系为$p = \lambda h$，$\lambda$为上覆岩层的容重，可取为27 kN/m³。μ、a、b、m四个参数取决于岩性，一般以组合形式出现，每种岩性对应一个参数组(μ、a、b、m)。由于缺乏不同岩性的试验数据，这里暂不研究岩性的影响，只就$h$、${r_0}$、${p_0}$三个因素对围岩拉应力的影响进行分析。如无特殊说明，模型参数与第1.3节部分相同。

2.2.1   埋深对围岩拉应力分布的影响

以隧道埋深$h$为单一变量，$h$以0.1为步长，在0.1 km~10 km范围内取值，绘制不同埋深的隧洞围岩拉应力径向分布曲线(见图5)。表2 为不同埋深的圆形隧洞围岩${\sigma _{{\text{tmax}}}}$和$d$统计表。由图表可以看出，隧道埋深$h$对围岩最大拉应力${\sigma _{{\text{tmax}}}}$、$d$影响显著；${\sigma _{{\text{tmax}}}}$随$h$增加而单调增加；$d$随$h$的变化呈现出明显的阶段性，可分为I、II、III三个阶段。在阶段I，$h$在0 km~1 km范围内，$d = 0$，不受$h$影响；在阶段II，$h$在1.0 km~2.8 km范围内，$d$随$h$增加而增加；在阶段III，$h$在2.8 km~10 km范围内，$d$随$h$增加而减小，在$h$=2.8 km时，$d$达到峰值(需要说明的是：阶段划分节点是当前计算条件下得到的，其他条件下可能不同)。

图  5  埋深对隧洞围岩拉应力分布的影响

Figure  5.  The influence of buried depth on the tensile stress distribution of the surrounding rock of the tunnel

下载: 全尺寸图片 幻灯片

表  2  不同埋深圆形隧洞围岩σ_tmax和d统计表(部分)

Table  2.  σ_tmax and d of surrounding rocks in circular tunnel under different buried depths

埋深 h/km	拉应力峰值 σtmax/MPa	拉应力峰值所在位置到 洞壁的距离d/m
0.1	0.49	0.00
0.6	1.59	0.00
1.1	2.28	0.28
1.6	3.27	0.84
2.1	4.52	1.03
2.6	6.01	1.08

下载: 导出CSV 

| 显示表格

2.2.2   隧道半径对围岩拉应力分布的影响

以隧道半径${r_0}$为单一变量，${r_0}$以0.1 m为步长，在0.1 m~5 m范围内取值，做出不同半径的隧洞围岩拉应力径向分布曲线(见图6)。表3 为不同曲率半径下圆形隧洞围岩${\sigma _{{\text{tmax}}}}$和$d$统计表。由图表可以看出，隧道半径${r_0}$对围岩最大拉应力${\sigma _{{\text{tmax}}}}$无影响，对$d$影响明显；$d$随${r_0}$增加而增加。${r_0}$与$d$的比值为定值，说明二者关系呈现为正比例关系，比例系数为3.98(需要说明的是：该比例系数是当前计算条件下得到的，其他条件下可能不同)。

图  6  曲率半径对隧洞围岩拉应力分布的影响

Figure  6.  The influence of radius of curvature on the tensile stress distribution of the surrounding rock of the tunnel

下载: 全尺寸图片 幻灯片

表  3  不同半径下圆形隧洞围岩σ_tmax和d统计表(部分)

Table  3.  σ_tmax and d of surrounding rocks in circular tunnel under different radius of curvature

曲率半径 r0/m	拉应力峰值 σtmax/MPa	拉应力峰值所在位置到 洞壁的距离d/m	隧洞 半径r0/d
1	4.25	0.25	3.98
2	4.25	0.50	3.98
3	4.25	0.75	3.98
4	4.25	1.00	3.98
5	4.25	1.26	3.98

下载: 导出CSV 

| 显示表格

2.2.3   支护应力对围岩拉应力分布的影响

以支护应力${p_0}$为单一变量，${p_0}$以0.2 MPa为步长，在0 MPa~20 MPa范围内取值，做出不同支护应力的隧洞围岩拉应力径向分布曲线(见图7)以及${\sigma _{{\text{tmax}}}}$和$d$的统计表(见表4)。由图表可以看出，${p_0}$对围岩${\sigma _{{\text{tmax}}}}$、$d$影响显著；${\sigma _{{\text{tmax}}}}$、$d$随${p_0}$的变化呈现出阶段性特征，可分为I、II两个阶段。在阶段I，${p_0}$介于0 MPa~8.8 MPa范围内，${\sigma _{{\text{tmax}}}}$不随${p_0}$变化，$d$随${p_0}$的增加而减小；在阶段II，${p_0}$>8.8 MPa，${\sigma _{{\text{tmax}}}}$随${p_0}$增加而减小，$d$不随${p_0}$变化。这说明，支护应力的作用是诱导拉应力集中向洞壁转移。

图  7  支护应力对隧洞围岩拉应力分布的影响

Figure  7.  Influence of support stress on the tensile stress distribution of tunnel surrounding rock

下载: 全尺寸图片 幻灯片

表  4  不同支护应力隧洞围岩σ_tmax和d统计表(部分)

Table  4.  σ_tmax and d of surrounding rocks in circular tunnel under different support stresses

支护应力 p0/MPa	拉应力峰值 σtmax/MPa	拉应力峰值所在位置到 洞壁的距离d/m
0	4.25	1.00
5	4.25	0.32
10	4.23	0.00
15	3.90	0.00
20	3.32	0.00

下载: 导出CSV 

| 显示表格

2.3   对认识深部围岩板裂现象的启示

板裂发生的条件：当围岩最大拉应力${\sigma _{{\text{tmax}}}}$大于等于围岩抗拉强度${\sigma _{{\text{ts}}}}$时板裂发生，即：

根据第2.2.1节，围岩${\sigma _{{\text{tmax}}}}$随埋深增加而增加的特征可以推测：当埋深比较浅时，围岩${\sigma _{{\text{tmax}}}}$<${\sigma _{{\text{ts}}}}$，难以发生板裂；当埋深比较深时，围岩${\sigma _{{\text{tmax}}}}$>${\sigma _{{\text{ts}}}}$，板裂发生。由此可见，围岩板裂只在达到一定埋深后才会发生板裂，存在临界埋深。这与围岩在浅部不发生板裂而在深部发生的特征是一致的。而且板裂形成的岩板的厚度随埋深先增加后减小。而根据$d$随埋深增加而呈现出的“不变-增加-减小”的特点可以推测，板裂最先所形成的岩板的厚度随着埋深增加，也会是先厚后薄，是否如此，仍需要试验或现场验证。

根据低第2.2.2节，围岩${\sigma _{{\text{tmax}}}}$不随曲率半径变化、$d$随曲率半径增加而增加可以推测：如果在某一埋深(决定了应力水平)下围岩能够发生板裂，则即使曲率半径减小，仍然能发生，但曲率半径越小，板裂形成的岩板厚度越小。这对板裂试验研究具有启示意义：当开展隧洞现场试验困难时，这一特点可为在实验室采用小尺寸的试件和孔洞来确定隧洞围岩板裂发生的应力条件提供理论支撑。在实验中，矩形孔洞的曲率半径比圆形孔洞大，形成的岩板厚度也更厚，也能够印证$d$与曲率半径的正相关性。

根据第2.2.3节围岩${\sigma _{{\text{tmax}}}}$、$d$随支护应力增加呈现为“不变-减小”“减小-不变”的特征可知，在支护应力增加的过程中，先后出现了两种作用，即先促使${\sigma _{{\text{tmax}}}}$的位置向洞壁转移，再使${\sigma _{{\text{tmax}}}}$减小。由此可以推测：当支护应力较小时，支护应力的作用是前者，难以抑制板裂的发生。这就意味着，隧道支护后仍然可能支护失效，发生围岩时滞型岩爆。一般地，隧洞开挖后，喷射混凝土层的厚度较小，很难提供较大的支护应力，这可能难以消除围岩板裂。

另外，围岩板裂与分区破裂都是深部岩体破坏形式，就破坏特征而言，从洞壁到围岩内部的过程中，破裂面-岩板交替出现的为围岩板裂特征，破裂区-非破裂区交替出现的分区破裂现象^[30]，二者具有一定相似性。考虑SDEM的厚壁筒应力半解析解对研究深部围岩分区破裂化现象也有类似的启示。本文只以弹性为假设，得到的最大拉应力集中条带更多反映的是起裂的位置，而起裂后如何扩展^[31]，相关的理论计算和数值分析技术仍有待开发。

3   讨 论

3.1   考虑SDEM的围岩塑性圈确定

在不考虑SDEM随围压变化的情况下，若出现塑性区，切向应力表现为随极半径增大先增大后减小^[32]，而图1中在弹性假设前提下，考虑弹模随围压变化，随着极半径增大，切向应力先增大后减小。关于在实际工程中应该如何划分塑性区的大小，是值得继续深入研究的问题。应该分为以下两种情况讨论。

1) 在浅部(小埋深)低应力条件下，围岩塑性圈的划分可以采用三圈层理论。首先，岩石力学认为，隧洞开挖后，如果应力超过强度，围岩将发生破坏，可被分为破碎圈、破裂圈和弹性圈。弹塑性围岩的应力推导，一般要用到静力平衡方程和摩尔-库仑准则；再利用均匀厚壁筒的应力公式，推导塑性圈半径。此时，一般认为围岩均匀(只有一个弹性模量)、围岩的破坏为剪切破坏。三圈层理论在浅部低应力条件下是适用的，而且通过现场声波探测等手段也能够探测到松动圈的存在。实际上，在浅部，隧洞开挖后，卸荷作用较小，围岩弹性模量虽然受围压影响，但此时围压变化不大，弹性模量径向变化也不大，弹性模量随围压变化的影响可以忽略，弹性模量可以近似为常数，围岩仍然可视为均质体。因此，在浅部，三圈层理论的推导是合理的。

2) 在深部(深埋)高应力条件下，围岩板裂、分区破裂和岩爆等破坏形式难以用三圈层理论解释。随着开挖深度增加，地应力越来越大，在高地应力条件下，硬脆性围岩的破坏形式发生了变化，常常以板裂(层裂)、岩爆、分区破裂等破坏形式出现。板裂表现为完整的岩石被破裂面切割，形成岩板，造成围岩呈现出“岩板-破裂面交替出现”的特点，分区破裂表现为“破裂区-非破裂区交替出现” 的特点，岩爆表现为“围岩破裂成板-岩板折断弹射-形成局部的V型凹坑”的特点。它们的特点和成因都难以用三圈层理论解释，原因是：① 在围岩劈裂成板的过程中，破裂的拉张性质明显，而三圈理论在推导时应用的是摩尔-库仑准则(默认为剪切破坏)；② 在深部，隧洞开挖后，卸荷作用较强，围岩弹性模量受围压影响较大，围压变化较大，弹性模量径向变化也较大，弹性模量随围压变化的影响不可忽略，弹性模量不再是常数，围岩不再是均质体。塑性圈半径在推导的过程中，均匀厚壁筒的应力公式也许不再适用。③ 深部围岩一般分层破坏，形成近似平行于轮廓线的岩板，并且是多期发生，但平行于轮廓线的裂纹不会无限产生，而是在围岩内部一定深度停止，这个深度叫“止裂深度”。止裂深度就类似于塑性圈的半径，在实际深部工程中对支护设计具有重要参考意义。但是，止裂深度的确定仍是尚未很好解决的难题，仍需继续研究。

3.2   基于MC准则的围岩应力分布

第2.2节部分的分析是基于格里菲斯强度准则展开的，下面结合Mohr-Coulomb(MC)准则进行考虑SDEM的隧道开挖后应力分布计算，围岩的应力重分布规律分析如下：

首先，利用微分方程法可以得到隧洞围岩的最大、最小主应力${\sigma _1}$和${\sigma _3}$。其次，结合格里菲斯强度准则，可以计算围岩各点拉应力，直接通过拉应力分布可以知道拉应力集中区域，即最容易发生张拉破坏的位置。这是因为围岩的抗拉强度是一个常数。但MC准则不同，它反映的岩石抗剪强度是一个关于正应力或者最小主应力的函数，不是一个定值^[33]。利用MC准则判断围岩稳定性，不仅要计算出最大主应力${\sigma _1}$，还要计算出不同位置围岩发生剪切破坏所要达到的最大主应力${\sigma }_1^\prime$。根据MC准则得到：

式中，$\varphi$、$C$分别为内摩擦角和内聚力；按表5取值。做出${\sigma _1}$和${\sigma} _1^\prime$的曲线如图8(a)所示。图8(b)为围岩最大主应力与发生剪切破坏所需最大主应力之差${\sigma} _1^\prime$-${\sigma _1}$的径向分布曲线。在${\sigma} _1^\prime$-${\sigma _1}$越小的位置，越容易发生剪切破坏。由图8可以看出，围岩的抗剪强度${\sigma} _1^\prime$在径向逐渐增加；而最大主应力${\sigma _1}$先增加，后减小；${\sigma }_1^\prime$-${\sigma _1}$的值在洞壁r=4的位置最小。这说明最容易在洞壁发生剪切破坏。

表  5  φ、C取值表

Table  5.  Values of φ and C

岩性	指标范围			指标均值
	内摩擦角φ/(°)	内聚力C/MPa		内摩擦角φ/(°)	内聚力C/MPa
砂岩	35~50	8~40		42.5	24

下载: 导出CSV 

| 显示表格

图  8  基于MC准则的围岩应力分布

Figure  8.  Stress distribution of surrounding rock based on MC criterion

下载: 全尺寸图片 幻灯片

3.3   非静水压力下非圆形隧洞围岩弹性模量分布

微分方程法推导的前提假设为静水压力条件，然而实际上很多情况下地下岩体并非如此，考虑非静水压力条件的隧洞弹性模量空间递变特征仍是一个值得探究的问题。

一般为非静水压力下非圆形隧洞围岩，这种条件下获得围岩的解析解和半解析解，可能要用到复变函数，求解过程较为复杂和困难。这时常常采用数值模拟的方法获得围岩应力的数值解。利用数值模拟方法绘制出直墙拱形隧洞围岩弹性模量的空间分布云图(侧压系数为0.5)，如图9(a)所示。由图可知，考虑SDEM的非静水压力非圆形隧洞围岩弹性模量在空间上依然是递变的，只是此时它不再是极半径的一元函数，同时也与极角有关。另外，通过对比考虑SDEM的直墙拱形隧洞围岩拉应力空间分布云图(见图9(b))与直墙拱形隧洞围岩破坏实验照片(见图9(c))也可以看出，拉应力集中条带与围岩两帮裂纹存在明显的对应关系。关于考虑SDEM的深埋隧洞围岩应力数值模拟方法，由于涉及到基于SDEM确定围岩弹性模量空间分布、弹性模量与主应力函数关系修正等关键难题，本文未加详细介绍。

图  9  考虑SDEM的直墙拱形隧洞围岩弹模和应力云图及破坏实验照片

Figure  9.  Elastic modulus and stress cloud charts of surrounding rock of straight wall arch tunnel considering SDEM and failure experiment photos

下载: 全尺寸图片 幻灯片

4   结 论

围绕考虑SDEM的厚壁筒应力分布和隧洞围岩拉应力分布及影响因素问题开展研究，通过理论分析和数值计算，主要得到以下结论：

(1) 推导出了一种考虑SDEM的厚壁筒应力半解析解的微分方程法；经过验证，它的计算结果与积分方程的一致；适用于静水压力下连续-各向同性-小变形-线弹性-弹性模量只与最小主应力相关-无限大厚壁筒的应力求解；可以用于近似求解考虑SDEM的静水压力下圆形隧洞围岩的应力。

(2) 结合格里菲斯强度准则，分析了埋深$h$、隧洞半径${r_0}$、支护应力${p_0}$对考虑SDEM的静水压力下圆形隧洞围岩拉应力分布(由围岩最大拉应力${\sigma _{{\text{tmax}}}}$及其到洞壁的距离$d$量化)的影响。结果表明：${\sigma _{{\text{tmax}}}}$与${r_0}$无关，随$h$增加而单调增加，随${p_0}$增加呈现出“不变-下降”的特征；$d$随${r_0}$增加而正比例增加，随$h$增加呈现出“不变-增加-减小”的特征，随${p_0}$的增加呈现出“下降-不变”的特征。分析结果对认识深部围岩板裂破坏现象具有一定的启示。

附录A：积分方程法

文献[1]给出的积分方程法求考虑SDEM的厚壁筒应力的过程介绍如下。

基本方程列出如下。

相容方程：

平衡方程：

物理方程：

联立式(A1)、式(A2)可得：

将式(A3)、式(A4)代入式(A5)得到：

将式(A6)等号两端同时积分得到：

需要说明的是，在文献[1]中如何由式(A6)推出式(A7)的过程不够详细。在静水压力下圆形隧洞的最小主应力等于径向应力：

岩石弹性模量与σ₃的关系^[1]：

将式(A8)、式(A9)代入式(A7)得到：

式中：

对式(A2)等号两侧积分得到积分方程：

式中：

求解积分方程式(A12)，即可得到$t$或${\sigma _r}$关于$r$的函数关系。由于形式复杂，其解析解难以求解，文献[1]采用高斯-勒让德数值积分方法求数值解。

为了实现高斯正交，令：

则式(A12)积分得到：

式中：

根据式(A16)来绘制应力分布曲线。