殷德胜，汪卫明，陈胜宏

（武汉大学 水资源与水电工程科学国家重点实验室，武汉 430072）

1 引 言

1959年法国Malpasset拱坝的溃坝事故，导致裂隙岩体渗流与应力的耦合问题成为岩石水力学的热门课题。由于该问题的边界条件和材料性质的复杂性，很难进行物理模型试验，人们普遍采用数值解法。

有限单元法是连续介质力学的代表，它有着完备的理论体系，在岩体渗流-应力耦合研究中的应用最为广泛。Noorishad等[1]利用土体固结分析理论，并运用有限元法中的节理单元，研究了非连续介质中的固液两相介质的耦合间题；Oda[2]将裂隙岩体渗透张量与弹性张量统一用裂隙网络几何张量表示，建立了较为严密的裂隙岩体渗流特性与应力特性之间的耦合关系；在国内，陈平和张有天[3]考虑了裂隙变形的非线性特性，以裂隙平均隙宽为计算参数，对立方定理作了修正，提出了渗流-应力耦合计算方法；王媛等[4]则建立的同时以渗流水压力和位移为未知值的4自由度全耦合法，这降低了耦合分析的难度。由于已有多种成熟的商业软件如ABAQUS、ADINA等能够模拟多场耦合问题，这将进一步加强有限元法在岩土工程渗流-应力耦合这一领域中的应用。

对于复杂的工程，已有的研究成果大多基于二维情况，由于三维裂隙岩体的几何复杂性，目前尚无一个强大的前处理器能够轻松地对其进行网格离散，这是限制有限元法在该领域进一步应用的主要原因。学者们只能退而求其次，大多采用等效连续介质模型进行耦合研究。然而，裂隙中的水主要通过裂隙等结构面流动，其渗流与传统的多孔介质渗流存在着本质上的差别。对裂隙岩体进行渗流分析，宜采用裂隙网络模型或双重介质模型[5]。

双重介质模型除裂隙网格外，将岩块视为渗透系数较小的连续渗透介质，且能考虑岩块孔隙与岩体裂隙之间水量的交换，是更接近于实际的模型。柴军瑞等[6]将各种裂隙和孔隙按规模和渗透性分为4级，建立起岩体渗流场与应力场耦合分析的多重裂隙网络模型。综合来说，同连续介质方法相比，双重介质模型优势明显，但由于继承了网格剖分的困难，实现难度相对更大。综合起来，不连续介质力学分析方法结合裂隙网络渗流模型则是个相对合理可行的思路。

不连续变形分析方法（discontinuous deformation analysis method，DDA）属于不连续介质力学范畴，十分适合处理裂隙岩体渗流、应力的计算问题。张国新等[7]应用DDA法，对日本某事故隧道进行了二维渗流-应力耦合分析，得到了边坡失稳坍塌的机制。应用DDA法进行三维渗流-应力的耦合分析，在计算中尚存在接触难以准确确定、罚函数如何施加等诸多困难，目前的研究成果很少。

块体单元法也属于不连续介质力学范畴，是一种能对节理裂隙岩体进行渗流、应力分析的数值方法[8]。在块体单元法的基础上，本文考虑了节理裂隙的刚度和导水系数随开度的变化关系，并将常规的块体单元法结构计算模块和渗流分析模块整合成一个系统，应用两场交叉迭代技术，完成了随机裂隙岩体的渗流-应力耦合研究。该方法前、后处理简单，计算过程也不复杂，能有效地应用于三维随机裂隙岩体的耦合分析中。

2 块体单元法原理简介

岩体的变形主要都集中于节理、裂隙等结构面上，相比而言，块体的变形可以忽略不计。在此假设条件下，块体单元法认为块体为刚体，只考虑结构面的变形和强度特性，以块体单元形心位移为基本未知量，通过建立块体单元的平衡方程、结构面变形和块体单元位移的几何相容方程、结构面的本构方程，综合推导出块体系统的整体平衡方程。

同时，断层、夹层、节理和裂隙等结构面不仅是岩体中的软弱面，且由于其渗透系数通常远大于完整岩石，故它们也是主要的渗流通道。假定岩块不透水，即在岩块中渗透系数k =0，渗流只发生在裂隙面上，根据裂隙网络交叉处的流量守恒和分片变分原理，建立了与块体单元法相配套的随机裂隙岩体的渗流分析方法[9]。

有关块体单元法的应力、渗流分析的详细过程请参阅文献[8]。

3 耦合机制

岩体的变形大部分由结构面产生，结构面隙宽的微小改变将引起渗流量的重大变化。同时，渗流场的变化会引起水荷载的变化，从而导致应力场的重分布。因而，裂隙岩体的结构场和渗流场之间存在着重要的耦合关系。

然而，在应力场和渗流场相互影响的定量描述上，目前尚无定论。由于裂隙面上的法向应力引起裂隙变形相对显著，为抓住主要矛盾，本文仅研究法向应力对裂隙变形特性、渗透特性的影响。

3.1 结构面变形对裂隙刚度系数的影响

结构面的相对变形影响着刚度系数，从而间接地影响应力场的分布。陈胜宏等[10]曾提出有关节理面变形和应力关系的“充填模型”，将节理裂隙视为一种均匀介质，并据此推导了有关节理面法向、切向刚度系数与应力之间简单而实用的关系式。

以单裂隙岩体说明裂隙渗流与应力的耦合关系。设图1中的裂隙面受法向应力σz和切向应力τzx,τzy作用，由于裂隙面的开度较小，故应变满足：

式中： εx、 εy分别为x、y向的正应变， γxy、 γyx为xy平面上的剪应变。

图1 单裂隙岩体Fig.1 Rock mass containing one fracture

“充填模型”假设裂隙面上，“充填介质”是各向同性体。据此，可得到裂隙面局部坐标下的弹性本构关系为

式中：λ、G为裂隙内充填物的拉梅常数。记e0为裂隙面应力为0时的开度，ux、 uy、 uz分别为x、y、z三向的相对位移，则裂隙开度满足：

于是式(2)可以变成：

式中： kn、 ks分别为裂隙面的法向、切向的刚度系数，可由充填物的弹性参数和开度确定：

由于开度e与法向变形uz相关，而uz又与应力有关联，故 kn、 ks都是应力的函数。

从式(4)中取第1行，可以得到：

结合公式(3)，上式可写成：

将式(7)积分后得到：

把式(9)代入式(5)，即可得到法向、切向刚度系数与法向应力的关系：

记 kn0= λ+ 2G /e0,ks0= G/e0，用它来表示初始刚度，即裂隙在无应力状态时的刚度系数，则式(10)可写成：

式(11)表明了裂隙刚度系数与法向应力之间存在着指数关系。

3.2 结构面变形对导水系数的影响

立方定律是研究裂隙岩体渗流的基础。然而，天然节理并非完全光滑平直，也不是纯粹的孔隙，因而立方定律的适用条件有限。近几十年来，许多考虑了节理面的粗糙性和凸起物高度分布情况的修正立方定律[11－12]相继被提出，其主要思想就是用等效水力隙宽代替几何平均隙宽。然而，等效水力隙宽和几何平均隙宽是完全不同的两个概念[13]，若将其直接用于耦合研究，不但理论说不通，且流量与隙宽的3次方成正比，会导致耦合分析中计算过于敏感和收敛困难。故立方定律在耦合研究中直接应用较少。

“充填模型” 同样适用于结构面的渗流特性的研究。该模型中假设裂隙中填充物的渗透系数不变，而导水系数随着开度的变化而变化，并结合式(9)，得到：

式中：kd为导水系数；kd0为初始导水系数，即裂隙面无应力时的导水系数，可通过试验获得，其他同式(10)。

式(11)、(12)能较好地模拟裂隙刚度系数和导水系数与应力之间的关系。本文运用这套关系式，结合块体单元法进行了裂隙岩体渗流-应力的耦合研究。

3.3 耦合步骤

根据结构面变形对应力场、渗流场的影响关系式，使用Fortran语言，编制了块体元渗流-应力的耦合分析程序。该程序中，充分考虑块体元应力、渗流分析的共同点，将2个分析模块整合在一个程序中，显著减少了公共计算部分的计算量，且相比传统的两场迭代，不需要2个程序相互调用时读写0时文件的时间开销，且降低了文件传递导致的累积误差，耦合计算的效率和精度均较高。

程序中，耦合分析的具体实施步骤如下：

(1)由初始刚度系数和初始导水系数，结合边界条件，计算得到初始应力场 {σ}(0)和初始渗流场{φ}(0)；(2)结合式(11)、(12)，根据{σ}(0)计算裂隙面新的法向、切向刚度系数和导水系数，并重新生成结构计算的刚度矩阵和渗流计算的导水矩阵；(3)由新的导水矩阵和渗流边界条件，重新计算渗流场{φ}(i)；(4)由新的渗流场 {φ}(i)计算出渗透荷载，并加上其它荷载，重新计算应力场 {σ}(i)；(5)计算应力、水头的增量： {Δ σ} ={σ}(i)-{σ}(i-1)和{Δφ}={φ}(i)-{φ}(i-1)，并取{Δσ}和{Δφ}的2范数。若和小于给定的误差值，则迭代结束，否则重新执行步骤(2)，直至渗流场和应力场都满足收敛精度为止。

4 算例验证

薛娈鸾[14]曾用基于双重介质模型的复合单元法研究过裂隙岩体的渗流-应力耦合问题。本文采用她文中的算例，从规律上对本文方法进行验证。如图2所示的长方体岩块试件，坐标系在左下角，试件的x方向长为1.0 m，y方向宽为0.1 m，z方向高为0.2 m。在z =0.1 m处有一条初始开度为1 mm的水平贯穿性裂隙。由于块体单元法假设块体为刚体，本文采用卓家寿等提出的界面元法[15－16]，用自定义的虚拟结构面将试件等分成若干块体，将变形累积在结构面上，以此来近似模拟。

对于均质材料，根据文[16]的建议，虚拟结构面的法向刚度、切向刚度按式(13)选取：

式中： h1、 h2为块体形心到结构面的距离；E为弹性模量；μ为泊松比。如图2所示，总共引入9条虚拟结构面，将试件分成20个块体。另外，虚拟结构面的导水系数的取值为实际裂隙导水系数的1%。所有结构面的参数见表1：

取岩块试件的左侧水位为 10 m，右侧水位为5 m，其它外表面为不透水边界。考虑水压力并在上表面正中间的 23、53两点处各施加集中力 F =0.1 MN，不考虑试件的自重。在除上表面外的其余5个侧面均施加法向约束。

图2 算例块体系统Fig.2 Block element system of the example

表1 裂隙面的参数Table 1 Parameters of fractures

整理耦合和不耦合两种工况下的成果，图3为结构面上的水头分布，图4为不耦合和耦合两种工况下y =0.05 m断面上的流速矢量图。

图3 结构面上的水头分布Fig.3 Hydraulic distribution on the structure

图4 流速矢量图Fig.4 Hydraulic velocity vector diagram

从图3中可以看出，不考虑耦合作用时，结构面上的水头呈线性分布；考虑耦合作用后，集中力作用点的左侧部位，结构面的水头趋于上游水位，右侧部位的水头则趋于下游水位。对比图4可以看出：当不考虑耦合时，整个裂隙面上的流速是相同的；考虑耦合后，集中力作用点底部压应力较大区域的导水系数减小，渗透坡降增加，流速增加，其他区域的流速减小。结果显示了本文耦合机制和程序的合理性。

5 工程算例

图5是贵州省乌江沙沱水电站非溢流坝段的典型剖面的示意图。其岩基上含有F88断层，j4、j5、j6三条夹层，并有两组倾角约为45°的正交裂隙。沿坝轴线方向取宽度为1 m的坝段参与计算，对于断层和夹层，采用确定性方式模拟；对于裂隙，采用蒙特卡罗法随机生成，并用块体识别方法[9]找出所有块体。块体系统如图6所示。

图5 重力坝剖面Fig.5 Profile diagram of gravity dam

图6 块体系统Fig.6 Block element system

对于地基部分，断层、夹层、裂隙等结构面的参数可以通过实验的方式获得；对于坝体部分，仍按卓家寿的方法，按5 m的间距获得虚拟结构面的参数。帷幕采用文献[17]的方法，即将帷幕块体周边结构面视为相对不透水面，其导水系数取为较小值，本文取值为1×10-8m/s。

断层、夹层的材料参数见表 2，裂隙参数见表3，坝体和基岩的参数见表4。

表2 结构面参数Table 2 Parameters of discontinuities

表3 裂隙参数Table 3 Parameters of fractures

表4 基岩和坝体混凝土参数Table 4 Parameters of rock foundation and dam concrete

同样，计算也分为不耦合和耦合两种工况。上下游水位见图5，考虑的荷载有坝体和基岩的自重、静水压力和渗透荷载作用。

整理了两种工况的计算成果对比。图7为裂隙网络的水头等势线，图8为水力坡降矢量图。

图7 裂隙网络等水头线（单位：m）Fig.7 Hydraulic potential isolines of fracture network (unit: m)

图8 水力坡降矢量图Fig.8 Hydraulic gradient vector diagrms

从图 7、8中可以看出：考虑耦合作用以后，坝前地基的水力坡降逐渐减小，水头普遍增大；坝后基础的变化规律与坝前类似；坝体下方地基的水平向主导的裂隙上的压应力增加，导水系数减小，坡降、流速均有所增加，进而导致建基面的扬压力增加，对坝体的稳定不利。

6 结 论

（1）块体单元法本身即是为解决节理、裂隙岩体的应力、渗流问题而提出，将其应用于裂隙岩体的渗流-应力耦合中，有着先天的优势。而且能考虑三维岩体中含有的复杂随机裂隙网络，这是其他模型所不能考虑或很难考虑的；

（2）法向应力对裂隙的渗透特性影响很大。在高压应力的区域，导水系数减小，渗透坡降增加，流速增大；其他区域流速减小。从算例中可以看出，重力坝分析中若不考虑耦合作用，则可能导致扬压力的计算值偏小。对实际工程会产生不利影响。

然而，就目前的计算机发展水平，采用块体单元法进行渗流-应力耦合分析的计算量仍显偏大，需要进一步优化数据结构、改进积分算法，甚至采用分布式并行计算，使该方法真正意义上用于工程的三维分析中。

[1]NOORISHAD J, AYATOLLAHI M S, WITHERSPOON P A. A finite-element method for coupled stress and fluid flow analysis in fractured rock masses[J]. International Journal of Rock Mechanics and Mining Sciences, 1982:185－193.

[2]ODA M. An equivalent continuum model for coupled stress and fluid flow analysis in jointed rock masses[J].Water Resources Research, 1986, 22(13): 1845－1856.

[3]陈平, 张有天. 裂隙岩体渗流与应力耦合分析[J]. 岩石力学与工程学报, 1994, 13(4): 299－308.CHEN Ping, ZHANG You-tian. Coupling analysis of seepage/stress for joint rock[J]. Chinese Journal of Rock Mechanics and Engineering, 1994, 13(4): 299－308.

[4]王媛, 徐志英, 速宝玉. 裂隙岩体渗流与应力耦合分析的四自由度全耦合法[J]. 水利学报, 1998, (7): 55－59.WANG Yuan, XU Zhi-ying, SU Bao-yu. Four-freedom complete method for the seepage-stress coupled analysis in fissured rock masses[J]. Journal of Hydraulic Engineering, 1998, (7): 55－59.

[5]张有天. 岩石水力学与工程[M]. 北京: 中国水利水电出版社, 2005.

[6]柴军瑞, 仵彦卿. 岩体渗流场与应力场耦合分析的多重裂隙网络模型[J]. 岩石力学与工程学报, 2000, 19(6):712－717.CHAI Jun-rui, WU Yan-qing. Research on multiple-level fracture network model for coupled seepage and stress fields in rock mass[J]. Chinese Journal of Rock Mechanics and Engineering, 2000, 19(6): 712－717.

[7]张国新, 武晓峰. 裂隙渗流对岩石边坡稳定的影响——渗流、变形耦合作用的DDA法[J]. 岩石力学与工程学报, 2003, 22(8): 1269－1275.ZHANG Guo-xin, WU Xiao-feng. Influence of seepage on the stability of rock slope coupling of seepage and deformation by DDA method[J]. Chinese Journal of Rock Mechanics and Engineering, 2003, 22(8): 1269－1275.

[8]陈胜宏. 计算岩体力学与工程[M]. 中国水利水电出版社, 2006.

[9]殷德胜, 汪卫明, 陈胜宏. 三维随机裂隙岩体渗流分析的块体单元法[J]. 岩土力学, 2009, 30(8): 1269－1275.YIN De-sheng, WANG Wei-ming, CHEN Sheng-hong.Block element method for seepage analysis in three dimensional random fracture network[J]. Rock and Soil Mechanics, 2009, 30(8): 1269－1275.

[10]陈胜宏, 王鸿儒, 熊文林. 节理面渗流性质的探讨[J].武汉水利电力学院学报, 1989, 22(1): 53－60.CHEN Sheng-hong, WANG Hong-ru, XIONG Wen-lin.Study on the seepage characteristics of joint surface[J].Journal of Wuhan University of Hydraulic and Electric Engineering, 1989, 22(1): 53－60.

[11]LOMIZE G M. Flow in fractured rocks[M]. Moscow:Gesenergoizdat, 1951.

[12]LOUIS C. Rock mechanics[M]. New York: Elsevier Science, 1974.

[13]王媛. 单裂隙面渗流与应力的耦合特性[J]. 岩石力学与工程学报, 2002, 21(1): 83－87.WANG Yuan. Coupling characteristic of stress and fluid flow within a single fracture[J]. Chinese Journal of Rock Mechanics and Engineering, 2002, 21(1): 83－87.

[14]薛娈鸾, 陈胜宏. 岩石裂隙渗流与法向应力耦合的复合单元模型[J]. 岩石力学与工程学报, 2007, 26(增刊):2613－2619.XUE Luan-luan, CHEN Sheng-hong. Composite element model of seepage-normal stress coupling for rock fractures[J]. Chinese Journal of Rock Mechanics and Engineering, 2007, 26(Supp.): 2613－2619.

[15]卓家寿, 赵宁. 不连续介质静、动力分析的刚体-弹簧元法[J]. 河海大学学报, 1993, 21(5): 34－43.ZHUO Jia-shou, ZHAO Ning. Piecewise rigid body-interface spring method for problems of discontinuous medium[J]. Journal of Hohai University,1993, 21(5): 34－43.

[16]方义琳, 卓家寿, 章青. 具有任意形状单元离散模型的界面元法[J]. 工程力学, 1998, 15(2): 27－37.FANG Yi-lin, ZHUO Jia-shou, ZHANG Qing. Interface stress element method for distinct model with elements of arbitrary shape[J]. Engineering Mechanics, 1998, 15(2):27－37.

[17]汪卫明, 徐明毅, 陈胜宏. 复杂边界条件下的岩体网络渗流分析[J]. 岩石力学与工程学报, 2001, 20(4): 473－476.WANG Wei-ming, XU Ming-yi, CHEN Sheng-hong.Analysis method of network seepage for blocky rock mass with complicated boundary condition[J]. Chinese Journal of Rock Mechanics and Engineering, 2001,20(4): 473－476.