渗流应力作用下层理岩体稳定性数值分析

刘远桢

(贵州省机场集团有限公司，贵阳 550005)

1 概 述

岩体强度在很大程度上主要由岩体间的结构面控制，然而层理是一种较为常见的结构面，当岩体中的层理面由于渗流有孔隙水压力作用时，这一作用将会对岩体的强度产生一定影响[1]。在岩体系统中，存在大量的地下水作用，通过物理、化学和力学等作用改变岩体的内部结构，这些作用将会给岩体施加静水压力和动水压力，这些渗流现象表现出岩体系统内的渗流场对应力场的影响作用[2]。因此，研究含层理岩体中渗流场对应力场之间的相互作用具有一定的理论意义。

此外，在进行层理岩体的变形和稳定分析时，对渗流影响的评估也非常重要，且存在较大困难。一方面的主要原因在于渗流产生渗透压力，减少层理面上的有效法向应力，从而影响岩体的变形和稳定；另一方面渗流又具有对应力和变形很强的依赖性，岩体在荷载荷载作用后，由于应力场的逐渐改变和层理法向及切向的应力逐渐增大，将会引起岩体中孔隙或层理张开度发生变化，从而岩体的渗流通道将发生变化，层理中渗流的流速和流体压力将重分布[3]，这种作用将会使变形和渗流存在较为复杂的耦合作用。岩体中的水对岩体具有物理化学和物理力学等各种作用，前者表现为水的存在降低了岩体的强度参数和变形参数，地下水的软化使岩土黏结力和摩擦力减小、变形参数改变、抗压强度降低等[4]；后者表现在地下水压力使有效正应力减小、增加渗流力及静水压力等。对于孔隙介质如表现为连续介质块体和很破碎的岩体，渗流力作为一种体积力作用在岩体上；对于层理岩体块体是相对不透水的情况下，水作为一种面力作用在结构面上，使岩体更容易沿结构面发生破坏。

鉴于此，本文基于等效连续介质模型和Louis经验公式，从宏观上建立能够表征层状边坡耦合特性的各向异性岩体渗流应力耦合模型，同时对本文的耦合方程模型进行数值计算，分析层理渗流应力耦合规律，研究结果可为类似工程应用提供一定的参考依据。

2 层理岩体渗流应力耦合模型

已有研究表明，在考虑层理作用时，假设层理岩体的变形为弹性变形；不考虑岩体损伤、破坏对渗流的影响[5]。因此，基于等效连续介质模型，通过等效弹性张量和等效渗透张量参数来表征岩体的各向异性特性，从而建立含有层理岩体的渗流应力耦合作用关系模型，具体情况如下。

2.1 应力控制方程

在分析层理岩体的受力情况时，假定为平面问题，则平衡方程为：

σij,j+Fj=0 (i,j=1,2)

(1)

式中：σij,j为j方向上的应力张量；Fj为体积力。

同时考虑孔隙水压力作用下岩体本构方程为：

(2)

几何方程为：

εij=(Ui,j+Uj,i)/2

(3)

同理由弹性力学可知，对于二维平面应变问题，其刚度矩阵[D]可由各向异性的工程弹性常数进行表示，公式如下：

(4)

(5)

式中：Ei(i=1,2)为岩体在i方向的等效弹性模量；νi(i=1,2)为岩体在i方向的等效泊松比；G12为等效剪切模量。

2.2 渗流连续性方程

在分析渗流问题时，近似看作平面二维渗流问题，分析时不考虑渗流水的可压缩性。因此得到稳定条件下的渗流连续性方程为:

(6)

2.3 耦合作用方程

研究表明，渗流应力耦合作用主要体现在表征单元体的渗透系数与应力(应变)间的相互关系，该关系是研究渗流应力耦合问题的核心[6-9]，但目前尚没有形成统一的理论。其中Louis经验公式仅含有一个耦合参数，该经验公式不仅形式简单，而且能够反映应力与渗透系数的非线性关系，是目前在工程中应用较为广泛的公式。

在平面应变条件下，研究各向异性岩体渗流应力耦合问题，三向应力状态下复杂程度较大。研究表明，层理岩体在三向应力作用下以垂直于层理的方向应力对渗流作用起主导控制[10]。考虑实际岩体的层理呈密集分布，且法向上的层理分布较为稀疏，与垂直平面的相比较，平面内结构面的法向应力成为影响层理渗流的主要因素。因此，本文基于Louis经验公式，将层理岩体的主渗透系数表示为矩阵形式，即：

(7)

式中:K11、K22分别为沿层理面和垂直于层理面方向的初始等效渗透系数；22、11分别为对应于K11、K22所在方向的法向应力，可通过应力矢量的坐标变化求得；α、β分别为Biot系数和渗流耦合作用参数，可通过室内试验或查阅相关文献确定，本文计算时，均取α=1,β=0.1。

2.4 层理岩体的稳定性判别式

为确定层理岩体在极限平衡状态下的稳定性，本节对层理岩体的稳定性判别式进行推导，具体情况如下。图1为含层理岩体单元的受力状况。

图1 计算简图

由图1可知，当层理面处于稳定状态和极限平衡状态时，由Mohr-Coulomb强度准则可得，岩体层理面上的剪应力应满足剪切强度要求，即：

|τ|≤cj+σtanφj

(8)

式中：cj、φj为层理面上的黏聚力、内摩擦角；当层理岩体处于极限平衡状态时，式(8)取等号。

式(9)、式(10)分别表示层理面上的正应力与切应力与最大主应力和最小主应力之间的应力状态关系式。

自相关矩阵Rxx的特征结构决定了LMS算法的收敛性，当特征值相差很大时，算法收敛速度很慢，同时受输入信号的功率变化的影响。

(9)

(10)

联立式(8)、式(9)、式(10)得：

σ1sinβcosβ-σ3sinβcosβ≤σ1sin2βtanφj+σ3cos2βtanφj+cj

(11)

将式(11)两边移项整理得：

σ1sinβ(cosβtanφj-sinβ)+σ3sinβ(cosβ+sinβtanβ)+cj≤0

(12)

运用三角函数运算法则简化得出无渗流状态下的稳定性判别式为：

σ1cosβsin(φj-β)+σ3sinβcos(φj-β)+cjcosφj≤0

(13)

当考虑层理面内含有孔隙水压力时，稳定性判别式为：

σ1cosβsin(φj-β)+σ3sinβcos(φj-β)+cjcosφj-pwsinφj≤0

(14)

式中的符合含义同上。当式(14)成立时，则层理岩体处于稳定状态；相反当式(14)不成立时，层理岩体将失去稳定性。

3 数值模型与验证

3.1 数值模型的建立

为研究本文公式的合理性与适用性，基于实际工程岩体芯样分布情况，本文选用含有30°层理倾角的灰岩进行验证，其岩体试样的基本参数如下：弹性模量E=5.272GPa,泊松比为0.27，重度为26.4 kN/m3，饱和重度为27 kN/m3，黏聚力为12.72 kPa，内摩擦角为49.5°，抗拉强度为7.41 kPa；层理倾角为30°，层理面上黏聚力为0 kPa，内摩擦角为40°。

基于上述层理岩体的基本物理指标，由式(7)得到模型在渗流应力场下的主渗透系数，采用Midas GTS NX数值模拟软件建立渗流应力耦合模型，得到渗透应力作用下的层理岩体模型，见图2。

图2 数值分析模型

3.2 试验结果分析

为研究渗流应力作用下，层理岩体的渗流场和应力场的分布情况，本文选取荷载应力为q=10 kN/m、q=50 kN/m的两种类似进行数值分析，得到相应的总水头、总位移以及剪应变云图，具体情况见图3-图5。

图3为层理岩体在荷载作用下的渗流场水头云图。由图3可知，在荷载作用下层理岩体的总水头场呈层状结构分布，上部水头最大，底部水头最小，由上至下逐渐减小。

图3 渗流场下的总水头

图4为荷载10和50 kN/m作用下层理岩体总位移场分布情况。由图4可知，在小荷载作用下的总位移在底部最大，在大荷载作用下的总位移在层理岩体的顶部。当荷载10 kN/m作用时，最大位移区达到总位移场的2.3%；当荷载50 kN/m作用时，最大位移区达到总位移场的0.8%。这说明随着荷载作用的增加，最大位移场的作用位置发生了变化，由底部转向顶部，在实际工程中应该注意这种变化带来的影响。

图4 不同荷载作用下层理岩体的总位移场

图5为在荷载10和50 kN/m作用下层理岩体剪切应变分布情况。由图5可知，在荷载作用下的剪应变分布主要集中在层理面的两侧，但均未出现最大剪应变区。这说明在这两种荷载作用下，层理岩体沿层理结构相对较少，从层理岩体的位移上看，该岩体处于稳定状态。

图5 不同荷载作用下层理岩体的剪应变场

为验证本文建立模型的适用性，结合式(14)与数值模拟进行拟合，计算出层理岩体的安全系数，计算结构见表1。

表1 安全系数计算结果

由表1可知，通过量取最大位移处的孔隙水压力代入公式进行计算并与数值模拟得出的安全系数值。结果表明：①30°层理岩体在渗流和不同大小应力作用下的变形位移是不一样的，随荷载作用的增大，位移逐渐由下部上移；②相对于数值模型结果，理论计算与模拟结果基本相同，理论结果偏安全。因此，本文建立的模型适用性较强，可应用于类似工程。

4 结 论

本文对含有层理岩体的渗流应力作用进行较为详细的分析，建立了相应的理论模型。采用Midas GTS NX数值软件建立数值模型进行验证，得到的主要结论如下：

1) 通过等效连续介质模型和Louis经验公式，建立了层理岩体渗流应力耦合模型，该模型能反映层理岩体在渗流应力作用下的变形特征，具有一定的理论参考价值。

2) 基于岩体极限平衡状态推导的30°层理岩体渗流应力耦合的稳定性情况判别式的计算结果与数值计算结果较好吻合，具有一定现实意义。

3) 根据实际工程层理岩体，建立了相应的数值模型，对现有模型进行验证，适用性较好。同时发现随着荷载作用增大，位移场作用发生转移，由层理岩体的底部向顶部转移。

4) 研究结果可为类似工程提供一定的参考依据，具有一定的工程应用价值。但受到计算方法和数值模拟的局限，具体的结果还有待作进一步应用和推广。