基于有限元法計(jì)算三江平原地下水開發(fā)潛力研究

宮興龍，付 強(qiáng)，李 衡，郭 佳(東北農(nóng)業(yè)大學(xué)水利與建筑學(xué)院，哈爾濱 150030)

0 引 言

從20世紀(jì)90年代以后隨著墾區(qū)興起打井種稻,水稻的種植技術(shù)取得了長(zhǎng)足的發(fā)展，產(chǎn)量穩(wěn)步提高，同時(shí)由于打井大量開采地下水，降低了地下水位，澇漬堿危害減小。但人們擔(dān)心,抽取的地下水水溫低，影響著影響水稻米質(zhì)與產(chǎn)量，在不斷擴(kuò)大水井稻的同時(shí)需要不斷加大抽取地下水，導(dǎo)致地下水位不斷下降，使這些地區(qū)可能面臨超采，經(jīng)濟(jì)資源環(huán)境發(fā)展難可持續(xù)[1-3]。故本文僅重點(diǎn)論證地下水資源開發(fā)能否支撐本農(nóng)場(chǎng)糧食生產(chǎn)、經(jīng)濟(jì)發(fā)展和生態(tài)環(huán)境可持續(xù)問題,即地下水資源開發(fā)潛力分析。很多人進(jìn)行了三江平原地下水開發(fā)潛力研究，但大部分都是基于水量平衡法[4-9]，基于此本文嘗試通過分析三江平原灌區(qū)典型區(qū)域853農(nóng)場(chǎng)的水文地質(zhì)條件，采用Galerkin法推導(dǎo)出該地區(qū)承壓水的有限元方程組，并采用了改進(jìn)平方根法求解該方程組，進(jìn)而計(jì)算出853農(nóng)場(chǎng)的地下水開發(fā)潛力。

1 地下水流數(shù)學(xué)模擬

承壓水地下水運(yùn)動(dòng)采用的數(shù)學(xué)模擬模型如下[8]:

(1)

H(x,y,0)=H0(x,y)t=0,(x,y)∈Q

(2)

H(x,y,t)|Γ1=φ1(x,y)t=0,(x,y)∈Q

(3)

(4)

式中:T為導(dǎo)水系數(shù)；S為貯水系數(shù)；H0為水頭初始值；φ1為Γ1(第一類邊界)的已知水位函數(shù)；q為Γ2(第二類邊界)上的側(cè)向補(bǔ)給量；n為邊界Γ1的外法線方向；Q為方程建立區(qū)域。

1.1 基于Galerkin法構(gòu)建地下水運(yùn)動(dòng)有限元方程

對(duì)式(1)～式(4)采用Galerkin法得到有限元方程為：

(5)

式中：φi為基函數(shù)；n為除了第一類邊界上的節(jié)點(diǎn)外的節(jié)點(diǎn)數(shù)。

應(yīng)用分布積分法(5)式得：

(6)

式(6)左端第一項(xiàng)應(yīng)用Green公式:

(7)

利用式(7)對(duì)式(6)整理：

(8)

以三角元離散式(8)：

(9)

在單元上式(9)變形為：

(10)

式中:e表示三角單元上的變量。

1.2 有限元方程的構(gòu)建

為了構(gòu)建有限元方程對(duì)(10)變形為：

(11)

式中:Δe為三角元面積；qai,qik為邊界段的流量；Lai,Li為邊界段長(zhǎng)度；a,b,c采用下式計(jì)算：

(12)

式中:(xi,yi)、(xj,yj)、(xk,yk)為三角元三個(gè)節(jié)點(diǎn)i,j,k的坐標(biāo)。

1.3 方程組求解

式(12)利用矩陣形式為：

(13)

式中：導(dǎo)水系數(shù)矩陣D、貯水系數(shù)矩陣P、常數(shù)項(xiàng)矢量矩陣F。上述方程是具有對(duì)稱，正定系數(shù)矩陣的線性方程，因此采用改進(jìn)平方根法求解[9,10]。

1.4 抽水井有限元方程

“853”農(nóng)場(chǎng)的取得地下水主要方式是利用抽水井進(jìn)行抽水，如果抽水井在三角單元節(jié)點(diǎn)上應(yīng)在三角單元方程(13)式上加一項(xiàng)：

(14)

式中：rwi為水井的半徑。

如抽水井在單元的其他段，單位時(shí)間，單位面積上的垂直流出(入)含水層的量為ω=-Q/Fw(Fw單元面積)，單元e上的節(jié)點(diǎn)i,j,k依次有：

(15)

2 計(jì)算實(shí)例

2.1 研究區(qū)地質(zhì)概況

853農(nóng)場(chǎng)位于三江平原東部，寶清縣境內(nèi)，屬于溫暖半濕潤(rùn)農(nóng)業(yè)氣候區(qū)，多年平均降水量557.2 mm。在地理位置上，東南部緊靠完達(dá)山，西北及北部以平原河流撓力河為界，西南的以生產(chǎn)蛤蟆的蛤蟆通河為邊界。本區(qū)有基巖裂隙孔隙含水層、第四系孔隙承壓含水層、第三系裂隙承壓含水層,第四系孔隙含水層系統(tǒng)在區(qū)內(nèi)最重要，含水層較厚供水量充沛,從山前向河谷含水層逐漸增厚，地下水埋深逐漸變淺。第三系裂隙孔隙水系統(tǒng)內(nèi)流動(dòng)較弱，承壓水的賦存、運(yùn)輸開采條件差于第四系，富水性差于第四系。該區(qū)農(nóng)業(yè)灌溉用水主要來自于第四系含水層，該層地下水易于補(bǔ)給和排泄，因該區(qū)內(nèi)地下水的開采還可能深入到第三系含水層,但開采量比較少，基于此本次模擬以第四系含水層為主。該區(qū)的主要補(bǔ)給方式是河流滲漏、降水入滲和側(cè)向地下徑流；人工抽水井開采是地下水排泄的主要方式,在開采量小的時(shí)間里還有一部分河流排泄和側(cè)向徑流排泄。

2.2 研究區(qū)域離散

為建立研究區(qū)的地下水流運(yùn)動(dòng)介質(zhì)中的水文地質(zhì)模型，基本步驟是要對(duì)實(shí)際水文地質(zhì)條件加以概化，建立水文地質(zhì)概念模型。根據(jù)上述水文地質(zhì)特殊條件,本文主要致力于研究第四系含水層，山區(qū)第四系覆蓋較薄，僅有少量孔隙裂隙潛水可供利用，其貯存于完達(dá)山麓，本次模擬不計(jì)。第三系裂隙孔隙承壓水儲(chǔ)存于巨厚的第四系含水層之下，在目前的開采條件和經(jīng)濟(jì)條件下尚難以對(duì)其進(jìn)行開采，且收集資料困難，所以本次模擬也不考慮。由于第四系含水層其覆蓋有10～35 m厚的亞黏土層,最厚達(dá)42 m，具有良好的隔水能力，地下水具有承壓性質(zhì)。故將含水層在垂向上概化為1層，即第四系孔隙承壓含水層，其上腹亞黏土，下部為不透水基巖。水平三條河流概化為流量邊界。采用有限元法數(shù)值求解式(7)對(duì)上述數(shù)學(xué)模型進(jìn)行求解，對(duì)比例尺為1：10萬的八五三農(nóng)場(chǎng)綜合水文地質(zhì)圖及地形圖進(jìn)行掃描，然后用數(shù)字化工具進(jìn)行矢量化處理為柵格文件,將柵格文件導(dǎo)入River Tool工具，提取出撓力河，蛤蟆通河流，七里沁河和圖1右下方的邊界。模型作為計(jì)算模擬區(qū)的剖分底圖計(jì)算區(qū)面積為1 182 km2，采用1 890個(gè)三角元進(jìn)行離散見圖1。圖1上邊的左邊邊線為撓力河，右邊的線為七里沁河，下邊的線為不透水邊界，左下邊的線為蛤蟆通河。

圖1 853農(nóng)場(chǎng)三角元離散圖Fig.1 Discrete map of Triangular element in 853 farm

2.3 水位觀測(cè)點(diǎn)

853農(nóng)場(chǎng)有長(zhǎng)期地下水觀測(cè)井，本文選取7個(gè)長(zhǎng)期觀測(cè)井，分別為853農(nóng)場(chǎng)1分場(chǎng)6小隊(duì)的觀井編號(hào)32號(hào)，2分場(chǎng)1小隊(duì)的觀井編號(hào)6號(hào)3分場(chǎng)7小隊(duì)的觀井編號(hào)10號(hào)，4農(nóng)場(chǎng)1小隊(duì)的觀井編號(hào)25號(hào)，5農(nóng)場(chǎng)2小隊(duì)的觀井編號(hào)20號(hào)，6農(nóng)場(chǎng)5小隊(duì)的觀井編號(hào)14號(hào)，6農(nóng)場(chǎng)9小隊(duì)的觀井編號(hào)13號(hào)。

2.4 降水量

降雨通過分析1995-2001年的月降水量數(shù)據(jù)，分析得到最大降水年份是1959年降水量為729.6m，多年平均降水量523.35。最小降水年份2004年458 mm。其過程線見圖2?？梢姳镜貐^(qū)降水量主要集中在6-8月份。本地區(qū)的降水入滲系數(shù)空間分布難于確定，因此本文根據(jù)清河灌區(qū)的設(shè)計(jì)值采用0.06。

圖2 月降水量數(shù)據(jù)圖Fig.2 Monthly rainfall data map

2.5 地下水開采控制水位判斷

文獻(xiàn)[4]還指出，根據(jù)水文地質(zhì)條件，開采技術(shù)條件和社會(huì)生產(chǎn)需要，利用地下水變化直接反映開采程度的特點(diǎn)，提出合理的“地下水開采控制水位”概念，防止超量開采地下水導(dǎo)致地下水枯竭。筆者暫提出“理想水位(埋深)”概念，853農(nóng)場(chǎng)承壓開采為主區(qū)的理想水位本文控制在10～15 m。

2.6 目前開采條件下地下水開發(fā)程度

采地下水一直存在爭(zhēng)議，本文以“853”農(nóng)場(chǎng)的1998年到2012的觀測(cè)井?dāng)?shù)據(jù)進(jìn)行分析。以“853”農(nóng)場(chǎng)的5農(nóng)場(chǎng)2小隊(duì)的觀井編號(hào)20，因?yàn)槠湓谇搴庸鄥^(qū)內(nèi)，存在開采地下水進(jìn)行水稻田灌溉的量大于其他各觀測(cè)井所在的地區(qū)。從圖可以看出地下水下降最大的一年是2012年，地下水為45.712 7 m。20號(hào)觀測(cè)井的地面高程為58.912 7，其差值為13.5，由于其為承壓區(qū)不超過15 m，因?yàn)楸疚徽J(rèn)為853農(nóng)場(chǎng)目前開采地下水不存在超采，而且還有一定開采的空間。由八五三農(nóng)場(chǎng)五分場(chǎng)2隊(duì)地下水位高程曲線圖(圖3)可以看出其地下水恢復(fù)地比較快，可見地下水很豐富,故有很大開發(fā)潛力。

圖3 5分場(chǎng)2隊(duì)地下水位高程曲線圖Fig.3 The underground water level elevation curve in fifth parvial field second team

2.7 基于數(shù)值法推求地下水開發(fā)潛力研究

“853”農(nóng)場(chǎng)從近年地下水動(dòng)態(tài)分析，每年4月恢復(fù)水位埋深大約3.0 m，尚小于理想埋深最高值(10 m)?！?53”農(nóng)場(chǎng)屬于以承壓水為主的地下水開采由。撓力河河道寬廣，一年時(shí)間里水深變化不大，在56.5 m波動(dòng)，本文據(jù)此給出撓力河水位為56.5 m。蛤蟆通河水位較淺，本文采用的水位為河道高程，同理對(duì)于七里沁河也采用同樣處理。對(duì)7個(gè)水位觀測(cè)點(diǎn)的多年1月份的地下水水位分析得出其多年平均初始水位，采用距離倒數(shù)法插值得到計(jì)算所用的初始水位(見圖4)。第四系承壓含水層由于資料較少，其參數(shù)的選取主要是根據(jù)含水層水文地質(zhì)特征,及前人的抽水試驗(yàn),參照經(jīng)驗(yàn)值給出初始值彈性釋水系數(shù)0.003 5，導(dǎo)水系數(shù)1 560 m2/d?？刂扑徊捎?0～15 m。為了計(jì)算地下水開發(fā)潛力本例按照承壓水理想埋深10～15 m計(jì)算地下水開發(fā)潛力，地下水為控制在水深10 m是計(jì)算的開發(fā)潛力為8 328.09萬m3/a、地下水水位控制在15 m時(shí)，可得其開發(fā)潛力為9 456.09萬m3/a。

圖4 853農(nóng)場(chǎng)多年春季平均水位圖Fig.4 The average water level diagram of spring for many years in 853 farm

2.8 激發(fā)出江河對(duì)地下水補(bǔ)給

“853”農(nóng)場(chǎng)周圍三面鄰河，北及北部以撓力河為界，西南的蛤蟆通河，西北的撓力河，東北的越嶺河，因此在控制開采地下水水位降到一定程度，就發(fā)生江河加大補(bǔ)給地下水[6]，基于此本位在控制承壓區(qū)10～15 m，計(jì)算其激發(fā)的河流補(bǔ)給量。由于蛤蟆通河和七里沁河水位相對(duì)于承壓水比較高，激發(fā)補(bǔ)給量不明顯。因此本文的激發(fā)江河補(bǔ)給主要考慮撓力河，由計(jì)算可得，在控制降深水位取15 m，可得其激發(fā)的補(bǔ)給量為256.09 萬m3/a。

3 結(jié) 語

綜上所述，由“853”農(nóng)場(chǎng)從1998-2013的觀測(cè)井?dāng)?shù)據(jù)及地下水水位控制數(shù)據(jù)分析得出，目前“853”農(nóng)場(chǎng)地下水開采不僅沒有超采，而且還有很大開發(fā)潛力，地下水位下降只是7、8月份，之后隨著地下水開采量減少，地下水水位迅速的恢復(fù)。由基于有限元法構(gòu)建的地下水運(yùn)動(dòng)數(shù)值模型計(jì)算得到“853”農(nóng)場(chǎng)可以開發(fā)9 456.09 萬m3/a水量，“853”農(nóng)場(chǎng)的北部靠近撓力河和西部的蛤蟆通的地區(qū)地下水開發(fā)量明顯大于東部的靠近七里沁河的開采量。在合理控制地下水情況下，可以激發(fā)出江河更多的補(bǔ)給，撓力河中部補(bǔ)給量為256.09 萬m3/a。

□

[1] 王立權(quán),馮建維.水稻熱的隱憂[J].水利天地,2004,(6):18-19.

[2] 張惠斌,于 東,姚章村.論“打井種稻”與“循環(huán)經(jīng)濟(jì)”[J].水利科技與經(jīng)濟(jì),2006,12(12):820-821.

[3] 王立權(quán),姚章村.高寒地區(qū)井灌水田井水增溫技術(shù)的探索與實(shí)踐[J].地下水,2006,28(3):46-49.

[4] 何 璉.中國三江平原[M].哈爾濱:黑龍江科學(xué)技術(shù)出版社,2000.

[5] 姚章村,安瑞強(qiáng),董經(jīng)財(cái),等.三江平原建三江地下水激發(fā)補(bǔ)給初步分析[J].黑龍江水專學(xué),2001,28(3):20-23.

[6] 李宏勛,陳 杰,李 偉,等.291農(nóng)場(chǎng)地下水資源開發(fā)潛力初步分析[J].水利科技與經(jīng)濟(jì),2008,14(1):49-51.

[7] 王貴玲,劉志明,劉花臺(tái),等.地下水潛力評(píng)價(jià)方法[J].水文地質(zhì)工程地質(zhì),2003,(1):63-66.

[8] 薛禹群,謝春紅.地下水?dāng)?shù)值模擬[M]. 北京：科學(xué)出版社,2007.

[9] 張德良.計(jì)算流體力學(xué)教程[M]. 北京：高等教育出版社,2009.

[10] 吳頌平,劉趙淼 譯.計(jì)算流體力學(xué)基礎(chǔ)及其應(yīng)用[M]. 北京：高等教育出版社,2013.