裂隙網絡無壓滲流分析的初流量法

姚 池，姜清輝， ，葉祖洋，周創兵，

（1. 武漢大學 水資源與水電工程科學國家重點實驗室，武漢 430072；2. 武漢大學 土木建筑工程學院，武漢 430072）

1 引 言

當前，描述裂隙巖體滲流的數學模型主要有3種方法，即等效連續介質模型、離散介質模型和雙重介質模型。等效連續介質模型將裂隙巖體等效為連續體[1－3]，假定地下水在巖體中的流動服從Darcy 定律，然后運用經典的連續滲流理論進行分析，目前應用較為廣泛，但該模型不能很好地反映空間結構面組合形成的裂隙網絡的導水作用及優勢流效應，在某些情況下數值模擬結果與實際出入很大。離散介質模型假定巖塊不透水[4－6]，滲流只發生在巖體裂隙內，地下水通過相互連通的裂隙網絡進行傳導，巖體滲透特性受裂隙的跡長、走向、傾角和隙寬等參數控制，這種方法抓住了巖體中滲流主要受裂隙等結構面控制這一本質特征，具有廣泛的應用前景，其主要缺點在于實際工程中裂隙網絡資料不易獲取，因此，常通過Monto-Carlo 法隨機生成具有相同統計規律的裂隙網絡來代替真實的裂隙網絡。雙重介質模型同時考慮了裂隙的強導水性、巖塊的儲水作用及巖塊和裂隙的水量交換[7－8]，但該模型對裂隙系統的幾何形態和空間分布存在較嚴格的假設，從而限制了其應用。

由于巖體本身賦存環境的復雜性，裂隙的發育極不規則，裂隙網絡的幾何構成往往非常復雜，導致滲流在巖體中的運動顯示出明顯的非均勻性和各向異性，很多情況下并不存在所謂的 REV（representative elementary volume）。此時，如用等效連續介質方法進行滲流分析實際上是不合適的，應該采用離散裂隙網絡的模擬方法。總體上，相對基于多孔連續介質滲流分析方法，目前基于離散裂隙網絡模型滲流問題研究相對較少，尤其對有自由面的裂隙網絡滲流問題，國內外只有少數的幾篇論文專門對其進行了探討。Jing 等[9]利用DDA（discontinuous deformation analysis）形成的裂隙網絡，采用剩余流量法求解裂隙網絡自由面。王恩 志[10]借鑒初流量法的思想，提出了一種計算裂隙網絡滲透自由面的迭代算法。柴軍瑞等[11]提出含自由面巖體裂隙網絡非線性滲流分析的傳導矩陣調整法。由于上述方法主要借鑒求解連續介質滲流自由面的一些常用的直覺化求解方法，這些方法不具有嚴格的理論的基礎，且因為沒有解決溢出點的奇異性，計算過程中往往會出現數值的不穩定，對于含有大規模的復雜裂隙網絡則難以收斂。

變分不等式方法是近年來求解連續介質無壓滲流分析理論上的一個重要進展[12－13]，本文把單條裂隙概化為平行板模型，通過將Darcy 定律延拓到整個計算區域，定義了在整個區域上的非線性邊值問題，并將潛在溢出面邊界條件歸納為Signorini型邊界條件，建立與PDE 提法等價的變分不等式(variational inequality， VI)提法，從而在理論上解決了出滲點的定位問題。通過采用有限元方法的線單元模擬裂隙，在有限元計算格式中，引入連續的罰Heaviside 函數，顯著改善了離散裂隙網絡滲流自由面迭代過程中的數值不穩定性。

2 有自由面裂隙網絡滲流問題的描述

對于一個任意的二維裂隙網絡，裂隙網絡的連通性及裂隙的跡長、水力開度、位置等參數決定了地下水滲流在巖體中的分布，由于孤立的裂隙、只與一條裂隙相交的枝狀裂隙及裂隙端部的獨頭不能作為滲流通道，對流場沒有貢獻，因此，在進行滲流分析前需對隨機生成的裂隙進行必要的刪減和截斷，經過凈化處理形成連通的裂隙網絡，凈化處理前后裂隙網絡的對比可參見圖1。

不失一般性，以圖1 所示裂隙巖體邊坡滲流問題為例，對于有自由面的無壓滲流問題，區域Ω（ Ω=Ωw∪ Ωd）上的滲流實際上僅在自由面Γf以下的濕區Ωw中運動。在多孔連續介質滲流理論中，自由面是一個由計算域上壓力水頭為0 且法向流速為0 的點組成的連續曲面。而對于裂隙巖體而言，當忽略巖塊滲流時，裂隙巖體的滲流是沿裂隙網絡定向進行的，自由面在裂隙處表現為潛水面，依次連接裂隙處的潛水面位置可形成裂隙巖體滲流的自由面邊界（見圖1）。

圖1 裂隙巖體邊坡滲流示意圖 Fig.1 Schematic diagram of seepage flow through a fractured rock slope

假定邊坡上游水位面位于A 點，下游水位面位于D 點，BC 為不透水邊界，E 為溢出點，AE 為由裂隙處的潛水位構成的自由面。假定巖塊不透水，地下水僅在自由面以下的裂隙網絡濕區 wΩ 內產生流動。

裂隙網絡濕區內任一點的總水頭φ 可以定義為

式中：z 為垂直坐標分量；p 為裂隙水壓力；wγ 為水的重度。

圖2 裂隙ij 的局部坐標系 Fig.2 Local coordinate system of the fracture ij

把巖體中的裂隙概化為光滑平行板模型，在二維裂隙網絡中，單個裂隙上的滲流是一維運動。如圖1 裂隙ij 的局部坐標所示，對任意裂隙ij 建立局部坐標系，根據Darcy 定律，在裂隙ij 內的流速可表示為

式中：ijk 為裂隙ij 的滲透系數。根據立方定理可得，

式中：g 為重力加速度；ijb 為裂隙ij 的隙寬；υ 為水的運動黏滯系數。

水在裂隙ij 內的流動應滿足連續性方程：

可知，ijv 在裂隙ij 上為常量。設有m 條裂隙匯聚于節點i，如圖3 所示。由于節點i 本身不具備儲存流體的能力，依據質量守恒原理，節點總流量的代數和為0，即

圖3 通過節點i 的滲流量 Fig.3 Schematic diagram for flow distribution at fracture intersection.

水在裂隙網絡濕區內的運動應滿足連續性方程（4）和（5）以及下列邊界條件： （1）水頭邊界條件

（2）流量邊界條件

式中：Γq為隔水邊界（規定流出為負）。

（3）自由面邊界條件

式中：p 為裂隙ij 與自由面的交點。

（4）出滲面邊界條件

至此，給出了帶自由面裂隙網絡滲流問題的基本方程和邊界條件，自由面fΓ 的位置事先未知，因此，確定其位置成了帶自由面滲流分析的主要內容。

3 裂隙網絡滲流全域上邊值問題的PDE(partial differential equation)提法

為了建立全區域Ω 上的邊值問題，借鑒初流量法的思想，將Darcy 定律延拓到整個區域

式中：0v 為初流速，初流速0v 的引入是為了抵消干區 dΩ 實際不存在的流速，其具體形式為

其中

顯然，(10)式中的 vij仍滿足連續性方程。為了回避Ωw的不確定性，將構造一個定義在全域上的邊值問題，該邊值問題的基本方程由連續性方程、質量守恒方程和廣義Darcy 定律構成，水頭邊界條件和流量邊界條件保持不變。剩下的邊界AGFD 為Signorini 型互補邊界條件，

綜上所述，定義在整個裂隙網絡上邊值問題的提法是：求函數φ，使其滿足控制方程（4）和（5）、外部邊界條件（6）、（7）和式(13)以及內部自由面邊界條件式（8）。

4 裂隙網絡滲流全域上邊值問題的VI 提法

上述PDE 提法中含有未知的內部自由邊界fΓ和潛在的溢出面Signorini型互補條件，數值求解中選取其試探函數非常困難。通過建立與PDE 提法等價的變分不等式(VI)提法，使得自由邊界fΓ 以及潛在溢出邊界上的互補條件轉化為自然邊界條件，從而大大減小數值求解的難度，并改善數值計算的收斂性和穩定性。

定義試探函數集VIΦ 為

則在數學上與上述PDE 提法完全等價的變分不等式(VI)提法可表述為：在ΦVI中求一函數φ，使得對?ψ ∈ΦVI，都有如下不等式成立

其中， ξ (φ ，ψ - φ)和 ζφ(ψ - φ)的形式為

為了完成兩種提法的等價性證明，首先利用分部積分展開式 ξ (φ ，ψ - φ ) - ζφ(ψ - φ)：

式中：n 為整個區域所有裂隙節點的總數；im 為匯聚于i 節點的裂隙條數。

4.1 PDE 提法?VI 提法

如果φ 是PDE 提法的解，將連續性方程式（4）和（5），外部邊界條件（6）、（7）和（13）代入式（18），并注意到條件式（14），得

因此，φ 也是VI 提法的解。

4.2 VI 提法?PDE 提法

假定φ 是變分(VI)提法的解， ?ψ ∈ΦVI都滿足式（18）。分別取ψ = φ + θ1和ψ = φ - θ1，其中θ1為在所有裂隙節點上等于0 的任意函數，可以得到連續性方程為

因此，式（18）可以簡化為

在式（21）中分別取ψ = φ + θ2，ψ = φ - θ2，其中θ2為在 Γφ、 Γq和 Γs邊界上裂隙節點等于0 的任意函數，可以得到內部任意節點上的質量守恒方程為

因而式（21）可以簡化為

對式（23）分別取ψ = φ + θ3，ψ = φ - θ3，其中 θ3為在 Γφ和 Γs邊界上裂隙節點等于0 的任意函數，可以得到 Γq上裂隙節點的流量邊界條件為

因此，式（23）可以進一步簡化為

考慮到上式等號右邊的第二項為一常數，要使不等式對于 ?ψ ∈ΦIV都成立，則等號右邊的第一項必須非負，即

考慮到式（14），很明顯可以得到

在式（25）中當然也可以取 zψ= ，如此就有

由式（27）、（28）以及在sΓ 上iizφ ≤ ，可得sΓ 上的互補條件為

為導出內部自由面fΓ 邊界條件，將式（18）中的Ω 表示成 wΩ 和 dΩ 上的積分之和為

利用上面已經導出的條件，可得對于?ψ ∈ΦIV，

分別取ψ = φ + θ4和ψ = φ - θ4，θ4為在 Γφ和 Γs邊界上裂隙節點等于0 的任意函數，如此可得：

至此，已經導出了PDE 提法中的全部方程及邊界條件。

5 有限元數值求解格式

5.1 自適應罰Heaviside 函數

盡管可以對式（15）進行離散來求其數值解，但由于其中的Heaviside 函數的不連續性，將導致數值不穩定性，主要原因是Heaviside 函數為非連續的階梯式函數，在對單元積分時會造成數值的跳躍。為避免這種數值跳躍，采用自適應罰Heaviside 函數來代替Heaviside 函數，其具體表達式為

式中：2λ 為自由面附近由干區到濕區的過渡層厚度，在計算中可以取為穿過自由面裂隙平均長度的1/10；? 為自適應放大系數，為獲得較好的收斂性，可適當調整其大小。

5.2 有限元求解格式

將裂隙網絡中的單條裂隙看成是有限元中的線單元，單元節點即為裂隙的兩個端點，對裂隙ij，其上任意一點的水頭可以用節點i、j 的水頭表示為

式中： Ni、 Nj為形函數， Ni= 1 -l /lij， Nj= l /lij。

經過有限元近似后得到的變分不等式提法為：尋找一水頭向量 φr+1∈I，對于 ?ψ ∈，都有以下不等式成立，

其中

式中：n 表示整個區域所有裂隙節點的總數；r 和r+ 1表示迭代次數。求解離散形式的變分不等式的算法有很多，本文選用Zheng 所建議的算法[14]。

5.3 迭代過程

Signorini 型變分不等式方法求解有自由面的滲流問題時，溢出邊界的判斷和自由面的迭代是同時進行的。為表述方便，將節點總集N 劃分為4 個子集：

式中：inN 為整個裂隙網絡區域內部節點集合。

基于上述離散型VI 提法的無壓滲流分析方法的有限元迭代算法如下：

步驟1：構造總體滲透傳導矩陣K ，設對邊界Γs上的 ?i，都有 i ∈。求解方程 Kφ1=0，得到φ1。

步驟2：令 1r=

（3）如果 |φr+1- φr||＜ ε|| φr||，則轉入步驟3，否則令 r = r+ 1，轉入（1）。

步驟3：令φ =φr，計算結束。

6 算 例

基于本文所建立的裂隙網絡無壓滲流分析的變分不等式提法，采用Visual C++語言和面向對象的程序設計方法，研制開發了相應的計算程序FracSeepage，并用多個算例對程序進行了考察，驗證了程序的有效性和魯棒性。

6.1 均質矩形壩的滲流分析

如圖4 所示，有一均質矩形壩，壩高為12 m，寬為10 m，底部邊界隔水，上游水位為10 m，下游水位為2 m。采用間距為0.2 m，有著常水力開度的兩組正交裂隙形成的裂隙網絡來模擬均質矩形壩的滲流，則裂隙的等效隙寬可以表示為

式中：k 為壩體材料的滲透系數，B 為裂隙間距。計算中選取均質壩體的滲透系數 k = 4.13× 10-3mm/s，則裂隙的等效隙寬 b= 0.1 mm。

對于該算例，圖4 給出了自由面的經驗解及裂隙網絡滲流數值計算結果，其中經驗解為

對于滲流溢出點，經驗解為z=4.47 m，數值解為z=4.27 m。表明數值解與經驗解較為吻合，但數值解的自由面降落略大于經驗解。

圖4 均質矩形壩滲流自由面位置的比較 Fig.4 Comparison of locations of the free surface from different approaches

根據Dupuit 假定，均質壩體滲流的總流量可以表示為

式中：1h 和2h 分別表示上下游水位，L 表示壩體寬度。

采用 Dupuit 公式得到壩體滲流量為1.98 × 10-5m2/s，基于裂隙網絡滲流分析的數值解為 1.95 × 10-5m2/s，兩者誤差僅為1.79%，表明數值解的計算結果是合理可靠的。

6.2 復雜裂隙網絡滲流分析

為了研究裂隙巖體的滲透特性和水力耦合特性， 國 際 合 作 項 目 Decovalex[15]第 5 期（Decovalex-2011）的Task 提供了一個含有7 797裂隙、29 814 個節點的復雜裂隙網絡模型即DFN（discrete fracture network），如圖5 所示。產生該DFN 模型的裂隙參數來自英國Sellafield 地區的現場地質測量[16]，發育有4 組裂隙，裂隙產狀服從Fisher 分布，裂隙跡長服從負指數分布，表1 給出了裂隙的幾何參數及統計分布規律。

圖5 裂隙網絡模型 Fig.5 Geometry of fracture system in the DFN model

表1 裂隙幾何參數 Table 1 Fracture parameters used for DFN generation

本文采用該DFN 模型進行裂隙網絡無壓滲流分析，假定模型左側水頭邊界為10 m，右側水頭邊界為20 m，底部為不透水邊界。為了分析裂隙的水力開度對裂隙網絡滲流的影響，計算考慮了兩種情況：（1）所有裂隙的隙寬為常數65 um；（2）裂隙的隙寬服從對數正態分布且裂隙跡長與隙寬相關，兩者的關系函數為[17]

式中：minl 和maxl 為最小和最大的裂隙跡長；D 為分形維數；g 為截斷誤差函數；imb 和inb 為裂隙隙寬的上、下限。

圖6給出了DFN模型在兩種不同水力開度分布條件下裂隙網絡內的流量分布及自由面位置，圖7給出了模型左邊界裂隙出口的流量分布。從圖6、7可以看出，對于所有裂隙隙寬為常水力開度的情況，裂隙網絡內部以及裂隙在下游邊界出口的流量分布相對比較均勻，而且自由面形狀變化平緩；對于裂隙隙寬服從對數正態分布且與跡長相關的情況，裂隙網絡內部及出口的流量分布主要由隙寬較大、跡長較長的裂隙所控制，優勢流的效應比較明顯，充分反映了長大裂隙的導水作用。同時，受裂隙隙寬隨機分布影響，自由面形狀在局部起伏較大。圖8還給出了裂隙網絡內部水頭等勢線的分布，可以看出，水頭自裂隙網絡右側向左側隨整體水流方向逐步遞減，與連續介質滲流規律基本一致。

圖6 裂隙網絡內部流量分布 Fig.6 Flow rates distribution in the DFN model

圖7 下游邊界出口裂隙流量分布 Fig.7 Flow rates distribution in each fracture intersecting the left vertical boundary of the model

圖8 裂隙網絡水頭等勢線分布 Fig.8 The water head contours in unit of m

7 結 論

（1）引入初流速來抵消在裂隙網絡干區實際不存在的流速，將Darcy 定律擴展到整個區域，定義了在整個區域上的非線性邊值問題，并將潛在溢出面邊界條件歸納為Signorini 型邊界條件，建立等價的變分不等式提法，從而在理論上解決了出滲點的定位問題。

（2）通過使用連續型的自適應Heaviside 函數 代替階梯型的Heaviside 函數，避免了在自由面附近由干區到濕區的過渡區域進行積分時形成的數值跳躍，從而顯著提高了在精確定位自由面和溢出點的數值求解過程中的穩定性和程序的魯棒性。

（3）來自均質壩體和Decovalex Task 的復雜裂隙網絡算例驗證了本文算法的有效性，同時表明，該算法對復雜裂隙網絡具有良好的適應性，并能反映裂隙網絡滲流的結構特征。

[1] 張有天， 陳平， 王鐳. 有自由面滲流分析的初流量法[J]. 水利學報， 1988，8(1)： 18－26. ZHANG You-tian， CHEN Ping， WANG Lei. Initial flow method for seepage analysis with free surface[J]. Journal of Hydraulic Engineering， 1988， 8(1)： 18－26.

[2] 毛昶熙. 滲流計算分析與控制[M]. 北京： 中國水利水電出版社， 2003.

[3] LONG J C S， REMER J S， WILSON C R， et al. Porous media equivalents for networks of discontinuous fractures[J]. Water Resources Research， 1982， 18(3)： 645－658.

[4] 毛昶熙， 陳平， 李祖貽， 等. 裂隙巖體滲流計算方法研究[J]. 巖土工程學報， 1991， 13(6)： 1－10. MAO Chang-xi， CHEN Ping， LI Zu-yi， et al. Research on fractured rock mass seepage calculation method[J]. Chinese Journal of Geotechnical Engineering， 1991， 13(6)： 1－10.

[5] LONG J C S， GILMOUR P， WITHERSPOON P A. A method for steady fluid flow in random three-dimensional networks of disc-shaped fractures[J]. Water Resources Research， 1985， 21(8)： 1105－1115.

[6] 張奇華， 鄔愛清. 三維任意裂隙網絡滲流及其解法[J]. 巖石力學與工程學報， 2010， 29(4)： 720－730. ZHANG Qi-hua， WU Ai-qing. Three-dimensional arbitrary fracture network seepage model and solution[J]. Chinese Journal of Rock Mechanics and Engineering， 2010， 29(4)： 720－730.

[7] 王嬡， 速寶玉. 裂隙巖體滲流模型綜述[J]. 水科學進展， 1996， 7(3)： 276－282. WANG Yuan， SU Bao-yu. Comment on the models of seepage flow in fractured rock masses[J]. Advances in Water Science， 1996， 7(3)： 276－282.

[8] BARENBLATT G I. Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks[J]. Journal of Applied Mathematical Mechanics， 1960， 24 (5)： 1286－1303.

[9] JING L， MA Y， FANG Z. Modeling of fluid flow and solid deformation for fractured rocks with discontinuous deformation analysis (DDA) method[J]. International Journal of Rock Mechanics and Mining Sciences， 2001， 38(3)： 343－355.

[10] 王恩志. 剖面二維裂隙網絡滲流計算方法[J]. 水文地質工程地質， 1993， 20(4)： 27－29. WANG En-zhi. Seepage calculation method in fissure networks on vertical section[J]. Hydrogeology & Engineering Geology， 1993， 20(4)： 27－29.

[11] 柴軍瑞， 仵彥卿. 巖體裂隙網絡滲流自由面的確定方法[J]. 工程勘察， 2000， 1： 23－24. CHAI Jun-rui， WU Yan-qing. The method for determination of the position of free surface in fractured rock masses[J]. Geotechnical Investigation & Surveying， 2000， 1： 23－24.

[12] ZHENG H. A new formulation of Signorini’s type for seepage problems with free surfaces[J]. International Journal for Numerical Methods in Engineering， 2005， 64(1)： 1－16.

[13] CHEN Y， ZHOU C， ZHENG H. A numerical solution to seepage problems with complex drainage systems[J]. Computers and Geotechnics， 2008， 35(3)： 383－393.

[14] ZHENG T S， LI L， XU Q Y. An iterative method for the discrete problems of a class of elliptical variational inequalities[J]. Applied Mathematics and Mechanics， 1995， 16(4)： 351－358.

[15] NIRE X. Evaluation of heterogeneity and scaling of fractures in the Borrowdale Volcanic Group in the Sellafield area[R]. UK： [s. n.]， 1997.

[16] JING L， TSANG C F， STEPHANSSON O. DECOVALEX－An international cooperative research project on mathematical models of coupled THM processes for safety analysis of radioactive waste repositories[J]. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts， 1995， 32(5)： 389－398.

[17] ALIREZA B， LANRU J. Hydraulic properties of fractured rock masses with correlated fracture length and aperture[J]. International Journal of Rock Mechanics and Mining Sciences， 2007， 44(5)： 704－719.