劉東澤, 譚 燕, 盧應發
(湖北工業大學土木建筑與環境學院, 武漢 430068)
2020年11月1日,水利部、國家發展改革委公布,三峽工程完成整體竣工驗收全部程序,標志著三峽工程正式完工,自三峽工程進行大規模蓄水以來,在防洪與發電等方面發揮了巨大作用,但也帶來了一系列問題。其中三峽庫區的滑坡與崩塌問題最為嚴重,其形成原因也十分復雜,大量學者對庫區的滑坡、崩塌形成機理進行研究,肖捷夫等[1-2]、江強強等[3]、朱元甲等[4]通過建立大型滑坡實驗對滑坡進行分析,認為水位的快速降落與降雨是影響滑坡變形的重要因素,閆國強等[5]、劉新榮等[6]、周劍等[7]對典型滑坡成因進行分析,認為強降雨與庫水位升降會使滑坡產生較大變形,王英珍等[8]、夏軍強等[9]對黃河、荊江河崩岸建立理論模型,分析多種因素對河岸的影響,還有大量學者利用數值模擬對滑坡的形成因素進行分析[10-16]。
顯然,在諸多因素中,庫水位變化與降雨對三峽庫區的滑坡崩塌形成有重要影響,因此現建立庫水位降落聯合降雨的耦合數值分析模型,以新田崩塌體為對象,建立二維、三維計算模型,分析崩塌體在不同環境下的變形規律,并建立數值預測模型,預測最危險情況下的變形特征。
新田崩塌體區位于武隆縣江口鎮,地處烏江右岸,崩塌體前緣為烏江主河道,崩塌體全貌見圖1。崩塌體地貌形態為一斜坡,坡體地形起伏較大,崩塌體地形上陡中緩下陡。據野外調查測繪,崩塌體后緣受地形特征及巖土界面控制,呈圈椅狀展開,崩塌體高程275 m,縱長250 m,橫寬400 m,平均厚度15 m,主變形方向157°。崩塌體兩側邊界的劃定主要以地形特征變化及變形特征綜合確定,西側邊界延伸走向170°,以土坎變形及地形特征、沖溝為界,其外圍未見變形現象,而崩塌體區內存在變形。斜坡中上部堆積有厚度較大的崩坡積物,前緣受烏江水不斷沖刷形成高陡的臨空面,為崩塌體變形提供了空間。并且,崩塌體物質結構較為松散,水流易于下滲,下伏為相對隔水的灰巖、頁巖,地下水沿上部崩坡積物與灰巖頁巖接觸面匯集,大大增加了其接觸面的變形性能。

圖1 新田崩塌體全貌
根據勘察資料,選取滑坡體的主剖面為計算剖面(圖2)。滑體的變形模量、凝聚力和摩擦角,類比周邊崩塌體的參數取值,綜合確定新田崩塌體數值計算參數具體見表1。

表1 新田崩塌體有限元計算物理力學參數取值表

圖2 新田崩塌體地質剖面圖
由于崩塌體前緣被庫水淹沒,則模型的水頭邊界為浸沒的崩塌前緣,崩塌體表面處取降雨引起入滲的流量邊界,當降雨強度大于坡面巖土體的入滲速度時,取巖土體的入滲速度值作為邊界流量值;當降雨強度小于坡面巖土體的入滲速度時,取降雨強度值作為邊界流量值,取降雨強度計算。模型底面和兩側為自由滲流邊界,在此由于基巖的滲透性很小,可認為是不透水邊界。各邊界條件的表達如圖3所示。

v為入滲速度,f為產流
2.3.1 二維計算模型
根據新田崩塌體的地質條件和地形地貌特征,利用GEO-SLOPE有限元軟件對新田崩塌體的主剖面進行穩定性計算,采用四邊形單元進行有限元網格剖分,節點數為3 806,單元數為3 775,網格圖如圖4所示。

圖4 二維計算網格模型圖
2.3.2 三維計算模型
為進一步分析崩塌體的變形特征,結合二維計算結果進行三維建模,選取新田崩塌體三維數值計算模型的范圍為:沿烏江水流方向為720 m,垂直烏江水流方向為720 m,模型底面高程為50 m。計算域包含崩塌體、與基巖接觸帶和基巖,整個計算域剖分了38 291個六面體單元,共計43 200個節點,計算模型中的斷層采用實體單元模擬,三維計算模型與網格如圖5所示。

圖5 三維計算網格模型圖
三峽工程建成后,水位呈周期性變化,庫區水位升級變化過程如圖6所示。三峽庫區在汛期來臨前,水位會降至145 m,在遇到洪水時,水位會達到162 m,洪水過后水位又會迅速下降,在蓄水時,水位能達到175 m,然后又會緩慢下降。

圖6 三峽水庫運行水位曲線圖
武隆縣庫區地處長江支流烏江兩岸,烏江庫水位的升降參考長江調度,依據武隆區降雨資料統計分析得到:重現期50 年一遇的汛期日降雨量q全=63.371 mm/d,非汛期日降雨量q枯=45.484 mm/d;重現期20 年一遇的汛期日降雨量q全=59.065 mm/d,非汛期日降雨量q枯=40.571 mm/d。在工況設計時,根據水位調動與降雨資料,設置9種工況對崩塌體的穩定性進行分析,工況組合見表2,同時為分析崩塌體破壞機理,設置5種工況進行三維變形分析,詳細組合見表3。

表2 二維計算工況及荷載組合

表3 三維計算工況及荷載組合
用Morgenstern-Prince法對新田崩塌體主剖面進行穩定性計算,當蓄水時,對崩塌體的穩定性無影響,即工況1條件下的穩定系數為1.094,工況5條件下的穩定系數為1.044。
計算中考慮飽和-非飽和滲流場與應力場耦合的有限元法,在不同的外部環境下得到新田崩塌體的穩定性系數,水位漲落耦合不同強度降雨得到耦合變化曲線如圖7所示。

圖7 不同工況下穩定系數變化圖
在水位緩慢下降過程中,由圖7(a)可知,在僅有水位變化的條件下,穩定系數變化不大,但在水位下降時遭遇3 d的持續降雨,穩定性系數開始急劇降低,主要原因是庫水位降落與降雨引起崩塌體滲流場的改變,同時在降雨早期階段,降雨強度對穩定性影響并不明顯。
在水位降落速度較快時,根據圖7(b)可知,降雨強度越大,早期對穩定性影響較大,但隨著時間的推移,不同強度降雨對穩定性影響逐漸減小,當在50 年一遇的暴雨時,崩塌體在第7 天穩定性系數就降至1.042,而在降雨強度稍弱的20 年一遇暴雨時,則要達到第9 天。
根據圖7(c)可知,在145 m上升至175 m過程中,僅有庫水位上漲時,穩定性基本不變,表明水位上漲對崩塌體影響不大,在遭遇暴雨時,早期穩定性出現了短暫的上升,隨后穩定性開始降低,總體而言,穩定性降低較小。
聯合分析,得到以下結論:工況6為新田崩塌體最危險工況,即水位快速下降與50年一遇暴雨同時出現的情況對崩塌體最為不利,同時在水位上升時可以短暫提高崩塌體穩定性。
在工況1(水庫蓄水至175 m水位)條件下得到孔壓計算結果,由圖8(a)所示的孔壓分布圖可知:崩塌體175 m水位以下區域的孔壓受水位的直接影響。
結合工況2(175 m水位緩降到145 m)與工況4(壩前162 m水位驟降到145 m水位)孔壓計算結果[圖8(b)、圖8(d)],分析可知:在沒有降雨影響時,崩塌體的孔隙水壓力主要受水位下降的影響,在整個崩塌體中,主要受影響的區域為涉水的滑坡前緣,隨著水位的降低,崩塌體前緣的地下滲流場發生較大變化,地下水位也隨之下降,崩塌體前緣的孔隙水壓力消散緩慢,表現為應力滯后現象,崩塌體的中后部區域的孔隙水壓力變化不大,可以忽略不計。

圖8 不同條件下孔壓分布圖
在水位下降時遭遇50年一遇降雨的情況下,得到工況3(175 m水位緩降到145 m+50年一遇降雨)與工況5(壩前162 m水位驟降到145 m+50年一遇降雨)計算結果[圖8(c)、圖8(e)],結果表明:崩塌體前緣地下水分布受庫水位變化和降雨聯合作用影響,崩塌體中后部的水壓力受降雨影響,但對于崩塌整體而言,孔隙水壓力變化不大。
綜合圖8所示5種工況下孔隙水壓力云圖可知,水位下降對涉水庫岸的孔隙水壓力影響較大,在不涉水區域影響不大,降雨主要影響區域為崩塌表層。
水庫蓄水至175 m水位后大部分區域拉、壓應力分布變化較小。根據工況1條件下的主應力場分布云圖(圖9) 可知,拉應力(第一主應力)主要分布在崩塌體前緣和崩塌體后緣;壓應力(第三主應力)主要分布在崩塌體上部、中部和下部(除崩塌體前緣和后緣)的區域,其絕對值的最大值為616 kPa。

圖9 工況1條件下應力場分布圖
在沒有強降雨的條件下,選取工況2與工況4第三主應力云圖[圖10(a)、圖10(c)],水位由175 m緩降到145 m時,壓應力的分布均有所改變,壓應力值的絕對值均有所增大,壓應力(第三主應力)最大值的絕對值由616 kPa增大為625 kPa;壓應力(第三主應力)主要分布在崩塌體上部、中部和下部。庫水位由壩前162 m驟降到145 m時,壓應力場整體變化較大,壓應力(第三主應力)的絕對值最大值為627 kPa;壓應力(第三主應力)主要分布在崩塌體上部、中部和下部。
結合工況3與工況5的應力云圖[圖10(b)、圖10(d)],在庫水位下降時遭遇50 年一遇降雨時,壓應力的分布均有所改變,壓應力值的絕對值均有所減小,庫水位由175 m緩降到145 m又遇到50 年一遇降雨時,壓應力(第三主應力)最大值的絕對值由625.6 kPa減小為625.2 kPa。庫水位由壩前162 m驟降到145 m又遇到50年一遇降雨時,壓應力(第三主應力)最大值的絕對值由627 kPa減小為625 kPa。

圖10 不同工況下第三應力場分布圖
結合孔隙水壓力的變化可知,水位變化與降雨作用會導致滲流場的改變,而滲流場與應力場會相互影響,因此在不同工況下應力場的改變主要是滲流場的改變。
圖11為工況1條件下水平、垂直位移分布圖,在蓄水狀態下崩塌體下部發生較大水平位移,最大位移值為0.005 m,崩塌體中上部水平位移相對較小,主要原因是崩塌體前緣涉水導致崩塌體下部產生較大變形。在崩塌體上部發生較大垂直位移,方向向下,最大位移值為0.032 m,崩塌體中下部垂直位移相對較小。這主要是由于崩塌體中上部較為陡峭,使得崩塌體在自重作用下產生較大的向下垂直位移,崩塌體下部較為平緩,向下垂直位移較小。

圖11 工況1條件下位移分布圖
對比工況2~5的水平位移分布圖(圖12),四個工況下崩塌體中上部均發生較大水平位移,在沒有降雨影響時,水位降落產生的水平位移均比蓄水條件下的位移大,即靜止水位下產生的水平位移最小,主要原因是水位下降改變滲流場,使其產生向外的滲透力,導致崩塌體的滑動力增加,由此表現出更大的水平位移,但水位下降速率對水平位移影響不大。在強降雨條件下,崩塌體的滲透力進一步增大,滑動力也隨著暴雨的持續而增大,導致工況3、5的位移比工況2、4所得水平位移較大,其中工況5產生的水平位移最大,工況5在崩塌體中上部發生較大水平位移,量值為0.014~0.019 m,最大位移值為0.019 m;崩塌體下部水平位移較小,量值為0~0.014 m。

圖12 不同工況下水平位移場分布圖
對比工況2~5的垂直位移分布圖(圖13),崩塌體上部和中部發生較大垂直位移,方向向下,崩塌體下部垂直位移較小,其中庫水位下降疊加強降雨條件下產生的垂直位移最為明顯,在工況5條件下崩塌體位移最大,主要是由于庫水下降引起崩塌體向下滑移,同時在自重和降雨作用下,使崩塌體產生向下的垂直位移增大。上部最大垂直位移值為0.077 m。在沒有降雨影響時,對比工況2、4垂直位移場分布,在工況2崩塌體上部和中部垂直位移量值為0.033~0.042 m,最大位移值為0.043 m,在工況4上部和中部垂直位移量值為0.032~0.048 m,最大位移值為0.048 m,結果表明水位下降速率對崩塌體中上部的垂直位移有影響。

圖13 不同工況下垂直位移場分布圖
取5個工況下主剖面的等效塑性應變分布計算結果如圖14所示,在蓄水狀態下,等效塑性應變較大值集中分布在與基巖接觸帶的中部和上部,等效塑性應變最大值為0.007。在強降雨的影響下,工況3等效塑性應變較大值分布范圍較工況2有所增大。等效塑性應變較大值由0.011增大為0.013。工況5等效塑性應變分布范圍較工況4有所增大。等效塑性應變較大值由0.011增大為0.015。
綜合工況1~工況5的數值計算分析結果,新田崩塌體在工況5(162 m水位驟降到145 m+50 年一遇降雨)條件下,庫水位下降與強降雨的影響造成滑體內滲流場發生很大的變化,使坡體內的滲透力(降雨入滲和坡體內外水位差引起)發生很大的改變,進而引起位移發生很大的變化;相比較其他4個工況,新田崩塌體在工況5(162 m水位驟降到145 m+50年一遇降雨)條件下變形增量最大,變形破壞主要發生在崩塌體后部和前部,后部主要發生大范圍垂直位移。因此,工況5(162 m水位驟降到145 m+50 年一遇降雨)為對新田崩塌體穩定性最不利工況。
由數值計算分析結果可知庫水位由162 m水位驟降至145 m水位遭遇50 年一遇暴雨為對新田崩塌體穩定性最不利工況。提取崩塌體主滑剖面上地表位移監測點XT01、XT02和XT03(圖15)在最不利工況條件下的位移變化及滑帶代表性位置點的等效塑性剪應變等信息(圖16),通過這些信息的數據擬合,建立崩塌體變形破壞的數值預測模型。

圖15 監測點在三維模型中的位置

圖16 各監測點位移變化趨勢
為了獲得崩塌體等效塑性應變的預測模型,在崩塌體與基巖接觸處分別選取的A、B、C三個代表點進行分析,A、B、C三點位置如圖17所示,等效塑性應變變化趨勢如圖18所示。

圖17 代表點位置分布圖

圖18 等效塑性應變變化趨勢
綜合監測位移趨勢與等效塑性應變變化趨勢可知工況5(庫水位由162 m水位驟降至145 m水位+50 年一遇暴雨)條件下,處于崩塌體中后部的地表變形監測點XT03和XT02的位移變化比處于崩塌體前部的地表變形監測點XT01明顯;由于崩塌體較陡,所有監測點的垂直位移要比水平位移變化明顯。
同時崩塌體與基巖接觸處的前中后部三個監測點A點、B點以及C點的等效塑性應變均大于0,說明崩塌體可能發生整體的破壞;C點和B點的等效塑性應變值大于A點,說明崩塌體中后部更容易破壞。
(1)根據二維穩定性分析可知,新田崩塌體最危險工況為水位從162 m水位快速降至145 m,同時遭遇50 年一遇暴雨,在該工況條件下,崩塌體的穩定系數為1.042,崩塌體處于欠穩定狀態。
(2)在庫水位上升的過程中,整體穩定性系數相較于水位下降過程有所提升,即新田崩塌體庫水位上升比水位下降時安全。
(3)根據三維變形破壞分析可知,新田崩塌體在工況5(162 m水位驟降到145 m+50 年一遇降雨)條件下,變形最大和等效塑性剪應變增量最大,工況5為對新田崩塌體穩定性最不利的工況。
(4)根據三維變形破壞分析可知,在工況5(162 m水位驟降到145 m+50 年一遇降雨)條件下,新田崩塌體在崩塌體中后部和前緣分別出現局部的拉剪和壓剪破壞,但未發生整體破壞。其中,崩塌體中后部滑帶較前緣滑帶陡,破壞的范圍較前緣大。
(5)綜合二維、三維計算分析可知,庫水位改變對新田崩塌體整體穩定性影響不大,但會導致涉水區域的局部破壞,強降雨對整體穩定性影響較大,在后期監測中,要重點關注暴雨天氣。