曾楊晨，李雙龍，王 鵬，譚文帥

(1.江西省南昌縣象湖聯(lián)圩管理站，江西 南昌 330200；2.長沙理工大學(xué)水利工程學(xué)院，湖南 長沙 410076；3.湖南省株洲市佳聯(lián)建設(shè)工程有限公司，湖南 株洲 410004)

0 引 言

壩基及近壩山體抬升變形現(xiàn)象在國內(nèi)外相對少見。國外的托克托古爾重力壩[1]、英古里重力壩[2]在蓄水后都發(fā)現(xiàn)水工結(jié)構(gòu)物有抬升變形現(xiàn)象；我國湖南省江埡大壩近壩山體最大抬升位移達21.8 mm[3，4]，銅街子水電站壩體右岸山體表面產(chǎn)生24.3～28.9 mm抬升位移[5]，金沙江向家壩左岸近壩山體抬升量超過13.0 mm[6]等。針對以上抬升變形問題，國內(nèi)外學(xué)者進行了研究。Xiaoli Liu[7]等認為銅街子大壩及右岸近壩山體抬升變形的主要原因是分布在壩址區(qū)承壓含水層作用和壩基強度參數(shù)的弱化以及構(gòu)造應(yīng)力的組合影響；王蘭生[8]等利用MTS剛性壓力試驗機對巖體擴容行為進行試驗?zāi)M，認為承壓含水層D2y內(nèi)的孔隙水壓力增高是造成壩基和近壩山體抬升的主要原因；伍法權(quán)[9]等結(jié)合監(jiān)測數(shù)據(jù)與江埡水庫地質(zhì)資料，認為熱水含水層中的有效應(yīng)力減小，導(dǎo)致巖體卸荷擴容；祁生文[10]等利用FLAC3D數(shù)值軟件對江埡大壩進行了抬升變形過程模擬，認為水庫蓄水是抬升變形的本質(zhì)原因。Fuzhang Yan[11]等通過三維數(shù)值模擬，得到了抬升變形的空間位移分布；蔣中明等[6]提出了水庫蓄水條件下的抬升變形水文地質(zhì)結(jié)構(gòu)模式，探討了不同壩基地質(zhì)條件對抬升變形的影響。

由此可見，針對壩基及近壩山體抬升變形問題，現(xiàn)有成果側(cè)重于外因研究，即抬升變形的本質(zhì)原因在于水庫蓄水，同時結(jié)合一定的特殊地質(zhì)條件導(dǎo)致巖體有效應(yīng)力降低，使巖體產(chǎn)生擴容變形，然而在抬升變形的內(nèi)因條件方面研究不足，如巖體參數(shù)、水文地質(zhì)條件等。本文擬基于應(yīng)力場-滲流場耦合理論，充分考慮蓄水條件下滲透力作用，采用三維數(shù)值模擬方法，研究壩基及近壩山體中巖體變形模量及滲透系數(shù)對抬升變形的影響規(guī)律，進一步揭示抬升變形的變形機制。

1 抬升變形數(shù)值理論

1.1 流固耦合基本方程

根據(jù)文獻[6]、[9]、[10]，水庫抬升變形的誘因是水庫蓄水，蓄水條件下壩基巖體在巨大滲透力作用下產(chǎn)生上抬變形，采用巖土數(shù)值軟件FLAC3D可以實現(xiàn)對抬升變形過程的模擬。本文擬基于應(yīng)力場-滲流場耦合理論，利用FLAC3D內(nèi)嵌fish語言[12]，模擬壩基及近壩山體不同巖體參數(shù)在蓄水條件下的抬升過程。FLAC3D在模擬巖土體流固耦合問題中，將巖土體都作為孔隙介質(zhì)來考慮，符合Darcy[13]定律，滿足Biot[14]方程。下面給出FLAC3D計算采用的流固耦合基本方程。

(1)平衡方程

假定小變形條件下，流體質(zhì)點平衡方程為：

-qi，j+qv=?ζ/?t

(1)

式中：qi為滲流速度，m/s；qv為流體源強度，1/s；ζ為單位體積流體體積變化量；t為時間，s。又有：

?ζ/?t=(1/M)×?p/?t+α(?ε/?t)-β(?T/?t)

(2)

式中：M為Biot模量，N/m2；p為孔隙水壓力，Pa；α為Biot系數(shù)；ε為體積應(yīng)變；β為熱膨脹系數(shù)(1/℃)；T為溫度，℃。此方程考慮了流固耦合作用中的溫度效應(yīng)，在本文研究中并沒有考慮地下水溫度的影響，故本文計算方程溫度項可刪除。

(2)運動方程

假定介質(zhì)為均質(zhì)、各向同性且密度不變，利用Darcy定律描述運動方程為：

qi=-k[p-ρfxjgj]

(3)

式中：k為介質(zhì)滲透系數(shù)，m2/Pa·s；ρf為流體密度，kg/m3；gi，(i=1，3)為重力加速度的三個分量，m/s2。

(3)本構(gòu)方程

流固耦合作用下，應(yīng)力通過改變孔隙介質(zhì)滲透屬性來影響滲流場，滲流場的改變導(dǎo)致滲透力對應(yīng)力場的改變，致使發(fā)生體積應(yīng)變，其本構(gòu)方程為：

Δσij～+αΔp=Hij*(σij，Δεij)

(4)

式中：Δσij～為應(yīng)力增量；Hij*為給定函數(shù)；εij為總應(yīng)變。

(4)相容方程

還需滿足相容方程，即：

εij*=(vi，j+vj，i)/2

(5)

式中：vi，j為介質(zhì)中某點的速度，m/s。

(5)連續(xù)性方程

聯(lián)立平衡方程(2)和本構(gòu)方程(4)可得到流體的連續(xù)性方程：

(1/M)(?p/?t)+(n/s)(?s/?t)=(1/s)(-qi，j+qv)-α(?ε/?t)

(6)

式中：s為節(jié)點飽和度，n為節(jié)點孔隙率，其他變量同前。

1.2 流固耦合數(shù)值實現(xiàn)

在FLAC3D數(shù)值計算中，將模型區(qū)域劃分為四面體，將孔隙水壓力和飽和度作為節(jié)點變量，依賴于流體連續(xù)性方程的有限差分節(jié)點表達式，以力學(xué)平衡作為起始狀態(tài)，通過若干流體步及力學(xué)步相互轉(zhuǎn)換循環(huán)計算，達到系統(tǒng)計算平衡；通過流體步得出孔隙水壓力增量，力學(xué)步得出體積應(yīng)變[12]。從而計算出計算模型整體滲流場及應(yīng)力應(yīng)變場分布。

2 數(shù)值模型

壩基及近壩山體的抬升變形現(xiàn)象實質(zhì)為巖土體的擴容行為。影響此種擴容行為的因素諸多，除現(xiàn)有研究成果[6～11]提到的外因(水庫蓄水導(dǎo)致壩基及近壩山體巖層滲流場變化而產(chǎn)生垂直向上的滲透力分量)及內(nèi)因(壩址區(qū)特定的水文地質(zhì)結(jié)構(gòu)模式)之外，還存在諸如土體力學(xué)及水力學(xué)參數(shù)、壩址區(qū)巖層地下水溫度以及壩體接觸空氣溫度等其他可能影響因素，這些因素都將可能引起壩體及近壩山體不同程度的抬升變形。

依據(jù)應(yīng)力場與滲流場耦合作用理論可知，壩基巖層滲透系數(shù)的改變，將引起壩基巖層整體滲流場的改變，進而改變應(yīng)力場，產(chǎn)生相應(yīng)變形。水庫蓄水，巖層孔隙(裂隙)充水，孔隙(裂隙)水壓力增大，有效應(yīng)力降低，由此產(chǎn)生回彈擴容變形。因此，可初步認為：巖土體變形模量影響擴容變形的大小，而滲透系數(shù)的大小對流固耦合作用效應(yīng)會產(chǎn)生不同程度的影響，進而影響抬升變形。

為驗證以上觀點，進一步揭示內(nèi)因?qū)μ冃蔚挠绊?，本文對壩基巖層變形模量及滲透系數(shù)對抬升變形的作用效應(yīng)展開研究，進一步深入探討抬升變形的形成機制以及變形規(guī)律，獲得水庫蓄水作用下庫區(qū)抬升變形的位移空間分布規(guī)律。以下對壩基及近壩山體的抬升變形進行三維數(shù)值分析，分別研究不同巖體變形模量、不同滲透系數(shù)對水庫樞紐區(qū)抬升位移分布的影響規(guī)律。

2.1 計算模型

由于實際壩址區(qū)巖層分布復(fù)雜，壩基及近壩山體巖層的巖體參數(shù)分布差異較大。根據(jù)蔣中明等[6]的研究成果，壩基及近壩山體抬升變形與水庫樞紐區(qū)的水文地質(zhì)結(jié)構(gòu)模式有重大關(guān)系，本文采用文獻[6]中數(shù)值計算模型。采用江埡水庫[8，9]的簡化地質(zhì)模型，將整體模型分為大壩、灌漿帷幕、滲透系數(shù)相對較低的隔水層和滲透系數(shù)較大的含水層。整體三維計算模型范圍取值如圖1所示，以壩軸線分別取上游、下游800.00 m，左岸、右岸取450.00 m作為四周邊界，底部高程-200.00 m作為底部邊界。上游模擬計算水位高程236.00 m，下游模擬計算水位高程125.00 m。大壩為混凝土重力壩，大壩建基面高程114.00 m，壩頂高程245.00 m，壩基底寬102.00 m。壩基灌漿帷幕最大深度處在相對隔水層中，以起到隔水效果，相對隔水層厚度90.00 m。

圖1 模型地質(zhì)結(jié)構(gòu)圖

圖2 計算模型網(wǎng)格圖

圖2為FLAC3D數(shù)值計算網(wǎng)格。采用四面體網(wǎng)格單元剖分，節(jié)點數(shù)共計16 719個，單元數(shù)共計88 097個。

2.2 初始條件

邊界條件：四周豎向邊界采用水平位移約束邊界，模型底部采用位移固定邊界；上游庫水位以下表面按線性變化施加深度122 m庫水壓力，下游河道水位以下表面按線性變化施加深度11.00 m庫水壓力。水庫庫底表面及下游河道表面采用水頭邊界；模型底部及四周豎直邊界采用不透水邊界。

巖土體本構(gòu)關(guān)系采用摩爾-庫倫本構(gòu)模型，滲流計算服從達西定律?；谖墨I[10]地質(zhì)參數(shù)，本文數(shù)值研究采用的地質(zhì)計算參數(shù)取值見表1。

表1 基本參數(shù)表

2.3 計算方案

2.3.1 變形模量計算方案

以表1參數(shù)作為基本參數(shù)，將含水層與隔水層變形模量分別按0.1、0.5、1.0、1.5、2.0、3.0倍的變換系數(shù)同比縮小或增大且其他參數(shù)保持不變建立計算方案，計算方案見表2。方案1壩基巖體含水層變形模量為1.10 GPa，相對隔水層變形模量為0.96 GPa。顯然，相比實際重力壩壩址區(qū)基巖變形模量參數(shù)，方案1中所采用的參數(shù)明顯偏小，實際勘測選址也不會將低參數(shù)巖層作為優(yōu)選壩址區(qū)，但作為抬升變形規(guī)律定性分析過程，故也將其納入計算方案中，為類似低變形模量壩址巖層的抬升變形分析提供參考。

2.3.2 滲透系數(shù)計算方案

將地質(zhì)巖層分為滲透系數(shù)相差較大的相對隔水層和含水層。建立研究計算方案：第一種，采用表1基本參數(shù)，將相對隔水層的滲透系數(shù)確定不變，變換含水層巖體的滲透系數(shù)，其他力學(xué)及水力學(xué)參數(shù)均不變，考察含水層滲透系數(shù)變化對抬升變形的影響；第二種，將含水層的滲透系數(shù)確定不變，變換相對隔水層巖體的滲透系數(shù)，其他力學(xué)及水力學(xué)參數(shù)均不變，考察相對隔水層滲透系數(shù)變化對抬升變形的影響。計算方案見表3。

表2 不同變形模量計算方案表

2.4 計算步驟

按照以上擬定計算方案參數(shù)，分別賦予模型進行計算。計算步驟如下：首先，施加重力荷載生成初始自重應(yīng)力場，接著，進行大壩一次性施工完成模擬，生成大壩施工完成后應(yīng)力場，位移清除；然后，進行生成初始地下滲流場；最后，考慮流固耦合作用，施加上游蓄水荷載，模擬壩基及近壩山體抬升過程。

表3 不同滲透系數(shù)計算方案表

3 抬升變形空間位移分布規(guī)律

3.1 巖體變形模量對抬升變形的影響

表4為計算模型含水層與相對隔水層在不同變形模量參數(shù)下，最大抬升位移量及其最大抬升位置統(tǒng)計表。由表可知，隨著含水層與相對隔水層變形模量不斷增大，其抬升位移量呈現(xiàn)減小趨勢。方案1中含水層變形模量為1.10 GPa，相對隔水層為0.96 GPa且傾向下游(傾角為38°)，計算最大抬升量為176.90 mm；當(dāng)含水層變形模量增加到33.00 GPa，相對隔水層為28.80 GPa，其最大抬升量僅為1.5 mm，這表明壩基及近壩山體巖層變形模量參數(shù)對抬升變形的作用影響較大。而從最大抬升部位來看，最大抬升位移位置均出現(xiàn)在下游河道近壩山體，基本在壩軸線往下游200 m范圍內(nèi)。當(dāng)變形模量足夠大時，也出現(xiàn)在下游河道較遠范圍，可能原因是巖層變形能力有限，加上壩體重力荷載，滲透力不足以引起近壩河道巖層與近壩山體巖體上抬變形，因其最大抬升部位離壩體較遠。

針對6種不同變形模量參數(shù)計算方案，圖3給出了各方案抬升位移空間分布云圖。由圖可知，壩基及近壩山體巖層的變形模量參數(shù)不同時，其抬升變形量及抬升空間分布亦不同。由于水庫上游庫底水壓力最大，造成最大沉降位移量出現(xiàn)在水庫底部位置。在上游水庫蓄水滲透力的作用下，壩體及下游近壩庫岸、河道部分區(qū)域出現(xiàn)抬升變形，且抬升變形量隨變形模量的增大而減小，而近壩山體的上抬變形很有可能由于邊坡表面傾角的增大造成局部滑坡。而離壩體較遠處的下游近壩山體、河床及高程較高山體依然產(chǎn)生沉降變形，由于不受水庫蓄水導(dǎo)致的滲透力影響，故沒有表現(xiàn)上抬趨勢。由于計算模型中相對隔水層傾向下游(傾角38°)，因而無論是沉降變形還是抬升變形都是關(guān)于河道中心剖面對稱分布。

表4 最大抬升位移統(tǒng)計表

圖3 變形模量敏感計算鉛直位移空間分布圖

圖4給出了方案1至方案6計算模型在河道中心剖面的抬升位移等值線圖。等值線云圖直觀的表達了河道中心剖面抬升位移分布情況。由圖可知，隨著變形模量的增大，河道中心剖面抬升位移減?。欢咸ё冃螀^(qū)域均發(fā)生在相對隔水層的上盤巖層，在抬升變形與沉降變形之間存有一位移變形為零的“分界面”。圖中位移等值線表明，隨著變形模量的增加，分界面逐漸減小，即抬升位移區(qū)域范圍逐漸減小。

圖4 河道中心剖面鉛直位移分布圖

3.2 巖體滲透系數(shù)對抬升變形的影響

3.2.1 含水層滲透系數(shù)變化對抬升變形的影響

圖5為計算模型相對隔水層滲透系數(shù)不變的情況下，含水層滲透系數(shù)變換計算的鉛直位移空間分布圖。由圖可知，左右兩岸鉛直位移整體分布關(guān)于河道中心剖面對稱，因庫水壓力荷載所致，各方案最大沉降位移均出現(xiàn)在水庫庫底。壩體及下游近壩河床、近壩山體出現(xiàn)不同程度抬升變形，最大抬升位移量均出現(xiàn)在下游近壩山體位置，其原因是：下游近壩河道由于承受最大河水壓力，其抬升量比近壩山體抬升量略??；而近壩山體雖然承受一定水位的河水壓力，但是在較大滲透力的作用下出現(xiàn)抬升，兩種效應(yīng)綜合導(dǎo)致最大抬升量位置出現(xiàn)在近壩山體(見表5)。圖中還表明，隨著含水層滲透系數(shù)的減小，其抬升量整體呈先增大后減小趨勢，抬升范圍有向下游移動的趨勢并且不斷增大。因此，在含水層滲透系數(shù)不斷變化的方案中，抬升變形的產(chǎn)生是含水層與隔水層滲透系數(shù)特定“組合”的結(jié)果，這種“組合”作用直接影響抬升變形的大小及范圍。

圖6為在河道中心剖面鉛直等值線位移圖。由圖可知，與上述計算成果規(guī)律類似，在抬升變形與沉降變形之間存有一位移變形為零的“分界面”，并且隨著滲透系數(shù)的減小，其河道下部巖體抬升范圍越大，抬升量整體呈增大趨勢。

圖5 含水層滲透系數(shù)變換計算鉛直位移空間分布圖

圖6 河道中心剖面鉛直位移分布圖

計算方案巖層滲透系數(shù)K/(cm/s)抬升位移概況含水層相對隔水層最大抬升位移量/mm最大抬升位移部位a11.0E-025.0E-063.1壩體壩段中間部位a21.0E-035.0E-069.1壩軸線往下游方向170 m河道近壩山體a31.0E-045.0E-0614.5壩軸線往下游方向247 m河道近壩山體a41.0E-055.0E-0613.2壩軸線往下游方向267 m河道近壩山體

3.2.2 相對隔水層滲透系數(shù)變化對抬升變形的影響

圖7給出了在含水層滲透系數(shù)不變的情況下，相對隔水層滲透系數(shù)逐漸減小的抬升位移空間分布圖。由圖可知，隨著相對隔水層滲透系數(shù)的逐漸減小，抬升位移量并沒有大幅度增加或減小，而是在9.0～10.0 mm范圍內(nèi)，抬升位移分布關(guān)于河道中心剖面對稱，并且抬升范圍基本一致(見圖8河道中心剖面鉛直位移等值線圖)。可見，當(dāng)含水層滲透系數(shù)處于某一定值時，相對隔水層的變化幅度不會太大，而是在某一范圍值內(nèi)變化。

表6為抬升位移統(tǒng)計表。由表可知，最大抬升位移量出現(xiàn)在下游近壩山體，并且與相對隔水層滲透系數(shù)變化不是單一相關(guān)關(guān)系，在下游近壩山體180 m附近變化。

圖7 相對隔水層滲透系數(shù)變換計算鉛直位移空間分布圖

圖8 河道中心剖面鉛直位移分布圖

計算方案巖層滲透系數(shù)K/(cm/s)抬升位移概況含水層相對隔水層最大抬升位移量/mm最大抬升位移部位b11.0E-035.0E-0510.2壩軸線往下游方向154 m河道近壩山體b21.0E-035.0E-069.1壩軸線往下游方向182 m河道近壩山體b31.0E-031.0E-0610.3壩軸線往下游方向170 m河道近壩山體b41.0E-031.0E-089.7壩軸線往下游方向190 m河道近壩山體

4 結(jié) 論

本文基于應(yīng)力場-滲流場耦合理論，采用三維數(shù)值模擬方法，針對巖體特性參數(shù)對抬升變形的影響展開研究，考察變形模量及滲透系數(shù)對抬升變形的敏感性影響規(guī)律，進一步揭示抬升變形的形成機制。得到以下結(jié)論：

(1)壩基巖體變形模量參數(shù)對抬升變形的影響較大，抬升變形量對變形模量大小較為敏感。隨著巖體變形模量的增大，抬升變形量呈減小趨勢；由于受水位差引起的滲透力及下游水位壓力的綜合作用，發(fā)生最大抬升部位均在下游近壩山體，而不是下游河床底部。上抬變形區(qū)域均發(fā)生在相對隔水層的上盤巖層，在抬升變形與沉降變形之間存有一位移變形為零的“分界面”，“分界面”的形成有可能導(dǎo)致巖土體及水工結(jié)構(gòu)物產(chǎn)生剪拉變形。

(2)改變滲透系數(shù)對抬升變形規(guī)律及其空間分布的影響較不敏感。在相對隔水層滲透系數(shù)不變情況下，隨著含水層滲透系數(shù)的減小，其抬升量整體呈增大趨勢，抬升范圍有向下游移動的趨勢并且范圍擴大；在含水層滲透系數(shù)不變情況下，隨著相對隔水層滲透系數(shù)的減小，其抬升量并沒有大幅增加或減小，而是穩(wěn)定在一個范圍值；空間位移分布形狀相似；最大抬升位移部位同樣發(fā)生在下游近壩山體。近壩山體抬升變形過大有可能造成邊坡局部滑坡。