長距離調水工程渠道裂隙巖質邊坡滲流模擬

黃兆虎，方攀博，韓延成，王 棟，鄧曉川，宋廣增，趙嘉誠

（1. 濟南大學，山東 濟南 250022； 2. 山東省水利勘測設計院有限公司，山東 濟南 250013；3. 山東省調水工程運行維護中心濰坊分中心，山東 濰坊 261000）

0 引 言

明渠輸水工程是實現我國水資源調度，保證水資源平衡的重要舉措，例如南水北調中、東線，引黃濟青工程、膠東調水工程等均采用明渠形式。然而，由于長距離調水工程距離長，經常會經過裂隙發育地區。裂隙發育地區的巖質邊坡，其滲流特性和與土質邊坡渠道不同［1］。裂隙的復雜分布以及地下水在裂隙巖體中的分布和滲流的不均勻性致使裂隙巖體滲流問題研究極其復雜。由于裂隙滲流的復雜性，以往長距離調水工程渠道滲流都是按孔隙模型進行滲流計算的，忽略了裂隙對滲流場的影響，而實際運行中滲透破壞時常發生，進而影響供水安全甚至發生嚴重的輸水事故。

滲流對大壩、渠道等水工建筑物邊坡穩定性、滲透破壞具有重要影響［2］，國內學者在渠道邊坡滲流方面進行了廣泛的研究。高丹盈等［3］用有限元方法進行了南水北調工程大沙河段渠道排水非穩定滲流計算分析；朱國勝等［4］對逆止式排水系統對南水北調中線渠道邊坡滲流控制效果進行了研究；鐘登華等［5］構建了南水北調中線工程某渠段的工程地質信息三維可視化模型，運用耦合方程，對研究區域進行滲流模擬，并采用流固耦合進行渠道穩定性計算，結果表明渠道滲流對其穩定性影響較大；左海鳳、王光謙等［6］建立了南水北調中線磁縣段三維水文地質結構模型，進行了地下水排水數值模擬；方攀博等［7］研究了渠道排水管的不同布設位置對渠道滲流及揚壓力的影響；陳思涵等［8］研究了高地下水位地區管、井相結合的復合排水方案及滲流場的模擬、減壓效果研究。

國內外學者在裂隙滲流方面也進行了廣泛的研究，在裂隙巖體渠道邊坡滲流方面，上世紀四五十年代前蘇聯的一些學者Волоэько 和Ломизе 等就對裂隙介質（巖體）的水力學問題進行了試驗研究，并最早提出了單裂隙巖體地下水運動的立方定律［9］。但在1963年法國Malpasset拱壩失事事故發生后，巖體裂隙滲流開始了更充分的研究。1965 年Snow［10］創立了巖體裂隙滲透率張量理論。李定方、萬力［11］采用有限元法分析滲流對巖體高邊坡的影響。張家發、李思慎、葉自桐［12］采用有限元分析法模擬分析了多年平均降雨條件下高邊坡巖體的穩定滲流狀態。在巖體飽和-非飽和流的滲透規律方面，Brooks-Corey 在1964 年提出了適用多孔介質的滲透特性的Brooks-Corey 模型。Mualem 于1976 年首次提出了適用于多孔介質的Mualem 模型，并驗證了此模型的適用性［13］。Van Genuchten［14］于1980 采用物理實驗法擬合孔隙介質模型得出非飽和裂隙滲流的重要模型Van Genuchten（VG）模型。黃濤等［15］將VG 模型和Mualem 模型結合在一起提出適用于兩相多孔介質的VG-Mualem 模型，經計算VG-M 模型與試驗結果的擬合效果最佳。平揚等［16］基于壓水試驗分析了不同自重應力以及巖性等因素對裂隙巖體滲透性的影響；魯志強等［17］進行了裂隙水作用下水利巖質邊坡抗傾覆動力穩定性研究，胡小昕等［18］進行了基于裂隙網絡有限元的地下洞室圍巖位移反演分析。

綜上所述，學者們分別對渠道邊坡滲流、裂隙滲流進行了廣泛的研究，目前在渠道設計時，滲流分析通常采用孔隙模型，在裂隙對渠道邊坡滲流的影響及孔隙模型與裂隙模型對渠道滲流影響對比方面缺乏研究。

山東省膠東調水工程是南水北調東線工程中山東“T”字型調水大動脈的重要組成部分，是緩解膠東地區水資源供需矛盾的重要水利基礎設施。但該工程萊州段部分渠段裂隙較多，所以通水以來每次停水檢查時均會出現大量的滲透破壞（見圖1），嚴重影響工程安全和供水安全。換言而知，在遠小于設計地下水位的情況下，仍然產生了大量的滲透破壞。破壞位置呈現點狀分布，開展渠道裂隙滲流的研究具有十分重要的現實意義。

圖1 渠道襯砌滲透破壞情況Fig.1 Seepage damage of channel lining

1 研究區概況及水文地質勘探

研究區位于膠東引黃調水工程萊州段，萊州市境西南部沿海低緩丘陵區，地處趴埠周家橋（樁號62+509）至后趴埠東交通橋段（樁號64+909）之間。東經119°49′33.213″～119°51′20.358″，北緯37°10′10.380″～37°11′19.691″。研究區西北為渤海，直線距離2.0 km。距離萊州市城區9 km，屬于溫暖帶季風區大陸性氣候，四季分明氣候溫和濕潤（圖2），多年平均氣溫12.4 ℃。1990-2019 年多年平均降水量為638.80 mm，最大年降水量為907.0 mm（2001 年）。研究區內地表水體不豐富，在研究區西側有朱流河經過，為季節性河流，常年水量短缺。研究區地層巖性簡單，主要為古元古界粉子山群崗箭組二云片巖、黑云片巖夾黑云變粒巖，巨屯組含石墨黑云斜長變粒巖及石墨大理巖等（圖3）。

圖2 萊州市1990-2019年降水量圖Fig.2 Precipitation map of Laizhou City from 1990 to 2019

圖3 研究區地質圖Fig.3 Geological map of study area

根據項目研究需要，在研究區布設了三眼勘測井（渠左南SK1 井、渠右南SK2 井、渠右北SK3 井），得到了巖芯，進行了滲水試驗［19］，安裝3 套水位自動觀測設備和1 套渠道水位觀測設備。在濟南大學進行了相關滲水實驗，滲水實驗所得滲透系數如表2所示。

2 模型計算原理

2.1 裂隙生成原理

目前直接獲取整體巖體邊坡的裂隙幾何分布是十分困難的，研究者通常采用鉆探得到的巖芯，采用統計方法來獲取裂隙幾何信息［20］。但這種方法所研究的范圍較小，覆蓋面不足，不能很好體現巖體內大范圍的裂隙幾何特征。本文通過鉆探獲取巖芯（圖4）統計裂隙幾何基本信息，基于蒙特卡羅法，利用MATLAB 語言編制二維離散化裂隙網絡程序生成渠道邊坡裂隙分布。

圖4 研究區域明渠兩岸鉆孔巖心Fig.4 Borehole cores on both sides of open channels in the study area

2.1.1 隨機裂隙網絡基本假定

對于要生成的裂隙網絡，進行如下5條假定［23］。

（1）裂隙網絡中裂隙的主要幾何參數有：跡長L、裂隙寬度、裂隙的分布密度和裂隙的傾角θ（即跡線與X軸方向的夾角）。

（2）對于生成裂隙，主要由裂隙的中心點（X0，Y0）、裂隙的跡長L、傾角θ、裂隙的兩端端點確定。裂隙兩端點的坐標表達式為：

端點1（X1，Y1）：

端點2（X2，Y2）：

（3）生成裂隙結構面的跡線為直線段，結構面的長度為裂隙的跡長、寬度為裂隙的張開度，結構面的產狀由裂隙的傾角確定。

（4）裂隙的數量由生成域的面積乘該組的裂隙密度得到。

（5）若有裂隙不位于和不完全位于研究邊界，位于研究域外的裂隙部分進行刪除。

2.1.2 裂隙網絡幾何參數確定

裂隙的生成需要對裂隙的各參數進行確定。具研究表明［21，22］，一般情況裂隙中心點服從均勻分布，裂隙跡長服從負指數分布或對數正態分布，裂隙傾角服從正態分布或對數正態分布，裂隙的密度服從負指數分布，裂隙的張開度服從負指數分布或對數正態分布。本文通過鉆孔調查、根據裂隙各參數分布規律進行計算分析，將研究區的裂隙組分為兩組，得到各組裂隙的初始幾何參數如下表1所示。

表1 裂隙初始幾何參數Tab.1 Initial geometrical parameters of fracture

2.2 飽和-非飽和滲流微分方程

由于渠道兩側降雨入滲對渠道邊坡滲流有影響。因此本文采用飽和-非飽和滲流微分方程進行滲流計算。

Darcy 定律是用來描述多孔介質中流體流動的方程，以孔隙水壓力為自變量，其微分形式如下［24］：

式中：K(θ)為土壤含水率函數；Hm為土水勢。

飽和滲流微分方程為：

式中：ε為孔隙率；ks為飽和滲透系數；p為孔隙水壓力；ρ為流體密度；g為重力加速度；D為縱坐標（方向同縱坐標z）。

非飽和滲流用Richards方程表示，以壓力水頭為自變量，其偏微分方程如下［25］：

式中：C為比水容量（dθ/dh）；Ss為儲水率也叫貯水率，1/m；θ為體積含水率；Se為有效飽和度；hp為壓力水頭（p）；Kr為非飽和滲透系數。其中非飽和水力參數C、Se、θ等用Van Genuchten 模型模擬。

Van Genuchten模型其原理表達式如下［26］：

式中：θ為體積含水率；h為負壓，m，取正值；θs，θr分別表示飽和體積含水率和殘余體積含水率；α、m、n為模型參數。

Kr的表達式為［27］：

式中：Se為有效飽和度；d為孔隙分布的分維。

3 研究區模型構建

為對比分析，本文分別建立考慮裂隙的模型（裂隙流模型）和不考慮裂隙的模型（孔隙模型）進行研究與對比。

3.1 研究區選取與模型生成

本文選取研究區域樁號62+809 位置處的局部渠道斷面作為此次的模擬計算區域，區域長高分別為55 和23 m。利用MATLAB 編制算法，以裂隙初始參數為依據，生成裂隙網絡模型如圖5所示。

圖5 隨機裂隙網絡圖Fig.5 Random fracture network diagram

將圖5 與渠道斷面疊加，并根據2.1.1 節所述假設，對不位于研究區域內的裂隙進行刪除和截?。徊捎萌切尉W絡剖分法進行有限元網格剖分（對研究重點區域進行網格加密），最終得到的裂隙流有限單元網格如圖6 所示，裂隙的開度則在COMSOL裂隙流模塊中設置。不考慮裂隙流的孔隙模型有限單元網格如圖7所示。

圖6 裂隙模型有限單元網格劃分Fig.6 Finite element meshing of fracture model

圖7 孔隙模型有限單元網格劃分Fig.7 Finite element meshing of pore model

3.2 模型相關參數設置

數值模擬采用COMSOL 軟件的Darcy 定律模塊進行數值模擬。在滲流計算中，采用Darcy 定律匹配Richards 方程作為滲流計算的基本方；非飽和水力參數采用Van Genuchten 模型計算；滲流相關參數采用野外實驗并參考原渠道設計滲流參數（見表2）。

表2 相關滲流計算參數Tab.2 Related seepage calculation parameters

3.3 模型邊界條件設置

在滲流模型中，上邊界與邊坡坡面為大氣降水邊界，左側設置水頭邊界。同等地下水位時渠道不輸水時揚壓力最大，為最不利工況，因此渠道底部設置為無流動邊界。右側按年庚乾［28］研究的滲流邊界設置（如圖8 所示），即當邊坡體內非飽和時，水不會向外滲出；當邊坡體內達到飽和時，水則不斷向外滲出，滲流邊界方程可表示為：

式中：n表示邊界外的法向量；u表示邊界的流速；ρ為液體密度；Ks為巖體飽和滲透系數；p為孔隙水壓力；L為耦合長度尺度（近似代表滲流邊界厚度）。

4 渠道邊坡滲流模擬及結果

4.1 不同降雨條件下模擬

在2013 年7 月，萊州連續降雨5 d，累計市平均降雨量達到500 mm 以上；在2022 年由于煙臺地區受第12 號臺風“梅花”影響，9 月14 日6 時-16 時，煙臺市發生一次強降水過程。據市水文中心統計，全市平均降水量208.7 mm，最大點為蓬萊區村里集鎮南官山站490.5 mm。因此本次模擬參考這兩次極端降雨的時段平均降雨量作為降雨過程。固定地下水水位（9.6 m），設置降雨強度分別為5 和20 mm/h，降雨時長分別為96 和24 h，并且在降雨停止后繼續模擬40 h的入滲時長（見表3），采用COMSOL 軟件分別構建孔隙模型和裂隙模型，采用Darcy 滲流模塊模擬，模擬方案如表3所示。

4.1.1 體積含水率隨降雨的變化過程

在自然界的渠道邊坡中，地下水位線以下的土體均處于飽和狀態，地下水位線以上的部分即為非飽和土體，并且含水率會隨著高程的增加而減少。對表3 中的4 種工況分別進行模擬，得到不同降雨強度條件下孔隙模型和裂隙模型在不同時刻的體積含水率分布圖。圖9 為初始時刻體積含水率圖，即深紅色部分為地下水位以下部分，圖10、圖11為雨強為20和5 mm/h不同時刻的體積含水率分布圖。

圖10 降雨強度20 mm/h下的模型體積含水率（地下水位9.6 m）Fig.10 Model volumetric moisture content at rainfall intensity of 20 mm/h（The water table is 9.6 m）

圖11 降雨強度5 mm/h下模型體積含水率（地下水位9.6 m）Fig.11 Moisture content of model volume under rainfall intensity of 5 mm/h（The water table is 9.6 m）

從模擬結果可以明顯看出，圖10 所示為降雨強度20 mm/h時，當降雨時長T=24 h 時，孔隙模型的濕潤峰最低點高程為12 m，而裂隙模型的濕潤峰最低點高程為11.5 m；此時孔隙模型入滲區體積含水率達到0.45 以上，裂隙模型入滲區體積含水率僅有0.35～0.4 左右。當降雨停止繼續入滲40 h 后，兩模型的入滲面積均達到最大，同時部分入滲水補充地下水，裂隙模型的入滲區體積含水率維持在0.25 左右，孔隙模型的入滲區體積含水率維持在0.3。

圖11 表明，降雨5 mm/h 的條件下，當T=96 h 時裂隙模型的濕潤峰運移速度比孔隙模型要快（裂隙模型濕潤峰最低點11 m，孔隙模型濕潤峰最低點11.6 m），此時裂隙模型入滲區體積含水率在0.3，而孔隙模型入滲區體積含水率在0.35；當降雨停止繼續入滲40 h 后，兩模型的入滲面積均達到最大，同時部分入滲水補充地下水，裂隙模型的入滲區體積含水率維持在0.25左右，孔隙模型的入滲區體積含水率維持在0.3。

對比孔隙模型和裂隙模型，可以看出，降雨強度一定情況下，兩種模型的濕潤峰運移速度有明顯的不同，裂隙模型的濕潤峰會比孔隙模型快，即時間相同時，濕潤峰深入深度更大，易造成地下水位的上漲。裂隙模型的體積含水率略小于孔隙模型；這是因為裂隙的存在，雨水進入裂隙后迅速下滲，入滲深度增加較快，同時裂隙內的入滲水量隨著裂隙下滲至地下水位，部分水分補充地下水導致入滲區內的體積含水率略小于孔隙模型。

模擬同時表明，降雨強度的改變對兩模型的入滲區域的體積含水率影響較大，入滲區體積含水率與雨強呈正比。

4.1.2 邊坡內水壓力模擬結果分析

在自然情況下的降雨入滲會造成邊坡內水壓力的變化，特別是坡腳處水壓力的變化容易造成渠道襯砌板鼓脹和脫落。本節選取表3 中工況1、2，針對裂隙的存在與否對邊坡內部水壓力變化的影響進行數值模擬，得到結果如圖12-14所示。

圖12 初始條件下邊坡內水壓力分布云圖（地下水位9.6 m）Fig.12 Distribution cloud diagram of pressure in slope under initial conditions（The water table is 9.6 m）

圖13 裂隙模型降雨20 mm/h邊坡內及坡腳處水壓力分布云圖（地下水位9.6 m）Fig.13 Cloud image of water pressure distribution in slope and foot of slope with 20 mm/h rainfall in fracture model （The water table is 9.6 m）

圖14 孔隙模型降雨20 mm/h邊坡內及坡腳處水壓力分布云圖（地下水位9.6 m）Fig.14 Cloud map of water pressure distribution in the slope and at the foot of the slope with rainfall of 20 mm/h for pore model（The water table is 9.6 m）

圖12所示為兩種模型初始時刻水壓力分布云圖，地下水位以下的壓力值呈現大于0 的數值。圖13（a）是采用裂隙模型在降雨強度為20 mm/h、降雨歷時為24 h 的下滲過程模擬40 h（降雨過程24 h、持續入滲40 h）得到的水壓力分布云圖；圖13（b）為坡腳處水壓力分布的局部放大圖。圖14（a）、圖14（b）采用孔隙模型模擬24 h+40 h 的邊坡水壓力分布云圖和坡腳處水壓力分布的局部放大圖。可以看出，在降雨下滲過程結束后，裂隙模型其壓力分布圖中部分區域的等值線呈現凸起，這說明此時地下水位已經上升至凸起部位（高程12 m處），而孔隙模型這種現象則不明顯。這是由于裂隙模型的濕潤鋒提前到達地下水位，補充地下水位，易導致地下水位提前上漲，而孔隙模型的濕潤峰運移速度比裂隙模型慢，地下水位上漲滯后于裂隙模型。但隨著地下水位的上升可能會導致襯砌板附近的壓力值會變大，從而對襯砌板產生一定的影響破壞。

4.2 襯砌板所受壓力分析

研究區此渠段自通水以來，渠道受揚壓力影響襯砌板遭到破壞，其破壞形式主要表現為塌板、鼓板以及脫落。襯砌的破壞導致渠道的輸水能力下降，嚴重影響工程安全運行。為研究渠道襯砌板受壓情況的變化情況，選取正常地下水水位情況（地下水位9.6 m），對孔隙和裂隙模型的襯砌板壓力情況進行對比分析。

若判斷襯砌板所受壓力是否達到安全壓力值，首先應先對安全壓力值進行計算判斷，襯砌板抗浮穩定性系數可表示為：

按照最不利情況進行計算：檢修期時假設襯砌板頂面水壓力為零，水壓力直接作用于襯砌板底面，此時采用混凝土襯砌板的飽和容重25 kN/m3，襯砌板厚度為0.06 m，查閱規范可知，當安全系數滿足Kf≥1.05 時，得出Δh≤0.14 m（即安全壓力最大值為1.4 kPa）。

取地下水位為9.6 m，根據COMSOL 模擬結果，可得到襯砌板所受壓力如圖15、圖16 和表4 所示。圖中橫坐標X=46.75 m至X=51 m 為渠道襯砌邊坡，X=51 m 至X=55 m 為渠道襯砌底板。

圖15 降雨強度20 mm/h下襯砌板壓力變化曲線Fig.15 Pressure curve of lining plate under rainfall intensity of 20 mm/h

圖16 降雨強度5 mm/h下襯砌板壓力變化曲線Fig.16 Pressure change curve of lining plate under rainfall intensity 5 mm/h

表4 襯砌板所受壓力分析表Tab.4 Pressure analysis table of lining plate

由圖表分析，可以看出對于存在裂隙的渠道，其裂隙的存在會對邊坡襯砌板的安全有影響。同一次降雨，裂隙模型的襯砌板壓力出現最大值處始終提前于孔隙模型，且最大值始終大于孔隙模型，這說明了裂隙的存在使得襯砌板的破壞不僅時間上提前也會發生空間上的擴大。當降雨條件改變時，可以看出高雨強條件下破壞位置比低雨強條件下偏左了0.6 m 左右，說明強降雨條件下襯砌板破壞的范圍更廣。

5 結 論

本文以膠東調水工程萊州段部分渠段為研究對象，基于蒙特卡洛法，生成了二維隨機裂隙網絡，并進行了不同降雨條件下孔隙模型和裂隙模型的滲流模擬，得到如下結論。

（1）在降雨條件不同時，兩模型入滲區域的體積含水率與雨強呈正比，裂隙模型的濕潤峰運移速度大于孔隙模型，與地下水位聯通提前，易造成邊坡處地下水位提前上漲；相同降雨強度下，裂隙模型入滲區域的體積含水率稍小于孔隙模型。

（2）在相同降雨條件、相同地下水位條件下，同一時刻裂隙模型的渠道坡腳處壓力值較孔隙模型的大，更易造成渠道襯砌板破壞。

（3）在相同降雨條件下，裂隙的存在使得襯砌板的破壞不僅時間上提前也會發生空間上的擴大。襯砌板壓力出現最大值處隨著降雨強度的增加其位置會向左偏移，即降雨強度越大，破壞面積越大。

因此，裂隙發育地區的渠道滲流計算應考慮裂隙對滲流、揚壓力的影響。若忽略裂隙滲流的計算，結果可能與實際不符，會給工程運行帶來潛在的安全風險。