有自由面滲流問題的等效裂隙網絡模型

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

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

作者簡介：

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

EditorialOfficeofYangtzeRiverThisisanopenaccessarticleundertheCCBY-NC-ND40license

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

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

摘要：

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

關鍵詞：等效裂隙網絡模型；自由面；穩定滲流；表征單元體

中圖法分類號：TV871

文獻標志碼：A" " " " " " " " " "DOI：1016232jcnki1001-4179202403025

0引言

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

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

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

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

1有自由面穩定滲流問題的描述

以圖1所示均質土壩穩定滲流為例，自由面把整個區域劃分為干區和濕區兩部分。干區為負壓區，濕區為正壓區。當自由面確定時，濕區即隨之確定。基于多孔介質滲流理論，濕區中地下水流動應滿足連續性方程：

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

根據Darcy定律，滲流速度可表示為

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

式（1）還應滿足如下定解條件：

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

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

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

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

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

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

2等效裂隙網絡模型

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

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

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

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

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

基于流量等效原則，多孔介質模型和等效裂隙網絡模型水平方向上的流量關系如下：

式中：dz為REV的垂直長度；Nx為REV內水平裂隙的數目，可以表示為

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

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

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

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

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

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

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

24等效裂隙網絡模型的一維形式

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

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

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

由于等效裂隙節點不能存儲水，根據質量守恒原理，流入和流出總和為零，即：

式中：mi為相交于節點i的裂隙總數。

等效裂隙網絡模型的邊界條件可以等效寫成：

（1）水頭邊界條件為

φi=φ（節點i∈AB+FG）（20）

（2）流量邊界條件為

qn=bijvij=q（節點i∈AG）（21）

（3）溢出面邊界條件為

φi=ziqij≤0（節點i∈EF）（22）

（4）自由面邊界條件為

φi=ziqij=0（節點i∈BE）（23）

至此，已經推導了等效裂隙網絡模型的滲流速度、連續性方程及邊界條件的等效形式，從而將有自由面穩定滲流的二維問題轉化為等效裂隙網絡滲流的一維問題。

3有限元數值求解格式

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

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

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

式中：λ為罰參數。

將式（24）代入式（18），等效裂隙網絡模型全區域的連續性方程為

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

根據變分原理，式（26）的泛函可以被表達為

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

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

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

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

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

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

由式（29）進一步推導得出

φ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）

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

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

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

4算例分析

41均質矩形壩

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

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

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

q=kh21-h222L（38）

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

42梯形壩

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

43土質斜心墻壩

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

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

5結論與展望

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

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

（2）通過引入連續型Heaviside函數，將濕區滲流流速擴展至整個研究區域，基于變分原理推導了等效裂隙網絡滲流分析的有限元求解格式，保證了自由面迭代過程的穩定性和收斂性。

（3）均質壩和梯形壩自由面對比分析，驗證了等效裂隙網絡模型求解有自由面穩定滲流問題的合理性和正確性，進一步證明了等效裂隙網絡模型與連續介質滲流模型滲流分析的等價性。

（4）土質斜心墻壩穩定滲流分析表明，穩定滲流條件下，上游壩殼堆石區自由面基本呈水平狀，隨著黏土心墻滲透系數的不斷減小，心墻內自由面降幅越大；下游壩殼堆石區自由面隨心墻滲透系數的減小逐漸下降。等效裂隙網絡模型可以很好地體現黏土心墻防滲體與上下游壩殼堆石區自由面的分布規律，并揭示出黏土心墻防滲體的阻水效應。

綜上所述，等效裂隙網絡模型可以準確有效求解有自由面的滲流問題，該方法簡單高效，對網格尺寸的敏感性較弱。可以預見，由于二維滲流問題和三維滲流問題具有一定的相似性，等效裂隙網絡模型還可以拓展至三維領域。但是，實際工程的滲流情況相較復雜，仍需進一步的工程應用研究來檢驗模型的精確度。

參考文獻：

［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］胡冉，陳益峰，周創兵，等非穩定滲流問題的變分不等式方法及工程應用［J］水動力學研究與進展，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］水利學報，1988（8）：20-28

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

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

［14］葉祖洋，姜清輝，姚池，等巖體裂隙網絡非穩定滲流分析與數值模擬［J］巖土力學，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］姚池，姜清輝，葉祖洋，等裂隙網絡無壓滲流分析的初流量法［J］巖土力學2012，6（33）：1896-1903

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

（編輯：鄭毅）