收稿日期：2023-02-16；接受日期：2023-05-07

基金項目：水資源與水電工程科學(xué)國家重點實驗室（武漢大學(xué)）開放研究基金項目（2019SGG04）

作者簡介：

闖喜宏，男，高級工程師，碩士，主要從事巖土工程建設(shè)管理相關(guān)工作。E-mail：384149491@qqcom

EditorialOfficeofYangtzeRiverThisisanopenaccessarticleundertheCCBY-NC-ND40license

文章編號：1001-4179（2024）03-0184-06

引用本文：闖喜宏，樊鑫成，秦慧凱有自由面滲流問題的等效裂隙網(wǎng)絡(luò)模型［J］人民長江，2024，55（3）：184-189

摘要：

為降低傳統(tǒng)連續(xù)介質(zhì)滲流模型的網(wǎng)格依賴性并提升數(shù)值計算穩(wěn)定性，基于連續(xù)介質(zhì)滲流模型，通過引入表征單元體將均勻多孔介質(zhì)水平和豎直孔隙概化為正交裂隙網(wǎng)絡(luò)，建立了有自由面滲流問題的等效裂隙網(wǎng)絡(luò)模型。推導(dǎo)了等效裂隙網(wǎng)絡(luò)模型的滲流速度、連續(xù)性方程及邊界條件的等效形式，將二維滲流問題轉(zhuǎn)化為一維等效裂隙網(wǎng)絡(luò)滲流問題。結(jié)合連續(xù)型罰Heaviside函數(shù)，將濕區(qū)滲流流速擴(kuò)展至整個研究區(qū)域，給出了等效裂隙網(wǎng)絡(luò)滲流分析的有限元求解格式。通過與均質(zhì)矩形壩和梯形壩其他數(shù)值解的對比分析，驗證了等效裂隙網(wǎng)絡(luò)模型的有效性。最后將該方法應(yīng)用于土質(zhì)斜心墻壩穩(wěn)定滲流分析，計算結(jié)果很好地反映了黏土心墻防滲體與上下游壩殼堆石區(qū)的自由面分布規(guī)律，揭示了黏土心墻防滲體的阻水效應(yīng)。

關(guān)鍵詞：等效裂隙網(wǎng)絡(luò)模型；自由面；穩(wěn)定滲流；表征單元體

中圖法分類號：TV871

文獻(xiàn)標(biāo)志碼：A" " " " " " " " " "DOI：1016232jcnki1001-4179202403025

0引言

巖土、水利、公路等工程普遍涉及巖土體自由面滲流問題，巖土體滲流分析通常采用連續(xù)介質(zhì)滲流模型［1-3］，通過引入表征單元體（RepresentativeElementaryVolume，REV）反映孔隙結(jié)構(gòu)滲透特性的統(tǒng)計平均規(guī)律，宏觀尺度上將巖土介質(zhì)作為一系列表征單元體組成的連續(xù)介質(zhì)體［4］。基于連續(xù)介質(zhì)滲流模型，假定滲透流體充滿包含固體介質(zhì)與孔隙的整個滲流區(qū)域，以Darcy定律為基礎(chǔ)建立其滲流控制方程。對于有自由面的穩(wěn)定滲流問題，由于自由面和出滲點的位置是未知的，進(jìn)一步增加了滲流問題的非線性程度與數(shù)值求解難度［5-7］。

目前，有限單元法是求解有自由面滲流問題較為成熟的數(shù)值計算方法，主要包括變網(wǎng)格法和固定網(wǎng)格法。變網(wǎng)格法［8］在每次自由面迭代過程中均需要重新劃分網(wǎng)格，對于復(fù)雜巖土體而言計算較為繁瑣，收斂性較差。固定網(wǎng)格法的核心是在計算過程中保持網(wǎng)格不變，代表性方法有Desai等提出的剩余流量法［9］、Bathe等提出的調(diào)整滲透系數(shù)法［10］、張有天等提出的初流量法［11］以及變分不等式法［12-14］。前三種方法都是通過自由面迭代盡量使干區(qū)的滲流量遠(yuǎn)遠(yuǎn)低于濕區(qū)的滲流量；變分不等式法具有嚴(yán)格的數(shù)學(xué)理論基礎(chǔ)，在理論上克服了溢出點的奇異性，將自由面邊界轉(zhuǎn)化自然邊界條件，從而降低滲流問題的數(shù)值求解難度。

與傳統(tǒng)的連續(xù)介質(zhì)滲流模型不同，等效管道模型將連續(xù)介質(zhì)劃分為三角形管道網(wǎng)絡(luò)［15-19］，假定三角形內(nèi)部基質(zhì)不透水，水只能沿著三角形各邊流動，三角形各邊則被一維管道網(wǎng)絡(luò)替代，通過滲透性等效簡化為一維滲流問題，從而提高數(shù)值求解的效率和精度。Xu等［15］首先提出了巖土介質(zhì)有壓滲流的拓?fù)涔芫W(wǎng)法，理論上推導(dǎo)了均勻各向同性介質(zhì)管網(wǎng)模型的等效水力參數(shù)。Ren等［16］進(jìn)一步將上述研究工作拓展至無壓滲流問題，同時考慮了多孔基質(zhì)與裂隙網(wǎng)絡(luò)的滲透特性，建立了裂隙巖體無壓滲流的等效管網(wǎng)滲流分析方法。Afzali等［17］和Abareshi等［18］采用等效管網(wǎng)模型求解了巖土介質(zhì)自由面滲流問題，探討了等效管網(wǎng)模型與連續(xù)介質(zhì)滲流模型有限元求解的等效性。Moosavian［19］提出了等效管網(wǎng)滲流分析顯式和隱式兩種方法，研究了線性和非線性流動對堆石壩滲流場的影響。然而，現(xiàn)有的等效管網(wǎng)模型僅適用于三角形網(wǎng)格劃分，等效水力參數(shù)計算依賴三角形單元尺寸，自由面穿過三角形單元需要重新劃分網(wǎng)格和計算等效水力參數(shù)，從而難以克服網(wǎng)格依賴性與降低數(shù)值求解難度。

為了求解有自由面的穩(wěn)定滲流問題，本文基于連續(xù)介質(zhì)滲流理論，引入表征單元體將均勻多孔介質(zhì)水平和豎直孔隙概化為正交裂隙網(wǎng)絡(luò)，建立了有自由面滲流的等效裂隙網(wǎng)絡(luò)模型，推導(dǎo)了滲流速度、連續(xù)性方程及邊界條件的等效一維形式，給出了相應(yīng)的有限元求解格式及對應(yīng)的迭代算法。通過對矩形壩、梯形壩以及土質(zhì)斜心墻壩等開展?jié)B流場分析，研究了裂隙網(wǎng)絡(luò)尺寸和滲透系數(shù)對自由面的影響，驗證了等效裂隙網(wǎng)絡(luò)模型求解有自由面滲流問題的正確性和有效性。

1有自由面穩(wěn)定滲流問題的描述

以圖1所示均質(zhì)土壩穩(wěn)定滲流為例，自由面把整個區(qū)域劃分為干區(qū)和濕區(qū)兩部分。干區(qū)為負(fù)壓區(qū)，濕區(qū)為正壓區(qū)。當(dāng)自由面確定時，濕區(qū)即隨之確定。基于多孔介質(zhì)滲流理論，濕區(qū)中地下水流動應(yīng)滿足連續(xù)性方程：

式中：vx和vz分別為x和z方向的滲流速度。

根據(jù)Darcy定律，滲流速度可表示為

式中：kx和kz為x和z方向的滲透系數(shù)；φ=z+p/ρg是總水頭，p為孔隙水壓力，z是垂直坐標(biāo)分量。

式（1）還應(yīng)滿足如下定解條件：

（1）水頭邊界條件（AB和FG）為

式中：為水頭邊界的已知水頭，在上、下游分別為1和2。

（2）流量邊界條件（AG）為

式中：q為流量邊界上的已知流量；nx和nz分別為邊界單位外法線向量在水平和垂直方向上的余弦。對隔水邊界，q=0。

（3）溢出面邊界條件（EF）為

（4）自由面邊界條件（BE）為

2等效裂隙網(wǎng)絡(luò)模型

多孔介質(zhì)中含有大量的不規(guī)則顆粒和孔隙，逐一研究單個孔隙結(jié)構(gòu)的流動特性非常困難，因此多孔介質(zhì)滲流分析往往采用連續(xù)介質(zhì)模型，通過引入表征單元體反映孔隙結(jié)構(gòu)滲透特性的統(tǒng)計平均規(guī)律。表征單元體（REV）是指多孔介質(zhì)滲透特性趨于基本穩(wěn)定時的最小體積，當(dāng)REV與滲流研究區(qū)域相比足夠小時，研究區(qū)域就可以看成由表征單元體組成的連續(xù)介質(zhì)體，從而建立多孔介質(zhì)滲流連續(xù)性方程。為簡單起見，均勻多孔介質(zhì)REV內(nèi)顆粒和孔隙的概化模型如圖2（a）所示，顆粒和孔隙分別呈方形和十字狀排列。如圖2（b）所示，由于水流只能在水平和豎直孔隙內(nèi)發(fā)生流動，因此可以將水平和豎直孔隙簡化為正交裂隙網(wǎng)絡(luò)，水平和豎直裂隙開度均勻且壁面光滑。

21水平與垂直裂隙的滲流速度

假定地下水僅在水平和垂直裂隙內(nèi)流動，根據(jù)立方定理，水平和垂直裂隙的滲流速度vf，x和vf，z可以表示為

式中：bx、bz分別為水平、垂直裂隙的等效開度；υ為水的運動黏滯系數(shù)。

22水平和垂直裂隙的等效開度

基于流量等效原則，多孔介質(zhì)模型和等效裂隙網(wǎng)絡(luò)模型水平方向上的流量關(guān)系如下：

式中：dz為REV的垂直長度；Nx為REV內(nèi)水平裂隙的數(shù)目，可以表示為

式中：Bx為水平裂隙間距。由于裂隙開度較小即bxBx，式（11）可以被簡化為

將上式代入式（10）并結(jié)合式（9）和式（10），可以得到：

同理，垂直裂隙的等效裂隙開度可以寫成：

式中：Bz是垂直裂隙間距。

23水平與垂直裂隙的等效連續(xù)性方程

對于水平裂隙滲流有vf，z=0，則式（1）可以被簡化為

同理，對于垂直裂隙滲流有vf，x=0，可以得到如下關(guān)系式：

24等效裂隙網(wǎng)絡(luò)模型的一維形式

如圖3所示，對任意等效裂隙ij建立一維局部坐標(biāo)系l，式（7）和式（8）可以被合并成：

式中：vij和bij分別為等效裂隙ij的滲流速度和等效開度。當(dāng)裂隙處于水平方向時即bij=bx，式（17）與式（7）完全一致，當(dāng)裂隙處于垂直方向時即bij=bz，式（17）與式（8）完全等效。

同理，式（15）和式（16）也可以被統(tǒng)一寫成一維形式：

由于等效裂隙節(jié)點不能存儲水，根據(jù)質(zhì)量守恒原理，流入和流出總和為零，即：

式中：mi為相交于節(jié)點i的裂隙總數(shù)。

等效裂隙網(wǎng)絡(luò)模型的邊界條件可以等效寫成：

（1）水頭邊界條件為

φi=φ（節(jié)點i∈AB+FG）（20）

（2）流量邊界條件為

qn=bijvij=q（節(jié)點i∈AG）（21）

（3）溢出面邊界條件為

φi=ziqij≤0（節(jié)點i∈EF）（22）

（4）自由面邊界條件為

φi=ziqij=0（節(jié)點i∈BE）（23）

至此，已經(jīng)推導(dǎo)了等效裂隙網(wǎng)絡(luò)模型的滲流速度、連續(xù)性方程及邊界條件的等效形式，從而將有自由面穩(wěn)定滲流的二維問題轉(zhuǎn)化為等效裂隙網(wǎng)絡(luò)滲流的一維問題。

3有限元數(shù)值求解格式

為了解決全區(qū)域（包含濕區(qū)與干區(qū)）自由面滲流問題，避免自由面迭代過程的數(shù)值不穩(wěn)定性和網(wǎng)格依賴性，通過引入連續(xù)型罰Heaviside函數(shù)，將濕區(qū)滲流流速擴(kuò)展至整個研究區(qū)域，等效裂隙滲流速度可以表示為

vij=Hλ（φ-z）gb2ij12υφl（24）

Hλ（φ-z）=1φ-z≥0（λ+φ-z）/λ0gt;φ-zgt;-λ0φ-z≤-λ（25）

式中：λ為罰參數(shù)。

將式（24）代入式（18），等效裂隙網(wǎng)絡(luò)模型全區(qū)域的連續(xù)性方程為

lHλ（φ-z）gb3ij12υφl=0（26）

根據(jù)變分原理，式（26）的泛函可以被表達(dá)為

I（φ）=Ω∫lij12Hλ（φ-z）gb3ij12υφl2dl+i∈Γqijφ（27）

泛函I（）對等效裂隙節(jié)點水頭求極值可得

lφi=Ω∫lijHλ（φ-z）gb3ij12υφlφiφldl+i∈Γbijqijφφi=0（28）

任意等效裂隙水頭插值函數(shù)可以簡單寫成

φ=Niφi+Njφj=［NiNj］φiφj=Nφe（29）

式中：φi和φj分別為等效裂隙端點i和j的總水頭；Ni和Nj是形函數(shù)，Ni=1-l/lij，Nj=l/lij，lij是裂隙長度；N=［NiNj］，φe=［φiφj］T。

由式（29）進(jìn)一步推導(dǎo)得出

φl=-1lijφi+1lijφj=-1lij1lijφiφj=Bφe（30）

因此，式（28）的矩陣形式可以寫為

Kφη+1=qη（31）

其中

K=Ω∫lijBTgb3ij12υBdl（32）

qη=Kλφη（33）

Kλ=ΩKeλ（34）

Keλ=∫lij1-Hλφη-zBTgb3ij12υBdl（35）

φ=φ1φ2…φN（36）

式中：η為迭代次數(shù)，N是整個區(qū)域所有裂隙節(jié)點總數(shù)。自由面迭代過程中，選用鄭宏等［13］建議的互補(bǔ)算法求解。算法的收斂條件定義如下：

‖φη+1-φη‖≤δ‖φη‖（37）

式中：δ為指定容差，本文取為0001。

4算例分析

41均質(zhì)矩形壩

首先利用姚池等［20］的均質(zhì)矩形壩模型驗證本文提出的等效裂隙網(wǎng)絡(luò)模型有效性。均質(zhì)矩形壩高12m，寬10m，底部邊界隔水，上游水位為10m，下游水位為2m。

采用5種等效裂隙網(wǎng)絡(luò)（Bx=Bz=01，02，025，05，10m）進(jìn)行均質(zhì)矩形壩滲流分析并確定自由面位置，根據(jù)式（13）和（14），相應(yīng)的等效裂隙開度分別為bx=bz=7985×10-5，1006×10-4，1084×10-4，1365×10-4，172×10-4m，其自由面的計算結(jié)果如圖4所示。比較可見，5種裂隙網(wǎng)格得到的自由面與姚池等［20］的計算結(jié)果幾乎完全一致，表明了采用等效裂隙網(wǎng)絡(luò)模型分析多孔介質(zhì)滲流的有效性且對網(wǎng)格依賴性較小。

均質(zhì)壩體滲流的總流量可以通過Dupuit公式計算：

q=kh21-h222L（38）

式中：h1和h2分別表示上下游水位，L表示壩體寬度。根據(jù)Dupuit公式得到均質(zhì)矩形壩的流量為198×10-5m2/s，5種裂隙網(wǎng)絡(luò)模型得到的流量和相對誤差如圖5所示，隨著網(wǎng)格尺寸的增大，流量逐漸減小，當(dāng)網(wǎng)格尺寸為025m時，相對誤差最小；其中最大誤差僅為556%，表明數(shù)值計算結(jié)果是合理可靠的。

42梯形壩

如圖6所示，Lacy和Prevost［21］的均質(zhì)梯形壩模型上底為2m，下底為7m，高為5m。底部為隔水邊界，上游水位和下游水位分別是5m和1m，滲透系數(shù)為01m/d。圖6給出了等效裂隙網(wǎng)絡(luò)模型滲流與Lacy和Prevost得到的自由面，3種等效裂隙網(wǎng)絡(luò)的自由面計算結(jié)果與Lacy和Prevost的數(shù)值解基本吻合，表明等效裂隙網(wǎng)絡(luò)模型適用于求解復(fù)雜幾何邊界模型的有自由面穩(wěn)定滲流問題。

43土質(zhì)斜心墻壩

如圖7所示，土質(zhì)斜心墻壩壩頂高程67m，壩頂寬6m，壩底寬180m。上游水位為60m，下游水位為0m，底部為不透水邊界。上下游壩殼堆石料的滲透系數(shù)k1=3×10-5m/s，黏土心墻防滲體的滲透系數(shù)k2=5×10-6m/s。等效裂隙網(wǎng)絡(luò)模型與Lacy和Prevost［21］得到的自由面對比如圖7所示，兩者吻合較好，表明等效裂隙網(wǎng)絡(luò)模型求解非均質(zhì)模型穩(wěn)定滲流問題具有可靠性。

為了研究黏土心墻滲透系數(shù)對壩體滲流場的影響，分別取黏土心墻的滲透系數(shù)為5×10-6，3×10-6，1×10-6m/s進(jìn)行滲流分析，自由面的計算結(jié)果如圖8所示。由圖可知，穩(wěn)定滲流條件下，黏土心墻的滲透系數(shù)由高到低自由面分布規(guī)律基本相同，上游壩殼堆石區(qū)自由面基本呈水平狀。但是，隨著黏土心墻的滲透系數(shù)不斷減小，心墻承擔(dān)的水頭損失逐步增大，即心墻內(nèi)自由面降落幅度更加顯著，說明黏土心墻的滲透系數(shù)越小，阻水效果越明顯；下游壩殼堆石區(qū)自由面隨心墻滲透系數(shù)的減小逐漸下降，心墻后自由面高程從2943m降至836m。

5結(jié)論與展望

針對有自由面的穩(wěn)定滲流問題，采用正交裂隙網(wǎng)絡(luò)替代連續(xù)介質(zhì)，建立了有自由面滲流的等效裂隙網(wǎng)絡(luò)模型，給出了相應(yīng)的有限元求解格式及對應(yīng)的迭代算法。通過矩形壩、梯形壩以及土質(zhì)斜心墻壩穩(wěn)定滲流分析，論證了等效裂隙網(wǎng)絡(luò)模型的正確性和有效性，主要結(jié)論如下：

（1）基于連續(xù)介質(zhì)滲流模型，通過引入表征單元體，將均勻多孔介質(zhì)的水平和豎直孔隙概化為正交裂隙網(wǎng)絡(luò)，在數(shù)學(xué)上嚴(yán)格推求了等效裂隙網(wǎng)絡(luò)模型的滲流速度、連續(xù)性方程及邊界條件的等效形式，將二維滲流問題簡化為一維滲流問題，從而克服了對三角形網(wǎng)格的依賴性并降低了數(shù)值求解難度。

（2）通過引入連續(xù)型Heaviside函數(shù)，將濕區(qū)滲流流速擴(kuò)展至整個研究區(qū)域，基于變分原理推導(dǎo)了等效裂隙網(wǎng)絡(luò)滲流分析的有限元求解格式，保證了自由面迭代過程的穩(wěn)定性和收斂性。

（3）均質(zhì)壩和梯形壩自由面對比分析，驗證了等效裂隙網(wǎng)絡(luò)模型求解有自由面穩(wěn)定滲流問題的合理性和正確性，進(jìn)一步證明了等效裂隙網(wǎng)絡(luò)模型與連續(xù)介質(zhì)滲流模型滲流分析的等價性。

（4）土質(zhì)斜心墻壩穩(wěn)定滲流分析表明，穩(wěn)定滲流條件下，上游壩殼堆石區(qū)自由面基本呈水平狀，隨著黏土心墻滲透系數(shù)的不斷減小，心墻內(nèi)自由面降幅越大；下游壩殼堆石區(qū)自由面隨心墻滲透系數(shù)的減小逐漸下降。等效裂隙網(wǎng)絡(luò)模型可以很好地體現(xiàn)黏土心墻防滲體與上下游壩殼堆石區(qū)自由面的分布規(guī)律，并揭示出黏土心墻防滲體的阻水效應(yīng)。

綜上所述，等效裂隙網(wǎng)絡(luò)模型可以準(zhǔn)確有效求解有自由面的滲流問題，該方法簡單高效，對網(wǎng)格尺寸的敏感性較弱?？梢灶A(yù)見，由于二維滲流問題和三維滲流問題具有一定的相似性，等效裂隙網(wǎng)絡(luò)模型還可以拓展至三維領(lǐng)域。但是，實際工程的滲流情況相較復(fù)雜，仍需進(jìn)一步的工程應(yīng)用研究來檢驗?zāi)Ｐ偷木_度。

參考文獻(xiàn)：

［1］徐穎，王偉，李艷玲，等高土石壩雙防滲墻滲流特性研究［J］人民江，2022，53（7）：181-186

［2］譚彬政，陶桂蘭，莊寧，等水位升降對荊江高灘岸坡滲流場影響分析［J］人民長江，2015，46（9）：36-40

［3］張友才，武雁剛，陸鳳，等堤防滲流值模擬與防滲方案研究［J］水利水電快報，2018，39（12）：43-48

［4］BEARJDynamicsoffluidsinporousmedia［M］NewYork：Elsevier，1972

［5］胡冉，陳益峰，周創(chuàng)兵，等非穩(wěn)定滲流問題的變分不等式方法及工程應(yīng)用［J］水動力學(xué)研究與進(jìn)展，2011（2）：239-251

［6］王慧，蘇永軍，王鳳瑞基于分形理論的巖石裂隙非線性滲流各向異性研究［J］人民長江，2019，50（2）：174-180，212

［7］ZHENGH，DAIH，LIUDAvariationalinequalityformulationforunconfinedseepageproblemsinporousmedia［J］AppliedMathematicalModelling，2009，33：437-450

［8］NEUMANSP，WITHERSPOONPAFiniteelementmethodofanalyzingsteadyseepagewithafreesurface［J］WaterResourcesResearch，1970，6（3）：889-897

［9］DESAICS，LIGCAresidualflowprocedureandapplicationforfreesurfaceinporousmedia［J］AdvancesinWaterResources，1983，6（1）：27-35

［10］BATHEK，KHOSHGOFTAARMFiniteelementfreesurfaceseepageanalysiswithoutmeshiteration［J］InternationalJournalforNumericalandAnalyticalMethodsinGeomechanics，1979，3（1）：13-22

［11］張有天，陳平，王鐳有自由面滲流分析的初流量法［J］水利學(xué)報，1988（8）：20-28

［12］張楓跨越堤防橋梁基礎(chǔ)對堤防滲透穩(wěn)定性影響分析［J］人民長江，2011，42（增2）：147-149

［13］鄭宏，戴會超，劉德富改進(jìn)的有自由面滲流問題的Bathe算法［J］巖土力學(xué)，2005（4）：505-512

［14］葉祖洋，姜清輝，姚池，等巖體裂隙網(wǎng)絡(luò)非穩(wěn)定滲流分析與數(shù)值模擬［J］巖土力學(xué)，2013（4）：1171-1178

［15］XUZ，MAG，LISAGraph-theoreticpipenetworkmethodforwaterflowsimulationinaporousmedium：GPNM［J］InternationalJournalofHeatandFluidFlow，2014，45：81-97

［16］RENF，MAG，YANGW，etalPipenetworkmodelforunconfinedseepageanalysisinfracturedrockmasses［J］InternationalJournalofRockMechanicsandMiningSciences，2016，88：183-196

［17］AFZALIS，MONADJEMIPSimulationofflowinporousmedia：Anexperimentalandmodelingstudy［J］JournalofPorousMedia，2014，17（6）：469-481

［18］ABARESHIS，HOSSEINIASEquivalentpipenetworkmodelforacoarseporousmediaanditsapplicationtoanalysisofflowthroughrockfillstructures［J］JournalofPorousMedia，2017，20（4）：303-324

［19］MOOSAVIANNPipenetworkmodellingforanalysisofflowinporousmedia［J］CanadianJournalofCivilEngineering，2019，46（12）：1151-1159

［20］姚池，姜清輝，葉祖洋，等裂隙網(wǎng)絡(luò)無壓滲流分析的初流量法［J］巖土力學(xué)2012，6（33）：1896-1903

［21］LACYS，PREVOSTJFlowthroughporousmedia：Aprocedureforlocatingthefreesurface［J］InternationalJournalforNumericalandAnalyticalMethodsinGeomechanics，1987，11（6）：585-601

（編輯：鄭毅）