劉昌軍
(中國水利水電科學研究院,北京 100038)
尾礦壩的壩體是由透水性強弱不等的粗細尾礦砂堆筑而成,而壩坡浸潤線是尾礦壩的生命線,它是直接影響壩體安全的一個非常重要的因素。我國多數尾礦壩采用上游法筑壩工藝,其顯著特點是隨著尾礦堆筑高度的不斷加大,尾礦壩壩體內浸潤線也不斷提高,導致壩面溢出水位升高,易引起壩體滲流破壞,危及壩體穩定性。因此需要對尾礦壩設計一定的排水系統以降低壩體浸潤線的高度。
在尾礦壩設計與施工中,虹吸井、導滲盲溝和輻射井等排滲措施都有廣泛的應用,但虹吸井和導滲盲溝往往徑向尺寸很小而縱向尺寸很大,數量眾多,且空間分布復雜。因此在用有限元模擬數量眾多的虹吸井和導滲盲溝作用下的滲流場分布就非常困難。目前堤壩滲控分析中對密集排水孔的模擬較多[1-5]。對尾礦壩虹吸井、導滲盲溝以及輻射井等排水措施的排水效果和滲流場進行數值模擬研究的文獻尚不多見。
為了解決難以精細模擬密集虹吸井和導滲盲溝的排水效果的難題,筆者提出一種適合求解虹吸井、導滲盲溝和輻射井等布置有復雜排水系統的尾礦壩三維滲流場的虹吸井子結構方法,該方法較好地解決了尾礦壩排水系統的精細模擬問題,可用于多種排水措施共同作用下多種排水方案的計算比較。以某尾礦壩為例,采用虹吸井子結構法研究了密集虹吸井和導滲盲溝作用下的滲流場分布,并對不同排滲方案的排滲效果進行對比分析,給出了該尾礦壩排水系統的合理布置方案。
非均質各向異性多孔介質的穩定飽和滲流問題的控制方程為


式中,kij為二階對稱的達西滲透系數矩陣;h為水頭;h1為已知水頭函數;ni為滲流邊界面外法線方向余弦,i=1,2,3;Γ1為第一類滲流邊界條件;Γ2為第二類滲流邊界條件;Γ3為滲流自由面;Γ4為滲流逸出面;qn為法向流量,以流出為正;x3為z方向的位置高程。
對式(1)~(5)可采用固定網格的有限單元法進行求解,根據變分原理,求解泛函數和支配方程為

式中:∏(h)為泛函數;K,P,F分別為飽和區的傳導矩陣、節點水頭列陣和已知節點壓力列陣;Ω為積分域。
文獻[6]提出了有自由面滲流場的求解方法——截止負壓法,文獻[7]在文獻[6]的基礎上對截止負壓法進行了改進。本文采用文獻[7]中改進的截止負壓法對滲流場進行求解。改進的截止負壓法的詳細公式推導見文獻[7]。
虹吸井排滲系統主要由虹吸井、觀測井、水封槽和虹吸管路組成[8],見圖1。其中,虹吸井是關鍵的降水設施,其成井質量要求較高。一般虹吸管內水頭差約8m左右。水封槽的作用在于保證虹吸系統的絕對真空度,即在虹吸系統啟動后應保證排水管出口始終位于一定的水位高度下。
虹吸井是靠虹吸作用實現自流排水。虹吸井滲流行為在算法上分為2種情況:一是虹吸井穿過自由面;二是虹吸井不與自由面相交。虹吸井穿過自由面又分為2種情況:一是虹吸井達到穩定狀態,井內水位和水封槽水位相同;二是虹吸井內水位大于水封槽水位,虹吸井內出水量較大,虹吸井以恒定流量出水。有自由面穿過時,虹吸井內邊界面上的滲流行為如圖2所示。

圖1 虹吸井排滲系統布置圖

圖2 虹吸井滲流行為
a.如果自由面水位低于cc′,則整個虹吸井起不到排水作用,虹吸井失效。因此可以在每個虹吸井水封槽水位高程處虛構一個數學開關器,詳見文獻[9]和文獻[10]。虹吸井計算時,先假設開關器打開,孔內邊界全部作為可能逸出邊界,每一步迭代后對孔內節點壓力進行甄別。
若虹吸井底面孔口節點壓力大于零,認為虹吸井全部位于非飽和區,此時虹吸井完全失效。其邊界條件數學表達式為

若虹吸井cc′處節點壓力大于零,認為虹吸井水封槽高程以下部分位于飽和區,此時邊界為定水頭邊界。其數學表達式為

式中:x3c為c點處z方向的位置高程。
若虹吸井cc′處節點壓力小于零,且底面節點壓力大于零,認為虹吸井自由面位于水封槽底部,此時虹吸井失效。
b.如果自由面水位位于bc或b′c′內,且逸出流量小于虹吸井設計恒定流量 Q,則 ac和a′c′面水頭為定水頭,水頭值為水封槽的高程,此時 bc或b′c′面為隔水邊界,因為虹吸井起動直至水位穩定,虹吸井內水位和水封槽高程保持一致。邊界條件數學表達式為

c.如果自由面水位位于bc或b′c′內,且逸出流量大于虹吸井恒定出水量,則虹吸井以恒定抽水量抽水,此時虹吸井邊界條件為流量邊界條件(但實際工程中,虹吸井設計恒定流量一般大于虹吸井內逸出流量)。其虹吸井內邊界條件數學表達式如下:

虹吸井在改進的截止負壓法計算中,作為二類邊界(可能逸出邊界)來處理,即初始迭代是把孔壁內節點高程作為邊界條件,然后逐步判別各節點是在飽和區還是在非飽和區。根據逸出點位置進行虹吸井孔壁逸出流量計算,計算逸出流量并與虹吸井設計恒定流量進行比較。然后根據上述情況轉化虹吸井邊界條件,進而實現虹吸井滲流行為的真實數值模擬。
對于導滲盲溝和輻射井的處理,同樣采用上述子結構方法,子結構內邊界按可能逸出邊界處理。
按照上述虹吸井滲流行為,結合改進截止負壓法原理,將子結構作為主網格整體結構的一部分一起來考慮[2,3,9,11],則虹吸井子結構總傳導矩陣及相應流量列陣寫成分塊形式,如式(12)所示。

式中:Kii為子結構內部未知節點水頭相互作用傳導矩陣;K ib,K bi分別為子結構內部未知節點水頭與出口未知節點水頭之間相互作用的傳導矩陣;K bb為子結構出口未知節點水頭相互作用傳導矩陣;p i,p b分別為子結構內部和子結構出口未知節點水頭列陣;Fi,Fb分別為子結構已知節點水頭對內部未知節點水頭和出口未知節點水頭貢獻流量(包括內部源匯項及非零流量邊界的貢獻流量)的列陣。
由式(13)~(14)可求得相應未知節點出口傳導矩陣及相應流量貢獻列陣K*bb,F*b如下:

用K*bb,F*b參與總平衡方程系數及常量矩陣的組裝,求解 pb后,回代求解 pi。

對于子結構內部已知水頭節點,其節點壓力為已知,因而其每步迭代求得的壓力增量為零,故可將滲透矩陣中該節點對應的對角元素乘以一個大數即可。
對于子結構內部可能逸出面上的節點,如果發現其節點壓力大于零,可以采用處理截止負壓的方法進行類似處理,即將滲透矩陣中該節點對應的對角線元素置為一個大數 λ,同時將其對應的節點不平衡力賦為-λp(p為該節點的壓力),因此可以迫使p=0。
在尾礦壩工程中,為精確反映虹吸井的三維尺寸效應,虹吸井通常是在特定區域內的母單元中進行剖分。虹吸井的布置方式和模式因工程特點和排水要求各不相同,在虹吸井子結構法中,虹吸井的子結構形式和子結構網格剖分模式參照文獻[2,9,11]中介紹的方法。筆者采用IDL語言編寫了有限元數值模擬軟件GWSS,該軟件提供了有限元網格剖分、虹吸井子結構網格剖分、滲流計算和滲流結果等值面(線)的顯示等功能。在虹吸井子結構網格剖分方面,該軟件實現了不同方向的虹吸井子結構單元網格的二次剖分和彎管子結構單元的剖分。其剖分原理和剖分步驟如下:①根據整體網格單元和事先控制的虹吸井和導滲盲溝的超單元,確定需要剖分子結構單元的母單元的節點信息和單元信息,同時搜索其相鄰子結構單元的母單元串。②找出最小單元截面周長l,虹吸井等效母單元邊長 a由公式a=πd/l(d為虹吸井的直徑)算出,進而計算母單元中子結構單元的局部節點坐標,形成局部節點和單元信息。③根據母單元新增單元和站點信息組裝子結構信息,進而得到整體網格和子結構網格的單元信息。
河北某尾礦壩,尾礦砂為中等顆粒,干密度為1.4 t/m3。該尾礦壩初期設計為碾壓土壩,壩高20m,壩頂寬度為3.5m,壩頂高程為450m,壩底高程為430m。外坡坡比為1∶2.5,內坡坡比為1∶2.0。現初期壩頂已被碎石覆蓋,碎石覆蓋后的壩頂高程為459.6m,外坡坡比為1∶1.36。初期壩壩坡多處有清水滲出。
可將尾礦壩地層分為3個部分,即天然地層、初期壩和尾礦堆積層。天然地層表面為殘坡積含礫粉質黏土,其下部為粉土及粉質黏土充填的碎石土層(強風化層)。初期壩為碾壓土壩。由于尾礦沉積的分選性不強,故砂性尾礦中普遍存在黏性尾礦夾層(3~10cm厚),且在同平面上不連續,呈透鏡體狀,但從整體上看仍存在一定規律:垂直層序上存在上粗下細的規律;水平層序上由堆積子壩向庫區尾礦砂由粗變細。根據室內篩分試驗判定該區尾礦砂大部分為尾礦中砂。

表1 尾礦壩各地層滲透系數 m/s
計算區域采用八節點六面體單元進行剖分,共剖分節點13120個,單元10500個,其三維網格見圖3,尾礦壩虹吸井和導滲盲溝等排水系統采用子結構進行剖分,其三維網格見圖4、圖5。
為節省篇幅,本文只給出了工況2和工況4計算結果。圖6和圖7分別為工況2和工況4順河向剖面水頭等值線,表2給出了不同工況設置虹吸井后虹吸井出水量。
從計算結果可以看出:①各工況滲流場的水頭分布規律合理,水頭等值線形態、走向和疏密程度都

圖4 虹吸井子結構網格(局部放大)

圖5 導滲盲溝子結構網格(局部放大)

表2 各工況虹吸井出水量 m3/d
尾礦壩各地層滲透系數見表1。工況1和工況2分別為正常蓄水位和設計洪水位下的滲流場分布,工況3和工況4分別為下游側虹吸井和上游側虹吸井失效的滲流場分布。準確反映了相應區域防滲和排水滲控措施特點、滲流特性和邊界條件,計算區域內虹吸井的排水效果得到了精細模擬,其滲控效果得到了準確反映。②采用虹吸井子結構方法可以精細模擬每個虹吸井的排水作用,可以計算單個虹吸井的排水量,如工況2,上排單井平均出水量為108.500m3/d,下排虹吸井平均出水量為7.590m3/d。③從工況3和工況4的計算結果看,下游側虹吸井完全失效后,地下水位逸出點高程為510m,而上游側虹吸井失效后,地下水位逸出點高程為523m。和工況2相比,下游側和上游側虹吸井失效后尾礦壩逸出點高程都相應增加,嚴重影響壩體穩定性,因此建議采用2排虹吸井設計方案,并確保虹吸井施工質量,防止虹吸井失效。④從工況1~4的計算結果可以看出,尾礦壩設置排水棱體的排水效果較差,只有最下層排水棱體起到排水作用,因此建議去掉上面2層排水棱體。

圖6 工況2順河向剖面水頭等值線(單位:m)

圖7 工況4順河向剖面水頭等值線(單位:m)
a.虹吸井排滲系統設計原理簡單,可操作性強,排滲效果良好,國內尾礦壩內使用較為廣泛,采用子結構技術模擬虹吸井的滲流行為效果良好,算法理論嚴密可靠。
b.在改進截止負壓法中引入虹吸井子結構算法,較好地解決了虹吸井作用下滲流場的求解問題,進而可研究尾礦壩密集虹吸井排水降壓效果。
c.工程算例的計算結果表明,采用虹吸井子結構法可以很好地模擬尾礦壩在虹吸井、導滲盲溝等排水措施作用下的滲流場分布,為尾礦壩滲流控制設計提供理論支持和技術指導。
[1]王鐳,劉中,張有天.有排水孔幕的滲流場分析[J].水利學報,1992(4):15-20.
[2]朱岳明,陳振雷.改進的排水子結構法求解地下廠房洞室群區的復雜滲流場[J].水利學報,1992(9):79-85.
[3]趙堅,沈振中.尾礦壩復雜排水系統滲流計算方法的改進[J].河海大學學報:自然科學版,1997(3),25(2):110-113.
[4]關錦荷,劉嘉忻,朱玉俠.用排水溝代替排水井列的有限元法分析[J].水利學報,1984(3):10-18.
[5]ZHAN Mei-li,SU Bao-yu.New method of simulating concentrated drain holes in seepage control analysis[J].Journal of Hydrodynamics:SerB,1999(3):27-35.
[6]速寶玉,沈振中,趙堅.用變分不等式理論求解有滲流問題的截止負壓法[J].水利學報,1996(3):22-29.
[7]張乾飛,吳中如.有自由面非穩定滲流分析的改進截止負壓法[J].巖土工程學報,2005,27(1):49-54.
[8]王鳳江.虹吸排滲系統在尾礦壩降水工程中應用[J].中國地質災害與防治學報,2005,16(1):102-105.
[9]朱岳明,陳建余,龔道勇,等.拱壩壩基滲流場的有限單元法精細求解[J].巖土工程學報,2003,25(3):326-330.
[10]崔皓東,朱岳明,吳世勇.有自由面滲流分析中密集排水孔幕的數值模擬[J].巖土工程學報,2008,30(1):440-445.
[11]陳建余.有密集排水孔幕的飽和-非飽和滲流場全精細數值模擬分析[J].中國農村水利水電,2004(11):53-60.