趙 樹 正， 田 宇 霖， 曹 山， 薛 宏 程

(重慶交通大學(xué) 河海學(xué)院，重慶 400016)

0 引 言

隨著國(guó)家長(zhǎng)江經(jīng)濟(jì)帶黃金水道的建設(shè)，沿江兩岸已探明的滑坡就有千余處[1]，若有巨型滑坡體高速滑入水中，其形成的涌浪不僅會(huì)影響過(guò)往船舶的航行安全，還會(huì)對(duì)沿岸居民的生命財(cái)產(chǎn)安全造成嚴(yán)重威脅[2]，因此有效預(yù)測(cè)滑坡涌浪這類次生災(zāi)害的影響范圍具有重要的科學(xué)意義。目前研究滑坡涌浪的方法主要有物理模型試驗(yàn)和數(shù)值模擬[3]，模型試驗(yàn)?zāi)軌蛴行У赜^測(cè)到涌浪生成、傳播和爬坡等過(guò)程，如岳書波等[4]和Heller等[5]通過(guò)物理模型試驗(yàn)，總結(jié)了二維水槽中涌浪的生成形態(tài)和初始波幅經(jīng)驗(yàn)公式；殷坤龍等[6]根據(jù)實(shí)際三維物理模型試驗(yàn)，發(fā)現(xiàn)涌浪在三維空間的波幅遠(yuǎn)小于二維情況。但是物理模型存在縮尺效應(yīng)顯著、數(shù)據(jù)測(cè)量較難和試驗(yàn)成本較高等問(wèn)題。因此，部分研究人員利用數(shù)值模擬方法對(duì)滑坡涌浪形態(tài)變化、傳播過(guò)程開展模擬研究，如孔增增[7]采用數(shù)值模擬等方法對(duì)V型河谷滑坡涌浪的最大波面高度以及涌浪的爬高進(jìn)行模擬，優(yōu)化了涌浪最大爬坡高度經(jīng)驗(yàn)公式；Huang等[8]用非線性Boussinesq水波模型研究了三峽庫(kù)區(qū)紅巖子滑坡涌浪。研究結(jié)果表明涌浪傳播過(guò)程受河谷地貌條件影響較大，滑坡對(duì)岸的城鎮(zhèn)是受涌浪影響的主要區(qū)域。數(shù)值模擬方法作為一種新興的研究手段，在解決地形復(fù)雜的水庫(kù)和河道滑坡涌浪問(wèn)題時(shí)具有明顯優(yōu)勢(shì)。但是滑坡涌浪屬于典型的流固耦合問(wèn)題，滑坡體高速滑入水中還會(huì)卷入大量空氣[9,10]，對(duì)計(jì)算區(qū)域網(wǎng)格劃分要求較高，尤其是進(jìn)行三維數(shù)值模擬難度較大。

本文利用CFD數(shù)值仿真技術(shù)，基于RNGκ-ε紊流模型與VOF方法，將碎土石滑坡體視為固體顆粒和液體的混合物，對(duì)已有模型試驗(yàn)研究成果進(jìn)行數(shù)值模擬，發(fā)現(xiàn)數(shù)值模擬結(jié)果與模型試驗(yàn)吻合較好。在此基礎(chǔ)上，模擬了某實(shí)際山區(qū)河流型水庫(kù)中碎土石滑坡涌浪的生成和傳播全過(guò)程，對(duì)涌浪的形態(tài)變化、波幅變化和傳播速度進(jìn)行了分析。研究結(jié)果可為實(shí)際工程中潛在的滑坡涌浪災(zāi)害預(yù)測(cè)提供科學(xué)依據(jù)。

1 數(shù)學(xué)模型

(1)

(2)

(3)

(4)

本文封閉連續(xù)方程(公式1)和運(yùn)動(dòng)方程(公式2)采用經(jīng)典的RNGκ-ε紊流模型(公式3和公式4)，并利用有限差分法將控制方程離散為代數(shù)方程，從而進(jìn)行數(shù)值求解。涌浪自由面的追蹤采用VOF方法，該方法不考慮氣體對(duì)于模型的影響，而只考慮純液體單元的影響，即在一個(gè)計(jì)算網(wǎng)格內(nèi)，所有流體相的體積分?jǐn)?shù)之和為1。另外，滑坡體用顆粒體來(lái)模擬，是固體顆粒和液體的混合物，這種混合物被視為不可壓縮流體，其邊界可以是自由表面。

2 數(shù)值模擬方法可靠性分析

根據(jù)已有模型試驗(yàn)研究成果建立滑坡涌浪計(jì)算模型，對(duì)散粒體滑坡入水產(chǎn)生涌浪的過(guò)程進(jìn)行數(shù)值模擬，以驗(yàn)證數(shù)值模擬方法的可靠性。

2.1 計(jì)算模型及網(wǎng)格劃分

在長(zhǎng)10 m，寬0.60 m，高1 m的矩形水槽內(nèi)開展了散粒體滑坡涌浪試驗(yàn)研究[11]。該試驗(yàn)的滑坡體材料為沙石顆粒(72%的顆粒粒徑為0.01～2 mm，18%的顆粒粒徑為2～5 mm，10%的顆粒粒徑為5～15 mm)，配比后的滑坡體密度為2 100 kg/m3，滑坡體的入水速度通過(guò)特制的滑動(dòng)裝置來(lái)控制，各項(xiàng)參數(shù)如圖1所示，其中h為水深，vs為滑坡體入水速度，ls為滑坡體長(zhǎng)度，s為滑坡體厚度，α為入水角度[11]。

圖1 模型試驗(yàn)中各參數(shù)示意圖

如圖2所示，本算例的計(jì)算區(qū)域與上述模型試驗(yàn)一致，計(jì)算區(qū)域全部為四邊形網(wǎng)格，網(wǎng)格間距為20 mm，網(wǎng)格總數(shù)約100萬(wàn)；邊界條件包括速度進(jìn)口、壓力出口和無(wú)滑移固壁邊界；初始相為水、空氣和滑坡體?，F(xiàn)以試驗(yàn)[11]中的一組試驗(yàn)參數(shù)值和數(shù)據(jù)作為數(shù)值模擬的依據(jù)，該組試驗(yàn)滑坡體下滑弗氏數(shù)F=1.65，滑坡體相對(duì)厚度S=0.23，水深h=0.3 m，滑坡體下滑速度vs=2.83 m/s，滑坡體長(zhǎng)度ls=0.75 m，入水角度α=33°?；麦w的平均顆粒粒徑d取1.8 mm，滑坡體緊密堆積體積分?jǐn)?shù)取0.63，自然休止角取34°。

圖a 計(jì)算區(qū)域示意圖 圖b 局部網(wǎng)格示意圖圖2 數(shù)值模擬計(jì)算區(qū)域及網(wǎng)格示意圖

2.2 計(jì)算結(jié)果對(duì)比與分析

圖3(a)～3(c)左側(cè)分別為模型試驗(yàn)在相對(duì)無(wú)量綱時(shí)間t*=0.46、t*=2.52和t*=3.43時(shí)的初始涌浪形態(tài)(以滑坡體與水面恰好接觸的時(shí)刻為t*=0)。而圖3(a)～3(c)右側(cè)分別選取了與 試驗(yàn)[11]時(shí)刻對(duì)應(yīng)的三個(gè)數(shù)值模擬涌浪形態(tài)，從圖3(a)(t*=0.46)可以看出數(shù)值模擬中的滑坡體沖擊水面后，水體的運(yùn)動(dòng)趨勢(shì)與試驗(yàn)得到的涌浪形態(tài)基本一致，但由于滑坡體沖擊水面形成的薄水舌和濺起的水花屬于強(qiáng)非線性水氣二相流，在模擬這種水氣間強(qiáng)烈混摻的情況時(shí)VOF法存在局限性，故未能模擬出濺射的水花部分；從圖3(b)(t*=2.52)中可以看出當(dāng)數(shù)值模擬的滑坡體完全滑入水中后，滑坡體的擴(kuò)散形態(tài)與文獻(xiàn)中的滑坡運(yùn)動(dòng)形態(tài)較為接近；由于y軸方向上數(shù)值模擬采用了對(duì)稱邊界條件，而文獻(xiàn)中水槽兩側(cè)的鋼化玻璃對(duì)涌浪傳播起到了限制作用，所以從圖3(c)(t*=3.43)中可以看到滑坡體堆積形態(tài)和最大波幅高度與模型試驗(yàn)相比基本吻合，但最大波幅出現(xiàn)的位置比試驗(yàn)傳播的稍遠(yuǎn)。總體上看，數(shù)值模擬得到的滑坡體初始堆積形態(tài)與試驗(yàn)結(jié)果吻合較好，波幅計(jì)算精度能夠滿足要求。

(a) t*=0.46

(b)t*=2.52

(c) t*=3.43圖3 不同時(shí)刻下的涌浪形態(tài)對(duì)比[11]

3 實(shí)際工程應(yīng)用

某水電工程位于瀾滄江一級(jí)支流的下游河段，是以發(fā)電為主的混合式開發(fā)水電站，屬于典型的山區(qū)河流型水庫(kù)。工程主要建筑物包括攔河大壩、泄水建筑物及引水發(fā)電建筑物等，水庫(kù)總庫(kù)容3 403萬(wàn)m3，電站總裝機(jī)容量100 MW，為中等規(guī)模Ⅲ等水電工程。工程勘測(cè)階段發(fā)現(xiàn)壩址上游左岸約200 m處有一古滑坡體(滑坡后緣高出正常蓄水位約60 m)，由于水庫(kù)蓄水后可能會(huì)引起兩岸山體地下水位抬升，導(dǎo)致岸坡巖土體強(qiáng)度降低，從而誘發(fā)滑坡。故本文針對(duì)該水電工程庫(kù)區(qū)蓄水后潛在的滑坡涌浪災(zāi)害開展了數(shù)值模擬研究。

圖4分別給出了時(shí)刻為1.5 s、3 s、4.5 s、7.5 s、10.5 s和13.5 s時(shí)庫(kù)區(qū)涌浪的波高變化過(guò)程(靜止水面高程設(shè)置為0)。當(dāng)滑坡體滑入水中，在近場(chǎng)處形成近30 m高的初始涌浪(t=4.5 s時(shí))，隨著涌浪向四周擴(kuò)散傳播，其沿滑坡體兩側(cè)河岸傳播時(shí)呈現(xiàn)出明顯的衰減趨勢(shì)，但主波區(qū)波高仍可達(dá)十余米。當(dāng)涌浪傳播到對(duì)岸山體時(shí)，波幅衰減到十米內(nèi)，但涌浪最大波峰離壩體較遠(yuǎn)。

圖4 庫(kù)區(qū)涌浪過(guò)程波高變化過(guò)程 (a)t=1.5 s; (b)t=3 s; (c)t=4.5 s; (d)t=7.5 s; (e)t=10.5 s; (f)t=13.5 s

圖5給出了不同時(shí)刻下涌浪傳播速度的變化過(guò)程。滑坡體在入水時(shí)具有較高的下滑速度，前緣入水時(shí)的速度可達(dá)18 m/s，后緣在7.5 s后也基本滑入水中，入水時(shí)最大速度可達(dá)24 m/s，具有較大的動(dòng)能。初始涌浪形成后，其傳播波速在庫(kù)區(qū)可達(dá)10～16 m/s，若有船舶經(jīng)過(guò)，海事風(fēng)險(xiǎn)極大。

4 結(jié) 語(yǔ)

采用RNGκ-ε紊流模型與VOF方法對(duì)散粒體滑坡涌浪進(jìn)行數(shù)值模擬。通過(guò)與已有文獻(xiàn)的試驗(yàn)結(jié)果對(duì)比分析，發(fā)現(xiàn)數(shù)值模擬得到的滑坡體水下初始堆積形態(tài)與試驗(yàn)結(jié)果吻合較好，且波幅計(jì)算精度較高?；诖耍⒘四成絽^(qū)河流型水電工程庫(kù)區(qū)的三維模型，模擬了碎土石滑坡涌浪在庫(kù)區(qū)內(nèi)的生成與傳播過(guò)程，并分析了涌浪的波幅和波速特征。研究結(jié)果表明，本文的數(shù)值模擬方法能夠較好地模擬山區(qū)河流型水庫(kù)內(nèi)的碎土石滑坡涌浪生成與傳播過(guò)程。通過(guò)數(shù)值模擬獲得的涌浪波幅和傳播速度等參數(shù)，可以為大壩安全評(píng)估和涌浪災(zāi)害影響范圍預(yù)測(cè)提供依據(jù)，還可為通航河流航行船舶的避險(xiǎn)范圍提供參考。值得說(shuō)明的是，將碎土石滑坡體視為固體顆粒和液體的混合物，在滑坡下滑的初始階段能夠較為真實(shí)地反映滑坡體在岸坡底部的堆積形態(tài)，但依然無(wú)法改變其流體的本質(zhì)，在庫(kù)區(qū)水面波動(dòng)趨于穩(wěn)定后，滑坡體會(huì)水平堆積于河道底部，但這并不影響前期涌浪生成和傳播數(shù)據(jù)的采集。另外，VOF方法在模擬強(qiáng)非線性水氣二相流時(shí)具有局限性，滑坡體入水時(shí)產(chǎn)生的水舌飛濺較難模擬，在實(shí)際工程應(yīng)用中可建立更加細(xì)致的三維模型和網(wǎng)格，盡可能確保滑坡體初始形態(tài)與實(shí)際情況吻合，從而進(jìn)一步提高模擬計(jì)算精度。

圖5 庫(kù)區(qū)涌浪過(guò)程水面速度變化過(guò)程 (a)t=1.5 s; (b)t=3 s; (c)t=4.5 s; (d)t=7.5 s; (e)t=10.5 s; (f)t=13.5 s