徐文杰

(清華大學(xué),清華大學(xué)水圈科學(xué)與水利工程全國重點(diǎn)實(shí)驗(yàn)室,北京 100084,中國)

0 引 言

滑坡涌浪是近水滑坡災(zāi)害鏈上的一個(gè)重要致災(zāi)方式,一直以來是國內(nèi)外工程地質(zhì)、水利工程領(lǐng)域等研究的難點(diǎn)和熱點(diǎn)。尤其,隨著大規(guī)模水電工程建設(shè)的不斷開展,庫區(qū)滑坡及涌浪災(zāi)害問題將直接影響水域航道、庫區(qū)沿岸居民及工程安全,甚至大壩和電站安全運(yùn)營。

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

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

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

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

圖2 滑坡涌浪過程示意圖

為了進(jìn)一步推動(dòng)滑坡涌浪災(zāi)害的分析方法的發(fā)展,本文系統(tǒng)地歸納總結(jié)了國內(nèi)外現(xiàn)有庫岸邊坡穩(wěn)定研究現(xiàn)狀,以及滑坡涌浪分析的3大類方法,即:物理模型試驗(yàn)法、解析法及數(shù)值模擬法。并對現(xiàn)有研究的成果和不足進(jìn)行了分析,在此基礎(chǔ)上探討了滑坡涌浪分析的關(guān)鍵問題及前沿方法,提出了基于FEM-SPH-DEM耦合的庫岸滑坡涌浪分析方法。

1 水動(dòng)力環(huán)境下庫岸邊坡穩(wěn)定分析

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

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

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

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

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

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

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

2 滑坡涌浪模型試驗(yàn)研究

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

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

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

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

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

4 滑坡涌浪數(shù)值模擬方法

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

4.1 基于單一流體動(dòng)力學(xué)的滑坡涌浪數(shù)值分析

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

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

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

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

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

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

圖6 CoSim中SPH-DEM耦合圖示

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

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

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

滑坡的范圍、變形及破壞過程等特征是涌浪形成的誘發(fā)因素。因此,滑坡這些特征的精細(xì)化模擬與分析是滑坡涌浪分析的關(guān)鍵。DEM雖然在模擬滑坡災(zāi)害動(dòng)力方面有較好的優(yōu)勢,但是由于其細(xì)觀參數(shù)選取方面存在一定的瓶頸,在實(shí)現(xiàn)滑坡體穩(wěn)定性定量分析時(shí)依然有不同的觀點(diǎn)。而相交于DEM而言,基于FEM的強(qiáng)度折減法則在滑坡穩(wěn)定性定量評價(jià)及滑動(dòng)面預(yù)測等方面有了很好的應(yīng)用。

因此,構(gòu)建多種方法的耦合(聯(lián)合)分析(圖8),充分利用不同分析方法的優(yōu)勢,實(shí)現(xiàn)滑坡涌浪的不同階段、不同過程的分析,是建立合理的滑坡涌浪分析和評價(jià)方法的關(guān)鍵。針對具體的滑坡體,基于FEM法開展降雨及庫水位變動(dòng)等水動(dòng)力環(huán)境下庫岸邊坡穩(wěn)定性的演化過程、失穩(wěn)范圍及失穩(wěn)機(jī)制的定量分析和評價(jià); 在此基礎(chǔ)上,基于DEM實(shí)現(xiàn)構(gòu)成滑坡體的顆粒接觸參數(shù)的反演與分析,從而實(shí)現(xiàn)災(zāi)害動(dòng)力學(xué)過程的模擬與預(yù)測; 最終,基于SPH-DEM方法實(shí)現(xiàn)滑坡涌浪、涌浪傳播等鏈生災(zāi)害動(dòng)力學(xué)全過程的模擬與分析。

圖8 FEM-SPH-DEM滑坡涌浪鏈生災(zāi)害動(dòng)力學(xué)分析思路

6 結(jié) 論

隨著世界各國水電工程建設(shè)的不斷發(fā)展,尤其是梯級水電工程群的建設(shè),滑坡涌浪鏈生災(zāi)害動(dòng)力學(xué)過程是影響甚至至于水電工程(群)穩(wěn)定性的重要因素。建立合理的分析與評價(jià)方法,將對于優(yōu)化工程建設(shè)、防災(zāi)減災(zāi)具有重要的理論意義和工程價(jià)值。

水動(dòng)力環(huán)境的劇烈變化是影響庫岸邊坡穩(wěn)定性的重要外部誘發(fā)因素。因此,構(gòu)建水-巖耦合的數(shù)值分析方法,是目前庫岸滑坡穩(wěn)定性演化及評價(jià)的發(fā)展趨勢。

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

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

基于先進(jìn)的數(shù)值計(jì)算分析方法,從滑坡涌浪災(zāi)害全過程的物理力學(xué)機(jī)制出發(fā),構(gòu)建多種方法的耦合(聯(lián)合)分析,充分利用不同分析方法的優(yōu)勢,實(shí)現(xiàn)滑坡涌浪的不同階段、不同過程的分析,對于深入揭示水動(dòng)力環(huán)境下庫岸滑坡動(dòng)力學(xué)、涌浪災(zāi)害機(jī)理及災(zāi)害鏈效應(yīng),同時(shí)對于發(fā)展現(xiàn)代巖土災(zāi)害分析理論方法,及指導(dǎo)滑坡涌浪災(zāi)害防治具有理論和工程意義。