壩基巖體裂隙滲流效應(yīng)數(shù)值模擬方法

蔣中明肖喆臻唐 棟

（1. 長(zhǎng)沙理工大學(xué) 水利工程學(xué)院，湖南 長(zhǎng)沙 410114；2. 長(zhǎng)沙理工大學(xué) 水沙科學(xué)與水災(zāi)害防治湖南省重點(diǎn)實(shí)驗(yàn)室，湖南 長(zhǎng)沙 410114；3. 長(zhǎng)沙理工大學(xué) 洞庭湖水環(huán)境治理與生態(tài)修復(fù)湖南省重點(diǎn)實(shí)驗(yàn)室， 湖南 長(zhǎng)沙 410114）

1 研究背景

1959年法國(guó)67 m高的Mal-passet薄拱壩與1963年意大利265 m高的Vajiont拱壩[1]均因忽視巖體中水的滲流效應(yīng)影響而引發(fā)了嚴(yán)重后果。在此之后，滲流效應(yīng)對(duì)大壩壩基、地下洞室以及邊坡等工程安全性的影響逐漸為大家所重視。為了分析巖體滲流對(duì)工程安全性的影響，眾多學(xué)者提出并發(fā)展了等效連續(xù)介質(zhì)模型、離散裂隙網(wǎng)絡(luò)模型和雙重介質(zhì)模型等裂隙巖體滲流模型。

等效連續(xù)介質(zhì)模型簡(jiǎn)單、理論完善，由Snow[2]提出，并由Oda、Neuman和Depner等學(xué)者[3-5]逐漸完善。該方法適合于裂隙密集的工程巖體滲流問(wèn)題分析。但是對(duì)于裂隙稀疏，或巖體滲透性樣本單元體積REV不存在的情況，等效連續(xù)介質(zhì)模型并不適用；周志芳等[6]研究了優(yōu)勢(shì)裂隙對(duì)滲流的影響，譚春等[7]提出使用灰色理論確定巖體表征單元體。此外，等效連續(xù)介質(zhì)模型只能大概反映區(qū)域流量的變化，并不能真實(shí)反映裂隙水的滲流路徑和小于REV區(qū)域滲流場(chǎng)的分布情況。

離散裂隙網(wǎng)絡(luò)模型假定滲流只發(fā)生在裂隙中，忽略了孔隙滲流。離散裂隙網(wǎng)絡(luò)模型由于忽略了巖塊的孔隙滲流，所以不適用于巖塊基質(zhì)滲透性較大的情況。另外，現(xiàn)階段三維裂隙網(wǎng)絡(luò)滲流的求解難度大，耗時(shí)長(zhǎng)，對(duì)解決工程尺度滲流問(wèn)題仍有相當(dāng)難度。裂隙滲流的求解方法也多種多樣，例如Dverstorp[8]假定巖體中滲流路徑是在裂隙中光滑的管道中進(jìn)行。N.Koudina[9]、李新強(qiáng)[10]則基于邊界元法求解三維多邊形裂隙網(wǎng)絡(luò)滲流問(wèn)題。何忱等[11]將巖石看作由不透水的四面體塊體組成的集合，滲流只在相鄰塊體之間的界面上發(fā)生，用界面網(wǎng)絡(luò)代替了巖石的孔隙結(jié)構(gòu)，裂隙也以界面形式顯式存在于巖體之中，并遵循Darcy定理。李海楓等[12]也采用Delaunay三角剖分技術(shù)生成離散裂隙網(wǎng)絡(luò)模型以研究裂隙巖體非飽和滲流問(wèn)題。

雙重介質(zhì)模型最初由Barenblatt等[13]提出，是另一種等效連續(xù)體模型。其基本思路是將孔隙介質(zhì)的連續(xù)體單元和裂隙網(wǎng)絡(luò)的離散裂隙視為同一空間內(nèi)可以相互交換流體的連續(xù)體。Aifantis[14]基于混合物理論提出并推導(dǎo)了流固耦合方程，但是其部分方程是基于物理意義得出而非數(shù)學(xué)推導(dǎo)。由于雙重介質(zhì)模型沒(méi)有從正面考慮離散裂隙與單元之間的連接問(wèn)題，所以導(dǎo)致節(jié)點(diǎn)上出現(xiàn)兩個(gè)水頭的情況，高海鷹和吳宏明[15]通過(guò)假設(shè)裂隙與單元相交的面水頭相同、流量不連續(xù)建立了混合模型，并對(duì)裂隙與單元連接處做了討論。黎水泉等[16]經(jīng)過(guò)數(shù)學(xué)推導(dǎo)，得到了雙重介質(zhì)模型的流固耦合方程。

為了綜合利用等效連續(xù)介質(zhì)模型和離散裂隙網(wǎng)絡(luò)模型二者各自的優(yōu)點(diǎn)，侯曉萍和徐青等[17]采用復(fù)合單元法研究了裂隙巖體非穩(wěn)定滲流分析方法，其主要思路是將處于連續(xù)介質(zhì)單元內(nèi)部的離散裂隙通過(guò)拓?fù)溆?jì)算，在單元內(nèi)部形成基于裂隙面的若干子單元，并通過(guò)建立平衡方程后統(tǒng)一求解；該方法對(duì)于密集裂隙情況下的滲流計(jì)算仍需要進(jìn)行發(fā)展。王臻等[18]用Oda的裂隙張量理論來(lái)獲得子區(qū)域的滲透張量，然后采用DC模型（離散連續(xù)模型）分析二維滲流問(wèn)題，取得了良好效果。近年來(lái)，國(guó)外部分學(xué)者提出了EFC模型（Embedded Fractured Continuum）[19-22]，即在連續(xù)介質(zhì)單元中考慮了裂隙影響來(lái)模擬裂隙巖體的滲流效應(yīng)。以上方法在本質(zhì)上還是屬于連續(xù)介質(zhì)滲流的分析方法。本文在上述方法的基礎(chǔ)上，通過(guò)對(duì)巖體中裂隙滲流特性分析，將單元中的裂隙滲透性按各向異性處理，并將之與塊體滲透性進(jìn)行疊加，從而實(shí)現(xiàn)用連續(xù)介質(zhì)模型模擬巖體裂隙滲流效應(yīng)的數(shù)值計(jì)算方法。

2 裂隙介質(zhì)EFE滲流模型

嵌入裂隙單元（Embedded Fractured Element）的滲透性由單元中的裂隙數(shù)量、各裂隙開(kāi)度以及裂隙周圍的巖塊滲透性共同確定。對(duì)大多數(shù)巖體而言，由于巖塊的滲透性相對(duì)較小，故其滲透性往往由巖體中包含的裂隙數(shù)量及其幾何和充填特征所控制。圖1所示的六面體單元中包含了2條水平貫通裂隙和1條鉛直貫通裂隙。如果巖塊和裂隙滲透性分別用Kr和Kf表示，單元在3個(gè)方向上的滲透系數(shù)分別用kxx、kyy和kzz表示，對(duì)于裂隙巖體，大多數(shù)情況下都滿足Kr ＜＜Kf條件，因此單元滲透性主要由Kf貢獻(xiàn)。具體地說(shuō)，圖1中的單元在z方向上的滲透性主要由f3控制；單元在x方向上的滲透性由f1和f2共同決定；單元在y方向上的滲透性則受到f1、f2和f3共同影響。假定圖1中的各裂隙開(kāi)度相同，則存在kxx =2 kzz和kyy =3 kzz關(guān)系。這表明裂隙的存在使得單元呈現(xiàn)出了顯著各向異性特點(diǎn)。含裂隙單元滲透性的各向異性特性可以采用各向異性連續(xù)介質(zhì)滲流模型來(lái)刻畫(huà)。

圖1 裂隙單元示意

大量試驗(yàn)研究表明，對(duì)于不可壓縮流體，巖塊（孔隙介質(zhì)）中流體運(yùn)動(dòng)服從達(dá)西定律[23]。當(dāng)縫寬較小，且水流流態(tài)為層流時(shí)，流體沿裂隙面流動(dòng)也服從達(dá)西定律。所以對(duì)于巖體來(lái)說(shuō)在恒定流情況下，滲流方向上的總流量Q由裂隙流量Qf和巖塊流量Qr組成：

式中：Q為滲流方向上的總流量；Qf為滲流方向上巖體裂隙的流量；Qr為滲流方向上巖體巖塊的流量。于是在水力梯度相同的情況下，可得整體等效滲透系數(shù)張量Kij為：

式中：A為滲流方向截面面積；Ar為滲流方向巖塊所占面積；Af為滲流方向裂隙橫截面面積；Kr為巖塊滲透系數(shù)；Kfij為裂隙滲透系數(shù)。由于裂隙開(kāi)度很小，一般情況下，存在以下關(guān)系：

于是式（2）可改寫(xiě)為：

式中：Lfi為單位面積上第i條裂隙的跡長(zhǎng)；bi為第i條裂隙的開(kāi)度；nf為滲流方向面積裂隙率，其中角度θ為裂隙面與滲流方向橫截面的夾角。

當(dāng)計(jì)算單元中沒(méi)有包含裂隙時(shí)，由式（6）可知計(jì)算單元滲透系數(shù)Kij =Kr，計(jì)算域內(nèi)巖體滲透性完全可以用孔隙連續(xù)介質(zhì)來(lái)描述。孔隙連續(xù)介質(zhì)可以各向同性材料，也可以是各向異性材料。工程分析時(shí)，由于計(jì)算單元的尺度遠(yuǎn)遠(yuǎn)大于單元內(nèi)的空隙尺度，因此絕大多數(shù)情況下采用各向同性的假定是合理的。

當(dāng)計(jì)算單元中包含裂隙時(shí)，單元滲透系數(shù)由Kr和Kf,ij共同決定。當(dāng)計(jì)算單元尺度大，包含裂隙數(shù)量少時(shí)，盡管Kf,ij一般相對(duì)較大，但由于式（6）中的nf相對(duì)較小，單元滲透性可能主要由巖塊滲透性決定；反之，當(dāng)計(jì)算單元尺度較小，nf相對(duì)較大，故單元滲透性主要由裂隙滲透性決定。

由于Kf,ij具有顯著的各向異性特點(diǎn)，因此計(jì)算域內(nèi)巖體的滲透性必然也將呈現(xiàn)出顯著各向異性特點(diǎn)，即含裂隙單元為完全各向異性材料。滲透張量可通過(guò)立方定律確定。采用Snow[2]的假設(shè)，即使不同方向裂隙組在裂隙網(wǎng)絡(luò)系統(tǒng)中相互連通，不同方向裂隙組內(nèi)的裂隙水流互不干擾，且裂隙貫穿整個(gè)單元，單元各個(gè)方向上的流量等于各裂隙在相應(yīng)方向上的流量和，即流量滿足疊加原理。于是，滲透張量Kf,ij的矩陣表達(dá)式如下[2]：

式中：Kf,ij為單元中包含的裂隙對(duì)計(jì)算單元滲透張量的貢獻(xiàn)；m為裂隙個(gè)數(shù)；g為重力加速度；b為裂隙開(kāi)度；υ為流體運(yùn)動(dòng)黏滯系數(shù)；αxi、αyi、αzi為第i個(gè)裂隙面單位法向向量在3個(gè)方向上的分量，可由裂隙面的傾向及傾角求出。

由此可見(jiàn)，式（6）既適用于單元中有裂隙的情況，也適用于單元中無(wú)裂隙分布情況。因此，采用式（6）對(duì)計(jì)算域內(nèi)不同部位的巖體分別采用不同屬性的滲透系數(shù)就可以真實(shí)地反映巖體中的裂隙滲流效應(yīng)（圖2中的裂隙等效單元，即EFE單元）以及塊體內(nèi)的孔隙滲流效應(yīng)（圖2中的基質(zhì)單元）。由于裂隙面寬度都很小，為了盡可能降低裂隙等效單元尺度過(guò)大對(duì)裂隙附近范圍的滲流形態(tài)的影響，可以采用如圖2中的網(wǎng)格加密技術(shù)對(duì)裂隙穿過(guò)的單元進(jìn)行網(wǎng)格加密，從而減少裂隙滲流單元的數(shù)量。

圖2 裂隙巖體網(wǎng)格剖分

3 EFE滲流模型的數(shù)值實(shí)現(xiàn)

對(duì)于工程尺度的滲流分析來(lái)說(shuō)，采用完全連續(xù)孔隙介質(zhì)模型或完全離散裂隙網(wǎng)絡(luò)模型都是不可接受的。在實(shí)際工程中，巖體中分布的裂隙尺度大小不均，如果全部納入到離散裂隙網(wǎng)絡(luò)模型中，其計(jì)算工作量將是不可接受的。

一種可行的方案是將一些小尺度的裂隙空間納入到孔隙空間范圍，并將其按孔隙介質(zhì)滲流（基質(zhì)單元/孔隙單元）進(jìn)行處理，而尺度大的裂隙及結(jié)構(gòu)面等按等效裂隙滲流單元（EFE）進(jìn)行處理是合理的。式（6）剛好具備同時(shí)描述裂隙滲流特性和孔隙滲流特性的能力。

利用式（6）對(duì)分析域同時(shí)進(jìn)行裂隙滲流特性和孔隙滲流特性分析的前提是計(jì)算軟件需要同時(shí)具備三維裂隙網(wǎng)絡(luò)生成及處理能力以及連續(xù)介質(zhì)滲流分析能力。在計(jì)算過(guò)程中，對(duì)所有單元根據(jù)試驗(yàn)值綜合確定其滲透系數(shù)。對(duì)存在裂隙的單元，根據(jù)裂隙的分布情況利用式（8）自動(dòng)計(jì)算其各向異性滲透系數(shù)，并對(duì)計(jì)算單元滲透系數(shù)進(jìn)行修正；對(duì)于不包括貫通裂隙的單元，則不需進(jìn)行滲透性修正。由此可見(jiàn)，上述利用裂隙等效單元獲取裂隙滲流效應(yīng)的方法只需在原有的滲流分析程序中添加包含裂隙的單元滲透系數(shù)處理模塊就可以實(shí)現(xiàn)。

圖3 計(jì)算流程圖

為實(shí)現(xiàn)含裂隙的大規(guī)模工程巖體滲流分析功能，基于FLAC3D5.01版本提供的DFNs（discrete fracture networks）功能，利用上述技術(shù)思路，二次開(kāi)發(fā)了能同時(shí)進(jìn)行滲流分析的連續(xù)孔隙單元和裂隙等效單元的FISH程序。需要說(shuō)明的是，盡管FLAC3D5.01版本提供了DFNs幾何信息生成功能，但該軟件平臺(tái)卻沒(méi)有提供如何利用DFNs進(jìn)行單元離散性質(zhì)模擬的相關(guān)模塊或命令，因此，需要進(jìn)行二次程序開(kāi)發(fā)。程序計(jì)算流程如圖3所示。（1）建立研究域內(nèi)的實(shí)體單元模型（孔隙單元）；（2）確定離散裂隙的生成域，在同一坐標(biāo)系下通過(guò)蒙特卡羅法生成離散裂隙網(wǎng)絡(luò)模型；（3）根據(jù)給定的裂隙等效單元精度參數(shù)和加密參數(shù)，對(duì)裂隙穿過(guò)單元網(wǎng)格進(jìn)行加密剖分；（4）將與裂隙接觸的單元標(biāo)記為裂隙等效單元，并儲(chǔ)存裂隙的幾何信息；（5）裂隙等效單元使用各向異性滲流模型，計(jì)算域內(nèi)其余的孔隙介質(zhì)單元使用各向同性滲流模型；（6）若單元為裂隙等效單元，根據(jù)單元中儲(chǔ)存的裂隙幾何信息按式（8）計(jì)算裂隙對(duì)滲透張量的貢獻(xiàn)，然后利用式（6）計(jì)算單元整體滲透張量；若單元為孔隙介質(zhì)單元，直接使用孔隙介質(zhì)的滲透系數(shù)。單元滲透性計(jì)算完成后對(duì)單元進(jìn)行滲透系數(shù)賦值；（7）最后施加滲流模型初始條件和邊界條件，進(jìn)行滲流計(jì)算。

需要說(shuō)明的是，在數(shù)值模擬過(guò)程中，往往需要同時(shí)考慮計(jì)算精度和計(jì)算效率的問(wèn)題，對(duì)于變化平緩的非重點(diǎn)區(qū)域需要布置稀疏的網(wǎng)格，而在變化劇烈的重點(diǎn)部位則需要布置密集的網(wǎng)格。因此，在程序分析過(guò)程增加了網(wǎng)格加密剖分環(huán)節(jié)。圖4為裂隙單元加密示意圖。圖中紅色圓盤為離散裂隙模型，灰色區(qū)域?yàn)榱严秵卧?，白色區(qū)域?yàn)闊o(wú)裂隙的孔隙介質(zhì)單元。由圖4（a）（d）可知，裂隙等效單元加密前的形狀為正方形，與離散裂隙的圓盤形不符；加密后的等效單元形狀更加接近于圓盤形。圖4（c）（f）則表明，隨著裂隙單元的加密，裂隙單元的形狀在走向方向上的形狀會(huì)更精確，所以其滲透張量誤差也會(huì)越小。由此可見(jiàn)，對(duì)裂隙等效單元進(jìn)行局部加密，可以有效地減小網(wǎng)格剖分帶來(lái)的誤差。

圖4 裂隙單元加密示意

4 程序及模型驗(yàn)證

為評(píng)價(jià)本文提出模型的正確性，分別以具有單一矩形貫穿裂隙巖體和文獻(xiàn)[9]中的算例來(lái)進(jìn)行裂隙巖體的滲透性和單一圓形非貫穿裂隙巖體滲透張量的驗(yàn)證。

圖5 算例1計(jì)算網(wǎng)格

算例1：研究域?yàn)橐粋€(gè)10 m×10 m×10 m的方形區(qū)域，巖塊（基質(zhì)）計(jì)算滲透系數(shù)取1×10-10 m/s，假設(shè)巖體中心有一個(gè)貫穿裂隙，隙寬為0.1 mm。如圖5，裂隙面與立方體區(qū)域中心重合，且裂隙面法向量垂直于x軸與z軸。令水力梯度為1，與水力梯度方向垂直的兩面為定水頭邊界，其余面為不透水邊界。

沿x方向上的總流量為裂隙流量和孔隙流量之和，其中裂隙流量Qf根據(jù)立方定律[2]求出：

式中：Lf為裂隙的跡長(zhǎng)；Jf為沿x方向上的水力梯度?？紫读髁縌r根據(jù)達(dá)西定律求出。計(jì)算得到沿x方向上的裂隙流量Qf為8.250825×10-2 m3/s，孔隙流量Qr為1×10-8 m3/s，所以沿x方向上的總流量Q為8.250826×10-2 m3/s。按照本文研究方法通過(guò)數(shù)值方法得到模型在x方向上的流量為8.250826×10-2 m3/s，在保留小數(shù)點(diǎn)后6位有效數(shù)字的情況下，數(shù)值計(jì)算結(jié)果和解析解結(jié)果相同?？梢钥闯鲈趲r塊透水性較差的情況下裂隙單元中孔隙流量對(duì)總流量的貢獻(xiàn)基本可以忽略。巖體滲流的優(yōu)勢(shì)通道為連通裂隙通道，因此，本文提出的考慮裂隙滲流效應(yīng)的巖體滲流分析方法較常規(guī)連續(xù)介質(zhì)滲流分析方法更加能反映裂隙滲流的優(yōu)勢(shì)效應(yīng)。

圖6 算例2幾何模型及邊界條件示意

算例2：假設(shè)邊長(zhǎng)為20 m立方形巖體中心有一個(gè)直徑為10 m，隙寬為1 mm的圓盤形裂隙，巖塊（基質(zhì)）計(jì)算滲透系數(shù)取1×10-8m/s。圓盤圓心始終與立方體區(qū)域中心重合，裂隙面法向量始終垂直于z軸，定義α角為y軸與裂隙面法向量形成的夾角，計(jì)算α角分別為0°、15°、30°、45°、60°、75°以及90°時(shí)的巖體整體滲透張量。令水力梯度為1，與水力梯度方向垂直的兩面為定水頭邊界，其余面上根據(jù)與兩定水頭邊界面的位置關(guān)系設(shè)置為線性變化的變水頭邊界（見(jiàn)圖6）。

根據(jù)文獻(xiàn)[11]可知該問(wèn)題中的滲透張量Kf為角度α的函數(shù)，且kzz不隨α角的變化而變化，Kf可以表示為：

根據(jù)文獻(xiàn)[1]可知滲透張量各滲透系數(shù)隨坐標(biāo)軸旋轉(zhuǎn)而變化，將坐標(biāo)系xoy順時(shí)針旋轉(zhuǎn)α角，其轉(zhuǎn)換關(guān)系為：

式（10）中的各個(gè)元素取值可通過(guò)式（11）得到。圖7為算例2所述問(wèn)題的解析解與數(shù)值解對(duì)比。依據(jù)流量等效和邊界相同的原則，求得各個(gè)面的滲流量后，研究域的滲透張量可參考文獻(xiàn)[24]方法確定。由圖可知，滲透張量的理論解和數(shù)值解在整體上吻合度很好，說(shuō)明本文采用的含裂隙介質(zhì)的滲流分析方法是合理可行的。

5 單元尺度大小對(duì)研究域等效滲透性影響分析

假設(shè)邊長(zhǎng)為20 m立方形巖體中心有一個(gè)直徑為10 m，隙寬為1 mm的圓盤形裂隙。圓盤圓心與立方體區(qū)域中心重合，裂隙面法向量與Y軸平行。單元未加密時(shí)裂隙單元邊長(zhǎng)為2 m，加密后使裂隙單元邊長(zhǎng)分別為1、0.5和0.25 m。巖塊（基質(zhì)）計(jì)算滲透系數(shù)取1×10-8 m/s。通過(guò)比較不同加密程度后研究區(qū)域的整體滲透性，來(lái)研究加裂隙附近單元的加密程度對(duì)滲透性的影響。

圖7 滲透張量數(shù)值解與理論解對(duì)比

圖8 剖面y=10m處單元尺度對(duì)比

在圖8中白色區(qū)域?yàn)榛|(zhì)單元，深灰色區(qū)域?yàn)榱严秵卧咨珗A圈為圓盤離散裂隙。隨著裂隙單元邊長(zhǎng)不斷縮小，裂隙位置處的灰色區(qū)域形狀逐漸從方形接近圓形，即裂隙影響的基質(zhì)單元范圍越來(lái)越小。表1給出了包含裂隙的單元尺寸變化對(duì)應(yīng)的研究域計(jì)算滲透系數(shù)變化情況，其中研究域的計(jì)算滲透系數(shù)是根據(jù)各個(gè)面的滲流量，通過(guò)達(dá)西定律計(jì)算得到的整個(gè)研究域的滲透系數(shù)。圖9為依據(jù)表1繪制的計(jì)算與滲透性與單元尺度變化關(guān)系。由表1和圖9可知，由于裂隙沒(méi)有穿透基質(zhì)巖石，故研究域綜合滲透性主要受基質(zhì)滲透性控制；但裂隙的存在使得水流在x和z方向的滲透性出現(xiàn)了一定程度的增加，但其計(jì)算滲透性的增加幅度與穿越裂隙單元的尺寸大小相關(guān)。隨著裂隙處單元邊長(zhǎng)不斷減小，計(jì)算滲透系數(shù)基本呈線性減少的趨勢(shì)。含裂隙單元尺寸越小，得到的研究域內(nèi)計(jì)算滲透系數(shù)約接近于真實(shí)值。單元尺度由0.25 m增加到2.0 m后，兩者間的計(jì)算誤差可達(dá)到19.1%。理論上，在裂隙未貫穿的情況下，研究域內(nèi)整體的滲透性均由基質(zhì)巖塊的滲透性控制。需要說(shuō)明的是，對(duì)于實(shí)際巖塊被裂隙貫穿的單元，其滲透張量計(jì)算值不受計(jì)算單元尺度大小的影響。因此，在計(jì)算中為盡量減少誤差，可只需對(duì)裂隙邊界單元采用較小的單元尺度即可。

圖9 滲透系數(shù)與單元尺度關(guān)系

表1 計(jì)算滲透系數(shù)與單元尺度關(guān)系表

6 壩基裂隙巖體滲流效應(yīng)分析

對(duì)于重力壩壩基而言，壩基巖體中分布的裂隙對(duì)于壩基揚(yáng)壓力與滲流路徑的空間分布規(guī)律有著重要的影響。為了解壩基巖體中的裂隙分布對(duì)滲流場(chǎng)的影響，擬定如圖10所示的三種計(jì)算模型進(jìn)行分析。模型一中的裂隙總體傾向下游；模型二壩基巖體中的裂隙總體傾向上游；模型三壩基巖體中不考慮裂隙影響。壩基巖體中裂隙統(tǒng)計(jì)數(shù)據(jù)如表2，隙寬均為1 mm。模型一和模型二中裂隙穿過(guò)的計(jì)算單元為各向異性滲流單元，其余為各向同性滲流單元；模型三中的計(jì)算單元均為各向同性滲流單元。含裂隙單元的滲透系數(shù)按式（6）和式（8）由程序自動(dòng)計(jì)算，壩基巖塊滲透系數(shù)取1×10-7m/s?？紤]到壩體混凝土滲透性相對(duì)較小，本算例分析時(shí)視為不透水介質(zhì)。三個(gè)計(jì)算模型的上游水頭均為40 m，下游水頭均為10 m，計(jì)算壩基的穩(wěn)定性滲流場(chǎng)。

表2 裂隙統(tǒng)計(jì)參數(shù)表

圖10 模型計(jì)算網(wǎng)格圖

圖11給出了邊界條件相同情況下壩基巖體中的部分流線分布圖。由圖可知，壩基中存在導(dǎo)水裂隙時(shí)，水流路徑主要沿裂隙通道分布。由此可見(jiàn)，考慮裂隙分布影響的ELE模型能夠反映壩基巖體滲流的真實(shí)路徑。裂隙產(chǎn)狀對(duì)壩基孔隙壓力的空間分布形態(tài)也存在顯著影響，如圖12所示。裂隙分布狀態(tài)對(duì)壩基揚(yáng)壓力分布規(guī)律也將產(chǎn)生不同程度的影響，如圖13所示，裂隙分布的隨機(jī)性是導(dǎo)致壩基揚(yáng)壓力分布的非規(guī)則性的直接因素?？紤]裂隙分布隨機(jī)性導(dǎo)致的壩基揚(yáng)壓力分布的非規(guī)則性與壩基揚(yáng)壓力的實(shí)測(cè)分布更為接近?？傊?，EFE模型可以很好地模擬離散裂隙的局部導(dǎo)水和疏水性能，從而真實(shí)有效地反應(yīng)壩基裂隙巖體的滲流狀態(tài)。

圖11 壩基巖體中流網(wǎng)圖

圖12 壩基巖體中水壓力圖（單位：1×104Pa）

圖13 壩基揚(yáng)壓力分布

7 結(jié)論

本文提出嵌入裂隙單元模型，并在FLAC3D軟件平臺(tái)上完成了程序?qū)崿F(xiàn)?；诹严稁r體滲流理論研究的基礎(chǔ)上，提出并實(shí)現(xiàn)了利用連續(xù)介質(zhì)各向異性滲流模型模擬裂隙滲流特性的數(shù)值分析方法，并驗(yàn)證了方法的合理性。（1）將離散裂隙網(wǎng)絡(luò)（DFNs）技術(shù)與連續(xù)介質(zhì)各向異性滲流分析技術(shù)結(jié)合，利用構(gòu)造的EFE滲流單元模擬裂隙導(dǎo)致的各向異性非均勻滲流效應(yīng)效果，能夠真實(shí)反映裂隙隨機(jī)分布所導(dǎo)致的各向異性滲透特性。（2）理論上使用EFE單元數(shù)值模型分析真實(shí)裂隙滲流問(wèn)題必然會(huì)存在一定的誤差，可以通過(guò)調(diào)整裂隙邊界處的EFE單元尺度大小來(lái)降低這種誤差。（3）本文提出的裂隙巖體滲流的連續(xù)介質(zhì)數(shù)值模擬方法可適用于工程尺度的大規(guī)模數(shù)值計(jì)算，其計(jì)算效率和精度均能滿足工程要求。