賀翠玲, 李鵬峰, 張 智, 荊海曉, 萬克誠,4, 段軍邦

(1.陜西省水生態(tài)環(huán)境工程技術(shù)研究中心, 陜西 西安 710065; 2.中國電建西北勘測設(shè)計(jì)研究院有限公司,陜西 西安 710065;3.西安理工大學(xué), 西北旱區(qū)生態(tài)水利國家重點(diǎn)實(shí)驗(yàn)室, 陜西 西安 710048; 4.國家能源水電工程技術(shù)研發(fā)中心高邊坡與地質(zhì)災(zāi)害研究治理分中心, 陜西 西安 710065; 5.黃河上游水電開發(fā)有限責(zé)任公司, 青海 西寧 810000)

1 研究背景

庫區(qū)滑坡涌浪是一種破壞性極高的次生災(zāi)害,高速失穩(wěn)體可在庫區(qū)產(chǎn)生巨大的涌浪,傳播至壩前可能會翻越壩體,甚至?xí)斐蓧误w潰決,對擋水建筑物及庫區(qū)居民造成嚴(yán)重威脅[1-3]。歷史上發(fā)生過較大的滑坡涌浪事故,如:1958年美國阿拉斯加州的Lituya Bay滑坡,涌浪波最大爬坡高度達(dá)到524 m[4];1963年意大利瓦伊昂水庫滑坡造成約2 000人死亡[5];我國三峽庫區(qū)也是滑坡涌浪災(zāi)害發(fā)生較多的區(qū)域,據(jù)不完全統(tǒng)計(jì),三峽庫區(qū)共監(jiān)測到滑坡3 028處,其中2003年7月13日千將坪滑坡產(chǎn)生的涌浪波在對岸爬高超過30 m,在庫區(qū)傳播距離超過10 km,造成了巨大的經(jīng)濟(jì)損失[6-8]。因此,滑坡涌浪生成過程、涌浪波在庫區(qū)演變及岸坡爬坡過程的精準(zhǔn)預(yù)測對防止水庫漫壩、降低滑坡涌浪災(zāi)害至關(guān)重要。

為了對滑坡涌浪災(zāi)害問題進(jìn)行預(yù)警和風(fēng)險(xiǎn)評估,眾多學(xué)者針對滑坡涌浪生成過程和浪高預(yù)測開展了細(xì)致的研究。黃智敏等[9]模擬了樂昌峽水庫鵝公帶滑坡體不同滑速工況對浪高、漫壩水量及動水壓強(qiáng)的影響;中國水利水電科學(xué)研究院提出浪高預(yù)測中的主導(dǎo)因素為滑坡體入水的速度和滑坡體的方量,并以此假設(shè)為基礎(chǔ)建立了浪高預(yù)測公式[10];潘家錚[11]根據(jù)條分滑坡體提出了計(jì)算滑坡涌浪初始浪高的方法,后續(xù)有學(xué)者用潘家錚方法和水科院經(jīng)驗(yàn)公式法論證了坡腳壓腳措施對降低涌浪高度的效果[12];Clous等[13]基于能量轉(zhuǎn)換,研究了散粒體滑坡產(chǎn)生的涌浪波波能問題;Evers等[14]基于理論分析研究了涌浪波在時(shí)間和空間上的脈沖性,探討了涌浪波在傳播過程中的震蕩危害;曹婷等[15]利用物理模型試驗(yàn)研究了滑坡體形狀對涌浪爬高的影響;肖莉麗等[16]針對三峽庫區(qū)滑坡涌浪問題建立了1∶200的物理試驗(yàn)?zāi)Ｐ?研究了近源區(qū)多因素對首浪高度的影響;黃錦林等[17]建立了1∶150的樂昌峽水庫物理模型,研究了滑坡體滑速對涌浪特性的影響,并將實(shí)驗(yàn)結(jié)果與多個(gè)經(jīng)驗(yàn)公式進(jìn)行了對比;馬鑫磊等[18]比選了在20多個(gè)典型滑坡涌浪實(shí)際災(zāi)害評價(jià)中用到的國內(nèi)外主流滑坡涌浪浪高預(yù)測方法,并指出了主流評價(jià)方法的不足;黃宇云等[19]應(yīng)用三維流體模擬軟件Fluent,假定長嶺皮水庫滑坡體為顆粒流,研究了滑坡體不同體積與壩址最大浪高的關(guān)系。

綜上所述,現(xiàn)有研究主要針對寬廣水域或無壩河段的滑坡涌浪問題,對于狹窄型河道水庫而言,其涌浪生成、傳播甚至于爬坡的影響因素更加復(fù)雜,因此寬廣水域滑坡涌浪規(guī)律未必能反映出近壩庫區(qū)狹窄河道的滑坡涌浪規(guī)律。此外,流體動力學(xué)三維數(shù)值模擬軟件Flow Hydro在滑坡涌浪問題研究中的應(yīng)用性能也需進(jìn)一步研究和評價(jià)。因此,本研究針對黃河羊曲水電站狹窄庫區(qū)1#變形體失穩(wěn)后不同滑速產(chǎn)生的涌浪災(zāi)害問題,基于Flow Hydro模擬軟件分析庫區(qū)滑坡涌浪的產(chǎn)生、發(fā)展和傳播過程,進(jìn)一步探討了滑坡入水速度對涌浪特征的影響以及對工程的危害,以期為大壩的安全提供技術(shù)支撐。

2 數(shù)值模型及驗(yàn)證

羊曲水電站位于青海省海南藏族自治州興?？h與貴南縣交界處的黃河干流上,屬黃河上游水系。工程規(guī)模為一等大(1)型工程,壩頂高程為2 721 m,防浪墻高度為1.2 m,水庫正常蓄水位為2 715 m[20]。通過地質(zhì)調(diào)查發(fā)現(xiàn),H1滑坡體位于距壩軸線約1.2～2.0 km的壩前左岸,1#變形體位于壩址上游左岸,其上游側(cè)邊界與H1滑坡相鄰,下游邊界距壩址約750 m,變形體平面上呈不規(guī)則扇形,前緣寬度為710 m,后部較窄部位寬度約為100 m,變形體平均厚度為25 m,體積約為500×104m3,該變形體具有明顯的強(qiáng)傾倒層和滑移拉裂層?；麦w所處位置庫區(qū)狹長,其上游和壩前河道寬廣,具有傳播方向上的二維特征和傳播變形方向上的三維特征。因此,本研究建立基于Flow Hydro的三維數(shù)學(xué)模型。

2.1 數(shù)值模型原理

滑坡涌浪問題涉及復(fù)雜的流體和固體強(qiáng)耦合過程以及流體飛濺和固體破碎變形。因此,流體的基礎(chǔ)三維特性通過求解三維Navier-Stokes方程獲得,在此基礎(chǔ)上進(jìn)一步耦合紊流RNG模型來實(shí)現(xiàn)流體復(fù)雜變形的模擬[21],流體表面識別采用改進(jìn)的VOF(volume of fuid)技術(shù)[22],其控制方程如下。

連續(xù)性方程:

(1)

動量方程:

(2)

式中:u、v、w分別為x、y、z方向的速度,m/s;Ax、Ay、Az分別為x、y、z方向的水流面積,m2;t為時(shí)間,s;ρ為水體的密度,kg/m3,模型材料屬性默認(rèn)加載值為1 000 kg/m3;p為壓強(qiáng),Pa;RSOR為源項(xiàng);VF為單元網(wǎng)格體中流體的體積分?jǐn)?shù),水氣交界面取值為0.5,(Gx,Gy,Gz)為作用在單位質(zhì)量流體微元的重力,N; (fx,fy,fz)為作用在單位質(zhì)量流體微元的黏性力,N;Uw=(uw,vw,ww)為源項(xiàng)的速度,m/s;us,vs,ws為相對速度,m/s。

紊流模型為RNG模型:

DiffT-εT

(3)

式中:DiffT為紊流擴(kuò)散項(xiàng);PT為紊流動能產(chǎn)生項(xiàng);GT為浮力產(chǎn)生項(xiàng);kT為紊流動能;εT為紊流耗散項(xiàng)。

對于滑坡與水體的流固耦合過程,只要能夠準(zhǔn)確描述滑坡體運(yùn)動過程,將該過程體現(xiàn)在網(wǎng)格上,通過VOF在網(wǎng)格上識別固體區(qū)就可實(shí)現(xiàn)滑坡體與水體的耦合。因此選用Flow Hydro中穩(wěn)定性和精度最為突出的GMO(general moving object)模型來描述滑坡運(yùn)動,由每個(gè)時(shí)刻網(wǎng)格的變化來逐步更新滑坡體的位置,采用VOF法捕捉每個(gè)時(shí)刻的水氣交界面、固液交界面及固氣交界面[23-24]。

2.2 數(shù)值模型建立

研究范圍的選取和模型網(wǎng)格的劃分需要考慮計(jì)算成本、計(jì)算精度及區(qū)域內(nèi)浪高干擾等多方面因素,本研究中涌浪模型下游側(cè)選取到羊曲水電站壩下游,只要不影響波浪在壩身的爬高即可;上游側(cè)選取到羊曲水電站壩上游地形開闊區(qū)域,該區(qū)域可以使得波浪盡可能衰減,減小反射波對觀測值的干擾;左、右岸范圍包含滑坡體即可。按照上述原則建立的滑坡涌浪三維數(shù)值模型區(qū)域及邊界條件如圖 1所示。圖 1中具體的模型范圍為:1#變形體上游側(cè)邊界向上游延伸1 km,下游至羊曲水電站大壩;左岸沿河地形模擬至3 000 m高程,右岸模擬至2 800 m高程。模型網(wǎng)格設(shè)計(jì)方案中,沿x方向(寬度方向)寬度為1 500 m,沿y方向(長度方向)長度為3 600 m,沿z方向(高度方向)高度為450 m;網(wǎng)格尺寸在x方向和y方向?yàn)?.0 m,z方向?yàn)?.5 m,網(wǎng)格總數(shù)為4 860×104。在邊界條件設(shè)置方面,左岸、右岸、上游及河床底部均設(shè)置為墻邊界(Wall),下游設(shè)置為開邊界,漫壩水體自由出流(Outflow),大壩頂部設(shè)置通量面,用來監(jiān)測漫壩水體通量,頂部邊界設(shè)置為大氣壓(Pressure=0)。

2.3 數(shù)值模型驗(yàn)證

為了對數(shù)值模型進(jìn)行驗(yàn)證,本研究按重力相似準(zhǔn)則設(shè)計(jì)建立幾何比尺為1∶200的物理模型,根據(jù)《水電水利工程滑坡涌浪模擬技術(shù)規(guī)程》(DL/T 5246—2010)相關(guān)要求[25],模型滑坡體滿足幾何相似和塊體的比重相似,河道和岸坡滿足阻力相似。圖2為本研究所建立的羊曲水電站滑坡涌浪試驗(yàn)?zāi)Ｐ?圖3為模型試驗(yàn)觀測設(shè)備布置圖?；麦w通過不同尺寸的混凝土塊組合模擬,滑車通過多鉸鏈的剛板組合實(shí)現(xiàn);涌浪時(shí)程變化數(shù)據(jù)的記錄采用CBG03智能浪高儀,沿程總共布置16個(gè)浪高儀(Ob1～Ob16);滑坡體下滑速度、滑坡區(qū)域涌浪爬高及建筑物前涌浪爬高采用高速攝像機(jī)采集,漫壩水量使用量筒測量。

在1#變形體滑坡體積為100×104m3、滑坡啟動高程為2 755 m(試驗(yàn)入水速度為14.05 m/s)、水庫正常運(yùn)行水位為2 715 m條件下,進(jìn)行數(shù)值模擬和模型試驗(yàn)觀測,為了實(shí)現(xiàn)數(shù)值模擬與模型試驗(yàn)中滑坡體的相似性,模擬時(shí)將滑坡體切分成小塊。選取滑坡處(浪高儀Ob8)和壩前(浪高儀Ob3)兩個(gè)測點(diǎn),通過涌浪浪高模擬數(shù)據(jù)與試驗(yàn)實(shí)測數(shù)據(jù)比較來驗(yàn)證數(shù)值模型的準(zhǔn)確性,圖4為該兩個(gè)測點(diǎn)涌浪浪高時(shí)程變化的模擬值與試驗(yàn)值對比。從圖4可以看出,兩個(gè)測點(diǎn)的涌浪浪高和相位基本一致,首波和次生波的模擬值與試驗(yàn)值吻合精度高,而尾波吻合精度較差。分析其原因:由于模型試驗(yàn)的實(shí)際入水速度為14.05 m/s,與目標(biāo)速度或者數(shù)值模擬給定的入水速度15 m/s有偏差,這兩種速度的差異會反映到后續(xù)涌浪浪高和波長上;此外,數(shù)值模型上游區(qū)域未做消波處理,而物理模型試驗(yàn)中在上游加設(shè)了溢流堰來防止反射波的影響,各因素疊加導(dǎo)致后續(xù)波的擬合效果較差。但后續(xù)涌浪浪高較小,對建筑物的影響相對于首波和次生波來說可以忽略。綜合分析認(rèn)為,該數(shù)值模型計(jì)算精度滿足要求,可用于后續(xù)研究。

3 結(jié)果分析與討論

在水庫位為正常運(yùn)行水位2 715 m條件下,分析討論1#變形體滑坡方量為100×104m3的矩形斷面分別以5、10、15、20、25、30及35 m/s的入水速度剛性整體下滑時(shí)產(chǎn)生的初始涌浪高度及對岸坡和建筑的影響?；掠坷藬?shù)值模擬工況見表1。

表1 羊曲水電站庫區(qū)滑坡涌浪數(shù)值模擬工況

圖5為工況4(入水速度為20 m/s)滑坡涌浪流場變化過程云圖。由圖 5可觀察到,1#變形體所處位置的庫區(qū)水面較為狹窄,滑坡體入水使得大部分水體被直接推向?qū)Π?僅有少部分水體形成首波向上、下游河道傳遞(圖5(b));滑坡體坍塌后的岸坡形成空缺,使水體開始后潰,后潰水體反射再次生成次生波(圖5(c));起初滑坡體推向?qū)Π兜乃w在重力作用下從對岸岸坡跌落,再次形成涌浪繼續(xù)向上、下游傳播(圖5(d))。李世貴[6]、鄧成進(jìn)等[26]通過研究發(fā)現(xiàn),滑坡入水產(chǎn)生涌浪后,波浪場呈半圓曲線向遠(yuǎn)端傳遞,但本研究中的涌浪生成過程具有明顯的狹窄庫區(qū)特性,地形嚴(yán)重影響了涌浪波的生成以及波場的反射和傳播,因而無明顯的圓形曲線傳播規(guī)律。

圖6為各工況滑坡處Ob8測點(diǎn)水面波浪高度(相對于庫水位2 715 m的浪高)時(shí)程線模擬結(jié)果。由圖 6可知,首波的波峰高度并不一定是最大的,某些情況下次生波的波峰高度最高,這是由于滑坡產(chǎn)生的波浪與對岸反射波疊加所致,在該過程中,涌浪波表現(xiàn)出明顯的震蕩性[27]。隨著滑坡體滑速的增大,首浪浪高呈減小的趨勢,這雖然與常規(guī)二維水槽試驗(yàn)的結(jié)果不同,但與模型試驗(yàn)觀測結(jié)果一致[28]。分析產(chǎn)生這一結(jié)果的原因可歸納為兩個(gè)方面:一方面隨著滑坡體速度的增大,能量傳遞到對岸岸坡導(dǎo)致涌浪爬坡高度增加,而縱向擴(kuò)散減小;另一方面如果滑坡體速度過大,則首波傳遞到測點(diǎn)時(shí)涌浪并未成形,還在發(fā)展過程中。由圖6還可看出,波速隨著滑坡體滑速的增大而增大,滑速為35 m/s時(shí)形成的滑坡涌浪最先傳播到測點(diǎn),滑速為5 m/s時(shí)的滑坡涌浪最后傳播到測點(diǎn)。各工況波峰高度均在壩頂高程附近,如果涌浪在傳播過程中衰減較小,就有可能發(fā)生漫壩。

圖7為工況4壩前區(qū)域9個(gè)時(shí)刻的波場云圖。由圖7可以看出,該工況下會發(fā)生多次波散射,波場復(fù)雜。由于波反射和疊加作用,某些位置的波高可能高于入射波。圖8為各工況壩前中部Ob3測點(diǎn)水面波浪高度(相對于庫水位2 715 m的浪高)時(shí)程線模擬結(jié)果。由圖8可知,首先,各滑速生成的涌浪波在從Ob8測點(diǎn)傳至Ob3測點(diǎn)過程中并沒有出現(xiàn)大幅度衰減(與圖 6對比分析)。其次,各工況的首波波高并非最大,最大波高為120～125 s出現(xiàn)的第4個(gè)波峰,這是由于壩前區(qū)域的多次波反射所致。此外,幾乎所有工況生成的涌浪波均會出現(xiàn)漫壩現(xiàn)象,漫壩位置主要在左、右壩肩處。

圖9為各工況溢洪道進(jìn)口(壩左)、壩中及電站進(jìn)水口(壩右)處的次生涌浪浪高變化規(guī)律。從圖 9中可以看出,當(dāng)滑坡體入水速度從5 m/s增加到10 m/s時(shí),壩前浪高明顯增大;當(dāng)入水速度在10～35 m/s之間時(shí),浪高保持緩慢增加趨勢。工況1～7建筑物前浪高大小排序如下:溢洪道進(jìn)口處(壩左)>電站進(jìn)水口處(壩右)>壩中,壩前水面呈兩端高中間低的態(tài)勢。

圖10為各工況滑坡涌浪引起的漫壩流量時(shí)程線。由圖10可以看出,所有工況都會發(fā)生多次漫壩,與前述壩前浪高分析一致。首波對漫壩流量的貢獻(xiàn)很小,漫壩流量的第一個(gè)峰值出現(xiàn)時(shí)刻在100 s左右,但是首波到達(dá)壩前的時(shí)刻約為40 s,其時(shí)間差距較大,這表明滑坡生成的涌浪波首波不一定是決定災(zāi)害等級的關(guān)鍵性指標(biāo),涌浪波后續(xù)的傳播和散射也非常重要,在實(shí)際工程中必須進(jìn)行準(zhǔn)確預(yù)測。圖11為各工況漫壩水量統(tǒng)計(jì)結(jié)果。圖11表明,漫壩水量隨著滑坡入水速度的增大而增加。

圖1 滑坡涌浪三維數(shù)值模型區(qū)域及邊界條件示意圖

圖2 羊曲水電站庫區(qū)滑坡涌浪試驗(yàn)?zāi)Ｐ?/p>

圖3 羊曲水電站庫區(qū)滑坡涌浪模型試驗(yàn)觀測設(shè)備布置

圖4 滑坡處和壩前兩測點(diǎn)涌浪浪高時(shí)程變化的模擬值與試驗(yàn)值對比

圖5 工況4(入水速度20 m/s)滑坡涌浪流場變化過程云圖

圖6 各工況滑坡處Ob8測點(diǎn)水面波浪高度時(shí)程線

圖7 工況4(入水速度20 m/s)壩前區(qū)域9個(gè)時(shí)刻的波場云圖

圖8 各工況壩前中部Ob3測點(diǎn)水面波浪高度時(shí)程線

圖9 各工況建筑物前次生涌浪浪高變化規(guī)律

圖10 各工況滑坡涌浪引起的漫壩流量時(shí)程線

圖11 各工況漫壩水量統(tǒng)計(jì)圖

4 結(jié) 論

在本研究中,使用三維數(shù)值模型模擬了黃河上游羊曲水電站庫區(qū)1#變形體滑坡涌浪的產(chǎn)生、傳播、爬高和漫頂全過程。通過1∶200比尺的物理模型進(jìn)行涌浪時(shí)程線的比較,驗(yàn)證了數(shù)值模型的有效性,模擬分析了滑坡入水速度對波高和漫壩過程的影響。得出以下結(jié)論:

(1)物理模型驗(yàn)證表明,本研究所建立的數(shù)值模型是準(zhǔn)確的,該數(shù)值模型對首波和次生波的生成、傳播及漫壩過程均能準(zhǔn)確刻畫。

(2)隨著滑坡入水速度的增加,首波波峰高度減小,而且首波波峰并不總是最高的,這與滑坡涌浪生成過程及庫區(qū)地形影響有關(guān)。

(3)在滑坡生成的涌浪向大壩傳播過程中,首波沒有發(fā)生明顯的衰減,傳播到壩址處時(shí),溢洪道進(jìn)口(壩左)波高最高,電站進(jìn)水口(壩右)波高次之,壩中波高最小。

(4)滑坡體入水速度在5～35 m/s時(shí),涌浪波傳播到壩前均會發(fā)生漫頂,但首波對漫壩總水量的貢獻(xiàn)很小。此外,漫壩流量峰值出現(xiàn)的時(shí)刻遠(yuǎn)滯后于首波到達(dá)壩前的時(shí)刻,這是由大壩附近波浪的多次散射所致,這一現(xiàn)象在實(shí)際工程中必須引起重視。