楊先艾，曾樹元，王 瑞，張 斌

（1.貴州烏江水電開發(fā)有限責(zé)任公司索風(fēng)營發(fā)電廠,貴州 貴陽 550215；2.中國電建集團(tuán)貴陽勘測設(shè)計(jì)研究院有限公司,貴州 貴陽 550081 ）

1 前言

降雨、水位變動(dòng)等作用通常會影響到河道或水庫岸坡的穩(wěn)定性,一旦失穩(wěn)將與水體發(fā)生強(qiáng)烈的相互作用及能量傳遞形成涌浪,并沿著河道傳播,一方面在對岸形成較大的涌浪爬坡,越岸及沖刷[1],另一方面,涌浪向上下游傳播,對航道航行構(gòu)成巨大威脅,對沿岸的水工建筑物、居民區(qū)造成破壞,對下游大壩產(chǎn)生巨大沖擊,甚至發(fā)生漫頂、沖毀等災(zāi)害,進(jìn)而造成水庫下游產(chǎn)生洪水災(zāi)害。與此同時(shí),滑坡體入水也會造成水庫有效庫容的降低,增加蓄水風(fēng)險(xiǎn)[2]。

目前滑坡涌浪災(zāi)害的主要研究手段有物理模型試驗(yàn),解析求解和數(shù)值模擬等。其中物理模型試驗(yàn),依賴于高精度的儀器設(shè)備,同時(shí)考慮到試驗(yàn)中關(guān)鍵因素如水深、滑坡體入水速度等的采集,以及相似準(zhǔn)則的影響,試驗(yàn)結(jié)果推廣到實(shí)際尺度時(shí)誤差較大；解析法基于一定的模型抽象及條件假設(shè),計(jì)算簡便,計(jì)算成本低,但在實(shí)際工程應(yīng)用中會涉及較多的經(jīng)驗(yàn)估算,且在解析可解的條件假設(shè)下實(shí)際偏差較大[3]。隨著數(shù)值技術(shù)方法的不斷發(fā)展,已逐漸成為重大地質(zhì)災(zāi)害研究、分析和評估的重要手段。針對滑坡涌浪問題,目前眾多國內(nèi)外學(xué)者已經(jīng)采用CFD、SPH-DEM耦合方法[4]、一維淺水方程、二維淺水方程[5]等數(shù)值方法開展了相關(guān)研究。其中二維非線性淺水方程在解決滑坡涌浪問題中具有表達(dá)方式簡單,配套數(shù)值求解格式精度高,計(jì)算規(guī)模低,計(jì)算效率高等特點(diǎn),在大范圍流域模擬中被廣泛采用。

本文基于二維淺水波方程理論,開展庫岸滑坡涌浪災(zāi)害鏈生動(dòng)力學(xué)分析,將滑坡體滑動(dòng)行為宏觀上考慮為均勻非牛頓流體,將滑坡體滑動(dòng)行為及庫水的涌浪傳播行為均采用基于二維淺水波理論進(jìn)行描述,采用有限體積法進(jìn)行數(shù)值求解。在此基礎(chǔ)上,以貴州省索風(fēng)營水利工程近壩庫區(qū)地質(zhì)災(zāi)害為例,開展庫區(qū)滑坡涌浪災(zāi)害動(dòng)力分析。

2 計(jì)算方法與原理

本文計(jì)算基于二維淺水波方程理論,對滑坡涌浪過程進(jìn)行計(jì)算模擬。淺水方程是一個(gè)滿足守恒定律的微分方程組,常用于描述流體在地表上的流動(dòng)過程。在采用淺水波理論進(jìn)行滑坡涌浪分析時(shí),滑坡視為一種流體顆?；旌衔镌谥亓︱?qū)動(dòng)下的運(yùn)動(dòng),其中運(yùn)動(dòng)過程中考慮了在準(zhǔn)靜態(tài)與慣性變形過程中的膨脹與收縮效應(yīng)。在計(jì)算求解中采用深度積分的連續(xù)介質(zhì)力學(xué)方程,Iverson 和 George融合臨界狀態(tài)土力學(xué)、流體力學(xué)與顆粒力學(xué)的相關(guān)概念,提供模型方程的綜合推導(dǎo)及求解方程的數(shù)值方法。

2.1 二維淺水波及深度積分法

淺水波方程是一個(gè)滿足守恒定律的微分方程組,常用于描述流體在地表上的流動(dòng)過程。淺水方程的基本微分形式為：

式中：U為滿足守恒定律的物理量,U=[h u/h v/F]T；h為水深；u為x方向流速；v為y方向流速；E、G分別為x方向和y方向的通量：

右端項(xiàng)為方程源項(xiàng),與實(shí)際求解問題相關(guān),對方程求解有重要影響,常包含重力和河床摩阻力,及其他外源項(xiàng)。

在進(jìn)行滑坡體動(dòng)力學(xué)計(jì)算時(shí),將滑坡體視為固體相與流體相的混合介質(zhì),通過給定固體分?jǐn)?shù)及固體參數(shù)表征滑坡體的特性,在淺水方程求解中常進(jìn)行有限體積差分的格式求解。

2.2 數(shù)值求解方法

對深度平均后的淺水方程控制微分方程組進(jìn)行求解時(shí),采取有限體積方法進(jìn)行數(shù)值求解,使用Godunov格式來求解相鄰網(wǎng)格單元界面處的黎曼解進(jìn)行數(shù)值更新,該方法對求解問題中的不連續(xù)問題求解具有更好的適應(yīng)性,如流動(dòng)中常出現(xiàn)的激波等問題,適于滑坡涌浪過程的模擬。在數(shù)值計(jì)算中基于結(jié)構(gòu)化的長方形網(wǎng)格進(jìn)行劃分網(wǎng)格單元,得到有限體積單元。在離散完成后,對單元網(wǎng)格之間需要沿單元各邊法向求解一維黎曼問題：

其中q=q（x,t）。當(dāng)容量函數(shù)k（x）給定后,基于Godunov求解格式對黎曼方程進(jìn)行變換得到求解方程：

其中Qin為當(dāng)前狀態(tài),A+、A-分別表示相鄰單元界面處的左右兩側(cè)。基于Godunov格式求解后根據(jù)實(shí)際求解條件,及方程數(shù)量求解得到單元處的黎曼解,通常選擇深度、速度等作為基本變量求解,通過迭代計(jì)算對運(yùn)動(dòng)過程進(jìn)行模擬。

3 卞家寨滑坡涌浪數(shù)值計(jì)算模型

3.1 滑坡區(qū)基本概況

卞家寨滑坡位于貴州省烏江流域索風(fēng)營水電站上游水庫庫區(qū)左岸,滑坡體位置距壩址直線距離約1.5 km,堆積體分布高程760 m～1100 m,堆積體平行河床平均長800 m,垂直河床方向平均長600 m,推測厚度50 m～130 m,總體積約3300×104m3。

基于滑坡及庫區(qū)地形信息及滑坡體地質(zhì)勘探資料,建立滑坡床、滑坡體表面地形三維幾何模型（圖1）。

圖1 卞家寨滑坡及庫區(qū)三維模型

3.2 三維數(shù)值計(jì)算模型

計(jì)算考慮水庫正常蓄水情況下蓄水位為837 m?；谏疃绕骄亩S有限體積網(wǎng)格選取為正方形網(wǎng)格,單元邊長選取為5 m,計(jì)算網(wǎng)格規(guī)模為980356。

表1 卞家寨滑坡涌浪模型材料參數(shù)值表

根據(jù)勘察報(bào)告中滑坡體材料力學(xué)特性參數(shù),選取模型計(jì)算參數(shù),根據(jù)二維淺水方程的深度平均模型,模擬計(jì)算中將滑坡體滑動(dòng)行為簡化為同時(shí)包含固體相與流體相的混合介質(zhì)的流動(dòng)行為,并通過分別給定混合介質(zhì)中固體相體積分?jǐn)?shù)及兩相各自材料參數(shù)來描述滑坡體的特性。其中模型中滑坡體混合介質(zhì)固體相密度取2700.00 kg/m3,流體相密度取1000 kg/m3,根據(jù)滑坡體流動(dòng)性流體相粘度取0.005 Pa·s,其他計(jì)算參數(shù)取值參見表 1。

計(jì)算采用重力加速度9.8 m/s3,初始時(shí)間按步長為1.0×10-6s,計(jì)算過程中時(shí)間步長根據(jù)收斂性要求自適應(yīng)調(diào)整,確保計(jì)算穩(wěn)定,計(jì)算區(qū)域在滑坡體位置上游位置處截?cái)?為防止邊界效應(yīng)影響,計(jì)算區(qū)域邊界均采用無反射邊界。

4 卞家寨滑坡涌浪計(jì)算結(jié)果分析

圖2 顯示了數(shù)值計(jì)算得到的不同時(shí)刻涌浪速度分布云圖。為了更好地分析,涌浪沿著滑坡運(yùn)動(dòng)方向及河道傳播特征,圖3、圖4分別顯示了不同時(shí)刻I-I'及II-II'斷面（圖1）庫水位演化。

圖2 卞家寨滑坡涌浪全過程

圖3 I-I'斷面涌浪水位隨傳播距離變化

圖4 沿河道中心線斷面涌浪水面變化曲線

從圖中可以看出,在滑坡初始階段,滑坡體在大約在 5 s時(shí)入水沖擊庫水（圖 3）,滑坡體動(dòng)能傳遞至水體,產(chǎn)生涌浪,形成的涌浪首波高度約6 m（圖3）,并向周圍傳播。此外,從滑坡體的失穩(wěn)特征上來看,位于上游側(cè)的滑坡體高程要相對下游側(cè)的滑坡體高,從而入水時(shí)速度（動(dòng)能）較高；并且從運(yùn)動(dòng)方向上來看,偏向于上游側(cè)。

在約8 s時(shí),滑坡滑動(dòng)速度逐漸降低；在涌浪傳播過程中,未達(dá)到對岸時(shí)隨傳播距離的增加,涌浪高度逐漸降低；在約17 s左右時(shí)涌浪抵達(dá)對岸,由于對岸地形較為陡立,形成強(qiáng)烈的反射,對應(yīng)20 s～40 s結(jié)果可看到滑坡附近水域低流速區(qū),為滑坡近場區(qū)域水面發(fā)生的初始涌浪水波與反射水波的疊加。與此同時(shí),由于滑坡體堆積于庫盆,導(dǎo)致臨近庫水位整抬升；在約40 s時(shí)較正常蓄水位837 m高出2.2 m左右（圖3）。

從河谷形態(tài)上來看,庫盆在卞家寨滑坡處,上游寬闊；而在滑坡體下游側(cè)河道急劇收縮,形成峽谷。庫盆在滑坡入水處的這種特殊的地形及滑坡的運(yùn)動(dòng)特征,使得產(chǎn)生的涌浪沿著庫盆向上游傳播較為強(qiáng)烈,而向下游傳播則相對較弱。從涌浪傳播特征上來看,受滑坡體上游側(cè)涌浪的影響,在滑坡體上游側(cè)的河灣處形成的涌浪較大,并產(chǎn)生較大的爬坡（圖2,20 s～60 s）。在滑坡體的下游側(cè),庫盆庫岸急劇收縮,然后又變開闊,從而在該部位形成一個(gè)明顯的“埡口”地貌。這種庫盆形態(tài),導(dǎo)致涌浪從上游傳播過來時(shí),由于庫盆變窄,使得涌浪流速急劇升高（圖2,20 s～60 s）,并上涌形成較高的涌浪。通過“埡口”后,由于河道變寬,涌浪從波速及波高上急劇降低,形成的涌浪高度約為4.5 m（841.5 m～837 m）,并涌向壩體（圖3）。

從圖4所示的不同時(shí)刻涌浪沿河道傳播可以看出,約60 s左右涌浪開始抵達(dá)壩址,80 s左右波峰抵達(dá)壩址,沖擊壩體同時(shí),產(chǎn)生爬坡,水位升高,約180 s左右時(shí)刻發(fā)生漫頂。

為分析涌浪漫頂產(chǎn)生的洪水特征,如圖5將壩體沿壩軸線均分為3段,圖6和圖7 分別顯示了各段的漫頂水深、漫頂流速以及漫頂流量變化曲線。從圖中可以看出：左壩肩漫頂水深峰值為0.71m,同時(shí)漫頂流速較小為2.0 m/s,最大漫頂流量為87.6 m3/s；右壩肩漫頂水深峰值為0.57 m,漫頂流速為1.8 m/s,最大漫頂流量為101.3 m3/s；壩體中部漫頂水深峰值0.61 m,漫頂流速最大為2.8 m/s,最大漫頂流量為65.2 m3/s。

圖5 大壩分段示意圖

圖6 壩址處沿壩軸線各部分涌浪漫頂水深及流速變化曲線

圖7 大壩分段示意圖

從漫頂洪水整體來看,平均漫頂水深為0.62 m,平均漫頂流速為1.99 m/s,最大頂流量為224 m3/s,漫頂全過程約30 s,且壩體中部漫頂流量最大。因此,從總體上來看,涌浪產(chǎn)生的漫頂不會對大壩產(chǎn)生重大危害,對下游影響也較小,但需要進(jìn)行應(yīng)急防控防止由于產(chǎn)生的洪水對大壩及附屬結(jié)構(gòu)的局部破壞。

5 結(jié)論

滑坡涌浪災(zāi)害是一個(gè)復(fù)雜的流固耦合過程,也一直以來是水庫滑坡研究的重點(diǎn)和難點(diǎn)之一。本文基于二維淺水波方程理論,開展滑坡涌浪問題分析。以貴州省索風(fēng)營電站近壩庫區(qū)滑坡為例,對其涌浪過程進(jìn)行了系統(tǒng)的分析。根據(jù)數(shù)值計(jì)算結(jié)果,滑坡產(chǎn)生的涌浪傳播受滑坡運(yùn)動(dòng)方向、河道地形控制?；麦w整體偏向上游運(yùn)動(dòng),使得涌浪主體方向偏向上游,并在滑坡上游側(cè)凹岸產(chǎn)生強(qiáng)大的涌浪；而下游由于河道存在一個(gè)“埡口”對涌浪產(chǎn)生一定的阻礙作用,從而可以有效底降低涌浪災(zāi)害對下游大壩的影響。

涌浪到達(dá)大壩后,會產(chǎn)生漫壩現(xiàn)象,但從漫壩過程來看,漫頂洪水在壩體中部產(chǎn)生的漫頂水深最大為0.61 m,形成的漫頂洪水流速也最大約為2.8 m/s；漫頂洪水洪峰流量為224 m3/s,漫頂洪水持續(xù)時(shí)間約30 s。因此,從總體上來看,涌浪產(chǎn)生的漫頂不會對大壩產(chǎn)生重大危害,對下游河道影響也較小,但需要進(jìn)行應(yīng)急防控,防止由于產(chǎn)生的洪水對大壩及附屬結(jié)構(gòu)的局部破壞的沖擊。