有限單元法庫區滲流場數值模擬及敏感性分析

王健平 蔡昌淵 劉 暉

(國網新源富春江水力發電廠，浙江 桐廬 311500)

·水利工程·

有限單元法庫區滲流場數值模擬及敏感性分析

王健平 蔡昌淵 劉 暉

(國網新源富春江水力發電廠，浙江 桐廬 311500)

采用FORTRAN程序，編寫了調用針對東部某抽水電站上庫區滲流場數值模擬的計算程序，研究了庫區滲流場的計算方法，并論證了該計算方法的正確性和實用性，分析結果證明該計算方法適應性強、精度高、應用方便、實用性好。

水文地質，有限單元法，滲流場，數值模擬，敏感性分析

庫區地下水滲流分析多以線性的達西定律為基礎，這是因為在大多數情況下，地下水滲流是滿足或近似滿足達西定律的，而且達西定律所表達的線性關系也使理論分析與數值分析較為簡潔，然而修建在巖溶地區的大壩基礎中可能存在強透水地質構造，滲流的雷諾系數很大，滲流速度與水力坡度不再遵從達西定律，此時必須考慮慣性力的影響。柴軍瑞(2001年)[1]、代群力(2000年)[2]認為有必要對巖溶區進行壩基地下水非線性滲流分析。由于非線性滲流的復雜性，至今仍沒有一個統一的公式來更好地描述這種流動，其中兩個常用的經驗公式模型為Forchheimer公式[3]，將其代入一般情況下三維空間滲流的連續性方程，即可得到非線性滲流公式。有限單元法是一種較為常用的數值方法。它是古典變分法與分塊多項式插值結合的產物，它的核心是對區域的離散化。1965年，津克維茨(Zienkiewiz)和張(Cheung)提出有限單元法適用于所有可按變分形式進行計算的場問題，為該方法在滲流分析中的應用提供了理論基礎。Leiws R W等(1987年)[4]將這一方法用于模擬多孔介質中的地下水。國內不少研究[5，6]認為有限單元法是一種先進有效的數值模擬法，用于滲流分析計算時可部分替代模型試驗，精度相對較高，可方便地模擬多種外部條件的特點。本文采用自己編寫的FORTRAN程序，調用了針對該具體工程條件編寫了有限元計算程序，并進行了參數反演分析，計算結果較為理想。

1 工程實例

東部某抽水蓄能電站裝機容量1 350 MW。電站樞紐由上水庫(見圖1)、下水庫、輸水系統、地下廠房等建筑物組成。上水庫面板堆石壩最大壩高183.5 m，總庫容1 703萬m3；下水庫面板堆石壩最大壩高33.4 m，總庫容1 676萬m3。工程區北鄰長江，南望太湖，屬寧鎮山脈的低山丘陵區，地勢總體是西高東低，地形切割不深，溝谷稍發育，以沖積堆積、剝蝕地貌為主，區內普遍可見數級剝夷面和3級～4級階地，以及相應的成層巖溶、深切河谷。

工程區從巖溶發育程度、地層巖性、地質構造及巖溶水賦存運移特征來劃為四類巖溶水文地質單元(區)，一類非可溶巖分布區；其中：A區——裸露型純碳酸鹽巖單斜巖溶發育強度弱～中等水文地質單元(區)，地表基巖大片裸露，巖性由白云巖類逐漸變為灰巖類；西側以F7斷層為界，北側以F4斷層為界，東側為侖山東面坡腳下水庫庫區F12-1斷層及其伴生的閃長玢巖脈和龍潭組(P2l)，南側以侖山南面坡腳奧陶系上統(O3)為界，位于工程區上水庫、輸水發電系統樞紐工程部位，為本文滲流場分析區；B區——覆蓋型純碳酸鹽巖與非可溶巖互層(或夾層)巖溶發育強度中等水文地質單元，分布在侖山坡腳，下水庫西側庫岸；C區——裸露型背斜巖溶發育強度弱～中等水文地質單元，位于上水庫的西側；非可溶巖由粉砂巖、粗面巖、閃長玢巖等組成，主要分布在下水庫的北側及侖山的南側坡腳侖山水庫一帶。

2 滲流場數值模擬及結果分析

2.1 數學模型的建立

由于枯水期地下水變化不大，因此可以將地下水運動概化為穩定流。該計算區域內的地下水運動可以概化為在一個非均質各向同性介質中的二維流動，其數學模型如下：

其中，Ω為計算區域；Γ1為第一類邊界；Γ2為第二類邊界；Hb為第一類邊界上的已知水頭；n為第二類邊界上的外法線方向；qb為第二類邊界上法向的單寬流量，流入為正，流出為負；W為降雨入滲補給量。將其簡寫為有限元方程式，即為：

[G]{H}={W}+{F}。

由此得到一個NN階線性方程組，解方程組就可得出NN個節點的水頭。

2.2 幾何模型的建立及參數的選取

根據地層巖性、巖溶水文地質單元，結合地形地貌特征，將研究區分為5個區(如圖2所示)。北部邊界沿著F3斷層，南部邊界至F9斷層，東部以姊妹橋F12斷層附近鉆孔25所在一線為邊界，西部以F7為邊界。

Ⅰ區：F4斷層與F3斷層之間的志留系墳頭組的砂巖、粉砂巖；Ⅱ區：F7斷層與F8斷層之間北至斷層F4，南至侖山山脊；Ⅲ區：由斷層F12與斷層F8與侖山山脊圍成的三角形區域；Ⅳ區：F7斷層與F8斷層之間，侖山山脊以南的部分；Ⅴ區：F8斷層以東，侖山山脊以南部分。將研究區域剖分為2 056個三角形單元，1 092個節點，見圖3。

為了充分利用觀測資料，將上庫區的鉆孔地下水水位作為上庫第一類邊界條件；下庫區泉水出露的點作為流量邊界。

參數選取采用鉆孔壓水試驗換算值和經驗值相結合的方法。根據鉆孔壓水試驗資料，滲透Ⅰ區滲透系數值取0.017 5 m/d；滲透Ⅱ區滲透系數取0.008 64 m/d；滲透Ⅲ區滲透系數值取0.019 5 m/d；滲透 Ⅳ區滲透系數值取0.021 m/d；滲透Ⅴ區滲透系數值取0.018 7 m/d。考慮庫盆中降雨入滲補給，根據侖山水庫降雨量資料，一月份平均降雨量為46.1 mm，即1.54 mm/d。在庫盆中地表水絕大部分滲入地下，因此降雨入滲補給系數取0.55。

2.3 數值模擬結果及分析

根據上述數學模型，采用FORTRAN編寫的二維有限元程序進行了模擬計算，庫區滲流場等水位線如圖4所示。

在上庫的主體部分，其西北、西及西南三側分水嶺均位于本亞區。從地下水等水位線圖中可以看出，地下水分水嶺與地表水分水嶺基本一致。地下水總體上由上庫向四周補給，由于南北兩側的非可溶巖(微透水層)的阻水作用，地下水在南側和北側以泉的形式排泄。侖山山脊一帶存在一個分水嶺，地下水在總體上向東流的同時，還向分水嶺的南、北兩側流動。在上庫南側，地下水向大哨泉以北的沖溝中排泄；東側地下水向正東排泄至姊妹橋方向；在北側，存在一近南北向的分水嶺，地下水向北流動的同時分別向東、西兩側排；在西側，地下水向鑰匙灣流動，這與地表水的流向基本一致。為了更好的表明計算的精度，分別在5個區各設立了2個地下水水位觀測點，共計10個地下水水位觀測點，現將計算值與實測值進行對比。計算結果與實測結果見表1及圖5。

表1 計算結果與實測結果

由表1及圖5不難看出，計算結果與實際測量值較為接近，其中最大誤差未超過5 m，分別出現在Ⅰ2和Ⅲ1處；而最小誤差為2 m，出現在Ⅴ2監測點，最大相對誤差為2.6%，可見模型具有一定的精度，同時證明該方法具有一定的正確性和實用性。

3 模型參數敏感性分析

由于研究區地質及水文地質研究處于進一步研究階段，許多參數還不確定，因此可以通過參數敏感性分析來研究不同參數時的地下水運動狀態，分析參數對滲流場的敏感性。

第Ⅳ滲透性分區中，巖溶發育較強，其滲透系數增加1倍，其余參數不變。計算結果如圖6所示。

第Ⅰ滲透性分區中，巖溶發育較弱，其滲透系數減小1/2倍，其余參數不變。結果見圖7。

第Ⅲ滲透性分區中，巖溶發育相對較弱，為了從保守的角度考慮，其滲透系數取值與第Ⅳ滲透性分區一致，其余參數不變。計算結果如圖8所示。

從圖6～圖8中可以看出，第Ⅳ滲透性分區滲透參數的變化對地下水滲流場變化有較大的影響，但總的趨勢與圖4一致。其他區域滲透參數變化時，地下水滲流場變化不明顯，即地下水滲流場對其他區域滲透參數變化不敏感。

通過對研究區滲流場進行計算分析可知，ZK1鉆孔以東侖山南北坡之間仍存在一分水嶺，地下水在向分水嶺南、北兩側流動的同時，總體上還向東流動。侖山南北坡由于非可溶巖的阻水作用，地下水在南側和北側以泉的形式排泄。本次計算中沒有單獨劃分這一高滲透性區域，因而計算結構沒有反映低水位槽的情況，滲流計算僅給出了研究區地下水運動的總體趨勢。

F8斷層以東的區域滲透性參數變化對地下水滲流場有較大的影響，其余區域滲透參數變化對地下水滲流場不敏感。

4 結語

1)根據實際工程應用的需要，采用FORTRAN程序，調用編寫的計算程序能夠真實模擬庫區滲流場情況，證明該方法具有一定的正確性和實用性。2)從地下水等水位線圖中可以看出，地下水分水嶺與地表水分水嶺基本一致。地下水總體上由上庫向四周補給，由于南北兩側非可溶巖的阻水作用，地下水在南側和北側以泉的形式排泄；東側地下水向正東排泄至姊妹橋方向；在西側，地下水向鑰匙灣流動。3)通過參數反演可知，第Ⅳ滲透性分區滲透參數的變化對地下水滲流場變化有較大的影響，其他區域滲透參數變化時對地下水滲流場影響不大。

[1] 柴軍瑞.壩基非達西滲流分析[J].水電能源科學,2001,19(4):1-3.

[2] 代群力.地下水非線性流動模擬[J].水文地質工程地質,2000(2):50-51.

[3] BEAR J.HYDRAULICS OF GROUND WATER[M].NEW YORK:MCGRAW-HILL,1979.

[4] LEIWS R W,SCHREFLER B A.THE FINITE ELEMENT METHOD IN THE DEFORMATION AND CONSOLIDATION OF POROUS MEDIA[M].NEW YORK:JOHN WILEY,1987.

[5] 李康宏，柴軍瑞.壩基滲流分析兩種數值方法的比較[J].紅水河，2003，22(4)：14-17.

[6] 閆澍旺，王瑞鋼，王學軍，等.白石水庫滲流場有限元分析[J].水力發電學報，2003(2)：53-61.

Numerical simulation of reservoir seepage field the finite element method and the parameter sensitivity analysis

WANG Jian-ping CAI Chang-yuan LIU Hui

(StateGridXinyuanFuchunjiangHydropowerPlant,Tonglu311500,China)

The paper adopts FORTRAN procedure, compiles the calculation program of the numerical simulation for some hydropower station reservoir in the east part of China, researches the calculation methods of the seepage field of the reservoir, indicates its accuracy and application, proves by the analysis result that the method is adoptable and practical with high accuracy and convenient application.

hydrogeology, finite element method, seepage field, numerical simulation, sensitivity analysis

1009-6825(2014)30-0232-03

2014-08-12

王健平(1964- )，男，高級工程師； 蔡昌淵(1978- )，男，工程師； 劉 暉(1987- )，男，助理工程師

TV139.1

A