臺蘭河地下水庫輻射井抽水過程的非穩(wěn)定滲流場的有限元分析

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

中國水利水電科學研究院，北京 100038

0 引言

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

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

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

1.1 輻射井結構和組成

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

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

1.2 輻射井的滲流行為

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

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

圖2 輻射管濾管內的滲流行為Fig.2 Drainage principle of radiation wells

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

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

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

式中，h0為輻射井豎井內的水位。

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

式中，z0為aa′或bb′逸出面上結點位置高程。

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

2.1 數(shù)學模型

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

式中：Q為源或匯項（1/s）；SS為貯水率（1/m），即單位體積的飽和土體內，當下降1個單位水頭時，由于土體壓縮和水的膨脹所釋放出來的貯存水量。

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

自由面邊界條件為

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

2.2 有限元法求解

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

在此基礎上，對時間域進行合理離散，可得到有限元時間隱式差分迭代格式：

式中：K為滲透系數(shù)矩陣；S為總質量矩陣；G為自由面貢獻矩陣；p為結點壓力列陣；F為等效結點流量列陣；Δt為時間步長。

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

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

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

2.4 預處理共軛梯度算法

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

PCG法假設矩陣A被分為2個矩陣，即：

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

2.5 輻射井子結構法

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

對于子結構內部邊界的處理與外部邊界相同即可，詳見參考文獻［11］。

3 工程應用

3.1 臺蘭河地下水庫工程概況

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

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

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

圖3 臺蘭河地下水庫工程平面布置圖Fig.3 Layout diagram of groundwater reservoir of Tailan River

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

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

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

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

3.4 抽水試驗過程

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

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

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

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

表1 抽水試驗期間停電時間Table 1 Time of power off during pumping test

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

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

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

圖6 初始滲流場各觀測孔計算值與觀值對比Fig.6 Contrast diagram of calculated and measured value of each aperture of sight of initial seepage field

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

圖7 計算時間步長內F2水位隨時間變化曲線Fig.7 Curve of F2water level during calculation time step

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

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

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

圖9 各觀測井的計算水位與觀測水位對比Fig.9 Contrast diagram of calculated and measured water level of each aperture of sight

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

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

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

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

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

4 結語

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

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

（References）：

［1］王國強，徐國斌，吳道祥，等.輻射井技術在貧水地區(qū)的應用［J］.水文地質工程地質，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］汪良鋒，朱君星，項紅海，等.尾礦庫排滲輻射井的設計探討［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］王鐳，劉中，張有天.有排水孔幕的滲流場分析［J］.水利學報，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］關錦荷，劉嘉忻，朱玉俠.用排水溝代替排水井列的有限元法分析［J］.水利學報，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］朱岳明，張燎軍.滲流場求解的改進排水子結構法［J］.巖土工程學報，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］劉昌軍.虹吸井子結構法在尾礦壩復雜滲流場求解中應用［J］.水利水電科技進展，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］劉昌軍，丁留謙，徐澤平，等.面板裂縫對積石峽面板堆石壩滲流影響的數(shù)值模擬分析［J］.水利水電科技進展，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］劉昌軍，丁留謙，高立東，等.文登市抽水蓄能電站三維滲流場有限元分析［J］.水電能源科學，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兩種迭代法在地下水數(shù)值計算中的應用對比［J］.西安科技學院學報，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］劉昌軍，丁留謙，楊凱虹，等.水位驟降條件下有密集排水孔大壩三維滲流場數(shù)值模擬研究［J］.水文地質工程地質，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］李濤，鄧銘江，王于寶，等.臺蘭河山前地下儲水構造及地下水庫可行性研究［J］.水文地質工程地質，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.