彎曲河道型庫岸滑坡涌浪傳播及其與大壩的相互作用

彭 輝，黃亞杰

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

1 研究背景

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

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

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

2 彎曲河道型庫岸滑坡涌浪的物理模型實驗簡介

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

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

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

3.1 基本理論

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

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

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

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

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

模型運動學(xué)邊界條件為

(1)

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

(2)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3.5 滑坡涌浪對壩體動水壓力作用分析

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

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

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

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

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

4 工程實例

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

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

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

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

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

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

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

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

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

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

5 結(jié) 論

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

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

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

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