庫岸滑坡涌浪鏈生災害動力學研究進展*

徐文杰

(清華大學,清華大學水圈科學與水利工程全國重點實驗室,北京 100084,中國)

0 引 言

滑坡涌浪是近水滑坡災害鏈上的一個重要致災方式,一直以來是國內外工程地質、水利工程領域等研究的難點和熱點。尤其,隨著大規模水電工程建設的不斷開展,庫區滑坡及涌浪災害問題將直接影響水域航道、庫區沿岸居民及工程安全,甚至大壩和電站安全運營。

滑坡涌浪災害在國內外時有發生。1958年美國Lituya海灣發生的大型滑坡(3060×104m3,巖質)(圖1),涌浪在對岸爬升達524m,波速達到150～200km·h-1(Fritz et al.,2009); 1959年意大利的龐特塞拱壩庫區發生滑坡(300×104m3,巖質)形成20m高的涌浪,使壩體漫頂; 1963年意大利瓦伊昂水庫左岸滑坡(近3×109m3,巖質)形成250m高的涌浪,震驚整個工程技術界(鐘立勛,1993); 2017年6月格陵蘭島海岸的Nuugaatsiaq村發生大型滑坡(5800×104m3,巖質),形成100余米高的涌浪(Quirin,2017)。

圖1 1958年美國Lituya海灣的滑坡涌浪(Fritz et al.,2009)

我國隨著水電工程建設的飛速發展,因庫區滑坡涌浪而帶來的災害損失更是屢見不鮮。1961年2月柘溪水庫蓄水期間,近壩庫區右岸發生滑坡(165×104m3)形成21m高的涌浪,3.6m高的涌浪翻過正在施工的壩體(杜伯輝,2006); 2007年6月水布埡庫區的大堰塘滑坡(300×104m3)涌浪在對岸爬坡高達50m 左右,下游20.8km 大壩處涌浪爬坡仍高4m左右(殷坤龍等,2008)。2009年7月小灣電站庫區壩上游的荒田滑坡(300×104m3),造成30余米高的涌浪(吳佳壕,2011)。此外,據統計1982～2001年期間三峽庫區發生地質災害70余次,其中規模較大的有40余次(國土資源部,2001),如:1982年7月的雞扒子滑坡(1500×104m3),將滑坡前緣長江岸坡向江中推進50余米,并形成巨大涌浪(李玉生,1984); 1985年6月的新灘滑坡(3000×104m3),形成涌浪在對岸最大爬坡高度達近50m(駱培云,1988); 2003年的千將坪滑坡(2400×104m3),形成近30m高的涌浪(王治華等,2003); 2008年龔家方滑坡(38×104m3),形成涌浪在對岸最大爬坡達13余米,影響范圍超過10km(Huang et al.,2012); 2015年6月巫山縣大寧河江東寺北岸發生滑坡(2.4×104m3)涌入庫區形成6m高涌浪(Huang et al.,2016)。

庫岸邊坡的失穩及涌浪是一個涉及滑坡動力學、巖土力學及流體力學等多學科的鏈生災害動力學問題,也是一個復雜的“流-固”耦合問題(圖2):一方面,庫水位波動、降雨等水力條件是影響庫區滑坡穩定性及災變的重要因素; 另一方面,形成涌浪的過程實質上是高速運動的滑坡體(固體)與水體(流體)相互作用的復雜過程; 其次,涌浪在傳播過程中遇到岸坡、船舶及大壩等工程結構時產生的沖擊、爬坡等災害效應也是“流體”與“固體”的相互耦合過程。隨著一些災難性滑坡涌浪事故的發生,庫岸滑坡涌浪問題的研究逐漸被國內外眾多學者所關注,并紛紛從不同角度開展了相關研究工作。

圖2 滑坡涌浪過程示意圖

為了進一步推動滑坡涌浪災害的分析方法的發展,本文系統地歸納總結了國內外現有庫岸邊坡穩定研究現狀,以及滑坡涌浪分析的3大類方法,即:物理模型試驗法、解析法及數值模擬法。并對現有研究的成果和不足進行了分析,在此基礎上探討了滑坡涌浪分析的關鍵問題及前沿方法,提出了基于FEM-SPH-DEM耦合的庫岸滑坡涌浪分析方法。

1 水動力環境下庫岸邊坡穩定分析

庫岸邊坡的失穩破壞是誘發涌浪的觸發因素。降雨、庫水位變動等水動力環境是誘發庫岸邊坡失穩的一個的重要外部因素(中村浩之等,1990; 薛果夫,等1998; Fiorillo et al.,2004; Mikos et al.,2004; 閔弘等,2004; Chang et al.,2005; Crosta et al.,2006)。

在降雨誘發土石混合體滑坡災方面,國內外眾多學者在該領域開展了大量的研究工作。胡明鑒等人(2002)，通過大型野外人工降雨試驗模擬研究了土石混合體斜坡暴雨、滑坡及泥石流的共生關系。俞伯汀等(2006)采用室內物理模型試驗研究了土石混合體邊坡內存在的地下水管道排泄系統形成過程及對斜坡穩定的影響。周中(2006)采用室內及現場試驗方法也研究了土石混合體的滲透破壞特性。Sassa(1998)在對日本發生的泥石流研究的基礎上,提出滑體的運動使得滑面處顆粒的破碎,從而導致滑面處土體液化,進而造成滑體的高速、遠程滑動現象。Wang et al.(2003)通過室內環剪試驗揭示了滑坡觸發泥石流演變模式及形成機理,指出了孔隙水壓力的產生是滑坡突變失穩的動因。Chae et al.(2012)指出降雨引起的滑坡體的體積含水率變化比孔隙水壓力波動更加明顯,提出將邊坡的體積含水率梯度作為滑坡預測的失穩閥值。汪丁建等(2016)基于Green-Ampt(GA)模型,同時考慮濕潤鋒以飽和帶滲流作用,推導了一種新的滑坡降雨入滲函數,彌補了GA 模型只適用于無限長坡的不足。

在庫水位變動作用下滑坡災害研究方面,王士天等(1997)通過大量的實際資料研究指出水庫滑坡有兩種:一種是庫水位達到敏感水位后滑體內孔隙水壓力分布達到新的平衡過程中產生的滑坡; 一種是發生在庫水位消落時的滑坡。Luciano(2007)研究表明水庫蓄水引起的孔隙水壓力變化控制著邊坡位移與力學變化行為。賀可強等(2017)針對庫水漲落誘發滑坡致災機理,提出了以庫水位下降引起的邊坡下滑動力增量作為邊坡的加載動力參數。胡修文等(2005)、Timpong et al.(2007)、賈官偉等(2009)、趙代鵬等(2013)、黃佩倫等(2020)、肖捷夫等(2020)分別采用模型試驗的方法,分析了庫水位變動作用下坡體內地下水位和有效應力變化情況,在此基礎上對滑坡失穩機理進行了探討。柴軍瑞等(2004)、鄭穎人等(2004)、劉新喜等(2005)及汪斌等(2007)分別采用有限元數值計算的方法系統研究了庫水位變動環境下坡體內地下水位變化及對邊坡穩定的影響。

水庫蓄水后庫岸巖土體的水-巖相互作用是一個復雜的物理-力學-化學耦合過程,將導致其物理力學狀態的劣化。尤其水庫長期運行過程中,消落帶范圍內的巖土體,一方面受強烈水動力作用的影響,導致內部細小礦物顆粒不斷帶走、掏空; 另一方面,巖體在強烈干-濕循環作用下,強度不斷降低。這兩種作用相互促進,加速巖土體的劣化發展。徐文杰等(2006)通過對天然及浸水后土石混合體的野外試驗研究表明,浸水后由于細粒被帶走導致黏聚力降低而摩擦角略有提升。此外,對水-巖作用下砂巖(劉新榮等,2008; 鄧華鋒等,2015)、千枚巖(崔凱等,2019)及紅層軟巖(謝小帥等,2019)等劣化規律進行了實驗研究。殷躍平等(2021)通過調查研究發現三峽庫區水位消落帶巖體劣化松動使得部分岸坡加速朝不穩定方向演化,給溶蝕巖體岸坡帶來了工程災變效應問題; 并通過對三峽庫區的灰巖浸泡-風干循環后力學參數研究表明,平均每循環下降率約為0.7%,隨著循環次數的增加巖樣強度下降可用指數函數形式擬合(黃波林等,2019)。

徐文杰等(2009)基于有限元及強度折減法對庫水位驟降過程中的流-固耦合及相應的穩定性變化特征進行了系統研究,表明庫水位作用下庫岸邊坡的穩定性演化是水動力環境與邊坡巖土體強度相互博弈的過程:庫水位驟降對邊坡穩定性影響較大,且其穩定系數下降隨著驟降幅度的增大而增大; 而庫水位上升過程中,邊坡的穩定性有所下降,當上升至某一臨界水位后邊坡穩定系數達到最低值,而后隨著庫水位的上升其穩定系數有所回彈(圖3)。當然這種演化規律還取決于具體的邊坡工程地質結構特征及與庫水位的空間關系。

圖3 邊坡穩定系數隨蓄水高程的變化曲線(徐文杰等,2009)

在近幾年來,隨著我國大型水電工程的建設越來越多的學者開始關注降雨與庫水位共同作用下滑坡災害問題。李曉等(2004)針對降雨和庫水位變動聯合作用下滑坡地下水動力場計算問題,提出了一種考慮降雨與庫水位逐日變化的庫岸地下非穩定滲流場的計算模式。羅先啟等(2005)、李卓等(2017)開展了降雨及庫水位共同作用下滑坡大型模型試驗研究,系統分析了降雨及庫水位共同作用庫岸滑坡的災害機理。在數值計算分析方面,劉新喜(2003)采用自編計算機程序開展了庫水位驟降與降雨耦合作用下滑坡體穩定性計算; 連志鵬等(2011)、鄭慧等(2012)、周永強等(2014)、王力(2014)及方四發(2016)等基于飽和-非飽和滲流理論分析了在降雨和庫水位變化的同時作用下,滑坡體內地下水滲流場及斜坡穩定系數的變化; 朱鵬譜(2008)、張桂榮等(2011)基于極限平衡法,對不同降雨強度及庫水位變動條件下的滑坡進行了穩定性分析,以研究最不利滑坡穩定的工況。

2 滑坡涌浪模型試驗研究

物理模型試驗是滑坡涌浪研究的重要手段,以揭示滑坡涌浪災害機理與規律,同時在此基礎上提出相應的解析計算公式。Wiegel et al.(1970)利用水動力模型試驗研究了滑坡涌浪與滑塊寬度、弗勞德數、水深及滑體速度等之間的關系。Kamphis et al.(1971)通過試驗的方法指出涌浪高度主要取決于滑體的體積及滑體的弗勞德數,并建立了相關關系。Enet et al.(2003),Fritz et al.(2003)、Altaie-Ashtiani et al.(2008)、Marcello et al.(2009)及Heller et al.(2015,2016)也分別采用不同形式的室內模型試驗研究了影響滑坡涌浪高度的因素、涌浪傳播特征及機制等; 尤其Ataie-Ashtiani et al.(2008)通過120組模型試驗分別從坡角、水深、入水速度、滑體幾何形狀(8種)、滑體密度及變形特性等影響因素著手開展研究,提出了涌浪波幅與周期的經驗公式。Heller et al.(2015)基于滑坡涌浪三維模型試驗研究了水體幾何特性對近場及遠場涌浪的影響,并提出了相應的經驗公式,同時基于SPH數值計算揭示了滑坡滑動及涌浪的動力學特性(Heller et al.,2016)。為研究松散滑坡涌浪特性,Ball(1970)、Davidson(1975)等率先開展了相關研究工作; Heller et al.(2011)通過200多組矩形水槽模型試驗,研究了靜水深度、沖擊速度、沖擊角度、滑坡厚度、體積、密度及粒度組成等對涌浪的影響,建立了滑坡涌浪波幅及周期與沖擊波參數的半經驗關系。Huber et al.(1997)開展了1000組2D和3D的散體顆粒涌浪模型試驗,指出涌浪波高主要取決于滑體體積與寬度之比。Fritz et al.(2004)和Mohammed(2010)開展了松散滑坡的2D和3D涌浪大型模型試驗,并將粒子圖像測速(PIV)技術、多重聲學傳感測量技術、激光測距儀、波高儀、高速攝像機等應用于試驗全過程量測,系統研究了滑坡體體積形狀、沖擊速度等對涌浪的影響,并根據試驗結果提出了首波波幅的衰減及波長和周期方程。Mulligan et al.(2017)在二維模型試驗的基礎上研究了散體滑坡入水過程中滑坡體與水體間的動量傳遞,并提出了相應的用于預測近場涌浪波幅的半經驗公式。Francesco et al.(2017,2021)研發了用于模擬高速散體滑坡涌浪的室內模型試驗裝置(圖4),并首次采用計算機視覺、PIV等技術實現涌浪的三維量測,通過40多個模型試驗提出了滑坡體參數與涌浪首波波幅、傳播特征及能量的關系; 通過對滑坡涌浪過程中能力轉換分析表明,52%消耗于滑坡體與滑動面的摩擦,42%消耗于湍流等過程,而只有6%傳遞給形成的涌浪(Francesco et al.,2021)。

圖4 Francesco et al.(2021)采用的滑坡涌浪模型試驗裝置圖

伴隨著大型水電工程建設,我國學者也開展了系列滑坡涌浪模型試驗研究。如,殷坤龍等(2012)采用室內大型模型實驗,研究了滑坡規模、入水速度、滑動面傾角、水深、岸坡坡角等因素對涌浪高度的影響; Wang et al.(2017),Huang et al.(2020)通過室內模型試驗對滑坡涌浪形成機制和傳播衰減效應開展了系統的研究; 韓林峰等(2018)基于淺水區和深水區的滑坡涌浪模型實驗,基于動量守恒原理,提出了最大近場撥付的預測模型。此外,針對水電工程庫區遇到的一些滑坡,也開展了針對具體庫區滑坡的涌浪模型試驗研究(如,小浪底水利樞紐(葉耀琪,1982)、李家峽電站(陶孝銓,1994)、三峽庫區(黃波林等,2013; 何慶,2014)、古水電站(丁軍浩等,2017)等),以研究庫區滑坡可能產生的涌浪災害。

3 基于解析法的滑坡涌浪研究

在滑坡涌浪早期的預測分析中主要為解析法。Kennard(1949)從流體力學的連續方程和運動方程出發提出了運動物體產生的表面波的解析工程,Noda(1970)在此基礎上考慮了滑坡體垂直下落和水平推動兩種極端情況,給出了運動物體引起涌浪高度與水深及運動速度的簡單關系式; 后續許多學者又基于Noda法進行了改進和發展,包括我國的潘家錚(1980)和陳學德(1984)等。Slingerland et al.(1979)、Frite et al.(2003)對滑坡涌浪及傳播函數等進行了探討。此外,基于實際工程案例及模型試驗成果,許多學者也提出了相應的經驗公式法,如:Ataie-Ashtiani et al.(2007)通過大量滑坡涌浪模型實驗數據,提出了大型滑坡涌浪波高的經驗關系; Wiegel et al.(1970)基于模型試驗提出了滑坡涌浪與滑體寬度、水深、Froude數及滑體速度等因素的關系; Valentin(2009)在總結部分經驗公式的基礎上,提出了不同類型滑坡涌浪的最大浪高、傳播浪浪高及爬高等的計算公式。黃種為等(1983)根據3個工程的模型試驗,初步探討了滑坡涌浪的計算方法,并指出滑速和滑體體積是影響滑坡涌浪高度的主要因素,并提出了較為簡單的估算方法; 汪洋等(2004)從明渠非恒定流出發,提出了滑坡失事點涌浪的衰減變化規律,并據此進行涌浪高度的疊加以得到滑坡失事點最終涌浪高度的方法。黃波林(2014)將水波動力學模型應用于滑坡涌浪災害模擬,并將其應用于三峽庫區龔家方滑坡涌浪災害模擬取得了較好的效果。

4 滑坡涌浪數值模擬方法

隨著數值計算技術的發展,在滑坡涌浪模擬方面取得了飛速的發展,尤其近幾年來許多學者開展了相關研究,大致可以分為基于單一流體動力學和基于流-固耦合的方法兩大類。

4.1 基于單一流體動力學的滑坡涌浪數值分析

基于單一流體力學計算方法,將滑體視為剛形體或非牛頓流體,從而大大簡化了滑坡涌浪數值分析方法。這方面的主要有:Davidson et al.(1975)及 Raney et al.(1976)利用有限差分法和實驗模擬計算對滑坡涌浪問題進行了研究,得出滑體最大速度與涌浪高度的關系; Pastor et al.(2009)采用歐拉算法將滑體和水分別視為兩種流體,對滑坡的涌浪問題進行了分析;Serrano-Pacheco et al.(2009),Basu et al.(2009),Abadie et al.(2010),Bosa et al.(2011)采用有限體積的方法對滑坡涌浪問題進行了模擬,Basu et al.(2009)和Bosa et al.(2011)分別將其應用于Lituya滑坡及瓦伊昂庫區滑坡涌浪災害; Cremonesi et al.(2011)采用拉格朗日有限元方法對滑坡涌浪問題進行了模擬;Zhao et al.(2016,2017)將滑坡體簡化為非牛頓流,基于有限元法提出了一種改進的水平集守恒法用于三維滑坡涌浪與傳播過程的三相流模擬(滑坡體、水與空氣)。Xiao et al.(2015)提出了一種用于模擬流體和類流體運動的Tsunami Squares方法,將其應用于三峽庫區龔家方滑坡運動、涌浪及傳播全過程的模擬。

近幾年來,光滑粒子流體動力法(SPH)在流體分析領域有了廣泛的應用和發展。Altaie-Ashtiani et al.(2008a,2008b)、Shi et al.(2016)、Shahab et al.(2017)將數值試驗與物理模型試驗進行了對比驗證取得了較好的效果,Vacondio et al.(2013)將其應用于瓦伊昂滑坡涌浪產生、傳播及沖擊壩體潰壩過程研究(圖5)。

圖5 Vacondio et al.(2013)基于SPH方法得到的Vajont滑坡涌浪速度(50s)

4.2 基于流-固耦合的滑坡涌浪數值分析

滑坡涌浪是一個復雜的流-固耦合過程。為了更好地描述滑坡體(固體)與水體(流體)間的相互作用,充分利用不同的計算方法各自的優勢,基于固體力學與流體力學耦合的計算分析方法成為一個新的研究方向。徐文杰在(2012)探討了耦合的歐拉與拉格朗日算法(CEL)在滑坡涌浪研究應用中的可行性。離散元法(DEM)、大變形不連續分析(DDA)在滑坡動力學模擬中有強大的優勢,并取得了很好成果。Mikola et al.(2014)、Wang et al.(2016)采用DDA及Tan et al.(2017)采用DEM來模擬滑坡體,并與SPH(模擬水體)進行耦合,模擬滑體入水涌浪,并與物理實驗進行了對比驗證,取得了較好的效果。此外,Zhao T et al.(2016)采用CFD-DEM耦合的方法分析了Vajont 滑坡涌浪過程。

Xu et al.(2021)研發了基于GPU加速并行的SPH-DEM耦合模擬程序CoSim-SPHDEM(圖6),并通過與Fahad(2010)的模型實驗對比,系統地驗證了SPH-DEM耦合法在滑坡涌浪分析中的可行性,并對體積相同的散體滑坡(Sg)、矩形塊體(Sc)、前緣為45°的楔形體(Sw)、兩個塊體(St)及4個塊體(Sf)等5種模型開展了數值實驗分析,研究了滑坡體的形態、類型等對滑坡涌浪過程中的能量轉化機制及涌浪特征的影響(圖7)。研究表明:滑坡體產生的涌浪首波傳播速度為Sf>St>Sc>Sw>Sg; 形成的涌浪首波波幅,為Sf>St>Sc≈Sw>Sg。

圖6 CoSim中SPH-DEM耦合圖示

圖7 相同體積不同類型滑坡涌浪分析(Xu et al.,2021)

5 基于FEM-DEM-SPH耦合的滑坡涌浪分析

縱觀目前對滑坡涌浪的分析方法可以看出,數值計算尤其是基于流體-固體耦合的數值計算方法(如,SPH-DEM,CFD-DEM等)可以較好的模擬滑坡體的變形破裂過程、又可以實現滑坡體與水體的強烈的非線性耦合作用,進而實現涌浪的流體動力學模擬,從而為滑坡涌浪分析提供了強有力的支撐。與基于有限體積(FVM)的流體動力學(CFD)數值計算方法相比,基于粒子的SPH數值方法可以較好地實現滑坡沖擊水體及涌浪過程中的“碎波”。

滑坡的范圍、變形及破壞過程等特征是涌浪形成的誘發因素。因此,滑坡這些特征的精細化模擬與分析是滑坡涌浪分析的關鍵。DEM雖然在模擬滑坡災害動力方面有較好的優勢,但是由于其細觀參數選取方面存在一定的瓶頸,在實現滑坡體穩定性定量分析時依然有不同的觀點。而相交于DEM而言,基于FEM的強度折減法則在滑坡穩定性定量評價及滑動面預測等方面有了很好的應用。

因此,構建多種方法的耦合(聯合)分析(圖8),充分利用不同分析方法的優勢,實現滑坡涌浪的不同階段、不同過程的分析,是建立合理的滑坡涌浪分析和評價方法的關鍵。針對具體的滑坡體,基于FEM法開展降雨及庫水位變動等水動力環境下庫岸邊坡穩定性的演化過程、失穩范圍及失穩機制的定量分析和評價; 在此基礎上,基于DEM實現構成滑坡體的顆粒接觸參數的反演與分析,從而實現災害動力學過程的模擬與預測; 最終,基于SPH-DEM方法實現滑坡涌浪、涌浪傳播等鏈生災害動力學全過程的模擬與分析。

圖8 FEM-SPH-DEM滑坡涌浪鏈生災害動力學分析思路

6 結 論

隨著世界各國水電工程建設的不斷發展,尤其是梯級水電工程群的建設,滑坡涌浪鏈生災害動力學過程是影響甚至至于水電工程(群)穩定性的重要因素。建立合理的分析與評價方法,將對于優化工程建設、防災減災具有重要的理論意義和工程價值。

水動力環境的劇烈變化是影響庫岸邊坡穩定性的重要外部誘發因素。因此,構建水-巖耦合的數值分析方法,是目前庫岸滑坡穩定性演化及評價的發展趨勢。

在滑坡涌浪災害方面的研究:解析法通常將滑坡體視為一個剛性體或剛性條塊集合,通過研究剛體質心的運動,將其視為移動的邊界,極大地忽略了滑體與水體間的強耦合作用; 物理模型試驗的方法固然能夠在解決滑坡涌浪問題方面具有一定的優勢,但是由于其試驗費用較高,而且模型的縮尺效應在一定程度上也會影響試驗的結果; 數值計算方法為深入研究滑坡涌浪問題開辟了新的方法,尤其是近些年來隨著計算機技術及數值算法的不斷發展,在滑坡涌浪災害分析方面體現了強大的優勢,當然由于涉及巖土體-水耦合、多相流及流體自由面等復雜的數值計算問題也給數值計算帶來了新的挑戰!

滑坡體的運動狀態對其涌浪起決定性作用,滑坡變形破壞機理是研究其運動狀態的重要基礎?；麦w的漸進破壞和啟動取決于自身的地質條件及外動力誘發因素; 運動過程滑坡將不斷破裂、解體,這種結構的改變進而又影響著外動力對其運動狀態的影響?；碌臑暮^程也是滑坡與外動力(如,降雨、地下水、地表水等)的強耦合過程。傳統的滑坡涌浪災害研究中,將滑坡體視為剛性體或非牛頓流體,忽略了滑坡災害機理,甚至忽略了滑坡體與庫水間的耦合作用,從而給滑坡體災害動力學特性研究帶來很大誤差甚至錯誤。運用滑坡災害動力學、流體力學及流-固耦合理論,從災害全過程角度揭示水動力環境下滑坡體漸進破壞、啟動、高速運動、到涌浪形成、傳播與爬坡災害鏈效應,是進行滑坡涌浪災害研究的關鍵。

基于先進的數值計算分析方法,從滑坡涌浪災害全過程的物理力學機制出發,構建多種方法的耦合(聯合)分析,充分利用不同分析方法的優勢,實現滑坡涌浪的不同階段、不同過程的分析,對于深入揭示水動力環境下庫岸滑坡動力學、涌浪災害機理及災害鏈效應,同時對于發展現代巖土災害分析理論方法,及指導滑坡涌浪災害防治具有理論和工程意義。