彭 輝，黃亞杰

(三峽大學(xué) 水利與環(huán)境學(xué)院，湖北 宜昌 443000)

1 研究背景

彎曲河道型水庫(kù)庫(kù)區(qū)因極端氣候、地震和不良地質(zhì)等條件的影響，易導(dǎo)致失穩(wěn)的變形體引起的巨大滑坡涌浪，不但會(huì)對(duì)沿岸的居民和過往的船只造成威脅，而且還會(huì)危及大壩的安全。因此研究彎曲河道型庫(kù)區(qū)滑坡涌浪傳播規(guī)律及其與大壩作用過程有著重要意義[1-2]。

許多復(fù)雜因素都可以影響彎曲河道型水庫(kù)庫(kù)區(qū)滑坡涌浪傳播，比如地形地貌、滑坡體大小等，因此經(jīng)驗(yàn)公式成為了一種簡(jiǎn)要的計(jì)算手段，常用的經(jīng)驗(yàn)公式有許多比如潘家錚法、水科院法、Noda E、美國(guó)土木工程學(xué)會(huì)法[3]。國(guó)內(nèi)外多位專家和學(xué)者，根據(jù)不同滑坡體的幾何形狀、材料、入射角度等特征因素對(duì)滑坡涌浪的影響開展了深入研究，研究手段有水工物理模型試驗(yàn)和數(shù)值仿真分析[4-5]。對(duì)于數(shù)值仿真分析滑坡涌浪有以下兩種方式簡(jiǎn)化：①簡(jiǎn)化流體力學(xué)的方程模型[6]；②簡(jiǎn)化計(jì)算模型[7]。這些簡(jiǎn)化方法的弊端在于沒有考慮地形的實(shí)際情況，而且很少考慮涌浪與大壩的相互作用和下游大壩對(duì)涌浪傳播的影響。對(duì)于水工物理模型試驗(yàn)需要大量的財(cái)力和人力支持，由于客觀條件又不能獲得準(zhǔn)確數(shù)據(jù)，限制了其推廣。

本文針對(duì)文獻(xiàn)[8]基于Flow-3D對(duì)其水工物理模型進(jìn)行三維數(shù)值模擬。通過對(duì)物理試驗(yàn)?zāi)Ｐ?∶1比例建模，充分考慮地形因素采用彎曲河道并在河道尾部增加重力壩。采用有限差分法進(jìn)行水流和水力學(xué)模擬分析滑坡涌浪在彎曲河道中的傳播規(guī)律[9]，通過三維數(shù)值模擬得到滑坡涌浪的傳播的完整的數(shù)據(jù)，進(jìn)而研究了涌浪與大壩的相互作用規(guī)律[10]。

2 彎曲河道型庫(kù)岸滑坡涌浪的物理模型實(shí)驗(yàn)簡(jiǎn)介

水工物理模型主要由滑坡體、滑動(dòng)控制板、河道滑坡區(qū)、水庫(kù)擴(kuò)大段和彎曲型河道構(gòu)成。物理模型試驗(yàn)采用不同滑坡體體積、水深、入水角度、落差等影響因素，研究了滑坡激起涌浪變化特點(diǎn)，并結(jié)合實(shí)驗(yàn)數(shù)據(jù)擬合出首浪高度計(jì)算公式[11]。具體模型如圖1所示。

圖1 彎曲河道模型實(shí)物圖Fig.1 Physical model of curved river channel

3 基于Flow-3D的滑坡涌浪數(shù)值模擬

3.1 基本理論

Flow-3D采用有限差分法來對(duì)控制方程進(jìn)行數(shù)值離散求解，利用極小殘差算法求解，該算法具有相對(duì)收斂快、精度高的特征，廣泛應(yīng)用于計(jì)算流體[12-13]。波浪在Flow-3D軟件中被視為不可壓縮的粘性流體[14]，所以流體運(yùn)動(dòng)控制方程為Navier-Stokes連續(xù)性方程。計(jì)算的流體視為牛頓流體?；麦w下滑時(shí)將會(huì)導(dǎo)致水面出現(xiàn)劇烈的抖動(dòng)與變形。因此將采用RNKκ-ε紊流模型來計(jì)算。RNKκ-ε湍流模型可以很好模擬流線彎曲的流動(dòng)[15-16]。通過面積分?jǐn)?shù)和體積分?jǐn)?shù)來定義有限元網(wǎng)格中的物體，有限元網(wǎng)格能夠通過流體體積(VOF)法精確地模擬出復(fù)雜的幾何形狀，還可以模擬三維流體的流場(chǎng)分布以及不同液體靜止或者運(yùn)動(dòng)復(fù)雜界面相互作用的液面變化。

它的計(jì)算原理就是定義一個(gè)單元流體函數(shù)，稱為F(單元內(nèi)所占流體體積除以該單元可容納流體體積)。F=0表示單元內(nèi)沒有流體；F=1表示流體完全充滿單元；0

3.2 三維滑坡涌浪數(shù)值模型

3.2.1 模型的邊界條件和初始條件的設(shè)定

三維數(shù)值模型中，本文將固體和流體接觸面的邊界條件設(shè)定為墻壁(wall)邊界條件，水槽的側(cè)壁和底部視為光滑無摩擦的自由液面，流體法向速度等于0。流體表面采用自由液面邊界條件，設(shè)置為一個(gè)標(biāo)準(zhǔn)大氣壓強(qiáng)。

模型運(yùn)動(dòng)學(xué)邊界條件為

(1)

模型動(dòng)力學(xué)邊界條件為

(2)

式中：η為波浪表面相對(duì)靜止水面高度；?φ為自由表面上所有速度分量變化率，?η為各方向水頭變化率，φ為各方向水流速度；g為重力加速度。

將試驗(yàn)滑塊簡(jiǎn)化為矩形剛性滑塊，模擬在無摩擦情況下整體滑入水中，滑塊僅受重力，不考慮空氣阻力。水槽中的水靜止，所有網(wǎng)格初始速度設(shè)定為0，滑塊以起始靜止?fàn)顟B(tài)下落。重力加速度g=9.81 m/s2，水的密度ρw=1 000 kg/m3。為了驗(yàn)證三維數(shù)值模型與水工物理模型試驗(yàn)結(jié)果的準(zhǔn)確性以及可行性，滑塊選擇了水工物理試驗(yàn)產(chǎn)生最大首浪高度所用試塊尺寸為0.4 m×0.3 m×0.3 m(長(zhǎng)×寬×高)，滑塊密度為2 400 kg/m3。

3.2.2 模型的尺寸及網(wǎng)格劃分

為了滿足與實(shí)際中滑坡的相似性，大多數(shù)河道的河流轉(zhuǎn)角處的彎折角度都在60°左右，故將河道模型彎折角設(shè)置為60°?；麦w入水方向與水流方向垂直。河道前部滑坡區(qū)域尺寸為1 m×1.3 m；彎曲河道段長(zhǎng)2.5 m，寬0.5 m，高1 m，河道截面為矩形。水庫(kù)段模型采用突變擴(kuò)大段鏈接。水庫(kù)區(qū)尺寸1.1 m×1.5 m，滑板長(zhǎng)1.5 m，寬度為0.5 m，滑動(dòng)面角度調(diào)節(jié)范圍為40°～90°。壩體為直立混凝土重力壩。三維數(shù)值模型具體尺寸如圖2所示。

圖2 彎曲型河道模型尺寸Fig.2 Dimensions of curved river model

在Flow-3D軟件中，網(wǎng)格劃分密度與模型計(jì)算精度、模擬效果呈正相關(guān)。因此，本文幾何模型計(jì)算區(qū)域采用均勻的矩形交錯(cuò)網(wǎng)格，將幾何模型導(dǎo)入Flow-3D軟件(STL格式)，共計(jì)2×106個(gè)。通過流固耦合模型(GMO)來模擬剛體的運(yùn)動(dòng)。計(jì)算時(shí)間步長(zhǎng)為0.062 5 s，計(jì)算時(shí)間為5 s。圖3顯示初始水位0.7 m、入射角60°、滑塊入水速度2.54 m/s的三維數(shù)值模擬示意圖。

圖3 三維數(shù)值模擬示意圖Fig.3 Three-dimensional numerical simulation

3.3 滑坡涌浪在彎曲河道中傳播過程分析

水工物理模型試驗(yàn)監(jiān)測(cè)點(diǎn)布置在圖2中M1—M4，這些監(jiān)測(cè)點(diǎn)的數(shù)據(jù)作為對(duì)比值是因?yàn)榛瑝K入水時(shí)產(chǎn)生涌浪效果顯著。以0.8 m為初始水位、入射角為60°、水面大氣壓為0.1 MPa，來分析滑坡涌浪在彎曲河道中的傳播規(guī)律。

滑塊在自重作用下，沿著斜面由靜止開始向下滑至水面并開始擠壓水體引起大尺度紊流，涌浪由落水點(diǎn)向四周轉(zhuǎn)播形成巨大浪花，附近水域水體迅速升高，最終形成了初始涌浪。初始涌浪在岸坡上爬至最高后水體下落與后續(xù)涌浪疊加傳入河道。河道左岸水域水面迅速抬升，在左岸爬升到最高處后迅速下落，并與后續(xù)傳來的波峰相互作用和震蕩。庫(kù)區(qū)水體也隨著落水點(diǎn)形成的涌浪傳播圈升起、降落，推動(dòng)著涌浪不斷向壩前傳播。在首浪駛出彎曲河道來到水庫(kù)擴(kuò)大段傳至壩前，由于壩體的阻擋，壩前水位持續(xù)升高，壩前水位達(dá)到最高。隨后壩前水位回落，涌浪在河道內(nèi)相互震蕩，涌浪波能損失，波高隨著能量的損失開始逐漸降低。雖然還有后續(xù)殘余的波浪向大壩方向傳播，但是水面逐漸趨于穩(wěn)定。

將相同工況下的監(jiān)測(cè)點(diǎn)同水工物理模型結(jié)果進(jìn)行對(duì)比，滑塊入水的瞬間作為開始計(jì)時(shí)點(diǎn)，數(shù)值模型與物理試驗(yàn)對(duì)比結(jié)果如下。

由數(shù)值模擬計(jì)算值與水工物理模型實(shí)測(cè)值對(duì)比(圖4)可以得出：①監(jiān)測(cè)點(diǎn)M2的首浪數(shù)據(jù)明顯高于其他監(jiān)測(cè)點(diǎn)數(shù)據(jù)，說明涌浪首浪對(duì)河道凹岸沖擊較大；②監(jiān)測(cè)點(diǎn)離落水點(diǎn)越遠(yuǎn)，誤差值越小；③各監(jiān)測(cè)點(diǎn)數(shù)值模擬計(jì)算值總體稍大于水工物理模型試驗(yàn)值，這主要是由于數(shù)值模擬中多采用理想條件計(jì)算，忽略水槽模型底部以及邊壁流體撞擊損失，忽略了摩擦阻力對(duì)于水流的影響。

圖4 各監(jiān)測(cè)點(diǎn)水工模型和數(shù)值模擬涌浪高度對(duì)比Fig.4 Comparison of surge height between hydraulic model and numerical simulation at each monitoring point

數(shù)值模擬可以直觀的顯示涌浪在彎曲河道中的傳播過程，且與水工物理模型試驗(yàn)的涌浪高度變化過程基本一致，數(shù)值上有一定誤差，可為彎曲河道傳播過程研究提供參考。

3.4 不同因素對(duì)壩前涌浪過程的影響分析

根據(jù)水工物理模型試驗(yàn)?zāi)M內(nèi)容：用控制變量法進(jìn)行多次物理試驗(yàn)后，得出該物理模型試驗(yàn)最主要的影響因素從大到小有：水深、滑塊入水角度等，則三維數(shù)值模擬主要從水深和滑塊入水角度兩個(gè)方面進(jìn)行驗(yàn)證水工物理模型的準(zhǔn)確性。

水工物理模型試驗(yàn)并沒有研究水深和滑塊入水角度對(duì)壩前的涌浪高度的變化規(guī)律。通過改變河道水深和滑塊下滑坡度來研究壩前涌浪變化規(guī)律。分析壩前水位變化得到滑坡涌浪對(duì)大壩最危險(xiǎn)的工況。

3.4.1 滑坡體角度對(duì)壩前涌浪過程影響分析

研究同水深不同入水角度對(duì)于壩前水位的影響，以初始水位為0.7 m，入射角度分別為40°、50°、60°、70°、80°來探究滑坡涌浪的水位變化規(guī)律。不同角度壩前水位隨時(shí)間變化如圖5所示(涌浪到達(dá)壩前為計(jì)時(shí)點(diǎn))。

圖5 不同角度壩前涌浪高度隨時(shí)間的變化線Fig.5 Variation line of swell height with time in front of the dam at different angles

通過對(duì)比分析不同角度壩前水位隨時(shí)間的變化線可以得出以下規(guī)律和結(jié)論：①可以發(fā)現(xiàn)不同滑塊下滑角度產(chǎn)生涌浪傳播到重力壩前，壩前水位線隨著時(shí)間變化的規(guī)律大體上是相似的。涌浪到達(dá)壩前，大量水體堆積處形成駐波，迅速抬高壩前水域高度，滑坡涌浪迅速達(dá)到峰值。在2 s時(shí)涌浪達(dá)到最大值。涌浪受到重力壩的阻礙在壩前引起爬高，受大壩反射，壩前水體與后續(xù)傳來的水體的動(dòng)能和沖擊波出現(xiàn)振蕩撞擊，使得壩前水面處的動(dòng)能和水面以下的沖擊波波能不斷發(fā)生轉(zhuǎn)化，當(dāng)后續(xù)能量消散后，水面開始恢復(fù)靜止。通過分析可知：初始水體的動(dòng)能和沖擊波能最大，涌浪的爬高最大。②壩前水位的變化與滑坡體下落角度之間不是線性關(guān)系，引起壩前水位最大的工況是滑塊下滑角度為60°，該工況下壩前涌浪的最大波高0.718 m，最深波谷為 0.678 m。

3.4.2 初始水位對(duì)壩前涌浪過程影響分析

由滑塊角度對(duì)壩前涌浪過程影響分析可知以60°為入射角所產(chǎn)生的壩前波高最大，則以60°為入射角，壩前水深為0.5、0.6、0.7、0.8 m不同水深條件下滑坡涌浪對(duì)于壩前水位的影響(涌浪到達(dá)壩前為記時(shí)點(diǎn))。

由不同水深壩前水位隨時(shí)間的變化線(圖6)可以看出壩前水位在這4種工況下，當(dāng)涌浪傳播到壩前時(shí)受到大壩的阻擋，然后隨之產(chǎn)生一個(gè)駐波，壩前涌浪爬高也最大，這幾種工況產(chǎn)生最大波峰的時(shí)間點(diǎn)比較接近。其中初始水位0.7 m時(shí)，壩前爬高最大，為0.016 m。在產(chǎn)生最大波峰后，后續(xù)水位變化規(guī)律就相差較大。產(chǎn)生這種情況的原因是當(dāng)水槽內(nèi)的初始水位較低時(shí)，滑坡體下落產(chǎn)生的機(jī)械能大部分會(huì)消耗在撞擊河床上，所以只有一小部分的機(jī)械能轉(zhuǎn)化成水體的能量，所以當(dāng)初始水位較低時(shí)，涌浪傳播到壩前產(chǎn)生的爬高較小。當(dāng)水槽內(nèi)的初始水位過高時(shí)，滑塊機(jī)械能將全部轉(zhuǎn)化為水體能量，水體體積過大將會(huì)稀釋滑塊機(jī)械能量，壩前波高又會(huì)降低。故最危險(xiǎn)工況為初始水位為0.7 m，入射角為60°。

圖6 不同水深壩前涌浪高度隨時(shí)間的變化線Fig.6 Variation line of swell height with time in front of dam with different water depths

3.5 滑坡涌浪對(duì)壩體動(dòng)水壓力作用分析

大壩承受的水壓力由靜水壓力和涌浪形成的動(dòng)水壓力組成。根據(jù)Flow-3D軟件數(shù)值模擬的計(jì)算結(jié)果，分析最危險(xiǎn)工況的動(dòng)水壓強(qiáng)隨著浪高變化的分布規(guī)律。當(dāng)涌浪傳播到壩前，造成壩前水面的大幅度波動(dòng)，從而引起壩前動(dòng)水壓力的不斷變化。通過觀察數(shù)據(jù)可知涌浪的動(dòng)水壓力ΔP與涌浪高度η的關(guān)系密切，所以將動(dòng)水壓力換算成動(dòng)水頭Δh=ΔPy，y=1/(ρg)。

由最危險(xiǎn)工況下壩前動(dòng)水頭與涌浪高度數(shù)值模擬結(jié)果隨時(shí)間變化(圖7)可以看出動(dòng)水頭和涌浪高度隨時(shí)間變化規(guī)律。動(dòng)水頭和涌浪高低的變化規(guī)律和波動(dòng)性基本一致，隨著壩前涌浪高度的增加，其動(dòng)水頭也逐漸增加，波峰時(shí)產(chǎn)生的動(dòng)水頭最大。在涌浪波峰時(shí)動(dòng)水頭高度小于涌浪高度，在波谷時(shí)動(dòng)水頭高度大于涌浪高度。

圖7 壩前動(dòng)水頭與涌浪高度數(shù)值模擬結(jié)果隨時(shí)間變化Fig.7 Time-histories of simulated dynamic water head in front of dam and surge height

從各波峰時(shí)段折減系數(shù)數(shù)值模擬結(jié)果隨時(shí)間的變化曲線(圖8)可以得到以下幾點(diǎn)結(jié)論：①這3個(gè)時(shí)段中t=5 s時(shí)，壩前涌浪的折減系數(shù)最大，t=7 s時(shí)折減系數(shù)最小，而t=2 s時(shí)的折減系數(shù)介于兩者之間，這說明在t=2 s時(shí)，涌浪剛傳播到壩前，此時(shí)涌浪相互作用相對(duì)較??；而t=5 s時(shí)，涌浪波動(dòng)作用較劇烈，而壩前水位和動(dòng)水頭差值較大；當(dāng)t=7 s時(shí)，折損系數(shù)最小，說明涌浪的能量在前一段時(shí)間的折減較大，此時(shí)涌浪的高度也相對(duì)較小，折減系數(shù)隨壩前水位降低而變小。② 3個(gè)時(shí)段的折損系數(shù)分布規(guī)律類似，都隨著壩前水位的減小先增加后減小，上述分析可知壩前動(dòng)水頭沿著水深分布有著一定程度的折減，所以當(dāng)直接采用涌浪水面的靜水頭作為大壩水壓力計(jì)算時(shí)，結(jié)果更加安全。

圖8 壩前折減系數(shù)數(shù)值模擬結(jié)果隨壩高分布Fig.8 Numerical results of change of reduction coefficient in front of dam with dam height

4 工程實(shí)例

某高壩庫(kù)區(qū)滑坡區(qū)位于大壩上游約1 300 m。岸坡為向河床突出的山梁，下游河流為60°折角型河流。岸坡巖體為花崗巖，滑坡范圍位于兩溝之間?；虑熬壐叱?06 m，后緣高程292 m，相對(duì)高差86 m，坡角在30°～60°之間。發(fā)生滑坡時(shí)，壩前為蓄水時(shí)期，水位高程206.2 m，河谷呈“U”型，平均水深70 m，河道寬度約為400 m?；麦w厚度為30～40 m，崩滑體長(zhǎng)度約為60 m?；滤俣燃s為15 m/s。

根據(jù)滑坡體與大壩位置示意圖(圖9)可以看出實(shí)際滑坡涌浪傳播與本模型極為相似，故為了讓數(shù)值模型更加符合實(shí)際模型，同時(shí)為了更好地觀察和分析數(shù)據(jù)，將采用實(shí)際尺寸的彎曲河道模型進(jìn)行三維數(shù)值模擬。由于某高壩實(shí)際初始水位和滑塊尺寸的數(shù)據(jù)均比水工物理模型試驗(yàn)對(duì)應(yīng)的三維數(shù)值模型數(shù)據(jù)擴(kuò)大了約100倍，故將三維數(shù)值模型進(jìn)行100倍擴(kuò)大。根據(jù)實(shí)際滑坡數(shù)據(jù)，研究滑塊體在下滑角度=50°的工況下，初始水位為70 m，滑坡體入水速度取15 m/s?；瑝K體的體積取60 m×30 m×30 m，根據(jù)實(shí)際山體表層平均材料密度為1 900 kg/m3，認(rèn)為滑坡體在下滑過程中為剛體運(yùn)動(dòng)且只受重力作用，不發(fā)生變形和破壞，重力加速度為9.81 m/s2，不考慮空氣阻力，時(shí)間步長(zhǎng)為0.1 s，計(jì)算時(shí)間為50 s，通過數(shù)值模擬可以獲得以初始水位70 m相對(duì)應(yīng)監(jiān)測(cè)點(diǎn)M′1、M′2、M′3、M′4(與圖2監(jiān)測(cè)點(diǎn)位置相同)位置的涌浪高度和水位變化情況。

圖9 滑坡體與大壩位置示意圖Fig.9 Schematic diagram of the location of landslide body and building

由圖10各監(jiān)測(cè)點(diǎn)M′1、M′2、M′3水位與試驗(yàn)數(shù)值模型初始水位為0.7 m、入射角為50°相對(duì)應(yīng)監(jiān)測(cè)點(diǎn)M1、M2、M3水位起伏變化規(guī)律基本一致，而監(jiān)測(cè)點(diǎn)M′4與M4略有不同，工程實(shí)例數(shù)值模擬監(jiān)測(cè)水位相比于試驗(yàn)數(shù)值模擬監(jiān)測(cè)水位起伏頻率偏大。

圖10 各監(jiān)測(cè)點(diǎn)涌浪高度數(shù)值模擬結(jié)果Fig.10 Simulated surge height at each monitoring point

由初始水位70 m相對(duì)位置的壩前涌浪高度及壩前動(dòng)水頭(圖11)可知，隨著涌浪高度的增大動(dòng)水頭也隨之增加，波峰時(shí)產(chǎn)生的動(dòng)水頭最大。壩前涌浪波動(dòng)效果與水工物理模型初始水位為0.7 m，入射角為50°的水位波動(dòng)效果略有不同，但動(dòng)水頭波動(dòng)規(guī)律基本一致。

圖11 壩前動(dòng)水頭與涌浪高度隨時(shí)間變化Fig.11 Variations of moving water head in front of dam and surge height with time

由數(shù)值模擬結(jié)果中獲取在波峰時(shí)(t=10、22、40 s)折減系數(shù)沿大壩高程情況如圖12所示。

圖12 壩前折減系數(shù)隨壩高分布Fig.12 Change of reduction coefficient in front of dam with dam height

由于水工物理模型尺寸過小，涌浪波動(dòng)頻率過高則折減系數(shù)明顯高于實(shí)際尺寸模型折減系數(shù)，但是兩次折減系數(shù)沿高程分布規(guī)律基本一致，折減系數(shù)都是隨高程的降低先增后減。

5 結(jié) 論

本文基于庫(kù)岸滑坡涌浪首浪高度試驗(yàn)研究，基于Flow-3D模擬技術(shù)，得出以下3點(diǎn)結(jié)論。

(1)數(shù)值模擬涌浪高度、水面高度起伏過程和物理模型試驗(yàn)結(jié)果相符。因?yàn)閿?shù)值模擬多采用理想條件計(jì)算，導(dǎo)致數(shù)值上有少許誤差，但是數(shù)值模擬還是能很好反映涌浪傳播過程與大壩相互作用的過程

(2)通過滑坡體下滑的角度變化以及模型水槽的水深這兩個(gè)變量的影響比較大，以0.7 m相同水深不同入射角和以相同入射角60°不同水位分析研究。最后分析得出：水深0.7 m、下滑角度60°的工況計(jì)算得出涌浪對(duì)壩體作用效應(yīng)最明顯。

(3)通過工程實(shí)例的驗(yàn)證分析得出：滑坡涌浪形成的動(dòng)水頭與涌浪高度和波動(dòng)頻率有關(guān)，涌浪高度和動(dòng)水頭的波動(dòng)性基本一致。波峰時(shí)形成的動(dòng)水頭最大，動(dòng)水頭沿大壩高程分布有一定的折減，折減系數(shù)隨壩高均為先增后減。因此采用最大涌浪爬高水頭壓力計(jì)算分析壩體穩(wěn)定應(yīng)力，其結(jié)果偏安全。