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

牟 宏，趙 巍

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

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

牟 宏，趙 巍

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

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

防滲體系；滲流控制；巖體裂隙；非線性問題；分析

0 前 言

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

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

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

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

(1)

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

2)Laplace方程式，滲流自由面當作流量補給，同時為下降流速條件，多用于結構穩(wěn)定且不易壓縮的堤壩。

(2)

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

(3)

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

(4)

式中：sS為貯存系數。

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

(5)

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

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

2 三維裂隙網絡滲流自由面解析分析

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

2.1 潛水滲流自由面方程

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

(6)

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

2.2 數值計算方法

在數值模擬分析計算時，定義面單元為裂隙單元，結點為裂隙之間的交叉點，在巖體三維裂隙網絡滲流機理的基礎上代入邊界條件，即可得到各個結點上水頭大小和流量值。

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

1)網格劃分整體的滲流區(qū)域，滲流計算第一步采用初始滲透系數作為整個滲流區(qū)域的滲透系數，然后建立整體滲透矩陣，在三維裂隙網絡滲流原理基礎上計算區(qū)域內各結點水頭大小和流量值。

2)第一步求解出滲流場各結點水頭值，與水流位置高度z進行比較，然后劃分整個滲流區(qū)域為三個子區(qū)域R1、R2和R3。定義R1區(qū)域在滲流自由面上，R2區(qū)域在自由面下，R3區(qū)域為交錯復合單元區(qū)。

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

4)子域R3區(qū)域的復合單元，采用插值法逐個進行判斷，通過對比分析高斯積分點的H和z計算結果，若H>z，則此高斯點認定為子域R2內，若H

5)對滲透矩陣進行重新組合，計算得出新的滲流場各個結點的水頭大小和流量值。前提是以收斂準則為依據，判斷精度要求，若滿足，結束迭代過程；若精度不滿足，重復運行過程(3)-(5)，最終使相鄰兩次迭代計算后的滲流自由面滿足精度要求，此時可得出最終的滲流自由面位置。

3 工程實例

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

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

3.1 計算模型及網格劃分

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

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

圖2 計算模型網格剖分

3.2 滲流計算參數與邊界條件

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

考慮到現場試驗滲透系數和實驗室試驗得出的滲透系數不同，滲流計算模型的滲透系數采用各大高校經反復演練并最佳擬合的物理參數，見表1。

表1 最佳擬合模型滲透系數

3.3 計算結果及分析

計算得出各鉆孔計算水位和絕對誤差如表2所示，(相對誤差等于絕對誤差和總水頭的比值)。各截面浸潤線實測值繪制的曲線圖見圖3-5所示。

表2 各鉆孔計算水位及誤差

圖3 浸潤面計算實測對比圖(ZK2-1-ZK8-1)

圖4 浸潤面計算實測對比圖(ZK2-2-ZK8-2)

圖5 浸潤面計算實測對比圖(ZK2-3-ZK8-3)

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

5 結 論

文章基于連續(xù)介質非穩(wěn)定滲流的原理，在滲流區(qū)域非線性問題的滲透矩陣和自由面三維裂隙網絡滲流分析的基礎上，利用有限元軟件和自編程序對滲流區(qū)域非線性問題進行了研究，認為在應用單元滲透矩陣調整法時，對計算模型選取和區(qū)域邊界條件的界定非常關鍵，特別在于滲流自由面的變化區(qū)域，有限元分析中的網格單元劃分進行越密，則計算精度越高。通過計算結果和實測資料進行對比得出，認為采用該自編程序進行滲流區(qū)域非線性問題的求解是可行的。

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

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

[3]梁業(yè)國，熊文林，周創(chuàng)兵.有自由面滲流分析的子單元法[J].水利學報，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