巖體滲流的一種改進(jìn)數(shù)學(xué)模型

張文捷,程榮蘭,詹美禮,盛金昌,王國(guó)進(jìn)

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

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

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

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

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

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

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

1.2 三維滲透突變單元計(jì)算理論與滲透機(jī)理分析

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

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

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

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

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

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

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

改進(jìn)模型的等效節(jié)點(diǎn)流量表達(dá)式為

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

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

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

2 算 例

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

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

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

2.1 模型驗(yàn)證

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

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

圖4 滲流場(chǎng)等值線及侵潤(rùn)面計(jì)算結(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種模型的流量?jī)H相差5.4%.這表明,滲透突變單元沒(méi)有改變整體滲流量,滲透突變單元模型是正確的.

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

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

2.2.1 K強(qiáng)/K弱=10000時(shí)

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

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

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

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

2.2.2 K強(qiáng)/K弱=1000,500,100,50,10時(shí)

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

圖5 K強(qiáng)=10000K弱條件下的滲流場(chǎng)等值線Fig.5 Contours of seepage field when K強(qiáng)=10000K弱

圖6 K強(qiáng)=1000K弱條件下的滲流場(chǎng)等值線Fig.6 Contours of seepage field when K強(qiáng)=1000K弱

圖7 不同K強(qiáng)/K弱條件下的滲流場(chǎng)等值線Fig.7 Contours of seepage field with different cases of K強(qiáng)/K弱

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

取K弱的數(shù)量級(jí)為m,突變單元的滲透性量級(jí)為n;n=m/10,即n小于m 1個(gè)數(shù)量級(jí).通過(guò)對(duì)強(qiáng)、弱透水介質(zhì)滲透性不同差異倍比條件下多組合計(jì)算分析,得出水頭突變單元連續(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é) 語(yǔ)

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

[1]祝云華,劉新榮,梁寧慧,等.裂隙巖體滲流模型研究現(xiàn)狀與展望[J].工程地質(zhì)學(xué)報(bào),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ì)水動(dòng)力學(xué)原理[M].北京:高等教育出版社,2007.

[3]王媛,速寶玉.裂隙巖體滲流模型綜述[J].水科學(xué)進(jìn)展,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]楊璐玲.巖體滲流場(chǎng)數(shù)值分析方法研究[D].天津:天津大學(xué),2003

[5]汝乃華.梅山連拱壩右壩座的錯(cuò)動(dòng)及其原因分析[J].水利學(xué)報(bào),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é)進(jìn)展,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]耿克勤,陳鳳翔,劉光廷,等.巖體裂隙滲流水力特性的實(shí)驗(yàn)研究[J].清華大學(xué)學(xué)報(bào):自然科學(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]魏進(jìn)兵,鄧建輝,高春玉,等.三峽庫(kù)區(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)亞,童富國(guó).泄洪霧雨區(qū)裂隙巖質(zhì)邊坡飽和-非飽和滲流場(chǎng)與應(yīng)力場(chǎ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ī)生成[J].巖石力學(xué)與工程學(xué)報(bào),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é)報(bào),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é)報(bào):自然科學(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ù)合土工膜選型及缺陷滲漏量試驗(yàn)研究[D].南京:河海大學(xué),2005.