巖體裂隙滲流區(qū)域非線性問(wèn)題的數(shù)值求解分析

牟 宏，趙 巍

(大慶松嫩工程管理處，黑龍江 大慶 163311)

巖體裂隙滲流區(qū)域非線性問(wèn)題的數(shù)值求解分析

牟 宏，趙 巍

(大慶松嫩工程管理處，黑龍江 大慶 163311)

防滲體系的構(gòu)建一直是壩工設(shè)計(jì)中的重中之重，滲流分析是評(píng)價(jià)大壩安全穩(wěn)定的主要指標(biāo)之一，對(duì)滲流控制方案可行性和邊坡抗滑穩(wěn)定性的評(píng)價(jià)非常重要。由于滲流控制方程是線性的，通常無(wú)壓滲流需指出滲流自由面的位置，而部分滲流區(qū)域邊界的未知性和確定的必要性，使得連續(xù)介質(zhì)的滲流分析轉(zhuǎn)化為非線性問(wèn)題。文章在三維裂隙網(wǎng)絡(luò)滲流理論基礎(chǔ)上，編制了求解非線性問(wèn)題的相應(yīng)程序，結(jié)合工程實(shí)例，通過(guò)對(duì)比數(shù)值模擬計(jì)算結(jié)果和工程實(shí)測(cè)結(jié)果，論證自編有限元程序的可行性。同時(shí)研究結(jié)果可以給類(lèi)似工程設(shè)計(jì)提供一定的參考價(jià)值。

防滲體系；滲流控制；巖體裂隙；非線性問(wèn)題；分析

0 前 言

滲流分析是大壩防滲設(shè)計(jì)的重要部分，在壩體出現(xiàn)安全事故中此方面是造成滲透破壞的重要原因之一[1]，問(wèn)題在于對(duì)大壩的滲流分析考慮不全面。大壩滲流部位主要是壩體、邊坡、基巖位置，還有許多由地下水運(yùn)動(dòng)、各類(lèi)繞壩滲流等引起，這些滲流問(wèn)題均是滲流區(qū)域的非線性問(wèn)題，必須確定滲流自由面。自由面一般是可以隨著降雨量和水位等外界因素的變化而變化，自由面上部為大氣壓。滲流區(qū)域非線性問(wèn)題計(jì)算分析的重要環(huán)節(jié)是對(duì)滲流區(qū)域邊界的確定，換言之是對(duì)滲流自由面和溢出面位置的確定，過(guò)程中確保迭代計(jì)算的穩(wěn)定性。邊界的確定是滲流分析研究的重要組成，同時(shí)滲流計(jì)算是評(píng)價(jià)大壩滲流安全穩(wěn)定的主要指標(biāo)，因此，對(duì)滲流區(qū)域非線性問(wèn)題的正確求解以及分析研究具有重要的意義。

1 連續(xù)介質(zhì)滲流非線性的滲透矩陣解析法

連續(xù)介質(zhì)滲流區(qū)域求解中數(shù)學(xué)模型的確定主要是滲流方程的定解條件。連續(xù)介質(zhì)滲流的非線性問(wèn)題分為兩種情形：飽和滲流和飽和-非飽和滲流，所以數(shù)學(xué)模型基于此分為兩種。通常對(duì)于前者的二維非穩(wěn)定流進(jìn)行分析時(shí)，定解方程[2]以下4種：

1)Boussinesq方程式，原理是滲流自由面作為流量補(bǔ)給，看做水平面滲流的支撐方程，通常多用于在地下水運(yùn)動(dòng)的研究。

(1)

式中：H為地下水深；ε為為蒸發(fā)強(qiáng)度；ω為入滲強(qiáng)度。

2)Laplace方程式，滲流自由面當(dāng)作流量補(bǔ)給，同時(shí)為下降流速條件，多用于結(jié)構(gòu)穩(wěn)定且不易壓縮的堤壩。

(2)

3)擴(kuò)散方程式，滲流自由面當(dāng)作流量補(bǔ)給條件，基本假定是以杜布依為基礎(chǔ)，故該方法多用于自由面變化程度小或者滲流坡降變化不大的土壩。

(3)

4)固結(jié)方程式，滲流自由面看成流量補(bǔ)給條件，對(duì)土體壓縮性進(jìn)行了定義，故多用于黏土筑壩的多種固結(jié)情形。

(4)

式中：sS為貯存系數(shù)。

在飽和-非飽和滲流分析中，常常使用Neuman方程式，此方程推導(dǎo)過(guò)程中不考慮滲流自由面的邊界條件，做出特定假設(shè)后多用于非均質(zhì)各向異性壩。所以該方程式轉(zhuǎn)化式如下。

(5)

式中：C(ψ)為容水度；ψ為壓力水頭。

對(duì)于連續(xù)介質(zhì)的非穩(wěn)定滲流分析，滲流分布利用飽和滲流分析中的固結(jié)方程和自由面邊界定解條件進(jìn)行計(jì)算是符合實(shí)際的。

2 三維裂隙網(wǎng)絡(luò)滲流自由面解析分析

多孔連續(xù)介質(zhì)滲流分析主要是為了明確滲流的區(qū)域位置和邊界條件，與裂隙巖體滲流分析的不同在于，不考慮巖塊滲流的基礎(chǔ)上，水流沿著裂隙方向流動(dòng)，滲流也是基于此定向進(jìn)行的。此時(shí)巖塊主要承受水荷載，潛水面僅產(chǎn)生于裂隙處。因此，巖體裂隙滲流中滲流自由面是承壓面和非連續(xù)的潛水面組合而成的曲面[3]。

2.1 潛水滲流自由面方程

巖體裂隙滲流分析的基本假定是巖塊無(wú)滲透能力，水流只存在于巖體的裂隙網(wǎng)絡(luò)中，各裂隙中水位高度的連接線組成的曲面定義為滲流自由面。故對(duì)各裂隙中水位的確定成為確定自由面的關(guān)鍵。此問(wèn)題的研究原理類(lèi)似連續(xù)介質(zhì)滲流，確定水位這是一非線性問(wèn)題，因?yàn)樗唤M成了巖體裂隙滲流的邊界條件，同時(shí)滿足既定位置上的水頭高度和流量與外界滲出量的相同。假定裂隙巖體水流流態(tài)平穩(wěn)，潛水自由面水量穩(wěn)定，則滲流自由面各裂隙處端頭水頭差為零，即流量變化值為零。則潛水面邊界方程見(jiàn)表達(dá)式(1)：

(6)

式中：Γ3為自由面邊界；z為水流位置高度。

2.2 數(shù)值計(jì)算方法

在數(shù)值模擬分析計(jì)算時(shí)，定義面單元為裂隙單元，結(jié)點(diǎn)為裂隙之間的交叉點(diǎn)，在巖體三維裂隙網(wǎng)絡(luò)滲流機(jī)理的基礎(chǔ)上代入邊界條件，即可得到各個(gè)結(jié)點(diǎn)上水頭大小和流量值。

三維裂隙網(wǎng)絡(luò)滲流數(shù)值模擬的基本單元采用的是二維滲流分析中的單元，本質(zhì)上和連續(xù)介質(zhì)二維滲流分析是相同的。對(duì)三維巖體裂隙網(wǎng)絡(luò)滲流的非線性問(wèn)題計(jì)算分析的過(guò)程如下：

1)網(wǎng)格劃分整體的滲流區(qū)域，滲流計(jì)算第一步采用初始滲透系數(shù)作為整個(gè)滲流區(qū)域的滲透系數(shù)，然后建立整體滲透矩陣，在三維裂隙網(wǎng)絡(luò)滲流原理基礎(chǔ)上計(jì)算區(qū)域內(nèi)各結(jié)點(diǎn)水頭大小和流量值。

2)第一步求解出滲流場(chǎng)各結(jié)點(diǎn)水頭值，與水流位置高度z進(jìn)行比較，然后劃分整個(gè)滲流區(qū)域?yàn)槿齻€(gè)子區(qū)域R1、R2和R3。定義R1區(qū)域在滲流自由面上，R2區(qū)域在自由面下，R3區(qū)域?yàn)榻诲e(cuò)復(fù)合單元區(qū)。

3)給定子域R1的滲透單元滲透系數(shù)為K/1000；子域R2的滲透單元滲透系數(shù)與原始值保持一致。

4)子域R3區(qū)域的復(fù)合單元，采用插值法逐個(gè)進(jìn)行判斷，通過(guò)對(duì)比分析高斯積分點(diǎn)的H和z計(jì)算結(jié)果，若H>z，則此高斯點(diǎn)認(rèn)定為子域R2內(nèi)，若H

5)對(duì)滲透矩陣進(jìn)行重新組合，計(jì)算得出新的滲流場(chǎng)各個(gè)結(jié)點(diǎn)的水頭大小和流量值。前提是以收斂準(zhǔn)則為依據(jù)，判斷精度要求，若滿足，結(jié)束迭代過(guò)程；若精度不滿足，重復(fù)運(yùn)行過(guò)程(3)-(5)，最終使相鄰兩次迭代計(jì)算后的滲流自由面滿足精度要求，此時(shí)可得出最終的滲流自由面位置。

3 工程實(shí)例

某尾礦庫(kù)工程位于秦嶺中部腹地山區(qū)，常年地貌是中高山侵蝕地貌。海拔位于1245-1948m之間，礦區(qū)地形南低北高，傾向?yàn)槟衔髯呦?。礦區(qū)為一級(jí)支溝，溝長(zhǎng)5.5km，底寬3-8m， “∨”型走向的河谷。初期壩和后期壩平面圖見(jiàn)圖1。工程防洪標(biāo)準(zhǔn)：初期校核洪水位為50-100a一遇；中、后期設(shè)計(jì)為200a一遇。河谷最大流量為160m3/h。工程所在庫(kù)內(nèi)布設(shè)壩面位移變形觀測(cè)點(diǎn)3個(gè)，水位觀測(cè)標(biāo)尺5個(gè)，浸潤(rùn)線觀測(cè)孔6個(gè)，定期觀測(cè)并進(jìn)行記錄。

圖1 初期壩和后期壩平面圖

3.1 計(jì)算模型及網(wǎng)格劃分

為更加符合工程實(shí)際觀測(cè)，使得有限元計(jì)算結(jié)果和實(shí)際更加貼合，計(jì)算模型選取庫(kù)區(qū)ZK2斷面和ZK8斷面之間的區(qū)域進(jìn)行滲流分析。從上到下各斷面的壩面線、浸潤(rùn)線、地面開(kāi)挖線首尾相連形成三條曲線來(lái)進(jìn)行計(jì)算分析。

尾礦庫(kù)滲流數(shù)值模擬計(jì)算模型的建立的依據(jù)為ZK2-ZK8橫斷面圖，計(jì)算區(qū)域是壩面線之下的大壩，區(qū)域不包括基巖下的覆蓋層，計(jì)算模型的網(wǎng)格單元剖分圖見(jiàn)2，采用軟件ADINA進(jìn)行的網(wǎng)格劃分，模型中單元采用八結(jié)點(diǎn)六面體等參單元，總計(jì)剖分單元1628個(gè)，節(jié)點(diǎn)2174個(gè)。

圖2 計(jì)算模型網(wǎng)格剖分

3.2 滲流計(jì)算參數(shù)與邊界條件

滲流分析方法為有限單元法，ADINA建模后利用編制的SEEP程序進(jìn)行計(jì)算。定義第一類(lèi)水頭邊界為ZK2斷面和ZK8斷面的各自三個(gè)橫斷面，水頭大小是各斷面上的鉆孔地下水位，其他表面均為第二類(lèi)邊界，通過(guò)第二類(lèi)邊界橫斷面上的各滲流量等于零。

考慮到現(xiàn)場(chǎng)試驗(yàn)滲透系數(shù)和實(shí)驗(yàn)室試驗(yàn)得出的滲透系數(shù)不同，滲流計(jì)算模型的滲透系數(shù)采用各大高校經(jīng)反復(fù)演練并最佳擬合的物理參數(shù)，見(jiàn)表1。

表1 最佳擬合模型滲透系數(shù)

3.3 計(jì)算結(jié)果及分析

計(jì)算得出各鉆孔計(jì)算水位和絕對(duì)誤差如表2所示，(相對(duì)誤差等于絕對(duì)誤差和總水頭的比值)。各截面浸潤(rùn)線實(shí)測(cè)值繪制的曲線圖見(jiàn)圖3-5所示。

表2 各鉆孔計(jì)算水位及誤差

圖3 浸潤(rùn)面計(jì)算實(shí)測(cè)對(duì)比圖(ZK2-1-ZK8-1)

圖4 浸潤(rùn)面計(jì)算實(shí)測(cè)對(duì)比圖(ZK2-2-ZK8-2)

圖5 浸潤(rùn)面計(jì)算實(shí)測(cè)對(duì)比圖(ZK2-3-ZK8-3)

綜合各鉆孔計(jì)算水位值和誤差的對(duì)比分析，各橫斷面鉆孔水位的實(shí)測(cè)值中間部位的計(jì)算誤差明顯>兩邊鉆孔的水位值，工程實(shí)測(cè)結(jié)果中自由面上、下數(shù)值大、中間部位數(shù)值小，分布規(guī)律基本一致。實(shí)測(cè)值均<有限元計(jì)算出的各鉆孔的計(jì)算水位值，采集點(diǎn)離鉆孔ZK2和ZK8越近，誤差越小，反之亦然。表格分析中最大相對(duì)誤差5.68%，在允許范圍內(nèi)，表明該計(jì)算模型有限元分析計(jì)算結(jié)果中的浸潤(rùn)面位置雖比實(shí)測(cè)值偏高，但數(shù)值變化規(guī)律基本一致，直接驗(yàn)證了該計(jì)算程序的正確性和可行性。一定程度上說(shuō)明，文章分析采用的最佳擬合模型和參照工程的各物理參數(shù)進(jìn)行計(jì)算是正確的，符合工程實(shí)際。

5 結(jié) 論

文章基于連續(xù)介質(zhì)非穩(wěn)定滲流的原理，在滲流區(qū)域非線性問(wèn)題的滲透矩陣和自由面三維裂隙網(wǎng)絡(luò)滲流分析的基礎(chǔ)上，利用有限元軟件和自編程序?qū)B流區(qū)域非線性問(wèn)題進(jìn)行了研究，認(rèn)為在應(yīng)用單元滲透矩陣調(diào)整法時(shí)，對(duì)計(jì)算模型選取和區(qū)域邊界條件的界定非常關(guān)鍵，特別在于滲流自由面的變化區(qū)域，有限元分析中的網(wǎng)格單元?jiǎng)澐诌M(jìn)行越密，則計(jì)算精度越高。通過(guò)計(jì)算結(jié)果和實(shí)測(cè)資料進(jìn)行對(duì)比得出，認(rèn)為采用該自編程序進(jìn)行滲流區(qū)域非線性問(wèn)題的求解是可行的。

[1]柴軍瑞.巖土體水力學(xué)非線性問(wèn)題[J].巖土力學(xué)，2003(24)：159-162.

[2]王均星，吳雅峰，白呈富.有自由面滲流分析的流形單元法[J].水電能源科學(xué)，2003，21(04)：23-25.

[3]梁業(yè)國(guó)，熊文林，周創(chuàng)兵.有自由面滲流分析的子單元法[J].水利學(xué)報(bào)，1997(08)：34-38.

Numerical Solution and Analysis for Nonlinear Problems about Rock-body Fissure Seepage Region

MU Hong and ZHAO Wei

(Daqing Songnen Project Management Administration, Daqing 163311, China)

Composition of the seepage protection system has been a priority for the dam project design and analysis of seepage is one of major index for evaluating the dam safety and stability, which is very important to evaluate the feasibility of seepage control scheme and slope stability against slip. Usually, the location of seepage free surface needs to be pointed out for unconfined seepage because the seepage control equation is linear, however, the unknown boundary of partial seepage region and necessity transforms the seepage analysis for continuous medium into the nonlinear problem. Based on the seepage theory of three-dimensional fissure net, this paper carried out a corresponding program to solve the nonlinear questions, combined with project cases, by comparing the calculated results with the observed results, it proves that the finite element grogram is feasible, simultaneously, the researched achievements may supply references for similar projects.

system of seepage protection; seepage control; fissure of rock body; nonlinear problem; analysis

1007-7596(2017)06-0006-04

2017-05-24

牟宏(1991-)，女，黑龍江大慶人，助理工程師；趙巍(1980-)，女，黑龍江大慶人，助理工程師。

TV

B