臺(tái)蘭河地下水庫(kù)輻射井抽水過(guò)程的非穩(wěn)定滲流場(chǎng)的有限元分析

劉昌軍，趙 華，張順福，丁留謙

中國(guó)水利水電科學(xué)研究院，北京 100038

0 引言

大口徑輻射井由于具有較好的排水效果，已在眾多水利工程中得到廣泛應(yīng)用。但輻射井因輻射管長(zhǎng)度、層數(shù)和數(shù)量的不同，其排水效果差異很大。如何進(jìn)行輻射井的合理設(shè)計(jì)，包括輻射井間距、輻射管的幾何尺寸（長(zhǎng)度、管徑和層數(shù)）等是應(yīng)用輻射管進(jìn)行排水設(shè)計(jì)的關(guān)鍵［1－2］。有限元分析是解決該問(wèn)題的重要手段，但輻射井的徑向尺寸很小而橫向尺寸很大、數(shù)量眾多，且空間分布復(fù)雜，給有限元模擬其精細(xì)滲流場(chǎng)帶來(lái)一定的難度。目前對(duì)堤壩滲控分析中密集排水孔的模擬方法較多［3－6］。但目前還沒(méi)有相關(guān)文獻(xiàn)對(duì)輻射井的排水效果和滲流場(chǎng)特性進(jìn)行數(shù)值模擬方面的研究。

為了解輻射井的抽排水效果，在充分研究了輻射井工作原理和滲流行為的基礎(chǔ)上，提出了輻射井子結(jié)構(gòu)法和子結(jié)構(gòu)滲流開關(guān)器等關(guān)鍵技術(shù)進(jìn)行輻射井滲流場(chǎng)的有限元分析方法，并聯(lián)合應(yīng)用迭代增量法和求解大型稀疏矩陣的預(yù)處理共軛梯度算法進(jìn)行非穩(wěn)定滲流場(chǎng)的分析。該方法既解決了輻射井附近滲流場(chǎng)的精細(xì)模擬，又解決了非穩(wěn)定滲流場(chǎng)自由面的迭代問(wèn)題，可用于多種排水措施和輻射井不同設(shè)計(jì)方案的計(jì)算比較。在此基礎(chǔ)上，采用IDL語(yǔ)言開發(fā)了三維可視化滲流有限元分析軟件 GWSS［7－9］。將該軟件應(yīng)用于臺(tái)蘭河地下水庫(kù)輻射井抽水過(guò)程和水位回復(fù)過(guò)程的的非穩(wěn)定滲流場(chǎng)的計(jì)算分析中，并利用各觀測(cè)孔的現(xiàn)場(chǎng)實(shí)測(cè)值對(duì)作者提出的輻射井子結(jié)構(gòu)法的計(jì)算理論和編制軟件的可靠性進(jìn)行驗(yàn)證。

1 輻射井的工作原理及其滲流行為

1.1 輻射井結(jié)構(gòu)和組成

輻射井排滲系統(tǒng)一般由豎井、輻射管等組成，見圖1。實(shí)際運(yùn)行中輻射井常常和虹吸井一起使用，即虹吸輻射井系統(tǒng)。虹吸輻射井系統(tǒng)除豎井和輻射管外還包括虹吸管和水封槽等。本文重點(diǎn)研究輻射井工作原理和方法，有關(guān)虹吸輻射井的研究成果將另文發(fā)表。輻射管是關(guān)鍵的降水設(shè)施，其成井質(zhì)量要求較高，另外輻射管的布置形式靈活，可以水平布置、也可以斜向布置，可根據(jù)地形地勢(shì)進(jìn)行選擇。

圖1 輻射井排滲系統(tǒng)組成Fig.1 Drainage system of radiation wells

1.2 輻射井的滲流行為

輻射井系統(tǒng)的工作模式有2種：一種是抽水，另一種是自流排水。對(duì)輻射井的滲流行為在算法上可分2種情況：一是自由面穿過(guò)輻射管；二是自由面不與輻射管相交。有自由面穿過(guò)的輻射管內(nèi)邊界面上的滲流行為如圖2所示。

1）如果自由面水位低于輻射管底部aa′，則該輻射管不起到排水作用，輻射管失效。因此可以在每個(gè)輻射管底部高程處虛構(gòu)一個(gè)數(shù)學(xué)開關(guān)器，詳細(xì)見文獻(xiàn)［6－7］。進(jìn)行輻射管計(jì)算時(shí)，先假設(shè)開關(guān)器打開，管內(nèi)邊界全部作為可能逸出邊界，每一步迭代后對(duì)孔內(nèi)結(jié)點(diǎn)壓力進(jìn)行甄別，記孔內(nèi)邊界各結(jié)點(diǎn)壓力為p。

圖2 輻射管濾管內(nèi)的滲流行為Fig.2 Drainage principle of radiation wells

若輻射管底面壓力p＜0，認(rèn)為輻射管全部位于非飽和區(qū)，此時(shí)輻射管完全失效。輻射管全部在非飽和區(qū)時(shí)，其邊界條件數(shù)學(xué)表達(dá)式為

式中：kij為滲透材料的滲透張量矩陣；h為水頭；xj為3個(gè)方向的坐標(biāo)；ni為內(nèi)邊界面外法線方向的方向余弦，i為3個(gè)坐標(biāo)主方向；x3為z方向坐標(biāo)。

若輻射管bb′面上的結(jié)點(diǎn)壓力p＞0，認(rèn)為輻射管全部位于飽和區(qū)，此時(shí)邊界為定水頭邊界，其數(shù)學(xué)表達(dá)式為

式中，h0為輻射井豎井內(nèi)的水位。

2）如果自由面水位位于aa′或bb′內(nèi)，則輻射管內(nèi)邊界為可能逸出邊界，其數(shù)學(xué)表達(dá)式為

式中，z0為aa′或bb′逸出面上結(jié)點(diǎn)位置高程。

2 非穩(wěn)定滲流計(jì)算理論

2.1 數(shù)學(xué)模型

非均質(zhì)各向異性有自由面穩(wěn)定滲流的基本微分方程如下：

式中：Q為源或匯項(xiàng)（1/s）；SS為貯水率（1/m），即單位體積的飽和土體內(nèi)，當(dāng)下降1個(gè)單位水頭時(shí)，由于土體壓縮和水的膨脹所釋放出來(lái)的貯存水量。

式（1）的定解邊界條件如下：

自由面邊界條件為

式中：h（x，y，z，t）為已知水頭函數(shù)；q為自由面上因自由變動(dòng)而引起的流量補(bǔ)給；μ為飽和差（自由面上升時(shí)）或給水度（自由面下降時(shí)），它表示自由面改變單位高度下，含水層單位截面積上吸收（自由面下降）或排出（自由面上升時(shí)）的水量，是無(wú)量綱數(shù)；θ為自由面外法線方向與垂線的交角；t為時(shí)間。

2.2 有限元法求解

對(duì)2.1節(jié)提出的數(shù)學(xué)模型采用Galerkin法和等參8結(jié)點(diǎn)六面體單元在空間域上進(jìn)行合理離散，將單元支配方程進(jìn)行集成，可得整體有限元支配方程：

在此基礎(chǔ)上，對(duì)時(shí)間域進(jìn)行合理離散，可得到有限元時(shí)間隱式差分迭代格式：

式中：K為滲透系數(shù)矩陣；S為總質(zhì)量矩陣；G為自由面貢獻(xiàn)矩陣；p為結(jié)點(diǎn)壓力列陣；F為等效結(jié)點(diǎn)流量列陣；Δt為時(shí)間步長(zhǎng)。

2.3 自由面處理的增量迭代法

由于自由面位置是未知的，式（3）無(wú)法一次求解，這里采用文獻(xiàn)［6－7］介紹的增量迭代法，可推導(dǎo)得如下適于計(jì)算的迭代格式：

式中：A為非線性方程的系數(shù)矩陣，

2.4 預(yù)處理共軛梯度算法

采用預(yù)處理共軛梯度算法（preconditioned conjugate gradients method，PCG）對(duì)線性方程組式（4）進(jìn)行聯(lián)立求解。

PCG法假設(shè)矩陣A被分為2個(gè)矩陣，即：

式中：M為預(yù)優(yōu)矩陣；N為剩余矩陣。M為A的預(yù)條件形式，稀疏性與A相同，總是對(duì)稱且正定的。有關(guān)預(yù)處理共軛梯度算法的具體求解方法見文獻(xiàn)［10］。

2.5 輻射井子結(jié)構(gòu)法

按照上述介紹的輻射井的滲流行為，結(jié)合有限元法求解原理，將子結(jié)構(gòu)作為主網(wǎng)格整體結(jié)構(gòu)的一部分一起來(lái)考慮［6－7，11］，則輻射井的非穩(wěn)定滲流子結(jié)構(gòu)總傳導(dǎo)矩陣及相應(yīng)流量列陣寫成分塊形式，如式（6）所示：

對(duì)于子結(jié)構(gòu)內(nèi)部邊界的處理與外部邊界相同即可，詳見參考文獻(xiàn)［11］。

3 工程應(yīng)用

3.1 臺(tái)蘭河地下水庫(kù)工程概況

臺(tái)蘭河地下水庫(kù)工程位于新疆維吾爾自治區(qū)阿克蘇地區(qū)溫宿縣境內(nèi)的臺(tái)蘭河灌區(qū)，距溫宿縣城33 km，距阿克蘇市46km，距烏魯木齊市955km，研究區(qū)位于314國(guó)道穿過(guò)2km范圍內(nèi)。臺(tái)蘭河地下水庫(kù)工程主要由滲庫(kù)工程、取水工程和輸水工程三大部分組成。滲庫(kù)工程是通過(guò)在上游修建簡(jiǎn)易人工分洪工程（蓄洪庫(kù)、擋水堤壩或坑、塘等），加大地表水地下水的回灌量，形成可灌的地下水庫(kù)；另外在下游選擇合適的地形地質(zhì)條件修建地下暗壩或天然暗壩，使上游地下松散地層內(nèi)一定空間的地下水位壅高，形成地下水庫(kù)。取水工程是在地下水庫(kù)取水點(diǎn)修建大口徑取水井截取第四紀(jì)松散地層中的地下水，采用低投入、高產(chǎn)出的方法提?。ɑ蜃粤鳎┑叵滤?，以起到盡量減少蒸發(fā)、以豐補(bǔ)歉、下游防洪等調(diào)節(jié)水庫(kù)的作用。

臺(tái)蘭河地下水庫(kù)工程是將臺(tái)蘭河洪水期的洪流和冬閑期河水轉(zhuǎn)換成地下水，再使地下水轉(zhuǎn)換為地表水資源，以達(dá)到增加臺(tái)蘭河流域灌溉水量的目的；利用地表水和地下水聯(lián)合調(diào)度臺(tái)蘭河流域需水過(guò)程，創(chuàng)新和優(yōu)化配置流域內(nèi)的水利工程布局。

目前取水井工程的一期已完成4眼大口徑輻射井及輸水管道的建設(shè)，可進(jìn)行自流式供水，流量達(dá)到0.3～0.5m3/s。臺(tái)蘭河地下水庫(kù)取水工程的二期完成“自流虹吸”輸水系統(tǒng)的建設(shè)。在一期建設(shè)的基礎(chǔ)上完成輸水廊道、集水廊道以及在輸水廊道距F1號(hào)井800m的下游布置F0號(hào)輻射井的施工建設(shè)。其中輸水廊道和集水廊道的最大降深均為10m，F(xiàn)0號(hào)輻射井和F1號(hào)輻射井的虹吸水位降深達(dá)到12 m，設(shè)計(jì)流量1.0～1.2m3/s；并且利用控制閘調(diào)控F0號(hào)井水位，使其在高水位狀態(tài)下將地下水引入高位水池，進(jìn)入自壓式滴灌系統(tǒng)［12］，示范區(qū)面積1萬(wàn)ha。整個(gè)取水工程的平面布置見圖3。

圖3 臺(tái)蘭河地下水庫(kù)工程平面布置圖Fig.3 Layout diagram of groundwater reservoir of Tailan River

3.2 研究區(qū)水文地質(zhì)條件

臺(tái)蘭河由北向南穿過(guò)洼地進(jìn)入古木別致低山丘陵地帶，海拔高程1 498～1 532m；又經(jīng)老龍口豁口下泄進(jìn)入洪積平原，在下游依次形成沖積礫質(zhì)傾斜平原，地形坡降10‰～16‰，東西向?qū)挾?0km；再向南為洪積細(xì)土平原，地面高程1 185～1 220m，地面坡降5‰～8‰。臺(tái)蘭河地下水庫(kù)示范區(qū)選擇在314國(guó)道以北山前沖洪積扇儲(chǔ)水構(gòu)造的小阿吉溝谷中，地下水位埋藏深度1.5～3.0m，含水層為砂卵礫石層。

3.3 取水工程系統(tǒng)

取水工程系統(tǒng)為“橫坎兒井”式的大口徑豎井和水平輻射井的集水系統(tǒng)［12］，垂直地下水流向布設(shè)3眼大口徑輻射井。由于地形坡度較大，在大口徑井下游約0.8km處布置1眼大口徑輻射集水井，采用“自流虹吸”原理，用鋼管將上述4眼井聯(lián)成一體，采用閘門進(jìn)行控制。垂直地下水流向的大口徑井井距為250m左右，井徑3.5m，井深30m，水平輻射井每井3層，每層8～10個(gè)輻射井。輻射井采用裸孔鋼管，管徑159mm×10mm，單根長(zhǎng)10～15m，最長(zhǎng)達(dá)24m。每個(gè)輻射井最大降深為10m，單井出水量大于1 200m3/h。在集水井下游修建引水管道進(jìn)行自流式引水。引水管道直徑1.2m，長(zhǎng)1.45 km，自壓式進(jìn)入農(nóng)田。

3.4 抽水試驗(yàn)過(guò)程

臺(tái)蘭河地下水庫(kù)目前取水工程共布置4眼大口徑輻射井（F0－F3）和12口觀測(cè)井（G01－G12），其平面位置分布見圖4。抽水試驗(yàn)選定F2輻射井為抽水井（其余輻射井作為觀測(cè)井），抽水時(shí)間從2010年6月14日14：30開始，一直持續(xù)到2010年7月26日4：00。其中經(jīng)歷了10次停電，具體停電時(shí)間見表1。

圖4 輻射井和觀測(cè)井平面布置圖Fig.4 Layout diagram of radiation well and observation well

3.5 抽水試驗(yàn)的數(shù)值模擬

計(jì)算模型范圍 本次重點(diǎn)研究的是F1號(hào)輻射井抽水過(guò)程和水位恢復(fù)的滲流場(chǎng)及示范工程附近地下水滲流場(chǎng)變化。從抽水試驗(yàn)中各觀測(cè)井的觀測(cè)資料可看出，F(xiàn)1號(hào)輻射井在連續(xù)抽水25d后，地下水影響范圍約為500m。因此，為減少邊界水位對(duì)地下水滲流場(chǎng)的影響，將計(jì)算模型平面范圍取為以F1號(hào)輻射井為中心，半徑1.5km。研究區(qū)含水層深度為300～500m，垂直向?yàn)楦叱?00m的平面至地表面。

表1 抽水試驗(yàn)期間停電時(shí)間Table 1 Time of power off during pumping test

邊界條件及初始條件 研究區(qū)位于臺(tái)蘭河洪積扇與細(xì)土平原交界帶上的洪水沖溝內(nèi)，地下水流沿洪水沖溝由北向南流，地面坡降5‰～8‰，地下水的水力坡降約為4.4‰。根據(jù)觀測(cè)孔和輻射井抽水前的水位擬合得到了研究區(qū)初始穩(wěn)定滲流場(chǎng)（圖5），根據(jù)初始滲流場(chǎng)確定計(jì)算范圍內(nèi)北部和南部邊界為定水頭邊界，水頭值分別為1 202.5m和1 187.0m。由于抽水試驗(yàn)歷時(shí)較短（約26d），在此期間沒(méi)有大的降雨，只下過(guò)2d小雨，且該地區(qū)蒸發(fā)量較大，降雨對(duì)地下水位幾乎沒(méi)有影響，因此計(jì)算模型的地表面邊界為隔水邊界，不考慮降雨入滲的影響。

為了驗(yàn)證初始滲流場(chǎng)的合理性，將邊界條件代入模型，利用穩(wěn)定滲流計(jì)算得到了研究區(qū)的初始滲流場(chǎng)，各鉆孔的計(jì)算水位和觀測(cè)水位擬合情況見圖6。由圖6可看出，穩(wěn)定滲流場(chǎng)計(jì)算得到的各觀測(cè)井的計(jì)算水位和觀測(cè)水位較為接近，且均勻分布在45°線上，由此可說(shuō)明模型上下游邊界選擇較為合理。

圖5 計(jì)算區(qū)域的初始滲流場(chǎng)Fig.5 Initial seepage field of calculation area

圖6 初始滲流場(chǎng)各觀測(cè)孔計(jì)算值與觀值對(duì)比Fig.6 Contrast diagram of calculated and measured value of each aperture of sight of initial seepage field

時(shí)間步長(zhǎng)選取 整個(gè)抽水試驗(yàn)過(guò)程共進(jìn)行了約40d，由于抽水試驗(yàn)各觀測(cè)井均采用了電子自動(dòng)記錄水位（每30min記錄一次），數(shù)據(jù)較大，為了能夠反映整個(gè)抽水過(guò)程的水位變化，計(jì)算模型的時(shí)間步長(zhǎng)選為1d，每個(gè)時(shí)間步長(zhǎng)末端時(shí)刻為12：00；為充分反映停電過(guò)程對(duì)地下水位的影響，每次停電時(shí)刻和停電后水泵開啟的時(shí)刻單獨(dú)作為一個(gè)時(shí)間步長(zhǎng)。計(jì)算時(shí)間從2010年6月14日14：30到2010年6月30日14：30結(jié)束，共40個(gè)時(shí)間步長(zhǎng)。計(jì)算時(shí)間步長(zhǎng)內(nèi)輻射井F2水位隨時(shí)間變化曲線見圖7。

圖7 計(jì)算時(shí)間步長(zhǎng)內(nèi)F2水位隨時(shí)間變化曲線Fig.7 Curve of F2water level during calculation time step

圖8 研究區(qū)及輻射井子結(jié)構(gòu)單元的三維網(wǎng)格Fig.8 Mesh diagram of study area and radiation well substructure element

研究區(qū)網(wǎng)格剖分及輻射井子結(jié)構(gòu)網(wǎng)格剖分 圖8給出了輻射井子結(jié)構(gòu)單元的三維網(wǎng)格圖。按照上述計(jì)算模型范圍，結(jié)合輻射井和觀測(cè)井的位置建立了研究區(qū)的有限元分析模型，對(duì)輻射井F2進(jìn)行子結(jié)構(gòu)模型處理，其他觀測(cè)井都布置在單元結(jié)點(diǎn)上，計(jì)算區(qū)共剖分單元12 000個(gè)，結(jié)點(diǎn)10 487個(gè)，剖分網(wǎng)格見圖8a。根據(jù)輻射井F2的結(jié)構(gòu)設(shè)計(jì)（圖1），確定了輻射井的子結(jié)構(gòu)母單元網(wǎng)格的平切圖見圖8b，剖分形成的輻射井的三維子結(jié)構(gòu)單元見圖8c。

滲流場(chǎng)計(jì)算及計(jì)算結(jié)果分析 圖9給出各觀測(cè)井計(jì)算水位和觀測(cè)水位對(duì)比圖。圖10給出了第1次停電時(shí)的第2層輻射管平面的地下水頭等值線圖和剖面水頭等值線圖。圖11給出了第10次停電后水泵啟動(dòng)時(shí)刻的第2層輻射管平面的地下水頭等值線圖和剖面水頭等值線圖。利用觀測(cè)井計(jì)算結(jié)果對(duì)模型參數(shù)進(jìn)行了校正分析，表2給出了采用改進(jìn)的潛水斜率解析法、Neuman算法和輻射井子結(jié)構(gòu)法的計(jì)算反演等方法計(jì)算得到研究區(qū)平均滲透系數(shù)和給水度值。從表2可看出，反演得到的滲透系數(shù)比前2種方法稍大，給水度較為接近。主要原因是前2種方法不適合于輻射井抽水試驗(yàn)滲透系數(shù)反演。另外，由于數(shù)值分析模型對(duì)地層進(jìn)行了概化，沒(méi)有反映地層的各向異性，反演值存在一定的誤差。

圖9 各觀測(cè)井的計(jì)算水位與觀測(cè)水位對(duì)比Fig.9 Contrast diagram of calculated and measured water level of each aperture of sight

從圖9結(jié)合計(jì)算結(jié)果可知：所有觀測(cè)井計(jì)算水位的平均絕對(duì)誤差為0.22m，單井水位的平均誤差最大值為0.40m，最小值為0.02m，各觀測(cè)井地下水位的模擬值與實(shí)測(cè)值的變化趨勢(shì)吻合較好，說(shuō)明模型擬合的地下水位變化趨勢(shì)與研究區(qū)的水文地質(zhì)情況一致；另外需要說(shuō)明的是G12觀測(cè)井的觀測(cè)水位比計(jì)算水位稍偏高，主要原因是G12的西側(cè)是小阿吉力克溝，洪水和抽水試驗(yàn)的排水對(duì)其有影響，造成了觀測(cè)水位高；G03，G06，G07，G10和F1觀測(cè)井的水位動(dòng)態(tài)變化較大，而G09和G12水位動(dòng)態(tài)變化較小，其原因是G03，G06，G07，G10和F1距離輻射井F2距離較近，G09和G12距離輻射井F2較遠(yuǎn)。

圖10 第1次停電時(shí)計(jì)算區(qū)平面和剖面地下水頭等值線圖Fig.10 Plane and profile groundwater head isoline in calculation area during the first time of power off

以上計(jì)算結(jié)果說(shuō)明作者提出的輻射井子結(jié)構(gòu)法在理論算法上具備嚴(yán)密性，軟件計(jì)算結(jié)果具有可靠性。從圖10和圖11滲流場(chǎng)分布來(lái)看，整個(gè)研究區(qū)的等值線形態(tài)和走向分布合理，特別是在輻射井周圍，水頭等值線呈梅花狀分布，能夠精細(xì)反映輻射井的三維滲流行為和工作狀態(tài)。

圖11 最后1次停電后計(jì)算區(qū)地下水等值線圖Fig.11 Plane and profile groundwater head isoline in calculation area after the last time of power off

表2 不同計(jì)算方法得到的研究區(qū)水文地質(zhì)參數(shù)Table 2 Hydrogeology parameters obtained from different calculation methods

4 結(jié)語(yǔ)

基于輻射井的滲流行為、子結(jié)構(gòu)法的基本原理和非穩(wěn)定滲流場(chǎng)的有限元計(jì)算原理，建立了基于輻射井子結(jié)構(gòu)法的地下水非穩(wěn)定滲流數(shù)值模擬的數(shù)學(xué)模型和相應(yīng)算法。文中利用臺(tái)蘭河抽水試驗(yàn)的實(shí)測(cè)數(shù)據(jù)驗(yàn)證了數(shù)值模擬結(jié)果的正確性，進(jìn)而說(shuō)明了用該方法分析輻射井的局部滲流行為和有輻射井影響的地下水非穩(wěn)定滲流場(chǎng)是一種可行、有效的新方法。該方法考慮了輻射井中每根輻射管的結(jié)構(gòu)、幾何尺寸和滲流特性，以及各輻射管的相互干擾對(duì)整體滲流場(chǎng)的影響，具有較好的工程應(yīng)用前景。

另外，對(duì)布置虹吸輻射井排水系統(tǒng)的滲流場(chǎng)精細(xì)模擬方法和滲透參數(shù)反演方法、虹吸輻射井的排水效果的理論分析以及含多種復(fù)雜排水措施的非穩(wěn)定滲流場(chǎng)的精細(xì)求解方法等尚需進(jìn)一步研究。

（References）：

［1］王國(guó)強(qiáng)，徐國(guó)斌，吳道祥，等.輻射井技術(shù)在貧水地區(qū)的應(yīng)用［J］.水文地質(zhì)工程地質(zhì)，2002（5）：57－59.Wang Guoqiang，Xu Guobin，Wu Daoxiang，et al.Application of Radial Well in Aquitarde of Granite－Gneiss Weathered Zone［J］.Hydrogeology and Engineering Geology，2002（5）：57－59.

［2］汪良鋒，朱君星，項(xiàng)紅海，等.尾礦庫(kù)排滲輻射井的設(shè)計(jì)探討［J］.金屬礦山，2010（2）：55－61.Wang Liangfeng，Zhu Junxing，Xiang Honghai，et al.Discussion on Design of Radiation Wells for Tailings Seepage Drainage［J］.Metal Mine，2010（2）：55－61.

［3］王鐳，劉中，張有天.有排水孔幕的滲流場(chǎng)分析［J］.水利學(xué)報(bào)，1992（4）：15－20.Wang Lei，Liu Zhong，Zhang Youtian.Analysis of Seepage Field Near a Drainage－Holes Curtain［J］.Journal of Hydraulic Engineering，1992（4）：15－20.

［4］關(guān)錦荷，劉嘉忻，朱玉俠.用排水溝代替排水井列的有限元法分析［J］.水利學(xué)報(bào)，1984（3）：10－18.Guan Jinhe，Liu Jiaxin，Zhu Yuxia.Finite Element Analysis for an Array of Wells Replaced by a Drainage Ditch［J］.Journal of Hydraulic Engineering，1984（3）：10－18.

［5］Zhan Meili，Su Baoyu.New Method of Simulating Concentrated Drain Holes in Seepage Control Analysis［J］.Journal of Hydrodynamics ：Series B，1999（3）：27－35.

［6］朱岳明，張燎軍.滲流場(chǎng)求解的改進(jìn)排水子結(jié)構(gòu)法［J］.巖土工程學(xué)報(bào)，1997，19（2）：69－76.Zhu Yueming，Zhang Liaojun.Solution to Seepage Field Problem with the Technique of Improved Drainage Substructure［J］.Chinese Journal of Geotechnical Engineering，1997，19（2）：69－76.

［7］劉昌軍.虹吸井子結(jié)構(gòu)法在尾礦壩復(fù)雜滲流場(chǎng)求解中應(yīng)用［J］.水利水電科技進(jìn)展，2011（3）：81－85.Liu Changjun.Application of Siphonic Wells Substructure Method in Solution to Complicated Seepage Field of Tailings Dam［J］.Advance in Science and Technology of Water Resources，2011（3）：81－85.

［8］劉昌軍，丁留謙，徐澤平，等.面板裂縫對(duì)積石峽面板堆石壩滲流影響的數(shù)值模擬分析［J］.水利水電科技進(jìn)展，2011（6）：50－54.Liu Changjun，Ding Liuqian，Xu Zeping，et al.Numerical Simulation and Analysis of the Influence of the Cracks on the Seepage Field of Jishixia Concrete Face Rock－Fill Dam［J］.Advance in Science and Technology of Water Resources，2011（6）：50－54.

［9］劉昌軍，丁留謙，高立東，等.文登市抽水蓄能電站三維滲流場(chǎng)有限元分析［J］.水電能源科學(xué)，2011（7）：57－60.Liu Changjun，Ding Liuqian，Gao Lidong，et al.Finite Element Analysis of Three－Dimensional Seepage Field of Wendeng Pumped Storage Power Station［J］.Water Resources and Power，2011（7）：57－60.

［10］馬馳.SIP和PCG2兩種迭代法在地下水?dāng)?shù)值計(jì)算中的應(yīng)用對(duì)比［J］.西安科技學(xué)院學(xué)報(bào)，2002，22（3）：59－62.Ma Chi.Contrasting the Application of the SIP and PCG2Methods for Groundwater Matrix Numerical Solution［J］.Journal of Xi’an University of Science and Technology，2002，22（3）：59－62.

［11］劉昌軍，丁留謙，楊凱虹，等.水位驟降條件下有密集排水孔大壩三維滲流場(chǎng)數(shù)值模擬研究［J］.水文地質(zhì)工程地質(zhì)，2011（3）：24－30.Liu Changjun，Ding Liuqian，Yang Kaihong，et al.Study on Three Dimensional Numerical Simulation of Seepage Field of Dam with Dense Drainage Holes Under Sudden Drop of Water Level［J］.Hydrogeology and Engineering Geology，2011（3）：24－30.

［12］李濤，鄧銘江，王于寶，等.臺(tái)蘭河山前地下儲(chǔ)水構(gòu)造及地下水庫(kù)可行性研究［J］.水文地質(zhì)工程地質(zhì)，2011（1）：30－34.Li Tao，Deng Mingjiang，Wang Yubao，et al.A Study of the Feasibility of the Groundwater Storage Structure and Groundwater Reservoir in the Piedmont of the Talian River［J］.Hydrogeology and Engineering Geology，2011（1）：30－34.