0 引言

地下水是影响边坡稳定的重要因素之一，一方面地下水在裂隙中渗流产生的渗透水压力可改变岩体的初始应力状态；另一方面岩体应力状态的改变引起岩体变形，并导致裂隙几何参数尤其是开度的变化，使岩体的渗透性发生改变，最终导致岩体渗流场的变化。统计资料^[1]显示，约90%以上的岩质边坡失稳破坏与地下水作用有关，据Hoek^[2]的研究，边坡饱水时的安全系数比边坡干燥时降低0.5~0.8。因此，研究渗流场与应力场耦合作用下地下水在岩质边坡中渗流规律，对边坡稳定性分析至关重要。

岩体中该耦合问题的边界条件和材料性质均非常复杂，很难对其开展物理模型试验，前人普遍采用数值方法求解边界条件，并分别在水电站大坝^[3-4]、土坡^[5]、岩体边坡^[6-7]等不同的工程中取得了较好的应用效果；但是这些研究仅基于介质的各向同性假设，并不能反映耦合情况下渗流场的各向异性特征。而层状岩体(沉积岩、部分变质岩等)具有一组占绝对优势的结构面，是复杂岩体的典型特征之一，其变形和渗透均具有各向异性特征^[8-9]，在这种岩体中，通常岩石坚硬不透水或弱透水，地下水主要赋存运动于节理裂隙中，节理的分布形式和渗透能力极大地影响着岩体中地下水流模式。即使对于土坡渗流，渗透各向异性也是一个不容忽视的影响因素^[10]。现场测试表明岩体渗透性具有一定的空间变化规律^[11]。周玉新等^[12]采用渗透张量虽然能够揭示层状岩体边坡裂隙控渗及渗流各向异性特征，但其只是对单一的渗流场进行分析，渗透各向异性系数始终是常数，难以反映耦合作用下边坡渗流场的各向异性特征以及渗透系数的空间变化规律。

本文针对层状边坡中更容易发生失稳破坏的顺倾向边坡，应用基于等效连续介质和Louis经验公式的各向异性岩体渗流应力耦合模型^[13]，以及COMSOL多物理场耦合软件对耦合作用下边坡渗流场进行数值模拟，以研究顺倾向层状边坡渗流场的各向异性特征及渗透系数的空间变化规律。

1 层状岩体各向异性弹性力学本构

层状岩体存在一组占绝对优势的结构面，是一种层合结构材料，因而其变形和强度都是各向异性的^[14]，考虑到岩体中定向节理系统的影响，将层状岩体简化为具有一个弹性对称面的各向异性连续介质，根据广义虎克定律，弹性阶段的应力应变关系可表示为

式中：σ_ij为应力张量；ε_kl为应变张量；D_ijkl为弹性张量。

岩体工程(如大坝、边坡等)在一定条件下可以简化为平面应变问题，那么D_ijkl可由各向异性工程弹性常数表示为

式中：E_i(i=1, 3) 为岩体等效弹性模量，1和3分别代表法向和切向，GPa；v_i(i=1, 3) 为等效泊松比；G₁₃为等效剪切模量，GPa。

孔隙水压力作用下，根据有效应力原理，岩体本构方程^[15]可表示为

式中：σ′_ij为有效应力张量；α为流固耦合参数；p为孔隙水压力，Pa；δ_ij为Kronecker常数。

2 层状岩体应力与渗透性关系

对于二维渗流问题，忽略水的可压缩性，稳态条件下渗流主方向上连续性方程为

式中：K_ij(σ)为与应力有关的主渗透张量；▽为Hamiltonia算子。

层状岩体存在一组占绝对优势的层理结构面，同时在其他方向也存在少量随机分布的结构面，如图 1所示。由于裂隙的控渗特征，层状岩体在沿着层理面方向和垂直层理面方向表现出强烈的渗透各向异性。基于文献[15]的岩体分类，将层状岩体简化为等效连续介质力学模型，如图 2和图 3所示。

研究^[16]表明，裂隙岩体在三向应力作用下以垂直于裂隙面方向的应力对渗流起主导作用；因此，鉴于层状边坡特有的结构面分布特征(层状结构面密集展布，法向稀疏)，二维条件下，基于Louis经验公式，应力影响下沿着和垂直层状结构面方向的渗透系数K₁₁和K₃₃可表示为：

式中：K₁，K₃分别为法向和切向层状结构面方向的初始等效渗透系数；β为耦合参数；σ_1n，σ_3n分别为对应于K₁，K₃所在方向的法向应力，可表示为

式(5)定义了应力变化对主渗透性的作用机制，体现了各向异性应力场作用下地下水渗流的非均匀性和各向异性，该模型在抚顺西露天矿南帮渗流场分析中取得了较好的应用^[17]。将层状岩体等效为连续介质模型时未考虑岩体的破裂损伤，适用于层状结构面产状、发育程度基本一致，且分布较均匀的岩体。一旦岩体破裂产生了优势渗流通道，地下水将沿着优势渗流通道运动，渗流极不均匀，此时采用等效连续介质模型会造成较大误差。

层状岩体渗流应力耦合控制方程由一系列非线性偏微分方程构成，其解是相互耦合的，即使赋予合适的边界条件，也难以得到解析解。本文采用数值模拟的方法，应用COMSOL软件，通过自定义弱解形式偏微分方程，对偏微分方程进行转换，形成统一的通式形式的微分方程组，并嵌入渗流场与应力场之间的交叉耦合项，采用不完全三角分解法(LU)分解的显式迭代对边坡应力场和渗流场进行数值求解，求解精度为10^-6。

3 数值模拟分析 3.1 层状边坡数值模型

按照平面应变问题建立二维层状边坡数值模型，对于顺倾向边坡，层状结构面倾角θ在0°~90°之间变化，边坡模型几何尺寸如图 4所示，坐标系原点O设在坡脚位置。计算区域被划分为3 869个三角形单元，选用的参数见表 1。

渗流边界条件为：模型左侧、右侧固定水头边界H分别为210，350 m；底部为相对隔水边界。

位移边界条件为：底边固定，左侧、右侧边界均约束法向位移，上部为自由边界。

应力边界条件为：荷载考虑自重和渗透压力。渗透压力分为静水压力和动水压力，求解时先将渗流场计算出的渗透压力转化为渗透体力，再将该渗透体力施加到对应的节点和单元上。

3.2 数值模拟结果

图 5给出了不同结构面倾角的边坡等水头线分布。由图 5可知，层状结构面的高渗透性及其在边坡结构面分布中的绝对优势，导致边坡渗流场表现出显著的结构面控渗特征。对比θ=0°和θ=90°时渗流场分布情况可知，θ=90°时，结构面垂直分布，垂直渗透性好，地下水更容易下渗，达到稳态后，潜水面较θ=0°低。θ=45°时，边坡层状结构面为顺倾，地下水容易沿着层状结构面排出，水位面降低，岩体内部某点的压力水头高于该位置的水位面。

在坡面和坡顶分别取测线AA′(x=200 m)和BB′(x=300 m)，监测边坡不同位置的水位降深与结构面倾角的关系，如图 6所示。由图 6可知，对于顺倾向边坡，水位降深随结构面倾角θ的增大先升高、后降低，呈现两头低、中间高的形态，倾角θ约为42°时水位降深最大，即潜水面最低；且越靠近溢出点，结构面倾角对水位降深影响越大。

地下水在边坡斜面上的溢出点位置见表 2，其中：Δz越小表明地下溢出水位越低，当θ≈42°时水位溢出点最低；对于该边坡模型，θ在30°~75°之间变化，地下水溢出点高度在坡高的1/4~1/3之间。

图 7所示为坡面AA′位置的渗透各向异性系数与埋深的关系。由图 7可见，随着埋深的增大，岩体渗透各向异性系数(K₁₁/K₃₃)不断减小；这是由于埋深越深，地应力越高，层状结构面趋于闭合，裂隙渗流通道受阻，裂隙控渗特征由显著到逐渐变弱，表现出向各向同性渗流过渡的趋势，该模拟结果与文献[18]不同钻孔的现场测试的结果相比，两者的变化趋势基本吻合。

图 8给出了坡面AA′位置上不同埋深处渗透各向异性系数与结构面倾角的关系。由图 8可知：对同一埋深位置，随着结构面倾角的不断增大，渗透各向异性系数也表现出先增大后减小的趋势，说明渗透各向异性系数与水位降深变化趋势有较好的一致性，其中结构面倾角约为42°时渗透各向异性系数最大，结构面水平分布(θ=0°)时渗透各向异性系数最小；对同一层状结构面产状的岩质边坡，埋深越深，岩体渗透各向异性系数减小越缓慢。

4 结论

应用基于等效连续介质模型和Louis经验公式建立的各向异性岩体渗流应力耦合模型，对顺倾向层状边坡的各向异性渗流规律进行了模拟分析，得到结论如下：

1) 层状结构面的高渗透性及其在层状边坡分布中的绝对优势，导致层状边坡呈现出强烈的各向异性渗流特征。对于顺倾向层状边坡，水位降深随结构面倾角θ的增大先升高、后降低，呈现两头低、中间高的形态，其中在倾角约为42°时水位降深最大，即潜水面最低。因此在层状边坡稳定性分析中应重视结构面倾角对地下水位的影响。

2) 渗流应力耦合作用下，层状边坡地下水为非均匀、各向异性渗流。埋深增大时，地应力升高，层状结构面趋于闭合，裂隙渗流通道受阻，岩体渗透各向异性系数随之降低，裂隙控渗特征由显著到逐渐变弱；表现出向各向同性渗流过渡的趋势。

3) 顺倾向层状边坡同一监测点的渗透各向异性系数随着结构面倾角的增大呈现先增大后减小的趋势，与水位降深的变化趋势相吻合。对于同一层状边坡，随着埋深不断增大，渗透各向异性系数衰减的幅度也会逐渐减弱。该模拟结果可用于层状边坡不同深度的渗透张量预测。