黃 棟，喬建平，張小剛

（1.長江科學院a.工程安全與災(zāi)害防治研究所；b.水利部水工程安全與病害防治工程技術(shù)研究中心；c.國家大壩安全工程技術(shù)研究中心，武漢 430010；2.中國科學院地表過程與山地災(zāi)害重點實驗室，

成都 610041；3.中國科學院水利部成都山地災(zāi)害與環(huán)境研究所，成都 610041）

庫水位變化對庫岸滑坡穩(wěn)定性影響的預(yù)測研究

黃 棟1a，1b，1c，2，3，喬建平2，3，張小剛2，3

（1.長江科學院a.工程安全與災(zāi)害防治研究所；b.水利部水工程安全與病害防治工程技術(shù)研究中心；c.國家大壩安全工程技術(shù)研究中心，武漢 430010；2.中國科學院地表過程與山地災(zāi)害重點實驗室，

成都 610041；3.中國科學院水利部成都山地災(zāi)害與環(huán)境研究所，成都 610041）

以三峽庫區(qū)曬網(wǎng)壩滑坡為研究對象，在充分考慮雙場耦合的條件下，結(jié)合極限平衡分析法，得到了庫水位變化過程中考慮位移增量和孔隙水壓力增量的庫岸滑坡穩(wěn)定性系數(shù)的計算方法，并對庫水位變化作用下曬網(wǎng)壩滑坡的長期穩(wěn)定性進行分析。結(jié)果顯示：隨著時間的變化，穩(wěn)定性系數(shù)逐步趨于一定值1.172，且坡體內(nèi)浸潤線也隨著庫水位的調(diào)度逐步保持穩(wěn)定，趨于不變。

水庫滑坡；滲流－應(yīng)力場耦合；極限平衡法；穩(wěn)定性分析

據(jù)不完全統(tǒng)計，自1982年以來，三峽庫區(qū)發(fā)生滑坡、崩塌、泥石流多達70多次，其中規(guī)模較大的就有40余次。大壩修建后所形成的狹長型水庫于2003年6月正式蓄水，隨著蓄水水位的不斷提升，水庫庫容不斷增大。水庫蓄水后，因水庫季節(jié)性或年際性水位變化，對庫岸進行反復沖刷、浸泡。古、老滑坡將會受到河流地質(zhì)作用和庫水動力作用，使得原已穩(wěn)定的滑坡再度復活。2010年10月26日，三峽工程正式達到175 m正常蓄水位，隨后庫水位將長期在145～175 m之間運行，水位變幅高達30 m。因此，對庫水位長期變化下庫岸滑坡穩(wěn)定性的研究顯得十分重要［1－2］。

目前，水庫運行條件下的滑坡穩(wěn)定性分析已有大量研究成果，理論和應(yīng)用都逐漸趨于成熟。但是，對庫水位長期變化作用下庫岸滑坡穩(wěn)定性的預(yù)測研究鮮有報導。因此，本文以三峽庫區(qū)曬網(wǎng)壩滑坡為研究背景，通過利用滲流－應(yīng)力場耦合的分析方法結(jié)合極限平衡法，對三峽庫區(qū)庫水位長期周期性調(diào)節(jié)作用下曬網(wǎng)壩滑坡穩(wěn)定性進行預(yù)測研究。

1 理論分析

1.1 應(yīng)力－滲流場耦合方程

1.1.1 滲流場有限元方程

對滲流的偏微分方程采用伽遼金的加權(quán)殘余方法，得到滲流方程有限元形式：

式中：［B］為梯度矩陣；［C］為單元水力傳導系數(shù)矩陣；｛H｝為節(jié)點水頭矢量；〈N〉為插值函數(shù)矢量；q為穿越單元邊的單位流量；τ為單元厚度；t為時間；δ為滲流方程瞬態(tài)分析時的存貯項；A為區(qū)域單元的綜合；L為單元邊的綜合。

滲流方程的有限元形式可簡寫為

式中：［K］為單元本征矩陣；［M］為單元質(zhì)量矩陣；｛Q｝為單元水流量。方程是進行瞬態(tài)滲流分析時的通用有限元方程，對于穩(wěn)定狀態(tài)分析，水頭不是時間的函數(shù)，項｛H｝，t消失，方程簡化為

此式是基本滲流方程的簡短有限元形式。

1.1.2 應(yīng)力場有限元方程

采用虛功原理獲得體系平衡方程［3－5］，其表達式為

式中：｛ε*｝為虛位移；｛δ*｝為虛應(yīng)變；｛Δσ｝為系統(tǒng)內(nèi)應(yīng)力；｛F｝為外部點荷載。

建立有限元方程：體模量；［B］為應(yīng)變矩陣；［D］為排水條件下有效的應(yīng)力－應(yīng)變本構(gòu)矩陣；［K］為剛度矩陣；［Ld］為耦合矩陣；｛Δδ｝為位移增量；｛Δuw｝為孔隙水壓力增量向量。

其方程的簡化形式可寫為

1.1.3 耦合方程

從以上推導可以得出用位移增量和孔隙水壓力增量為方程域變量的飽和－非飽和巖土體耦合固結(jié)控制方程，即

由于滲流－應(yīng)力場存在空間上的差異性和在時間上的非穩(wěn)態(tài)性，進行耦合場有限元分析就比一般彈塑性有限元分析要困難得多。在進行各矩陣計算的過程中，要區(qū)別飽和區(qū)和非飽和區(qū)，即要考慮單元與自由面的關(guān)系。如單元位于水下，則要考慮地下水的滲流和靜水壓力等問題；如單元位于水上，處于非飽和狀態(tài)，則可認為單元內(nèi)各點滲透流速為零，不存在滲流力和靜水壓力；另外，在巖土體體積改變時，其體積改變量并不等于單元體中水量的改變量。

1.2 穩(wěn)定性理論

Morgenstern-Price法既滿足力平衡又滿足力矩平衡條件，該法由Morgenstern和Price（1965）［6］以微分方程組形式提出。

圖1 條塊受力分析Fig.1 Forces on the stick

將邊坡劃分為n個條塊，取第j條塊為隔離體（圖1）。該條塊所受到的作用力有：重力Wj；水壓力合力Uj＝γuWj，γu為孔隙水壓力比；外力Qj，與豎直方向成ωj角；滑動面上有效法向力Nj′；調(diào)用的抗剪強度Tj／Fs，F(xiàn)s為安全系數(shù)；條塊間法向作用力Ej和Ej＋1，與底面的垂直距離分別是zj和zj＋1；條塊間的剪切力λfjEj和λfj＋1Ej＋1。于是有

其中：hj和bj為條塊高度和寬度；αj為底面傾角。

求解上述方程組應(yīng)滿足邊界條件：

式中：Ea和Eb為端部條間力；Ma和Mb為力矩。

這樣，式（8）和式（10）組成的方程組可簡化為

其中En＋1和Mn＋1分別稱為剩余力和剩余力矩，可遞推求得。

另外，式（8）和式（9）條間力函數(shù)f（s）為

式中：A為表征條間力函數(shù)形狀的系數(shù)；B（s）為表征峰值點位置的函數(shù)。

其中：s為滑體歸一化后的水平坐標；BBF為峰值點的位置。于是，通過Newton-Raphson法迭代求解方程組（11），便可求得λ和Fs。

2 數(shù)值計算與分析

2.1 模型的建立

為了分析庫水位的變動對曬網(wǎng)壩滑坡體穩(wěn)定性的影響，結(jié)合三峽庫區(qū)曬網(wǎng)壩滑坡勘察報告，采取如圖2所示的計算剖面（曬網(wǎng)壩滑坡剖面Ⅲ）和計算模型。由于計算滲流場的瞬態(tài)變化計算點上不能一次達到平衡，因此，如果計算中某點的計算水頭大于該點坐標，則下一次迭代中定義該點水頭為豎直方向坐標值，直到水頭值與坐標值相符為止。有限元網(wǎng)格劃分為四邊形法，劃分為7個區(qū)域，共剖分為1 544個結(jié)點，1 346個單元。計算模型的應(yīng)力場邊界條件為：水平方向固定，豎直向自由；底邊界豎直向和水平向均固定。

圖2 穩(wěn)定性計算模型Fig.2 Stability calculation model

2.2 參數(shù)的選取及工況

仍采用曬網(wǎng)壩滑坡穩(wěn)定性分析模型，巖土體參數(shù)也與之相同，模型及邊界條件設(shè)置見圖2。模型僅考慮庫水位對曬網(wǎng)壩滑坡的影響，不考慮巖土體物理力學參數(shù)的變化。巖土體物理力學參數(shù)取值如表1。

表1 物理力學參數(shù)Table1 Physical and m echanical param eters

根據(jù)三峽庫區(qū)調(diào)度（圖3），確定模型的計算工況：模擬三峽庫區(qū)5個水文年的調(diào)度過程，其中水位緩慢升、降速度為0.5，0.25 m／d，快速升、降速度為1，1.5 m／d，時間步長2 d一步，共609步，具體見圖4。

圖3 三峽庫區(qū)庫水位調(diào)度圖（長江水利委員會）Fig.3 Water level fluctuation in the Three-Gorges reservoir area（provided by Changjiang Water Resources Comm ission）

2.3 結(jié)果分析

根據(jù)文獻［7－9］最小穩(wěn)定性系數(shù)一般發(fā)生在庫水位下降時。因此，我們只需對5個水文年的庫水位下降時的穩(wěn)定性系數(shù)進行分析就可以得出最小穩(wěn)定性系數(shù)的變化規(guī)律。5個水文年庫水位下降時的穩(wěn)定性系數(shù)如圖5。

圖5 庫水位下降時穩(wěn)定性系數(shù)（5個水文年）Fig.5 Landslide stability factor during water level decrease in the reservoir（in five hydrological years）

圖4 庫水位調(diào)度圖（5個水文年）Fig.4 W ater level fluctuation（in five hydrological years）

對圖5分析可知，最小穩(wěn)定性系數(shù)變化趨勢符合一階指數(shù)衰減趨勢y＝y(tǒng)0＋A1·e－（x－x0）／t1。因此，對最小穩(wěn)定性系數(shù)進行非線性擬合得到曬網(wǎng)壩滑坡最小穩(wěn)定性系數(shù)非線性擬合曲線如圖6。

圖6 最小穩(wěn)定性系數(shù)擬合圖（5個水文年）Fig.6 Fitted curve ofm inimum landslide stability factor（in five hydrological years）

曬網(wǎng)壩滑坡最小穩(wěn)定性系數(shù)一階指數(shù)曲線為

式中：x為時間（d）；y為最小穩(wěn)定性系數(shù)。

對式（14）進行分析可得，當x趨向于無窮大（＋∞）時，最小穩(wěn)定性系數(shù)為定值1.172 45。因此，根據(jù)擬合曲線可知，隨著時間的變化穩(wěn)定性系數(shù)逐步趨于一定值1.172 45。

5個水文年庫水位緩慢下降時的浸潤線研究表明，坡體內(nèi)浸潤線隨著庫水位的調(diào)度逐步保持穩(wěn)定，趨于不變。

3 結(jié) 論

（1）以多孔介質(zhì)滲流理論、Biot固結(jié)方程、質(zhì)量守恒定律等方法為基礎(chǔ)，從流固耦合理論出發(fā)，推導了基于虛功原理的滲流場和應(yīng)力場的有限元空間域及時間域的離散方程，得到了用位移增量和孔隙水壓力增量為方程域變量的飽和－非飽和土體耦合固結(jié)控制方程。

（2）在充分考慮雙場耦合的條件下，結(jié)合極限平衡分析法，對曬網(wǎng)壩滑坡穩(wěn)定性進行了數(shù)值計算，得到了庫水位變化過程中考慮位移增量和孔隙水壓力增量的庫岸滑坡穩(wěn)定性系數(shù)。

（3）對曬網(wǎng)壩滑坡5個水文年作用下的最小穩(wěn)定性系數(shù)進行擬合得到曬網(wǎng)壩滑坡最小穩(wěn)定性系數(shù)非線性擬合曲線公式，公式顯示曬網(wǎng)壩滑坡在蓄水初期2～3 a將會產(chǎn)生變形破壞，隨后將保持穩(wěn)定。

［1］ 喬建平.滑坡減災(zāi)理論與實踐［M］.成都：科學技術(shù)出版社，1997.（QIAO Jian-ping.Theory and Practice of Landslide Disaster Reduction［M］.Chengdu：Science and Technology Press，1997.（in Chinese））

［2］ 中村浩之，王恭先.論水庫滑坡［J］.水土保持通報，1990，10（1）：53－64.（NAKAMURA H，WANG Gongxian.On Reservoir Landslide［J］.Bulletin of Soil and Water Conservation，1990，10（1）：53－64.（in Chinese））

［3］ 毛昶熙.滲流計算分析與控制［M］.北京：中國水利水電出版社，2003.（MAOChang-xi.Seepage Calculation Analysis and Control［M］.Beijing：China Water Power Press，2003.（in Chinese））

［4］ MATSUI T，SAN K C.Finite Element Slope Stability A-nalysis by Shear Strength Reduction Technique［J］.Soils and Foundations，1992，32（1）：59－70.

［5］ KULHAWY F H.Finite Element Analysis of the Behavior of Embankments［D］.California：University of California，1969.

［6］ MORGENSTERN N R，PRICE V E.The Analysis of the Stability of General Slip Surfaces［J］.Geotechnique，1965，15（1）：79－93.

［7］ BERILGENM M.Investigation of Stability of Slopes under Drawdown Conditions［J］.Computers and Geotechnics，2007，34（2）：81－91.

［8］ 劉新喜，夏元友，張顯書，等.庫水位下降對滑坡穩(wěn)定性的影響［J］.巖石力學與工程學報，2005，24（8）：1399－1444.（LIU Xin-xi，XIA Yuan-you，ZHANG Xian-shu，et al.Effects of Drawdown of ReservoirWater Level on Landslide Stability［J］.Chinese Journal of Rock Mechanics and Engineering，2005，24（8）：1399－1444.（in Chinese））

［9］ 廖紅建，盛 謙，高石夯，等.庫水位下降對滑坡體穩(wěn)定性的影響［J］.巖石力學與工程學報，2005，24（19）：3454－3458.（LIAO Hong-jian，SHENG Qian，GAO Shi-ben，et al.Influence of Drawdown of Reservoir Water Level on Landslide Stability［J］.Chinese Journal of Rock Mechanics and Engineering，2005，24（19）：3454－3458.（in Chinese） ）

（編輯：曾小漢）

2.2 錨桿應(yīng)力監(jiān)測

地下廠房共布置6個錨桿應(yīng)力監(jiān)測斷面，每個斷面內(nèi)設(shè)27～34支1點式或3點式錨桿應(yīng)力計對支護錨桿應(yīng)力和巖體應(yīng)力進行監(jiān)測，見圖3。

圖2 地下廠房典型變形監(jiān)測斷面圖Fig.2 Typical section of deformation m onitoring for the underground powerhouse

圖3 地下廠房典型錨桿應(yīng)力監(jiān)測斷面圖Fig.3 Typical section of anchor bar stressmonitoring for the underground powerhouse

2.3 錨索錨固力監(jiān)測

地下廠房（F8斷層范圍內(nèi)）上游側(cè)巖錨吊車梁上方設(shè)置2臺錨索測力計，下游側(cè)巖錨吊車梁上方設(shè)置1臺錨索測力計，對預(yù)應(yīng)力錨索錨固力進行監(jiān)測。

圖4 地下廠房典型鋼筋應(yīng)力監(jiān)測斷面圖Fig.4 Typical section of steel-bar stressmonitoring for the lattice beam of the underground powerhouse

2.4 格構(gòu)梁的鋼筋應(yīng)力監(jiān)測

地下廠房（F8斷層范圍內(nèi)）頂拱設(shè)置3個格構(gòu)梁鋼筋應(yīng)力監(jiān)測斷面，每個斷面內(nèi)布置10支鋼筋計對格構(gòu)梁鋼筋應(yīng)力進行監(jiān)測，見圖4。

3 監(jiān)測成果與圍巖穩(wěn)定性分析

3.1 圍巖變形

地下廠房圍巖表面變形最大值為15.79 mm，發(fā)生在F8斷層上游側(cè)邊墻47.63 m高程處。圍巖表面變形平均值為4.37 mm。巖體質(zhì)量較好部位的圍巖變形均＜2 mm。F8斷層較為破碎的部位，圍巖表面變形基本在5 mm以上。

地下廠房F8斷層典型部位的圍巖變形隨時間變化過程如圖5和圖6所示。地下廠房圍巖變形主要發(fā)生在開挖施工期間。開挖完成后，圍巖進入蠕動變形期，變形呈緩慢增加的趨勢，2009年12月以后地下廠房各部位的圍巖變形趨于穩(wěn)定，圍巖變形已收斂。

圖5 地下廠房圍巖頂拱變形過程線圖Fig.5 History of rock deformation in the arch roof of the underground powerhouse

圖6 地下廠房圍巖邊墻變形過程線圖Fig.6 H istory of rock deformation in the side-wall of the underground powerhouse

圍巖表面及內(nèi)部變形空間分布如圖7所示。樁號廠0＋120和廠房西端等部位的巖體較為完整，變形規(guī)律符合一般巖石開挖變形規(guī)律，變形主要發(fā)生在距圍巖表面2～3 m的松動圈范圍內(nèi)，表面變形最大，沿孔深方向變形越來越小。F8斷層穿越的部位，巖體質(zhì)量稍差，該區(qū)段內(nèi)頂拱部位的圍巖變形發(fā)生在距圍巖表面約5 m的范圍內(nèi)，基本呈現(xiàn)從孔口至孔底逐漸減弱的趨勢。

圖7 地下廠房圍巖變形空間分布圖（2011年6月30日）Fig.7 Spatial distribution of surrounding rockmass deformation of the underground powerhouse（Jun.30，2011）

地下廠房圍巖邊墻表面變形大于頂拱表面變形，這與地下洞室圍巖變形的基本規(guī)律相符。圍巖頂拱變形主要發(fā)生在開挖初期（I—IV層開挖），第IV層開挖完成后，圍巖頂拱變形趨于平穩(wěn)，第V—VIII層開挖對頂拱變形影響較小，見圖8。

圖8 開挖初期頂拱圍巖變形過程線圖Fig.8 H istory of rock deformation in the arch roof during initial stage of the excavation

3.2 支護錨桿應(yīng)力

錨桿最大拉應(yīng)力為150.5 MPa，發(fā)生在廠0＋120斷面上游側(cè)42.0 m高程處。地下廠房平均支護錨桿應(yīng)力為22.84 MPa。經(jīng)統(tǒng)計，地下廠房79.3%的觀測錨桿其應(yīng)力測值小于50.0 MPa，約5.4%的錨桿應(yīng)力位于50.0～100.0 MPa區(qū)間，約10.0%的錨桿應(yīng)力大于100.0 MPa。錨桿應(yīng)力增量主要發(fā)生在地下廠房開挖期間，之后逐漸趨于平穩(wěn)，見圖9。2009年12月后地下廠房內(nèi)的錨桿應(yīng)力總體上穩(wěn)定，僅隨圍巖溫度變化而變化，與溫度變化趨勢相反。錨桿應(yīng)力的變動范圍較小，平均變幅為5 MPa左右。

圖9 地下廠房錨桿應(yīng)力過程線圖Fig.9 History of anchor bar stress in the underground powerhouse

3.3 錨索錨固力

錨索錨固力最大值為1 143.29 kN。地下廠房開挖期間，錨索張力隨圍巖變形有少量增加，增加幅度約為50 kN。3臺錨索測力計均有少許松馳現(xiàn)象，松弛應(yīng)力約為15 kN。目前各臺錨索測力計的測值均已趨于穩(wěn)定，見圖10。

與圍巖變形和錨桿應(yīng)力相比，錨索錨固力變化較小，這主要是因為錨索位于地下廠房拱腳部位，其安裝時間為2006年3月，此時地下廠房頂拱的圍巖變形已基本趨于穩(wěn)定（見圖5）。

圖10 地下廠房錨索錨固力過程線圖Fig.10 History of anchor force in the underground powerhouse

3.4 格構(gòu)梁鋼筋應(yīng)力

地下廠房格構(gòu)梁最大鋼筋拉應(yīng)力63.76 MPa，發(fā)生在廠0＋95樁號頂拱下層鋼筋上；最大鋼筋壓應(yīng)力75.56 MPa，發(fā)生在廠0＋104樁號頂拱與拱腳之間的下層鋼筋上。平均鋼筋拉應(yīng)力為16.02 MPa，平均壓應(yīng)力為25.08 MPa。

格構(gòu)梁鋼筋應(yīng)力變化過程線如圖11所示。由圖中可見，鋼筋應(yīng)力變化主要發(fā)生在地下廠房開挖期間，當頂拱巖體發(fā)生變形時，巖體對格構(gòu)梁拱圈產(chǎn)生壓應(yīng)力，從而導致鋼筋應(yīng)力改變。開挖完成后，圍巖變形趨于穩(wěn)定，鋼筋應(yīng)力也趨于穩(wěn)定，但隨廠房內(nèi)溫度變化而呈正弦規(guī)律變化，在一定范圍內(nèi)波動。

圖11 地下廠房格構(gòu)梁鋼筋應(yīng)力過程線圖Fig.11 H istory of stress on the steel-bar of lattice beam of the underground powerhouse

從鋼筋應(yīng)力過程線來看，目前格構(gòu)梁鋼筋應(yīng)力主要受洞室氣溫的影響。當氣溫升高時，格構(gòu)梁混凝土溫度升高快于巖體溫度，混凝土膨脹，導致鋼筋壓應(yīng)力增加，拉應(yīng)力減少；當氣溫降低時，格構(gòu)梁混凝土溫度降低快于巖體溫度，混凝土收縮，導致鋼筋壓應(yīng)力減小，拉應(yīng)力增加。

3.5 圍巖穩(wěn)定性評價

從圍巖變形、支護錨桿應(yīng)力、錨索錨固力及格構(gòu)梁鋼筋應(yīng)力監(jiān)測的成果來看，地下廠房圍巖各項監(jiān)測效應(yīng)量正常，均小于設(shè)計指標，且進入電站試運行期后一直處于比較穩(wěn)定的狀態(tài)。

變形、錨桿應(yīng)力及鋼筋應(yīng)力的變化過程也反映出了三者之間的相互影響關(guān)系：當圍巖變形增加時，

Prediction of Landslide Stability of Reservoir Bank under the Influence of W ater Level Fluctuation

HUANG Dong1，2，3，QIAO Jian-ping2，3，ZHANG Xiao-gang2，3
（1.Yangtze River Scientific Research Institute，Wuhan 430010，China；2.Key Laboratory of Mountain Hazards and Surface Process，Chengdu 610041，China；3.Institute of Mountain Hazards and Environment，Chinese Academy of Sciences，Chengdu 610041，China）

Landslide stability affected by long-term reservoir water level fluctuation is predicted by finite equilibrium method considering the coupling between seepage field and stress field.Shaiwangba landslide in the Three Gorges Reservoir area is taken as a case study.Through numerical calculation，the stability factor considering the displacement incremental variable and pore water pressure incremental variable is acquired.By analyzing the stability factor during the decrease of reservoirwater level in 5 hydrological years，the nonlinear fitted curve ofminimum stability factor is obtained.Results show that the stability factor tends to be steady at1.172，and the landslide phreatic line tends to be constant with the variation of reservoir water level.It’s concluded that deformation will take place at the Shaiwangba landslide during the first2-3 years after the impoundment，and will remain stable in subsequence.

reservoir landslide；coupling between seepage field and stress filed；limit equilibrium method；stability analysis形進行監(jiān)測，見圖2。

TU457

A

1001－5485（2012）09－0103－04

10.3969／j.issn.1001－5485.2012.09.024

2012－07－05

中國科學院知識創(chuàng)新工程重要方向項目（KZCX2－YW－Q03－5）；國家自然科學基金（41001007）；國家科技支撐計劃課題（2012BAK10B04）

黃 棟（1980－），男，湖北黃石人，助理研究員，博士，從事滑坡機理研究，（電話）027－82829878（電子信箱）hd1017＠sohu.com。