山區河流型水庫碎土石滑坡涌浪三維數值模擬研究

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

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

0 引 言

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

本文利用CFD數值仿真技術，基于RNGκ-ε紊流模型與VOF方法，將碎土石滑坡體視為固體顆粒和液體的混合物，對已有模型試驗研究成果進行數值模擬，發現數值模擬結果與模型試驗吻合較好。在此基礎上，模擬了某實際山區河流型水庫中碎土石滑坡涌浪的生成和傳播全過程，對涌浪的形態變化、波幅變化和傳播速度進行了分析。研究結果可為實際工程中潛在的滑坡涌浪災害預測提供科學依據。

1 數學模型

(1)

(2)

(3)

(4)

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

2 數值模擬方法可靠性分析

根據已有模型試驗研究成果建立滑坡涌浪計算模型，對散粒體滑坡入水產生涌浪的過程進行數值模擬，以驗證數值模擬方法的可靠性。

2.1 計算模型及網格劃分

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

圖1 模型試驗中各參數示意圖

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

圖a 計算區域示意圖 圖b 局部網格示意圖圖2 數值模擬計算區域及網格示意圖

2.2 計算結果對比與分析

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

(a) t*=0.46

(b)t*=2.52

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

3 實際工程應用

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

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

圖4 庫區涌浪過程波高變化過程 (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給出了不同時刻下涌浪傳播速度的變化過程。滑坡體在入水時具有較高的下滑速度，前緣入水時的速度可達18 m/s，后緣在7.5 s后也基本滑入水中，入水時最大速度可達24 m/s，具有較大的動能。初始涌浪形成后，其傳播波速在庫區可達10～16 m/s，若有船舶經過，海事風險極大。

4 結 語

采用RNGκ-ε紊流模型與VOF方法對散粒體滑坡涌浪進行數值模擬。通過與已有文獻的試驗結果對比分析，發現數值模擬得到的滑坡體水下初始堆積形態與試驗結果吻合較好，且波幅計算精度較高。基于此，建立了某山區河流型水電工程庫區的三維模型，模擬了碎土石滑坡涌浪在庫區內的生成與傳播過程，并分析了涌浪的波幅和波速特征。研究結果表明，本文的數值模擬方法能夠較好地模擬山區河流型水庫內的碎土石滑坡涌浪生成與傳播過程。通過數值模擬獲得的涌浪波幅和傳播速度等參數，可以為大壩安全評估和涌浪災害影響范圍預測提供依據，還可為通航河流航行船舶的避險范圍提供參考。值得說明的是，將碎土石滑坡體視為固體顆粒和液體的混合物，在滑坡下滑的初始階段能夠較為真實地反映滑坡體在岸坡底部的堆積形態，但依然無法改變其流體的本質，在庫區水面波動趨于穩定后，滑坡體會水平堆積于河道底部，但這并不影響前期涌浪生成和傳播數據的采集。另外，VOF方法在模擬強非線性水氣二相流時具有局限性，滑坡體入水時產生的水舌飛濺較難模擬，在實際工程應用中可建立更加細致的三維模型和網格，盡可能確保滑坡體初始形態與實際情況吻合，從而進一步提高模擬計算精度。

圖5 庫區涌浪過程水面速度變化過程 (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