一種裂隙-連續(xù)介質(zhì)的耦合滲流計(jì)算方法

羅璨璨，宋 俊，尚高增，楊 錦，張建海

(1.四川大學(xué) 水力學(xué)與山區(qū)河流開(kāi)發(fā)保護(hù)國(guó)家重點(diǎn)實(shí)驗(yàn)室， 四川 成都 610065; 2.四川大學(xué) 水利水電學(xué)院， 四川 成都 610065)

E-mail：lccaakk@163.com

巖體中裂隙分布復(fù)雜，目前在數(shù)學(xué)上處理裂隙巖體滲流主要有：等效連續(xù)介質(zhì)模型、離散裂隙網(wǎng)絡(luò)模型和雙重介質(zhì)滲流模型。

等效連續(xù)介質(zhì)模型是由Snow[1]提出的，Long等[2]、Oda[3]、張有天[4]、田開(kāi)銘等[5-6]都對(duì)等效連續(xù)介質(zhì)模型進(jìn)行過(guò)研究。該模型對(duì)于解決裂隙密度較大的大尺度研究域問(wèn)題有較好的效果，簡(jiǎn)單實(shí)用，能滿(mǎn)足大多數(shù)工程的需要。

離散裂隙網(wǎng)絡(luò)模型認(rèn)為基巖為不透水介質(zhì)，水流只在裂隙網(wǎng)絡(luò)中運(yùn)動(dòng)，基于裂隙網(wǎng)絡(luò)幾何模型，以單裂隙立方定律為基礎(chǔ)，各裂隙交叉點(diǎn)單位時(shí)間內(nèi)流量平衡為準(zhǔn)則建立方程組求解，從而確定裂隙網(wǎng)絡(luò)中流體的真實(shí)流動(dòng)狀態(tài)。該模型近十幾年來(lái)發(fā)展迅速，其主要研究有：Wittke[7]首先提出網(wǎng)絡(luò)線(xiàn)素法，之后Wilson和Witherspoon等[8]分別采用三角形單元和線(xiàn)單元對(duì)二維裂隙網(wǎng)絡(luò)進(jìn)行模擬，提出了兩種有限元技術(shù)；對(duì)于三維問(wèn)題，Long[9]提出了三維圓盤(pán)裂隙網(wǎng)絡(luò)模型；萬(wàn)力等[10]將有限元方法引入到裂隙網(wǎng)絡(luò)滲流計(jì)算中，進(jìn)一步提出了三維裂隙網(wǎng)絡(luò)的多邊形單元滲流模型。

雙重介質(zhì)滲流模型由前蘇聯(lián)學(xué)者Barenblatt G I首先提出，Warren J E等對(duì)其進(jìn)一步發(fā)展。該模型可以分為兩種：一是擬穩(wěn)態(tài)流模型，研究代表人物有Barenblatt G I和Rott等；二是非穩(wěn)態(tài)流模型，研究代表人物有Zimmerman[11]。雙重介質(zhì)模型有許多優(yōu)點(diǎn)，如：在一定程度上描述了優(yōu)先流現(xiàn)象；考慮了裂隙系統(tǒng)和基巖系統(tǒng)的水流交換，符合實(shí)際。但是該模型也有不可忽視的缺點(diǎn)，如：基巖被裂隙分割成具有相同的大小和形狀的巖塊，過(guò)于簡(jiǎn)化；物質(zhì)交換系數(shù)難以確定，在應(yīng)用時(shí)也需要對(duì)其有效性進(jìn)行判定[12-14]。

隨著對(duì)裂隙介質(zhì)的深入研究，王月英等[15]提出了在研究區(qū)域上同時(shí)應(yīng)用多種流動(dòng)模型的耦合滲流計(jì)算方法的構(gòu)想。比如ECM-DFN橫向耦合模型，即某區(qū)域采用等效連續(xù)介質(zhì)滲流模型，相鄰區(qū)域采用離散裂隙網(wǎng)絡(luò)模型。但文獻(xiàn)[15]僅提出了該模型的設(shè)想，并未深入研究。張奇華等[16]深入研究了雙重介質(zhì)滲流模型，提出整體求解法和流量交換法兩種求解方法。胡井友等[17]驗(yàn)證了流固耦合數(shù)值模擬的有效性。徐建國(guó)等[18]探討了流固耦合效應(yīng)對(duì)圍巖穩(wěn)定性的影響。

由此可知，對(duì)于既有孔隙介質(zhì)，又有裂隙介質(zhì)的滲流區(qū)域，目前缺乏對(duì)于相鄰區(qū)域的耦合進(jìn)行分析的方法。針對(duì)裂隙巖基上修筑土石壩這種特殊滲流區(qū)域，本文對(duì)巖體裂隙區(qū)域和土體孔隙區(qū)域兩種不同介質(zhì)分別采用裂隙滲流和連續(xù)介質(zhì)滲透方法進(jìn)行模擬，進(jìn)而通過(guò)交界區(qū)節(jié)點(diǎn)水頭一致性條件，對(duì)兩個(gè)區(qū)域的滲流控制方程進(jìn)行耦合，形成整體滲透矩陣方程，建立了裂隙-連續(xù)介質(zhì)耦合滲流計(jì)算方法，并編制了裂隙-連續(xù)介質(zhì)耦合滲流分析程序，開(kāi)展了算例分析，印證了本文方法的正確性。

1 裂隙-連續(xù)介質(zhì)耦合滲流理論及程序

1.1 裂隙滲流理論

在裂隙網(wǎng)絡(luò)中，線(xiàn)單元及其兩端節(jié)點(diǎn)為基本單元，如果節(jié)點(diǎn)i為線(xiàn)單元j的一個(gè)端點(diǎn)，則稱(chēng)j單元銜接于i節(jié)點(diǎn)，節(jié)點(diǎn)i的度數(shù)為與其相銜接的線(xiàn)單元的個(gè)數(shù)。構(gòu)成閉合路徑的裂隙段集合稱(chēng)為該裂隙網(wǎng)絡(luò)的回路，其中所含裂隙段數(shù)目為回路的維數(shù)。

離散裂隙網(wǎng)絡(luò)滲流模型認(rèn)為巖塊的滲透性遠(yuǎn)遠(yuǎn)小于裂隙的滲透性，可以忽略不計(jì)。水只在裂隙網(wǎng)絡(luò)中流動(dòng)。假定選取裂隙研究域ABCD如圖1所示，其中有裂隙交叉點(diǎn)i，以i點(diǎn)為中心，作一通過(guò)各銜接線(xiàn)元中點(diǎn)的閉合曲線(xiàn)形成表征單元域。

圖1巖體裂隙網(wǎng)絡(luò)示意圖

根據(jù)質(zhì)量守恒原理，對(duì)于穩(wěn)定滲流，表征單元域內(nèi)水流運(yùn)動(dòng)方程為：

(1)

式中：qj(j=1,2,…N′)為j線(xiàn)元流進(jìn)或流出節(jié)點(diǎn)i的流量；N′為點(diǎn)i的度數(shù)，即交于i點(diǎn)的線(xiàn)元的總數(shù)；Qi為點(diǎn)i處的源匯項(xiàng)。

若裂隙研究區(qū)域內(nèi)有N個(gè)節(jié)點(diǎn)(裂隙交叉點(diǎn))，M個(gè)線(xiàn)單元(交叉點(diǎn)之間的裂隙段)，則可以組成N個(gè)形式為(1)的方程，寫(xiě)成矩陣形式為：

Aq+Q=0

(2)

式中：Q=(Q1，Q2，…，QN)T；q=(q1，q2，…，qN)T；A={aij}N×M。

A為裂隙網(wǎng)絡(luò)中反映線(xiàn)元與節(jié)點(diǎn)銜接關(guān)系的N×M階銜接矩陣[19]，矩陣中元素aij可表述為：

根據(jù)本文對(duì)裂隙回路的定義可知，已知裂隙研究區(qū)域內(nèi)有N個(gè)節(jié)點(diǎn)，M個(gè)線(xiàn)單元，則裂隙網(wǎng)絡(luò)中的回路個(gè)數(shù)至少為M-N+1，回路的個(gè)數(shù)與裂隙切割并圍繞的巖塊的個(gè)數(shù)相等。定義裂隙網(wǎng)絡(luò)中回路矩陣C為：

C={ckj}L×M

(3)

式中：L為回路總數(shù)；ckj為矩陣C中的元素，可表述為：

回路矩陣C的方向可以提前約定，其反映了裂隙段之間的關(guān)系，有以下性質(zhì)：每一行中非零元素的個(gè)數(shù)等于該行代表的回路的維數(shù)。

由銜接矩陣A與回路矩陣C的性質(zhì)可證，銜接矩陣與回路矩陣正交，故有：

ACT=0或CAT=0

(4)

由上式可以得到：

ΔH=ATH

(5)

式中：H為N×1階節(jié)點(diǎn)水頭向量；ΔH為線(xiàn)元兩端的水頭差。

根據(jù)單裂隙內(nèi)水流運(yùn)動(dòng)方程[11，20]，可推導(dǎo)單裂隙的單寬流量公式為：

(6)

流經(jīng)該裂隙的水流平均速度為：

(7)

式中：Vx，Vy為流體運(yùn)動(dòng)速度矢量；μ為流體的動(dòng)力黏滯系數(shù)；q為單寬流量；Jf為裂隙水流的水力梯度；γ為地下水的比重；b為裂隙寬度；Kf為裂隙的滲透系數(shù)。

由單裂隙滲流式(6)和式(7)可知，第j線(xiàn)元的流量為：

(8)

式中：bj,lj分別為線(xiàn)元的寬度和長(zhǎng)度；Kf為裂隙的滲透系數(shù)；ΔHj為j線(xiàn)元兩端的水頭差。

式中裂隙段j的水力傳導(dǎo)系數(shù)：

(9)

寫(xiě)成矩陣形式為：

q=Tl·ΔH

(10)

式中：

q=(q1，q2，…,qM)T； ΔH=(ΔH1，ΔH2，…,ΔHM)T；Tl=diag(Tl1，Tl2，…，TlM)

式(2)可以寫(xiě)成

(ATlAT)H=Q

(11)

1.2 連續(xù)介質(zhì)滲流理論及程序

本文主要研究二維滲流問(wèn)題，根據(jù)文獻(xiàn)[16]，二維穩(wěn)定滲流達(dá)西定律可以表達(dá)為：

(12)

式中：h=h(x,y,t)為水頭函數(shù)；kx、ky為x、y方向的滲透系數(shù)。

對(duì)于未知水頭值的邊界Γ1，如果已知該邊界上單位面積流入或流出的流量，則有流量邊界條件：

(13)

式中：n為Γ1的外法線(xiàn)方向；q(x,y)為Γ1上的已知流量。若邊界不透水，則q=0。

基于變分原理，式(12)可以等價(jià)為下列泛函取極小值：

(14)

式中：μs=ρg(α+nβ)；ρ為滲透水的密度；α為土體壓縮系數(shù)；n為土體的孔隙率；β為水的壓縮系數(shù)。

將滲流區(qū)域Ω離散化，剖分為m個(gè)互不相交的單元體e，單元含有M個(gè)節(jié)點(diǎn)，取插值函數(shù)為Ni，則單元體內(nèi)任一點(diǎn)的水頭表達(dá)式為：

(15)

依次對(duì)式(14)中各項(xiàng)求導(dǎo)并求最小值，再令所有單元的泛函的微分為0，最終將求解方程組變?yōu)椋?/p>

[K]{h}={F}

(16)

式中：{F}為已知常數(shù)項(xiàng)；[K]表示形式如下：

(17)

整體滲透矩陣中元素和單元滲透矩陣中元素的對(duì)應(yīng)關(guān)系如下：

(18)

式中：mi是共用節(jié)點(diǎn)i,j的單元的個(gè)數(shù)，由于每個(gè)節(jié)點(diǎn)的相關(guān)節(jié)點(diǎn)不多，因此，整體滲透矩陣是一個(gè)具有對(duì)稱(chēng)性的高度稀疏矩陣。

式(16)中等號(hào)右端的向量同樣由單元右端向量疊加得到，單元右端向量中各元素計(jì)算公式為：

(19)

式中：w為入滲或蒸發(fā)水量；q為邊界單位面積上流入的流量。

1.3 裂隙-連續(xù)介質(zhì)耦合方法及程序

1.3.1 耦合原理

為了解決連續(xù)介質(zhì)和裂隙介質(zhì)相鄰存在的區(qū)域的二維穩(wěn)定滲流問(wèn)題，本文建立了裂隙-連續(xù)介質(zhì)滲流法，其核心為：首先分別對(duì)整個(gè)滲流分析區(qū)域內(nèi)的裂隙介質(zhì)系統(tǒng)和連續(xù)介質(zhì)系統(tǒng)建立相應(yīng)的整體滲透矩陣，然后基于兩類(lèi)介質(zhì)共有節(jié)點(diǎn)的水頭相等以及流量平衡原則，形成整個(gè)滲流區(qū)域的整體滲透矩陣，進(jìn)行有限元分析，實(shí)際上是一種整體求解法。

由1.1節(jié)中基于圖論表示的裂隙網(wǎng)絡(luò)滲流理論可知，二維裂隙網(wǎng)絡(luò)穩(wěn)定滲流方程可寫(xiě)為：

[A][Tl][A]T{h}={Q}

(20)

式中：[A]為銜接矩陣；[Tl]為裂隙段水力傳導(dǎo)系數(shù)矩陣；{h}為節(jié)點(diǎn)水頭；{Q}為源匯項(xiàng)。

基于里茲有限單元法分析裂隙滲流，過(guò)程同1.2節(jié)中連續(xù)介質(zhì)滲流理論，將每一條裂隙段看作一個(gè)線(xiàn)單元對(duì)滲流區(qū)域進(jìn)行離散化，若j單元兩端分別為i節(jié)點(diǎn)和i+1節(jié)點(diǎn)，則有單元滲透矩陣為：

(21)

(22)

對(duì)于j單元?jiǎng)t有方程組為：

(23)

式中：Q和Qi+1分別為i節(jié)點(diǎn)和i+1節(jié)點(diǎn)的已知流量。

由所有的單元滲透矩陣疊加得到整體滲透矩陣為[K]，最終得到有限元方程組為：

[K]{h}={Q}

(24)

并且有：

[K]=[A][Tl][A]T

(25)

因此，在研究裂隙系統(tǒng)和連續(xù)介質(zhì)系統(tǒng)橫向耦合時(shí)，可以用線(xiàn)單元(裂隙網(wǎng)絡(luò)中兩個(gè)裂隙交叉點(diǎn)之間裂隙段)對(duì)裂隙介質(zhì)進(jìn)行離散化，用平面四邊形單元對(duì)連續(xù)介質(zhì)進(jìn)行離散化，對(duì)所有節(jié)點(diǎn)進(jìn)行統(tǒng)一編號(hào)，求得各個(gè)單元的滲透矩陣，按照節(jié)點(diǎn)連接關(guān)系組裝在一起形成總體滲透矩陣，進(jìn)行整體求解。

1.3.2 程序編制

由滲透矩陣疊加原則可以看出，整體滲透矩陣中的元素kij僅僅與i節(jié)點(diǎn)和j節(jié)點(diǎn)在含有該節(jié)點(diǎn)單元的局部位置有關(guān)系，與疊加時(shí)選擇單元的先后順序沒(méi)有關(guān)系。因此，本文在編制程序時(shí)，首先分別對(duì)裂隙系統(tǒng)和孔隙系統(tǒng)進(jìn)行整體滲透矩陣的組裝，最終對(duì)所有節(jié)點(diǎn)進(jìn)行統(tǒng)一編號(hào)，組裝成整個(gè)滲流系統(tǒng)的整體滲透矩陣。

2 程序應(yīng)用

2.1 工程概況

某土質(zhì)心墻堆石壩，建基面高程為2 930.00 m，壩高150 m，頂寬14 m，壩體內(nèi)設(shè)置黏土心墻防滲體，坡度均1∶0.25，頂寬6 m。大壩上下游坡面均在3 040.00 m高程處設(shè)4 m寬?cǎi)R道，上游馬道以上坡比為1∶2.5，以下為1∶2.25，下游坡比均為1∶2.0。壩基以下為深度為500 m左右的深厚覆蓋層，大致分為4層，在心墻底部設(shè)1.3 m厚、150 m深的混凝土防滲墻作為壩基防滲措施，插入心墻深度15 m，另在上游壩腳設(shè)置1.3 m厚、50 m深的副防滲墻，上游壩基面設(shè)水平防滲層。上游正常蓄水位3 070.00 m，對(duì)應(yīng)的尾水位為2 940.00 m。

2.2 模型參數(shù)

表1中給出了該土質(zhì)心墻堆石壩各個(gè)分區(qū)材料的滲透系數(shù)，本次計(jì)算采用已知水頭邊界，上游壩坡面按正常蓄水位施加，下游壩坡面按相應(yīng)的尾水位施加。

假設(shè)壩基為裂隙巖體，其中分布著大量裂隙，采用裂隙連通路徑程序進(jìn)行通路搜索后，得到裂隙連通網(wǎng)絡(luò)(見(jiàn)圖3)。

表1 材料滲透系數(shù)

圖3裂隙連通網(wǎng)絡(luò)模型概化圖

2.3 計(jì)算方案

本次計(jì)算采用表2中的6種方案，垂直防滲墻深度由150 m減小到110 m。主要計(jì)算內(nèi)容有通過(guò)壩體和壩基的單寬總滲流量和主防滲墻的最大滲透比降。

表2 計(jì)算方案

2.4 結(jié)果分析

圖4顯示，主防滲墻深度在150 m到120 m之間變化時(shí)，流經(jīng)土石壩壩基和壩體的單寬總滲流量幾乎保持不變，而當(dāng)主防滲墻深度減小至110 m時(shí)，單寬總滲流量從223 m3/年增長(zhǎng)至230 m3/年，增長(zhǎng)率為3.14%。其原因是，主防滲墻深度在120 m到150 m的這段位于巖塊內(nèi)部，根據(jù)圖10可知，縮短過(guò)程并沒(méi)有對(duì)裂隙網(wǎng)絡(luò)造成影響；而當(dāng)從120 m縮減到110 m時(shí)，本來(lái)由防滲墻阻斷的兩條裂隙，重新形成通路，連接到裂隙網(wǎng)絡(luò)中，形成了新的滲流路徑，導(dǎo)致滲流量的突然增加。由圖5可見(jiàn)，在主防滲墻深度的減小的同時(shí)，主防滲墻的最大滲透比降基本沒(méi)有變化。

圖4 單寬總滲流量隨防滲墻深度的變化

圖5 防滲墻最大滲透比降隨其深度變化

圖6防滲墻深度從150 m到110 m變化圖

3 結(jié) 論

本文針對(duì)裂隙巖基上修筑土石壩這種特殊滲流區(qū)域，根據(jù)裂隙系統(tǒng)與連續(xù)介質(zhì)系統(tǒng)交界處節(jié)點(diǎn)水頭一致的原則，將裂隙系統(tǒng)和連續(xù)介質(zhì)系統(tǒng)的單元傳導(dǎo)系數(shù)矩陣組集在一起，形成裂隙-連續(xù)介質(zhì)整體耦合滲流剛度矩陣，然后對(duì)整體方程求解滲流場(chǎng)。編制了裂隙-連續(xù)介質(zhì)耦合滲流分析程序，算例驗(yàn)證了程序的可靠性。

研究表明，程序可以有效描述裂隙巖體滲流特征和連續(xù)介質(zhì)滲流規(guī)律，顯示出良好的工程應(yīng)用前景，為裂隙與連續(xù)介質(zhì)的耦合滲流計(jì)算提供了新的方法。