裂隙巖質(zhì)邊坡非飽和降雨入滲特征分析

楊秀竹, 康 鏡, 劉正夫, 雷金山

(1.中南大學(xué) 土木工程學(xué)院, 湖南 長沙 410075; 2.江蘇省交通科學(xué)研究院, 江蘇 南京 210004)

?

裂隙巖質(zhì)邊坡非飽和降雨入滲特征分析

楊秀竹1, 康 鏡1, 劉正夫2, 雷金山1

(1.中南大學(xué) 土木工程學(xué)院, 湖南 長沙 410075; 2.江蘇省交通科學(xué)研究院, 江蘇 南京 210004)

[目的] 研究降雨入滲過程中，裂隙巖質(zhì)邊坡滲流場和暫態(tài)飽和區(qū)的變化規(guī)律，為裂隙巖體降雨入滲規(guī)律、工程技術(shù)及水土保持方面的研究提供支持。 [方法] 基于非連續(xù)裂隙模型，采用Matlab編制裂隙巖體滲透張量的計算程序；基于等效連續(xù)介質(zhì)滲流模型，采用FLAC3D內(nèi)置的FISH語言編制程序進(jìn)行非飽和滲流分析。 [結(jié)果] 依托實際工程分析得出了裂隙巖體的滲透張量，在降雨過程中隨著降雨強(qiáng)度的增加或者降雨時間的延長，邊坡表層區(qū)域零孔隙水壓力面不斷向內(nèi)部發(fā)展，暫態(tài)飽和區(qū)的范圍不斷擴(kuò)大，且暫態(tài)飽和區(qū)的出現(xiàn)發(fā)生在降雨強(qiáng)度大于邊坡的入滲率的情況。 [結(jié)論] 采用非連續(xù)裂隙網(wǎng)絡(luò)模型和等效連續(xù)介質(zhì)滲流模型對裂隙巖質(zhì)邊坡降雨入滲分析能較好反映裂隙巖體各向異性和非飽和滲流的特點。

邊坡; 降雨; 滲透張量; 等效連續(xù)介質(zhì); 非飽和滲流

文獻(xiàn)參數(shù)： 楊秀竹, 康鏡, 劉正夫, 等.裂隙巖質(zhì)邊坡非飽和降雨入滲特征分析[J].水土保持通報，2016,36(4)：143-147.DOI:10.13961/j.cnki.stbctb.2016.04.026

當(dāng)巖體破碎時，巖質(zhì)邊坡易受風(fēng)化降雨作用的影響。隨著降雨過程的發(fā)展，雨水入滲會引起邊坡非飽和區(qū)的負(fù)孔隙水壓力發(fā)生較大變化，該部分巖體的自重增加，并且使巖體軟弱結(jié)構(gòu)面強(qiáng)度下降，引起邊坡的失穩(wěn)。從而導(dǎo)致引起邊坡的失穩(wěn)的可能性增加，為工程及群眾生命財產(chǎn)安全造成較大威脅。因此，對降雨時裂隙巖體中的滲流特征以及滲流場的變化規(guī)律進(jìn)行研究，對降低安全風(fēng)險，防止滑坡災(zāi)害發(fā)生有較為重要的意義。滲流介質(zhì)是影響邊坡降雨入滲重要因素之一，劉子振[1]研究了降雨對非飽和黏土邊坡的入滲規(guī)律，但與研究土質(zhì)邊坡降雨入滲不同，裂隙巖體滲流具有的顯著各向異性特點是不能忽略的。滲透系數(shù)影響分析結(jié)果的正確性，張升堂等[2]用年降雨總量與產(chǎn)流降雨總量的統(tǒng)計推求降雨入滲系數(shù)，認(rèn)為產(chǎn)流降雨總量能比較清晰地反映流域治理程度對降雨入滲的影響，高鵬等[3]采用人工模擬降雨的試驗方法研究了不同降雨強(qiáng)度對土壤入滲的影響，建立了土壤水分入滲經(jīng)驗?zāi)Ｐ汀５€須從滲流介質(zhì)及其結(jié)構(gòu)特點的角度進(jìn)一步研究。分析裂隙巖質(zhì)邊坡降雨入滲的模型主要有：等效連續(xù)介質(zhì)模型、離散裂隙網(wǎng)絡(luò)模型和雙重介質(zhì)模型[4]。離散裂隙網(wǎng)絡(luò)模型認(rèn)為巖塊不能導(dǎo)水而忽略巖塊的導(dǎo)水性，但當(dāng)考慮裂隙巖體非飽和滲流時，巖塊的導(dǎo)水性是不能忽略的；雙重介質(zhì)模型中裂隙網(wǎng)絡(luò)和連續(xù)介質(zhì)巖塊耦合模型中的水交換量公式直接影響模型計算的準(zhǔn)確性，但水交換量難以準(zhǔn)確確定；當(dāng)巖體裂隙節(jié)理等結(jié)構(gòu)面充分發(fā)育時，等效連續(xù)介質(zhì)模型采用滲透張量矩陣能較好地模擬裂隙巖體滲流各向異性的特點，并且該模型只需知道裂隙的幾何信息的統(tǒng)計值，而對其準(zhǔn)確位置并不敏感。因此，本研究采用非連續(xù)裂隙網(wǎng)絡(luò)模型確定裂隙巖體滲透張量并通過等效連續(xù)介質(zhì)滲流模型對裂隙巖質(zhì)邊坡非飽和降雨入滲進(jìn)行分析，其研究成果可為裂隙巖體降雨入滲規(guī)律研究，工程技術(shù)及水土保持研究提供支持和參考。

1　等效連續(xù)介質(zhì)滲流原理

(1)

需要指出的是對于巖體而言，將裂隙等效為連續(xù)介質(zhì)模型必須對裂隙巖體進(jìn)行分析，判斷其是否滿足等效條件。Long[5]指出了等效裂隙巖體滿足：表征單元體REV存在；存在等效滲透張量，且裂隙巖體方向滲透率在極坐標(biāo)下能形成橢圓，即該等效滲透張量是對稱，才可以運用等效連續(xù)介質(zhì)模型。

2　滲透張量的確定

滲透張量是等效連續(xù)介質(zhì)滲流計算模型中的一個關(guān)鍵參數(shù)，它直接反映了裂隙巖體非均質(zhì)各向異性的特點。采用Matlab編制基于裂隙網(wǎng)格模型的滲透張量計算程序包括：運用蒙特卡洛法生成裂隙，并對裂隙進(jìn)行篩選，最終得到有效裂隙；計算不同大小分析區(qū)域中的滲透張量，并給出最小REV和對應(yīng)的滲透張量[6]。

2.1裂隙生成

生成裂隙采用的蒙特卡洛方法是一種依據(jù)大數(shù)定律的統(tǒng)計試驗方法，它采用隨機(jī)數(shù)進(jìn)行統(tǒng)計試驗，產(chǎn)生符合某種隨機(jī)分布函數(shù)的隨機(jī)數(shù)。在裂隙巖體的滲流分析中，可以根據(jù)裂隙幾何參數(shù)的概率模型[7]，在滿足給定的空間分布規(guī)律的情況下，隨機(jī)生成裂隙的長度、方向角、中心點位置。采用蒙特卡洛法生成裂隙的具體流程如圖1所示。

圖1　采用蒙特卡洛法生成裂隙流程圖

2.2滲透張量計算

在生成有效裂隙網(wǎng)格的基礎(chǔ)上，滲透張量的計算流程如圖2所示。

圖2　滲透張量計算流程圖

2.3滲透張量計算驗證

某巖體的生成域為30 m×30 m，分析域為15 m×15 m，共有2組裂隙，裂隙的中心點服從均勻分布，方向角和跡長均服從正態(tài)分布，隙寬為常數(shù)；具體裂隙參數(shù)詳見表1。

表1　試驗巖體裂隙參數(shù)

對裂隙進(jìn)行篩選之后的有效裂隙表征單元體尺寸為 ，在該尺寸內(nèi)可認(rèn)為裂隙是貫通的，該裂隙滲透張量系數(shù)為：

當(dāng)采用裂隙測量法時，計算出滲透系數(shù)張量：

對比兩種算法計算出來的滲透系數(shù)張量，發(fā)現(xiàn)兩者很接近，說明本研究對滲透張量的計算模塊是合理可靠的。

3　FLAC3D飽和-非飽和降雨入滲分析功能的二次開發(fā)

FLAC3D是目前國內(nèi)外巖土工程界應(yīng)用廣泛的有限差分軟件。但由于其在非飽和滲流計算中會將負(fù)孔隙水壓區(qū)的飽和度強(qiáng)制設(shè)為1[8]，這與實際情況不符，部分學(xué)者改進(jìn)FLAC3D飽和非飽和滲流分析模塊[9]，取得了一些進(jìn)展。因此，為了利用FLAC3D軟件平臺進(jìn)行裂隙巖體非飽和滲流分析，利用內(nèi)置的FISH語言對FLAC3D滲流模塊的二次開發(fā)，在迭代過程中根據(jù)飽和度不斷地修正非飽和區(qū)的滲透張量，求出比傳統(tǒng)滲流分析更準(zhǔn)確的滲流場。對FLAC3D的二次開發(fā)可分為兩個部分，即飽和—非飽和的滲流分析功能，以及降雨入滲及出滲邊界功能的FISH函數(shù)編制[10]。

進(jìn)行飽和—非飽和滲流分析時，在FLAC3D的計算過程中，自動根據(jù)上一計算增量時間步的結(jié)果，調(diào)整非飽和區(qū)的滲透張量，從而實現(xiàn)FLAC3D的非飽和滲流計算。本研究采用Genuchten[11]提出單裂隙毛管壓力—飽和度關(guān)系曲線和相對滲透率—飽和度關(guān)系曲線：

(7)

(8)

式中：Pc——毛管壓力； P0——進(jìn)氣值； Se——水相的有效飽和度； m,n——經(jīng)驗參數(shù)； kr——相對滲透率。

進(jìn)行降雨入滲和出滲邊界分析時，降雨入滲率取飽和滲透系數(shù)和降雨強(qiáng)度中的較小值。降雨入滲為動態(tài)的邊界條件，在降雨過程的分析中，需要在每一個時間步對邊坡表層節(jié)點的孔隙水壓進(jìn)行判斷，若孔隙水壓小于0，節(jié)點應(yīng)采用流量邊界；而當(dāng)孔隙水壓大于或等于0時，則需將該節(jié)點由流量邊界修改為已知水頭邊界。

4　算例分析

4.1工程概況

長沙大冰雪世界礦坑酒店項目位于長沙市岳麓區(qū)坪塘鎮(zhèn)，邊坡主要為巖質(zhì)邊坡，巖石為古生代泥盆系中統(tǒng)棋梓橋組淺灰色—深灰色薄層—厚層灰?guī)r，厚度大于210m。邊坡大部分基巖裸露，邊坡巖體表面裂隙性溶蝕風(fēng)化現(xiàn)象嚴(yán)重，沿斷層、裂隙及層面等結(jié)構(gòu)面溶蝕風(fēng)化現(xiàn)象較普遍，風(fēng)化裂隙發(fā)育。據(jù)1960—2013年長沙市氣象站資料統(tǒng)計：平均降雨量1 394.6mm，最大年降雨量1 751.2mm(1998年)，最小年降雨量708.8mm(1953年)，最大月降雨量515.3mm，最小月降雨量1.2mm，最大日降雨量192.5mm。

4.2滲透張量求解

研究區(qū)域邊坡巖體主要發(fā)育4組節(jié)理裂隙，將真實的裂隙節(jié)理進(jìn)行轉(zhuǎn)換，得到其在特征截面上的跡線長度和方向角，運用蒙特卡洛法隨機(jī)生成二維裂隙網(wǎng)絡(luò)所需的裂隙幾何參數(shù)統(tǒng)計詳見表2。

將以上數(shù)據(jù)輸入Matlab程序中，便可在 的裂隙生成域內(nèi)生成二維裂隙網(wǎng)絡(luò)模型(裂隙網(wǎng)絡(luò)圖略)。

表2　巖體裂隙幾何參數(shù)統(tǒng)計

為確定裂隙巖體的表征單元體(REV)尺寸，在上述全局域裂隙網(wǎng)絡(luò)的基礎(chǔ)上，對不同分析域尺寸下的滲透張量和對應(yīng)的均方誤差(RMS)進(jìn)行試算。分析域每次旋轉(zhuǎn)30°，共旋轉(zhuǎn)12次，并賦予相同的邊界條件；通過程序計算，現(xiàn)將不同尺寸下的RMS整理如圖3所示。

Johan[12]根據(jù)計算的RMS判斷擬合的橢圓質(zhì)量，認(rèn)為只有當(dāng) 時，各個方向滲透系數(shù)擬合的橢圓質(zhì)量較好，采用等效連續(xù)介質(zhì)模型計算時不會出現(xiàn)大的誤差。此時，該區(qū)域的尺寸可以視為該裂隙巖體的表征單元體(REV)。當(dāng)分析域尺寸小于5 m時，相應(yīng)的RMS大于0.2，其相應(yīng)擬合的橢圓質(zhì)量較差；當(dāng)分析域尺寸在5 m以上時，RMS小于0.2(穩(wěn)定在0.16附近)，擬合的橢圓質(zhì)量較好，表征單元體存在。參照

圖3，表征單元體尺寸取為6 m×6 m。

通過裂隙網(wǎng)絡(luò)水力學(xué)的方法求出6 m×6 m分析域內(nèi)裂隙巖體的沿梯度方向的滲透系數(shù)詳見表3。

圖3　不同分析域尺寸下的均方誤差值(RMS)變化

角度/(°)滲透系數(shù)/(m·s-1)橢圓半徑(1-m/s)角度/(°)滲透系數(shù)/(m·s-1)橢圓半徑(1/m/s)05.51E-06426.191805.51E-06426.19307.07E-06376.152107.07E-06376.15609.44E-06325.472409.44E-06325.47901.36E-05270.752701.36E-05270.751201.08E-05304.683001.08E-05304.681506.55E-06390.753306.55E-06390.75

通過最小二乘法將上表中的數(shù)據(jù)在極坐標(biāo)下擬合成滲透橢圓。擬合的橢圓的長短軸計算出滲透張量主值K1=1.33×10-5m/s，K2=5.54×10-6m/s，擬合的滲透橢圓的長軸方向與極坐標(biāo)的0°方向夾角為93.8°, 小于0.2，故可用這個滲透橢圓求出滲透張量來對裂隙網(wǎng)絡(luò)模型的滲透性進(jìn)行描述。考慮到裂隙部分填充、局部非連通以及裂隙測量開度變化等因素影響，計算并修正后的滲透張量(m/s)為：

4.3降雨滲流數(shù)值計算

根據(jù)實際工程地質(zhì)條件，建立邊坡降雨滲流計算模型，X坐標(biāo)在坡腳處沿水平方向指向坡內(nèi)，Z方向在坡腳處沿豎直方向向上為正，計算模型地下水位左側(cè)距模型底部為18 m，右側(cè)為9 m。模型中的人工填土、殘積粉質(zhì)黏土、灰?guī)rⅡ均按地勘報告中滲透系數(shù)選取，灰?guī)rⅥ滲透性由上文分析確定。計算模型網(wǎng)格數(shù)為869，節(jié)點個數(shù)為1 789個。計算模型如圖4所示。

4.3.1降雨持續(xù)時間對暫態(tài)飽和區(qū)的影響為分析降雨持續(xù)時間對零孔隙水壓面的影響，計算邊坡在降雨強(qiáng)度為300 mm/d情況下，降雨持續(xù)時間為1～4 d零孔隙水壓位置的變化圖，將數(shù)據(jù)導(dǎo)入TECPLOT，繪制不同降雨時間下零孔隙水壓面。分析得出，在降雨強(qiáng)度為300 mm/d時，隨著降雨時間的延長，灰?guī)rⅥ表層區(qū)域和斷層的零孔隙水壓面向內(nèi)部發(fā)展的范圍不斷擴(kuò)大，暫態(tài)飽和區(qū)的范圍不斷擴(kuò)大，且相比于灰?guī)rⅥ表層區(qū)域斷層處受到降雨持續(xù)時間的影響更大。但上覆土層的零孔隙水壓面增加區(qū)域有限。降雨時間在1～4 d范圍內(nèi)，對邊坡內(nèi)部的浸潤面的影響不大。

圖4　邊坡降雨滲流計算模型

4.3.2降雨強(qiáng)度對暫態(tài)飽和區(qū)的影響為了分析不同降雨強(qiáng)度對孔隙水壓的影響，分別對降雨持續(xù)4 d降雨強(qiáng)度為100，200，300 mm/d這3個工況下的邊坡滲流進(jìn)行計算，將計算結(jié)果導(dǎo)入TECPLOT，繪制出的不同情況下零孔隙水壓面。分析得出，降雨時間為4 d時，隨著降雨強(qiáng)度增大，灰?guī)rⅥ表層區(qū)域的零孔隙水壓面向內(nèi)部發(fā)展的范圍不斷擴(kuò)大，暫態(tài)飽和區(qū)的范圍不斷擴(kuò)大，但上覆土層的零孔隙水壓面增加區(qū)域有限。當(dāng)降雨強(qiáng)度為100 mm/d時，邊坡表面剛開始出現(xiàn)零孔壓，若降雨強(qiáng)度小于邊坡入滲率，雨水將全部入滲。當(dāng)降雨強(qiáng)度為200 mm/d時，斷層處還未出現(xiàn)零孔隙水壓面，說明斷層處仍處于非飽和狀態(tài)，但隨著降雨強(qiáng)度的增大，當(dāng)降雨強(qiáng)度為300 mm/d時，斷層F1區(qū)域內(nèi)出現(xiàn)零孔隙水壓，達(dá)到了飽和且影響范圍遠(yuǎn)大于灰?guī)rⅥ。在給定降雨時間為4 d時，降雨強(qiáng)度在100～300 mm/d范圍內(nèi)，對邊坡內(nèi)部的浸潤面的影響不大。

5　結(jié) 論

(1) 利用MATLAB軟件，編寫了基于蒙特卡洛方法的隨機(jī)生成裂隙程序和滲透張量計算程序。基于實際工程情況，通過模擬計算，得到了其裂隙分布，并計算出了表征單元體尺寸和滲透張量。

(2) 在模擬裂隙巖體的降雨入滲過程中引入滲透張量，并利用FISH語言編寫了非飽和入滲分析的修正模塊，完善了FLAC3D的各向異性非飽和滲流分析功能，為利用FLAC3D進(jìn)行復(fù)雜三維裂隙巖質(zhì)邊坡降雨入滲研究提供了基礎(chǔ)。

(3) 通過實際工程的研究發(fā)現(xiàn)，降雨強(qiáng)度大于邊坡的入滲率，可能導(dǎo)致暫態(tài)飽和區(qū)的出現(xiàn)。隨著降雨強(qiáng)度的增加或者降雨時間的延長，邊坡表層區(qū)域零孔隙水壓力面不斷向邊坡內(nèi)部發(fā)展，暫態(tài)飽和區(qū)的范圍不斷擴(kuò)大，但降雨過程中對邊坡內(nèi)部的浸潤面的影響較小。

[1]劉子振.持續(xù)降雨入滲非飽和黏土邊坡失穩(wěn)機(jī)理及其應(yīng)用研究[D].甘肅 蘭州:蘭州大學(xué),2014.

[2]張升堂,拜存有,萬三強(qiáng).黃土高原水土保持措施強(qiáng)化降雨入滲分析及灰色預(yù)測[J].水土保持通報,2004,24(2)：29-33.

[3]高鵬,穆興民,劉普靈,等.降雨強(qiáng)度對黃土區(qū)不同土地利用類型入滲影響的試驗研究[J].水土保持通報,2006,26(3)：1-5.

[4]張有天.巖石水力學(xué)與工程[M].北京：中國水利水電出版社,2005.

[5]Long J C S, Remer J S, Wilson C R, et al. Porous media equivalents for networks of discontinuous fractures[J]. Water Resources Research, 1982,18(3):645-658.

[6]劉海軍.基于蒙特卡羅法的巖體裂隙網(wǎng)絡(luò)模型及滲透張量的研究[D].黑龍江 哈爾濱：哈爾濱工業(yè)大學(xué),2011.

[7]劉曉麗,王恩志,王思敬.裂隙巖體表征方法及巖體水力學(xué)特性研究[J].巖石力學(xué)與工程學(xué)報,2008,27(9):1814-1821.

[8]陳育民. FLAC/FLAC3D基礎(chǔ)與工程實例[M].北京：中國水利水電出版社,2009.

[9]李毅,伍嘉,李坤.基于FLAC3D的飽和—非飽和滲流分析[J].巖土力學(xué),2012,33(2):617-622.

[10]蔣忠明,熊小虎,曾鈴.基于FLAC3D平臺的邊坡非飽和降雨入滲分析[J].巖土力學(xué),2014,35(3):855-861.

[11]Genuchten M T V. A Closed-form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils[J]. Soil Science Society of America Journal, 1980,44(5):892-898.

[12]柴軍瑞,徐維生.大壩工程滲流非線性問題[M].北京：中國水利水電出版社,2010.

Unsaturated Seepage Analysis of Fractured Rock Slope Under Rainfall Condition

YANG Xiuzhu1, KANG Jing1, LIU Zhengfu2, LEI Jinshan1

(1.SchoolofCivilEngineering，CentralSouthUniversity,Changsha,Hu’nan410075,China; 2.JiangsuTransportationInstituteCo.Ltd.,Nanjing,Jiangsu210017,China)

[Objective] Exploring the variation process of seepage field and transient saturated zone of fractured rock slope under the influence of rainfall to provide support for the study of the law of rainfall infiltration, engineering technology and water conservation. [Methods] Based on discontinuous fracture network model, a Matlab program was designed to calculate the permeability tensor. Based on equivalent continuous medium model, a program using FISH language was designed for unsaturated seepage analysis. [Results] Permeability tensor of fractured rock mass was obtained through the engineering case. With the increase of rainfall intensity and rainfall duration, the surface of zero pore water pressure will extend to the inside of the slope and the scope of transient saturated zone will expand. Moreover, the appearance of transient saturated zone will occur in the case when rainfall intensity is greater than the the infiltration rate of the slope. [Conclusion] The characteristics of anisotropic and unsaturated seepage in fractured rock mass can be better elucidated by using discontinuous fracture network model and equivalent continuous medium model and the models can be used to analyze the rainfall infiltration of fractured rock slope.

slope; rainfall; permeability tensor; equivalent continuous medium; unsaturated seepage

2015-02-22

2016-03-19

湖南省發(fā)改委項目“湖南城際軌道交通沿線隱伏型巖溶注漿充填材料應(yīng)用研究”(湘發(fā)改高技[2012]1493號)

楊秀竹(1972—),女(漢族),山東省萊州市人,博士,副教授,主要從事巖石力學(xué)、隧道與地下工程的教學(xué)科研工作。E-mail:xzyang1661@126.com，134812270@csu.edu.cn。

雷金山(1973—),男(漢族),湖南省湘鄉(xiāng)市人,博士,高級工程師,主要從事地質(zhì)災(zāi)害治理的研究工作。E-mail:jslei8522@126.com。

A

1000-288X(2016)04-0143-05

TV223.6