王子捷, 肖 明

(1.武漢大學(xué)水資源與水電工程科學(xué)國家重點實驗室，湖北 武漢 430072；2.武漢大學(xué)水工巖石力學(xué)教育部重點實驗室，湖北 武漢 430072)

在水電站地下洞室群施工及運行過程中，排水措施往往會由于滲流結(jié)晶作用、水流沖蝕作用等各種原因堵塞失效[1-2]，造成地下水位逐漸抬升，易使主要洞室群發(fā)生滲透破壞，故采用數(shù)值模擬方法對比分析排水失效前后滲流場及滲流量的變化在工程計算分析中具有重要意義。目前，基于有限單元法的排水孔模擬方法有多種，其中較為成熟的方法有空氣單元法[3]、子結(jié)構(gòu)法[4]、隱式復(fù)合單元法[5]等，但這些方法多數(shù)是考慮穩(wěn)定滲流的情況，對于排水孔的非穩(wěn)定滲流效應(yīng)研究較少。此外，目前多數(shù)研究僅針對于排水失效前后地下穩(wěn)定滲流場的情況，而實際上，排水措施失效過程及失效后滲流場的變化都是隨時間變化的非穩(wěn)定過程，僅僅采用穩(wěn)定滲流分析存在一定局限性。

本文基于排水孔模擬的隱式復(fù)合單元法[5]，考慮排水孔的給水度、單位貯存量等非穩(wěn)定滲流參數(shù)的等效方法，提出非穩(wěn)定滲流情況下排水孔的模擬方法，同時，結(jié)合非穩(wěn)定滲流的拋物型變分不等式提法[6]，探究地下滲流場在排水措施失效后隨時間的變化規(guī)律。最后，采用簡單的模型算例，驗證方法的合理性，并探究了排水失效后地下水隨時間的變化規(guī)律。

1 排水孔失效后非穩(wěn)定滲流場模擬方法

1.1 非穩(wěn)定滲流的拋物型變分不等式提法

由陳益峰等提出的非穩(wěn)定滲流的拋物型變分不等式提法[6]，將滲流溢出邊界上的Signorini型邊界條件及內(nèi)部自由面邊界條件轉(zhuǎn)化為自然邊界條件，從而極大程度地改善了數(shù)值計算的收斂性及穩(wěn)定性。

定義試探函數(shù)ΦPVI為

(1)

則該提法可表述為：在ΦPVI中求一隨時間變化的函數(shù)φ，使對?ψ∈ΦPVI，都有：

(2)

上述即為非穩(wěn)定滲流的拋物型變分不等式提法。

1.2 排水孔模擬的復(fù)合單元法

巖土材料的滲透特性一般是各向異性且非均勻的，而排水孔的滲透特征亦不同于巖體，它們的滲透系數(shù)、給水度、單位貯存量也各不相同，此時可用復(fù)合單元法來確定單元的滲流參數(shù)[5]。

排水孔可視為內(nèi)含在巖體中的空氣單元，而空氣可以視作一種強透水介質(zhì)。當(dāng)排水孔內(nèi)含在巖體單元或穿過巖體單元時，則可將其視為巖體單元中的一個子域，此時包含排水孔的巖體單元就成為了復(fù)合單元。

在模型進行離散，建立模型網(wǎng)格時，首先建立不考慮排水孔的整體有限元網(wǎng)格。模型網(wǎng)格建成后，再組織排水孔的起、終點單元號、坐標(biāo)、直徑等信息數(shù)據(jù)。在后續(xù)數(shù)值計算中，只需根據(jù)實體單元與排水孔之間的關(guān)系，修正包含排水孔或有排水孔穿越的實體單元的滲透傳導(dǎo)矩陣及單位貯存量、給水度等參數(shù)，這樣就可以實現(xiàn)對排水孔的數(shù)值模擬。在該方法中，排水孔隱含于巖體中，只需建立整體有限元網(wǎng)格，而不需要對排水孔進行實際單元模擬。故該方法有建模簡單，便于應(yīng)用的優(yōu)點。

由于排水孔的強排水作用，沿其軸向的滲透性顯著增強，垂直于軸向一定范圍的滲透性也有所加強，因此排水孔的走向與傾角會造成單元的滲透傳導(dǎo)矩陣有所不同。此時可將排水孔視做空間圓柱體，求出其軸向方向向量即為其主滲透方向(見圖1)。假定沿排水孔軸向向下為其局部坐標(biāo)系x′軸，其x′軸向主滲系數(shù)為kx′。為便于計算，令ky′的主滲透方向取在局部坐標(biāo)系x′軸與整體坐標(biāo)系y軸所形成的平面內(nèi)，可以求出y′軸的方向向量。取x′軸與y′軸的方向向量叉積為z′軸的方向向量。

圖1 排水孔滲透方向示意圖

如圖2所示的單元體內(nèi)，有兩根排水孔穿過，排水孔A與單元表面交點為2，排水孔B與單元表面的交點為1和3。根據(jù)其交點結(jié)合排水孔設(shè)計求出其軸向滲透系數(shù)kx′A及kx′B，將其平移至相交的平面內(nèi)疊加得到等效排水孔的滲透系數(shù)kx′(如圖3)。同理，也可以得到另外兩個方向的滲透系數(shù)ky′和kz′。

圖2 單元內(nèi)排水孔位置示意圖

圖3 多排水孔空間向量疊加圖

考慮每個排水孔在復(fù)合單元中的狀態(tài)，對單元內(nèi)的所有排水孔進行等效，然后計算得到每個主滲透系數(shù)與坐標(biāo)軸之間的方向余弦lx，ly，lz，組成坐標(biāo)轉(zhuǎn)換矩陣[Rp]：

(3)

則可在整體坐標(biāo)下得到排水孔的等效傳導(dǎo)矩陣：

[Kp]=[Rp]T[K′][Rp]

(4)

式中，[Kp]為排水孔的等效滲透張量。

根據(jù)單元滲透傳導(dǎo)矩陣計算原理，可將含排水孔復(fù)合單元的滲透傳導(dǎo)矩陣表示為：

(5)

式中：[A]為復(fù)合單元滲透傳導(dǎo)矩陣；vp為單元內(nèi)排水孔所占體積；vr為巖體體積。

將含排水孔復(fù)合單元視作一個整體，設(shè)其等效滲透矩陣為[Kd]，則整個復(fù)合單元的滲透傳導(dǎo)矩陣又可以表示為：

(6)

聯(lián)立式(5)、(6)，則可以求得復(fù)合單元的等效滲透傳導(dǎo)矩陣。

在非穩(wěn)定滲流計算中，除了滲透系數(shù)之外，還需考慮復(fù)合單元單位貯存量及給水度的等效。

對于單位貯存量Ss而言，根據(jù)其定義，可以認為其是各項同性參數(shù)。在工程實際中，排水孔一般不承擔(dān)貯水的作用，可以認為其正常工作時單位貯存量為0。由此，可直接采用體積等效的形式獲取復(fù)合單元整體的單位貯存量取值。如圖2所示的復(fù)合單元，所有排水孔在單元內(nèi)的體積vp，單元內(nèi)巖體體積為vr，單元體積為v，則有：

vr=v-vp

(7)

Ssrvr+Sspvp=Ssdv

(8)

式(8)中：Ssr、Ssp分別為巖體及排水孔的單位貯存量；Ssd為復(fù)合單元等效的單位貯存量。

求解式(8)，即可獲得復(fù)合單元的等效貯水率。

對于如圖4所示有自由面穿越的復(fù)合單元，還需求取其等效給水度。

圖4 有自由面穿越的復(fù)合單元示意圖

對于該單元體而言，將復(fù)合單元視作整體，其系數(shù)矩陣中與給水度有關(guān)的元素可記為：

(9)

將排水孔和巖體分開，復(fù)合單元的自由面邊界的系數(shù)矩陣中的元素又可寫成：

(10)

(11)

聯(lián)立式(10)、(11)即可求取有自由面穿越的復(fù)合單元的等效給水度。

1.3 排水孔逐漸失效過程模擬

排水系統(tǒng)在實際運行中，由于結(jié)晶作用，會隨著運行時間增長逐漸堵塞失效，其是一個涉及到溫度、pH、CO2、排水孔內(nèi)水流流速等多種因素影響的一個復(fù)雜過程，其堵塞失效的規(guī)律及其何時完全失效亟需系統(tǒng)性研究，本文為簡化計算，將排水孔由于結(jié)晶作用逐漸堵塞失效的過程視為排水孔孔徑隨時間增加勻速減小的過程，即排水孔的截面積及體積不斷減小，當(dāng)孔徑減小到0時，排水孔即等效于巖體，此時視為排水完全失效。其具體計算形式如下：

(12)

式中：D0為排水孔正常運作時的孔徑；T0為假定的排水孔完全失效的時間；t為已運行時長；Δt為非穩(wěn)定滲流計算的時間步長。

結(jié)合其他勘測資料及邊界信息，首先采用穩(wěn)定滲流方法獲取有排水措施情況下的地下滲流場分布，再以此為非穩(wěn)定滲流計算的初始條件，隨時間變化逐漸改變排水孔的孔徑大小，探究排水逐漸失效后地下滲流場的變化規(guī)律，具體流程見圖5。

圖5 排水措施失效過程非穩(wěn)定滲流計算流程圖

2 算例分析

2.1 計算模型及計算參數(shù)

為探究排水孔失效后滲流場變化規(guī)律，設(shè)置一個含排水孔的地下廠房簡化模型。模型范圍為250 m×250 m×250 m(x×y×z)，主廠房的尺寸為30 m×80 m×50 m(長×寬×高)，廠房底部距模型底部60 m，位于研究區(qū)域正中心，排水孔布置在廠房周圍，間距5 m，直徑80 mm，距離廠房邊界20 m，頂部排水孔呈45°斜向上。研究范圍內(nèi)整體有限元模型及廠房、排水孔布置見圖6及圖7。

圖6 整體有限元模型圖

圖7 地下廠房及隱式排水孔布置圖

計算條件設(shè)置如下：研究區(qū)域巖體的滲透系數(shù)取為1.0×10-7m/s，給水度取0.1，單位貯存量取1.0×10-6m-1，排水孔等效滲流參數(shù)按文獻[7]中建議取值。沿x方向左右兩側(cè)設(shè)置為定水位邊界，水位均為200 m，模型底面及沿y方向前后兩側(cè)面取為隔水邊界，主廠房洞室邊界取為潛在溢出型邊界，排水孔按前文所述的復(fù)合單元法進行模擬。

計算首先考慮排水孔正常作用的情況，然后以此情況下的計算結(jié)果作為初始條件，并假定排水系統(tǒng)在運行的第30 d后完全失效，按式(12)模擬其逐漸失效的過程，研究排水孔失效過程中地下滲流場的變化規(guī)律、主廠房滲流量的變化情況，并預(yù)測滲流場達到最危險情況的時間。

非穩(wěn)定滲流計算的總時長取為2個月，計算時步長取為2 d。

2.2 結(jié)果分析

不同時刻滲流場計算結(jié)果如圖8所示。當(dāng)?shù)叵聫S房處于完備滲控措施作用下時(t=0)，地下水自由面有效降低，地下水滲漏點僅出現(xiàn)在廠房底部，說明復(fù)合單元法能有效體現(xiàn)排水孔列陣的排水功能。當(dāng)排水孔逐漸失效后，自由面逐漸升高，且水位變化呈現(xiàn)出前期抬升速度較快、后期抬升速度較慢的規(guī)律，這是由于滲透壓力在前期排水孔失效過程中變化較大。在排水孔逐漸失效過程中，由于其孔徑減小，復(fù)合單元的滲流參數(shù)不斷變化，導(dǎo)致滲流場變化較大。30 d后，排水措施完全失效，此時滲流場雖未達到穩(wěn)定狀態(tài)，但此時包含排水孔的復(fù)合單元已視為普通的巖體單元進行計算，滲流參數(shù)不再改變，滲流場的變化只是隨著時間的增加逐漸趨于穩(wěn)定狀態(tài)。對于具體時刻的滲流場變化而言，在失效14 d后，地下水將排水系統(tǒng)及主廠房全部淹沒。當(dāng)?shù)竭_第54 d后，自由面趨于穩(wěn)定，不再變化。此時，自由面最低位置高出廠房頂拱50多米，廠房處于極其危險的狀態(tài)，在實際工程中應(yīng)盡量避免排水孔完全失效情況的發(fā)生。

圖8 排水失效后滲流自由面變化圖

圖9為廠房洞壁滲流量隨時間變化示意圖。由圖可知，當(dāng)排水孔完全生效時(0 d時)，主廠房滲流量僅有30.4 m3/d，滲流量較小，廠房運行安全。當(dāng)排水孔逐漸失效后，主廠房滲流量變化規(guī)律與自由面變化規(guī)律相一致，在前期迅速增加，一段時間之后趨于穩(wěn)定。地下水穩(wěn)定后，主廠房滲流量達到152.3 m3/d，是排水孔完全生效時的5倍多，極易發(fā)生滲透破壞，說明排水措施對廠房安全運行起到至關(guān)重要的作用。

總體而言，本文提出的方法能夠體現(xiàn)出排水失效后地下水非穩(wěn)定變化規(guī)律，能為后續(xù)的工程實際計算提供一定參考。

圖9 廠房滲流量隨時間變化圖

3 結(jié) 語

本文基于非穩(wěn)定滲流的拋物型變分不等式法及排水孔模擬的隱式復(fù)合單元法，提出排水孔模擬的非穩(wěn)定滲流方法，并考慮了排水孔逐漸失效的過程，結(jié)合一個簡易的地下廠房模型探究了排水措施失效后地下水的非穩(wěn)定滲流規(guī)律，可得到以下結(jié)論：①在排水系統(tǒng)逐漸失效的過程中，地下水自由面迅速升高，呈現(xiàn)出前期升高速度快、后期升高速度慢的規(guī)律，最后趨于穩(wěn)定；②滲流量變化規(guī)律與自由面變化規(guī)律相近，達到穩(wěn)定后廠房流量大小是有排水時的5倍，實際工程中應(yīng)盡量避免排水失效情況的發(fā)生。