張文捷,程榮蘭,詹美禮,盛金昌,王國進

(1.河海大學(xué)水利水電工程學(xué)院,江蘇 南京 210098;2.江西省水利廳,江西 南昌 330009;3.云南省電力設(shè)計院,云南 昆明 650051;4.昆明勘測設(shè)計研究院,云南 昆明 650051)

巖體滲流是影響高壩等工程安全的關(guān)鍵要素之一.巖體滲流不同于多孔介質(zhì)滲流[1-7],巖體由裂隙和被裂隙切割的巖塊構(gòu)成,巖塊的滲透性極弱,裂隙的滲透性極強,裂隙滲流不服從達西定律.國內(nèi)外眾多學(xué)者對巖體滲流進行了大量的研究,雖然在理論和數(shù)值模擬方面取得了許多成果,但是巖體滲流研究理論還不十分成熟.常用的數(shù)學(xué)模型主要有4種:等效連續(xù)介質(zhì)模型[8-9];離散裂隙網(wǎng)絡(luò)模型[10];雙重介質(zhì)模型(或多重介質(zhì)模型)[11];離散-連續(xù)耦合模型[12].這4種模型的不足之處是:(a)擬真性較差.離散裂隙網(wǎng)絡(luò)模型中裂隙分布的離散、無規(guī)律性決定了目前探測技術(shù)無法較完全地統(tǒng)計出計算所需的所有裂隙集合產(chǎn)狀,故使得數(shù)學(xué)模型擬真性降低;連續(xù)介質(zhì)模型不考慮裂隙與孔隙之間的水交換,將裂隙巖體也作為連續(xù)介質(zhì)處理,裂隙中的水流量被等效平均到整個巖體中,故該模型擬真性也較差.(b)水量交換項難以確定(如離散裂隙網(wǎng)絡(luò)模型、雙重介質(zhì)模型和離散-連續(xù)耦合模型).(c)無法反映巖體中真實存在的水頭間斷現(xiàn)象,更無法定量反映出水頭間斷程度.

為了探求一種合理反映巖體滲流的數(shù)學(xué)模型,本文從裂隙巖體離散-連續(xù)介質(zhì)滲流機理出發(fā),在2種強弱透水接觸面上引入無厚度虛擬滲透突變單元(或稱水頭突變單元),根據(jù)強弱透水介質(zhì)接觸面處節(jié)點流量相等建立改進的離散-連續(xù)介質(zhì)耦合模型,提出了反映水頭間斷程度的連續(xù)系數(shù)概念,并分析了合理的連續(xù)系數(shù)取值規(guī)律.

1 三維滲透突變單元有限元模型

1.1 滲透突變單元在巖體滲流數(shù)值模擬中的意義

在裂隙巖體滲流中,巖體塊體滲透性很小,基本只起儲水作用,裂隙是巖體的主要透水通道.在強導(dǎo)水的大裂隙壁面附近往往有很大的水頭跌落,形成較大的水頭間斷.

為了真實地反映巖體中滲流場分布的間斷突變性,本文借鑒位移突變單元的思想.這既能夠解決連續(xù)介質(zhì)模型中忽略巖塊與裂隙之間的水量交換的不足,又能夠彌補雙重介質(zhì)模型中水量交換量難確定的缺陷,可以比離散網(wǎng)絡(luò)模型更能夠高效、真實地模擬巖體滲流.引入滲透突變單元的基本思想是:在主干裂隙離散域(或斷層破碎帶)與巖體域接觸面上加入8節(jié)點無厚度虛擬單元.鑒于此無厚度單元在巖體滲流模擬中的意義,本文稱其為滲透突變單元.滲透突變單元是常規(guī)有限元的一種特殊單元模式.用滲透突變單元模擬水頭間斷,其前提條件是不改變整體滲流場流量,在等效節(jié)點流量相等原則下,通過在2種介質(zhì)接觸面上引入滲透突變單元來實現(xiàn)巖塊與裂隙(斷層)之間的水量交換,并模擬2種介質(zhì)交界面處水頭突變現(xiàn)象.

1.2 三維滲透突變單元計算理論與滲透機理分析

本文借鑒非連續(xù)變形Goodman接觸面的處理方式,在裂隙和巖體的交界面上引入滲透突變單元.該單元為8節(jié)點無厚度單元(圖1),其法向厚度為0,故只考慮法向滲透性,而忽略切向滲透性.

從數(shù)學(xué)間斷角度出發(fā),圖1中滲透突變單元左側(cè)面節(jié)點物理量 Φ-和右側(cè)面對應(yīng)節(jié)點物理量 Φ+不相等,即體現(xiàn)為同一點上物理量間斷不連續(xù).引入具體的滲流場分析,滲透突變單元左右側(cè)面上的不連續(xù)的物理量即為水頭H.在流量連續(xù)原則上,保證不改變整體滲流場流量,將圖1中的滲透突變單元運用到滲流分析中,當滲透突變單元完全不透水時,在滲透突變單元左、右側(cè)面上水頭絕對間斷;隨著滲透突變單元不同程度破損,水頭間斷程度隨之減小,當達到一定破損度后,滲透突變單元失效,左右兩側(cè)水頭值大致相等.

根據(jù)有限元基本理論,引入插值函數(shù)Ni,構(gòu)造函數(shù),令Mi=-Ni,Mi+4=-Ni+4,i=1,2,3,4,得到單元內(nèi)各節(jié)點的水頭差

圖1 無厚度滲透突變單元Fig.1 Discontinuous element with zero thickness

式中:ΔH——水頭差;H+,H-——滲透突變單元左、右側(cè)面的分布水頭;Hi——i節(jié)點水頭.

滲透突變單元法向厚度為0,故只考慮法向滲透性.借鑒土工膜防水性能優(yōu)劣判別思想[13],即當防滲土工膜的厚度比較小時,無法僅從滲透系數(shù)來衡量透水性,故將滲透系數(shù)和薄膜厚度一起考慮,常用透水率來衡量土工膜的防水優(yōu)劣程度.由于本文引入的滲透突變單元是一種虛擬的無厚度單元,故此也用虛擬滲透率來衡量滲透突變單元滲透性強弱.根據(jù)達西定律,當L=0時,流速只與虛擬滲透系數(shù)和滲透突變單元的水頭差有關(guān),即

式中:V——滲透流速;Kn——滲透突變單元的法向虛擬滲透率.

改進模型的等效節(jié)點流量表達式為

式中:K——滲透突變單元局部坐標系下的滲透矩陣;Qi——節(jié)點 i處的流量;Ωe——單元區(qū)域;s,t,n——3個方向的坐標.

滲透突變單元在本文數(shù)學(xué)模型中的作用之一是通過調(diào)整滲透率來恰當反映不同程度的水頭突變間斷現(xiàn)象,故結(jié)合間斷系數(shù)的數(shù)學(xué)、物理意義及在本文數(shù)學(xué)模型中的作用,這里將滲透突變單元虛擬滲透率Kn定義為連續(xù)系數(shù).Kn的合理取值是水量交換模擬合理正確的關(guān)鍵,也是正確反映水頭間斷程度的關(guān)鍵因素.

上述推導(dǎo)是在局部坐標系下進行的,為了集成整體滲透矩陣,需要將局部坐標轉(zhuǎn)化為整體坐標.最終得到的整體坐標下的節(jié)點水頭與節(jié)點流量關(guān)系式為

2 算 例

選擇一理想的水文地質(zhì)體進行研究.模型尺寸長×寬×高為100m×50m×50m,三維計算網(wǎng)格如圖2所示.沿地質(zhì)體長度方向(x方向)取左側(cè)的50m為弱透水介質(zhì),右側(cè)50m為強透水破碎帶.在強弱透水介質(zhì)接觸面x=50m的位置加入滲透突變單元.根據(jù)計算需求,對模型網(wǎng)格局部加密剖分,模型共劃分為20400個單元,塊體單元20000個,滲透突變單元400個,節(jié)點22932個.邊界條件:上、下游水位分別為49m和5m,兩者作為已知水頭邊界(圖3),位于下游水位以上的右側(cè)面邊界,視為可能出滲邊界.

圖2 滲流場計算網(wǎng)格模型Fig.2 Meshes for calculation of seepage field

圖3 透水介質(zhì)分布和邊界條件Fig.3 Distribution of permeable media and boundary conditions

2.1 模型驗證

擬定弱透水介質(zhì)的滲透系數(shù)K弱為1.0×10-6cm/s,強透水介質(zhì)的滲透系數(shù)K強為1.0×10-2cm/s.本文數(shù)學(xué)模型的驗證通過與常規(guī)數(shù)學(xué)模型計算結(jié)果的對比來實現(xiàn).

a.滲流場計算結(jié)果比較分析.沿 x方向取一典型剖面,將常規(guī)數(shù)學(xué)模型與本文數(shù)學(xué)模型計算得到的滲流場等值線以及侵潤面繪于同一張圖上,如圖4所示.從圖4可以看出,常規(guī)數(shù)學(xué)模型基本能反映該巖體中的滲流場分布,但是,在強弱透水介質(zhì)的接觸面(x=50m)處,水頭未在此處急劇下跌,而是在接觸面左側(cè)3～5m范圍內(nèi)緩慢降低,如浸潤線AB段,浸潤面凹至弱透水介質(zhì)中,這是計算的不合理之處,且從計算結(jié)果來看,常規(guī)數(shù)學(xué)模型未能反映強弱透水接觸面處的水頭突變現(xiàn)象.運用本文改進的離散-連續(xù)介質(zhì)耦合模型,通過調(diào)整連續(xù)系數(shù)取值,使得浸潤線從強弱透水介質(zhì)接觸面位置C點跌至B點,形成急劇變化的水頭間斷.從水頭間斷吻合度來看,BC段能較好地反映此水文地質(zhì)體中接觸面處的水頭間斷現(xiàn)象.

圖4 滲流場等值線及侵潤面計算結(jié)果比較Fig.4 Comparison of contours and phreatic faces of seepage field of conventional and proposed models

b.流量分析.在x=10m處截取一剖面,常規(guī)數(shù)學(xué)模型得到的流量為2.02m3/d,而本文數(shù)學(xué)模型得到的流量為1.90m3/d,2種模型的流量僅相差5.4%.這表明,滲透突變單元沒有改變整體滲流量,滲透突變單元模型是正確的.

上述分析表明,改進的離散-連續(xù)介質(zhì)耦合模型能反映復(fù)雜巖體結(jié)構(gòu)中的滲流場分布,可以通過改變滲透突變單元連續(xù)系數(shù)來控制水量交換,比雙重介質(zhì)模型確定水量交換項要快捷,且有常規(guī)數(shù)學(xué)模型流量和水頭間斷程度作為連續(xù)系數(shù)取值衡量指標,確保了水量交換的模擬正確性.此外,由于滲透突變單元特殊的單元性質(zhì),通過合理的連續(xù)系數(shù)取值,滲透突變單元又能定性、定量地反映出強弱透水介質(zhì)接觸面上的水頭突變間斷程度,這一點也是幾種常規(guī)數(shù)學(xué)模型無法模擬的.

2.2 水頭間斷程度與連續(xù)系數(shù)取值規(guī)律

2.2.1 K強/K弱=10000時

強弱透水介質(zhì)滲透系數(shù)取值分別為1.0×10-2cm/s和1.0×10-6cm/s,在x=10m處截取一平面,分別計算有無滲透突變單元時通過該截面的流量,統(tǒng)計結(jié)果見表1,滲流場分布如圖5示.

從圖5可以看出:當Kn=1.0×10-9s-1時,滲透突變單元相當于不透水,起到很大的阻水作用,滲流場與加入滲透突變單元前反映的滲流場嚴重不符;當Kn=1.0×10-7s-1時,浸潤線4的 CB段能解決AB段傾入弱透水介質(zhì)中的不合理之處,同時,水頭在CB段陡然下跌,反映了水在2種滲透性差異大的介質(zhì)交界面水頭不連續(xù)現(xiàn)象;當Kn=1.0×10-5s-1時,滲透突變單元的作用已體現(xiàn)不出.因此,Kn=1.0×10-7s-1,即小于K弱的1個數(shù)量級時,改進的離散-連續(xù)介質(zhì)耦合模型能合理反映K強/K弱=10000時的水頭間斷規(guī)律.

表1 不同Kn時加滲透突變單元前、后斷面流量比較Table1 Comparison of flux amounts with and without discontinuous elements with different Kn

從流量可以看出:當Kn=1.0×10-9s-1時,有無滲透突變單元時通過斷面的流量相差近92%,此時完全不符合不改變整體滲流場流量的條件;當Kn=1.0×10-7s-1時,二者流量僅相差5.4%;當Kn=1.0×10-6s-1時,可看作突變單元大部分破損,此時通過斷面流量基本等于加滲透突變單元前流量.這表明,當Kn取小于K弱1個數(shù)量級時,滲流場計算結(jié)果合理.由此可見,在K強/K弱=10000,且Kn小于K弱的1個數(shù)量級的條件下,改進的離散-連續(xù)介質(zhì)耦合模型能合理反映滲流場水頭間斷現(xiàn)象.

2.2.2 K強/K弱=1000,500,100,50,10時

K弱=1.0×10-6cm/s,K強根據(jù)不同要求相應(yīng)改變.各種連續(xù)系數(shù)取值條件下的滲流等值線如圖6和圖7所示.從流量和滲流場的比較可以看出:當K強/K弱=1000時,Kn的合理取值范圍是1.0×10-7s-1≤Kn≤3.0×10-7s-1;當K強/K弱=500,100,50時,Kn均取K弱數(shù)量級的0.3～0.5倍,即在本例中取 3.0×10-7s-1≤Kn≤5.0×10-7s-1;當K強/K弱=10時,Kn和K弱為同一個數(shù)量級,即在本例中取Kn=1.0×10-6s-1.

圖5 K強=10000K弱條件下的滲流場等值線Fig.5 Contours of seepage field when K強=10000K弱

圖6 K強=1000K弱條件下的滲流場等值線Fig.6 Contours of seepage field when K強=1000K弱

圖7 不同K強/K弱條件下的滲流場等值線Fig.7 Contours of seepage field with different cases of K強/K弱

2.3 連續(xù)系數(shù)取值規(guī)律

取K弱的數(shù)量級為m,突變單元的滲透性量級為n;n=m/10,即n小于m 1個數(shù)量級.通過對強、弱透水介質(zhì)滲透性不同差異倍比條件下多組合計算分析,得出水頭突變單元連續(xù)系數(shù)Kn取值大小和突變介質(zhì)滲透性倍比之間的關(guān)系,如表2所示.

表2 突變介質(zhì)滲透性倍比與水頭突變單元連續(xù)系數(shù)關(guān)系Table2 Relationship between permeability ratio of discontinuous media and continuous coefficient of head discontinuous elements

3 結(jié) 語

本文從離散-連續(xù)介質(zhì)滲透機理出發(fā),在接觸面處引入無厚度滲透突變單元,通過這種特殊模式的單元與常規(guī)單元的耦合,改進和完善了常規(guī)巖體滲流數(shù)學(xué)模型,改進后的巖體滲流數(shù)學(xué)模型(稱為改進數(shù)學(xué)模型)可以更加真實地反映巖體滲流特性;采用常規(guī)數(shù)學(xué)模型和改進數(shù)學(xué)模型對同一地質(zhì)體進行了滲流計算,計算結(jié)果驗證了改進數(shù)學(xué)模型的正確性和合理性;提出了反映水頭間斷程度的連續(xù)系數(shù)概念,通過不斷調(diào)整滲透突變單元連續(xù)系數(shù),能較容易地解決常規(guī)數(shù)學(xué)模型中水量交換難確定的問題;在不改變整體滲流場流量的前提下,用不同水文地質(zhì)條件下的水頭間斷程度作為衡量指標,總結(jié)出了能恰當反映不同間斷程度的水頭間斷現(xiàn)象對應(yīng)的連續(xù)系數(shù)取值規(guī)律,為滲透突變有限元中單元間斷強弱參數(shù)的選取提供了具體確定的物理背景、方法和取值依據(jù).

[1]祝云華,劉新榮,梁寧慧,等.裂隙巖體滲流模型研究現(xiàn)狀與展望[J].工程地質(zhì)學(xué)報,2008,16(2):178-183.(ZHU Yun-hua,LIUXin-rong,LIANG Ning-hui,et al.Current research and prospects in modeling seepage field in fractured rock mass[J].Journal of Engineering Geology,2008,16(2):178-183.(in Chinese))

[2]周志芳.裂隙介質(zhì)水動力學(xué)原理[M].北京:高等教育出版社,2007.

[3]王媛,速寶玉.裂隙巖體滲流模型綜述[J].水科學(xué)進展,1996,7(3):276-282.(WANG Yuan,SU Bao-yu.Comment on themodels of seepage flow in fractured rock masses[J].Advances in Water Science,1996,7(3):276-282.(in Chinese))

[4]楊璐玲.巖體滲流場數(shù)值分析方法研究[D].天津:天津大學(xué),2003

[5]汝乃華.梅山連拱壩右壩座的錯動及其原因分析[J].水利學(xué)報,1988(9):28-35.(RU Nai-hua.The slip of the right abutment of Meishan multiple-arch dam and its analysis[J].Journal of Hydraulic Engineering,1988(9):28-35.(in Chinese))

[6]王媛,速寶玉.單裂隙面滲流特性及等效水力隙寬[J].水科學(xué)進展,2002,13(1):61-67.(WANG Yuan,SU Bao-yu.Research on the behavior of fluid flow in a singlefracture and its equivalent hydraulic aperture[J].Advances in Water Science,2002,13(1):61-67.(in Chinese))

[7]耿克勤,陳鳳翔,劉光廷,等.巖體裂隙滲流水力特性的實驗研究[J].清華大學(xué)學(xué)報:自然科學(xué)版,1996,36(1):102-106.(GENGKe-qin,CHEN Feng-xiang,LIUGuang-ting,et al.Experimental research of hydraulic properties of seepage flow in fracture[J].Journal of Tsinghua University:Science and Technology,1996,36(1):102-106.(in Chinese))

[8]魏進兵,鄧建輝,高春玉,等.三峽庫區(qū)泄灘滑坡非飽和滲流分析及滲透系數(shù)反演[J].巖土力學(xué),2008,29(8):2262-2266.(WEI Jin-bing,DENG Jian-hui,GAOChun-yu,et al.Unsaturated seepage analysis and back analysisof permeability coefficientfor Xietan landslidein Three Gorges Reservoir area[J].Rock and Soil Mechanics,2008,29(8):2262-2266.(in Chinese))

[9]王環(huán)玲,徐衛(wèi)亞,童富國.泄洪霧雨區(qū)裂隙巖質(zhì)邊坡飽和-非飽和滲流場與應(yīng)力場耦合分析[J].巖土力學(xué),2008,29(9):2397-2403.(WANG Huan-ling,XU Wei-ya,TONG Fu-guo.Coupled analysis of fracture rock mass slope saturated-unsaturated seepage field and stress field in flood discharge atomized rain area[J].Rock and Soil Mechanics,2008,29(9):2397-2403.(in Chinese))

[10]宋曉晨,徐衛(wèi)亞.裂隙巖體滲流模擬的三維離散裂隙網(wǎng)格數(shù)值模型(I):裂隙網(wǎng)格的隨機生成[J].巖石力學(xué)與工程學(xué)報,2004,23(12):2015-2020.(SONGXiao-chen,XU Wei-ya.Numerical model of three-dimensional discretefracturenetwork for seepage in fractured rocks(I):generation of fracture network[J].Chinese Journal of Rock Mechanics and Engineering,2004,23(12):2015-2020.(in Chinese))

[11]盧剛,周志芳.降雨入滲下互層狀裂隙巖體非飽和滲流分析[J].巖土工程學(xué)報,2008,30(9):1399-1403.(LU Gang,ZHOU Zhi-fang.Analysis of unsaturated seepage in inter-bed layered fractured rock mass with rainfall infiltration[J].Chinese Journal of Geotechnical Engineering,2008,30(9):1399-1403.(in Chinese))

[12]黃勇,陳靜.巖體裂隙介質(zhì)各向異性耦合模型研究[J].河海大學(xué)學(xué)報:自然科學(xué)版,2009,37(2):171-174.(HUANG Yong,CHENJing.Anisotropic coupled model for fractured rock medium[J].Journal of Hohai University:Natural Sciences,2009,37(2):171-174.(in Chinese))

[13]劉鳳茹.復(fù)合土工膜選型及缺陷滲漏量試驗研究[D].南京:河海大學(xué),2005.