王健平 蔡昌淵 劉 暉

(國網(wǎng)新源富春江水力發(fā)電廠，浙江 桐廬 311500)

·水利工程·

有限單元法庫區(qū)滲流場數(shù)值模擬及敏感性分析

王健平 蔡昌淵 劉 暉

(國網(wǎng)新源富春江水力發(fā)電廠，浙江 桐廬 311500)

采用FORTRAN程序，編寫了調(diào)用針對東部某抽水電站上庫區(qū)滲流場數(shù)值模擬的計(jì)算程序，研究了庫區(qū)滲流場的計(jì)算方法，并論證了該計(jì)算方法的正確性和實(shí)用性，分析結(jié)果證明該計(jì)算方法適應(yīng)性強(qiáng)、精度高、應(yīng)用方便、實(shí)用性好。

水文地質(zhì)，有限單元法，滲流場，數(shù)值模擬，敏感性分析

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

1 工程實(shí)例

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

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

2 滲流場數(shù)值模擬及結(jié)果分析

2.1 數(shù)學(xué)模型的建立

由于枯水期地下水變化不大，因此可以將地下水運(yùn)動概化為穩(wěn)定流。該計(jì)算區(qū)域內(nèi)的地下水運(yùn)動可以概化為在一個(gè)非均質(zhì)各向同性介質(zhì)中的二維流動，其數(shù)學(xué)模型如下：

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

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

由此得到一個(gè)NN階線性方程組，解方程組就可得出NN個(gè)節(jié)點(diǎn)的水頭。

2.2 幾何模型的建立及參數(shù)的選取

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

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

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

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

2.3 數(shù)值模擬結(jié)果及分析

根據(jù)上述數(shù)學(xué)模型，采用FORTRAN編寫的二維有限元程序進(jìn)行了模擬計(jì)算，庫區(qū)滲流場等水位線如圖4所示。

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

表1 計(jì)算結(jié)果與實(shí)測結(jié)果

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

3 模型參數(shù)敏感性分析

由于研究區(qū)地質(zhì)及水文地質(zhì)研究處于進(jìn)一步研究階段，許多參數(shù)還不確定，因此可以通過參數(shù)敏感性分析來研究不同參數(shù)時(shí)的地下水運(yùn)動狀態(tài)，分析參數(shù)對滲流場的敏感性。

第Ⅳ滲透性分區(qū)中，巖溶發(fā)育較強(qiáng)，其滲透系數(shù)增加1倍，其余參數(shù)不變。計(jì)算結(jié)果如圖6所示。

第Ⅰ滲透性分區(qū)中，巖溶發(fā)育較弱，其滲透系數(shù)減小1/2倍，其余參數(shù)不變。結(jié)果見圖7。

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

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

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

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

4 結(jié)語

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

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

[2] 代群力.地下水非線性流動模擬[J].水文地質(zhì)工程地質(zhì),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] 李康宏，柴軍瑞.壩基滲流分析兩種數(shù)值方法的比較[J].紅水河，2003，22(4)：14-17.

[6] 閆澍旺，王瑞鋼，王學(xué)軍，等.白石水庫滲流場有限元分析[J].水力發(fā)電學(xué)報(bào)，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