鄧成進，黨發(fā)寧，陳興周

（1.西安理工大學(xué) 巖土工程研究所，陜西 西安 710048；2.中國電建集團 西北勘測設(shè)計研究院有限公司，陜西 西安 710065；3.西安科技大學(xué)，陜西 西安 710065）

1 研究背景

水庫庫區(qū)存在的滑坡、變形體、崩滑體等不良地質(zhì)體，在不利因素影響下可能失穩(wěn)高速滑入水庫并產(chǎn)生巨大涌浪，危及大壩及下游安全，因此研究庫區(qū)滑坡涌浪與大壩作用過程及其機理有著重要意義［1-2］。

庫區(qū)滑坡滑動過程的滑速和涌浪計算受多種復(fù)雜因素影響，經(jīng)驗公式計算成為一種主要的、簡便的手段，常用的經(jīng)驗公式有Noda E、潘家錚法、水科院法、美國土木工程學(xué)會法等［3］。任興偉在潘家錚法的基礎(chǔ)上提出了改進的經(jīng)驗公式，并用于計算新灘滑坡的初始涌浪高度［4］。水工物理模型試驗也是分析庫區(qū)滑坡涌浪及其作用機理的常用方法，比如李家峽庫區(qū)1 號滑坡、龍羊峽庫區(qū)近壩庫岸、三峽庫區(qū)龔家方崩滑體等工程均采用了大尺寸的水工物理模型試驗研究滑坡涌浪對大壩安全的影響［5-7］。龐昌俊通過二維滑坡涌浪試驗研究和探索了斜坡滑坡涌浪與水平、垂直滑坡涌浪的差異和關(guān)系［8］；岳書波通過試驗研究了滑塊在水槽中形成涌浪的形態(tài)及衰減規(guī)律［9］；這些試驗總結(jié)了一些對滑坡涌浪問題的基本認(rèn)識。大尺寸水工模型試驗要獲得完整的試驗數(shù)據(jù)需付出較大的代價，尤其滑體高速滑動所帶來涌浪劇烈變化而難以獲取準(zhǔn)確數(shù)據(jù)，限制了其推廣使用。

由于流體力學(xué)Navier-Stokes 方程式求解復(fù)雜，當(dāng)前一般從2 個方面來簡化實現(xiàn)滑坡涌浪過程的數(shù)值模擬。（1）簡化流體力學(xué)的方程模型。比如：鐘登華等直接采用前述水科院經(jīng)驗公式，用三維可視化動態(tài)模擬技術(shù)分析了某水庫滑坡失穩(wěn)入庫、激起涌浪的全過程［10］；周桂云通過自編程序GEO-FLOW 將淺水控制方程應(yīng)用于模擬滑坡涌浪的傳播動態(tài)過程［11］；黃波林采用基于水波動力學(xué)的FAST 軟件，通過簡化初始涌浪形成過程及涌浪源模型，對三峽庫區(qū)茅草坡滑坡涌浪災(zāi)害、龔家方崩滑體滑坡涌浪進行了模擬和預(yù)測［12-13］。（2）采用簡化的計算模型，即在二維剖面及小范圍的三維水槽中來研究滑坡涌浪?？娂獋惖冗\用SPH 方法實現(xiàn)了模擬涌浪產(chǎn)生的過程［14-15］；Viroulet 等采用二維SPH 法和有限體積法分別計算二維楔形滑坡涌浪規(guī)律，并采用水工物理模型試驗驗證［16］；李靜等采用SPH 法研究了小范圍三維矩形規(guī)則水庫中涌浪傳播過程以及涌浪對大壩壩面的沖擊力，分析了滑坡體寬度、入水速度以及庫區(qū)水深對壩面沖擊力的影響［17］。這些簡化方法對于研究滑坡涌浪形成、傳播過程及機理起到一定的促進作用；但無法考慮高壩大庫的實際河谷地形，無法真實模擬滑坡涌浪形成及傳播的全過程；尤其對于大范圍的高壩庫區(qū)數(shù)值模擬，模擬結(jié)果往往存在較大的誤差。且很少有考慮下游大壩對庫區(qū)滑坡涌浪傳播的影響，以及滑坡涌浪與大壩相互作用的過程。

隨著計算技術(shù)的發(fā)展，采用有限差分法進行水流、水力學(xué)模擬和分析在國內(nèi)外得到廣泛應(yīng)用［18-19］，但滑坡涌浪三維數(shù)值模擬的應(yīng)用仍在探索階段，其模擬的精度及可信度一直是關(guān)注的重點問題。本文建立了某高壩庫區(qū)3#變形體的三維滑坡涌浪數(shù)值模型，并結(jié)合水工物理模型試驗對其模擬的精度和有效性進行驗證，為高壩庫區(qū)滑坡涌浪數(shù)值模擬的應(yīng)用提供科學(xué)依據(jù)；并通過數(shù)值模擬獲得滑坡涌浪形成和傳播過程的完整數(shù)據(jù)，分析和研究滑坡涌浪與大壩相互作用機理。

2 庫區(qū)滑坡涌浪的水工物理模型試驗

某高壩庫區(qū)3#變形體位于大壩上游約1300 m（圖1）。岸坡呈向河床凸出的弧形山梁，走向約NE30°，下游河流流向NE45°左右。岸坡范圍為花崗巖體，溝梁相間，溝間為從岸頂?shù)狡履_連續(xù)延伸的山梁，3#變形體即位于兩溝之間，溝底一般有基巖出露。

圖1 3#變形體與建筑物位置示意圖

自水庫蓄水后3#變形體傾倒變形大，變形條件及成因機制復(fù)雜，有必要分析極端情況下整體失穩(wěn)下滑所產(chǎn)生的涌浪對大壩及下游安全的影響。為分析和研究滑坡涌浪對大壩安全的影響，建立了3#變形體庫區(qū)滑坡涌浪的水工物理試驗?zāi)Ｐ鸵妶D2。

圖2 水工物理模型試驗

試驗過程參照《滑坡涌浪模擬技術(shù)規(guī)程》［20］，模型滿足幾何相似、水流運動相似和動力相似，遵循佛勞德相似準(zhǔn)則。滑體滿足幾何相似和塊體的比重相似，河道和岸坡應(yīng)滿足阻力相似，幾何比尺為λL=280，相應(yīng)流量比尺為壓力比尺為λp=λL=280，時間比尺為庫區(qū)上游有足夠的長度，以避免上游反射波較快折回與原生波向下游傳遞時產(chǎn)生疊加。

由于邊坡表部巖體破碎，為方便觀測滑體落入水中后的堆積形狀，滑塊根據(jù)試驗方量由若干個標(biāo)準(zhǔn)塊體堆積組成?；麦w置于自制滑車上，滑車制成多節(jié)鉸式，置于滑床上面，滑坡體可隨滑車附貼著滑床下滑并取得較高的滑速；試驗時由起吊、擒縱裝置控制起吊及下放滑車，放置不同高度下滑可獲得相應(yīng)的入水速度。涌浪高度采用中國水科院水力學(xué)所研制的DJ800 型水工數(shù)據(jù)采集系統(tǒng)收集。水工模型試驗開展了3#變形體在不同方量、不同滑速、不同水位條件的涌浪試驗。試驗在沿岸及大壩共布置了6 個涌浪測點（圖1），1#—4#測點布置在左岸，基本代表涌浪沿岸傳播規(guī)律及過程；4#—6#測點布置于壩前，代表涌浪與大壩相互作用過程。在正常蓄水位條件下，當(dāng)滑體全部入水，滑速為24.6 m/s 時，各測點涌浪高度隨時間變化曲線見圖3所示。

圖3 各測點浪高隨時間變化曲線

由圖3所示，各測點涌浪高度隨時間不斷動蕩、起伏，在t=160 s 以后開始逐漸衰減趨于平緩。1#測點離落水點最近，其首浪高度最大；由于大壩的阻礙作用，4#測點處首浪高度比2#和3#測點的大，表明首浪高度隨著離落水點的距離先減小后增加。壩前4#、5#和6#測點處水面起伏過程基本一致，河床中部（5#測點）首浪高度比兩岸壩肩處（4#和6#測點）略小。壩前首浪高度在t=40～60 s 之間將大于壩頂超高，發(fā)生漫壩；隨后水面反復(fù)震蕩，在t=150～160 s 之間再次漫壩，此時兩岸壩肩處（4#和6#測點）形成的涌浪最高。

3 基于流體動力學(xué)的滑坡涌浪數(shù)值模擬

3.1 基本理論將水流運動視為不可壓縮的黏性流體運動，基于流體動力學(xué)采用離散完整的Navier-Stokes 方程式，通過體積分?jǐn)?shù)和面積分?jǐn)?shù)來定義網(wǎng)格中的障礙物，從而矩形網(wǎng)格能夠精確地模擬復(fù)雜的幾何形狀。控制方程在描述上與經(jīng)典方程有所不同，控制方程中含有體積分?jǐn)?shù)和面積分?jǐn)?shù)參數(shù)［21］。

Flow3D 采用有限差分法進行數(shù)值離散，將模型離散成空間矩形網(wǎng)格，而非常規(guī)的四面體或六面體網(wǎng)格，速度和面積分?jǐn)?shù)參數(shù)位于矩形網(wǎng)格邊界面的中心點上。求解方法采用廣義的極小殘差算法，其計算收斂速度相對較快，計算精度高，該方法在計算流體力學(xué)中得到了廣泛的應(yīng)用。

計算采用VOF 法追蹤自由水面［22］，引入了流體體積函數(shù)αq的概念，用于定義單元內(nèi)第q 相流體所占體積與單元體積的比值，空氣體積函數(shù)αq=0，流體內(nèi)部體積函數(shù)αq=1，自由液面的體積函數(shù)為αq=0.5，實現(xiàn)了對自由面的追蹤。若0＜αq＜1，則表示該單元為交接面單元。

RNG k-ε湍流模型可以很好地模擬高應(yīng)變率和流線彎曲程度較大的流動?；w入庫及庫區(qū)內(nèi)涌浪相互作用時會導(dǎo)致水面出現(xiàn)巨大的變形、破碎，適合采用RNG k-ε湍流模型計算模擬滑坡涌浪的產(chǎn)生及其傳播、作用過程。

3.2 庫區(qū)三維滑坡涌浪數(shù)值模型將3#變形體簡化為矩形剛性滑體，模擬在極端情況下整體滑動入水?；w整體模型順河向?qū)?20 m，厚度為26 m，落水點水面寬度約為530 m，庫區(qū)滑坡涌浪三維模型見圖4所示。三維模型上游、左右岸及下游大壩采用固壁邊界條件。模擬計算時間步△t=1.0 s，總計算時間步為160 s。

滑體滑動和堆積形態(tài)受河谷地形、以及水阻力的作用影響，且相互作用機理較為復(fù)雜，目前滑坡運動過程和堆積形態(tài)的模擬研究尚未有成熟辦法［23］。結(jié)合水工物理模型試驗的滑塊速度成果，在FLOW3D 中自定義滑體滑動速度與時間的關(guān)系，以控制滑體運動過程，模擬滑體沖擊水面及初始涌浪形成過程。

初始涌浪形成及涌浪爬坡過程中，水面發(fā)生劇烈變化，翻卷形成巨大的浪花，網(wǎng)格劃分的疏密程度對水面變化劇烈區(qū)域的計算精度影響較大，但網(wǎng)格劃分過于精細(xì)會造成數(shù)值模擬成本巨大，可能導(dǎo)致計算無法收斂［24］。對于高壩大庫滑坡涌浪的三維數(shù)值模擬計算，為了提高計算精度，大幅減少網(wǎng)格數(shù)量，在水體劇烈運動的區(qū)域采用相對精細(xì)的網(wǎng)格，與其他區(qū)域的各網(wǎng)格塊進行搭接。計算模型沿高程及河道分為9 個網(wǎng)格塊相鏈接，滑體附近的河段及中部（正常水位附近范圍內(nèi)）網(wǎng)格適當(dāng)精細(xì)，間距為2.5 m，其他網(wǎng)格間距5.0 m，模型網(wǎng)格化后共有單元約3325 萬個，有效單元約1108 萬個。

3.3 滑坡涌浪的形成及傳播過程分析通過對庫區(qū)滑坡涌浪整個過程的分析，可將滑坡涌浪形成及傳播過程分為三個階段：①初始涌浪形成及向?qū)Π秱鞑ルA段；②沿岸傳播階段；③涌浪與大壩相互作用階段。初始涌浪形成過程落水點附近水面高程云圖見圖5所示，在t=21～57 s 各時刻庫區(qū)水面高程等值線圖見圖6。

圖4 庫區(qū)滑坡涌浪計算三維模型

（1）初始涌浪形成及向?qū)Π秱鞑ルA段。隨著滑體滑入水庫，滑體擠壓水體形成初始涌浪，涌浪由落水點向四周水域傳播。在t=5～7 s 時由于滑體高速沖擊水面，附近的水體受擠壓后水面翻卷，迅速形成巨大浪花，表現(xiàn)為躍沖形態(tài)，形成了初始浪高。隨著初始涌浪形成，在水面形成了一個以落水點為中心的波峰圈向四周擴散、傳播。

（2）沿岸傳播階段。涌浪傳播至對岸后，在岸壁阻擋和涌浪推進的雙重擠壓作用下，左岸附近水域水面呈整體抬升，在岸坡上爬升至最高后水體回落，并與再次形成的波峰相互作用、震蕩。在落水點附近水面形成一系列涌浪傳播圈，庫區(qū)水體也隨著升起、降落，推動著涌浪繼續(xù)向壩前傳播。在t=33 s 時涌浪傳播至壩前水域，壩前水面開始緩慢上升。

（3）涌浪與大壩相互作用階段。涌浪傳播至壩前水域后，由于大壩及庫岸阻礙作用，壩前水域的水面繼續(xù)抬升，并逐漸在兩岸壩肩位置形成較高的水面，在t=44 s 時兩岸壩肩處水面高程達(dá)到壩頂高程2460.00 m；在t=52 s 時左岸壩肩處水面最大高程為2464.10 m，涌浪高度達(dá)12.10 m，超過壩頂約4.10 m，水體將漫過壩頂。隨后壩前水域水面反復(fù)震蕩，水體與大壩相互作用、疊加將形成更高的涌浪。

3.4 與水工物理模型試驗結(jié)果對比分析將各測點處涌浪高度變化與水工物理模型試驗結(jié)果進行對比見圖7。由圖7可得出以下規(guī)律：①1#—6#測點水工物理模型試驗和數(shù)值模擬的涌浪高度變化及起伏過程基本一致。其中1#—4#測點處水工物理模型試驗和數(shù)值模擬的水面起伏變化過程吻合程度較高；而5#、6#測點誤差較大，兩者吻合程度略有降低，但兩者波動相似性依然存在。②各測點水工物理模型試驗的首浪高度比數(shù)值模擬結(jié)果的大，計算存在一定誤差。③測點離落水點的距離越遠(yuǎn)，兩種方法得出涌浪高度的差別越??；隨著涌浪的傳播及能量消耗，后期兩種方法得出涌浪高度的差別也越小。

數(shù)值模擬采用VOF 法追蹤自由水面，網(wǎng)格劃分精細(xì)程度對模擬涌浪浪花翻滾的清晰度影響較大。當(dāng)數(shù)值模擬網(wǎng)格精細(xì)程度越高，就越能清晰的模擬出涌浪浪花翻滾、濺飛的微觀形態(tài)，其計算的涌浪高度就大，即越接近試驗或真實的水面高度。高壩大庫三維數(shù)值模擬的網(wǎng)格卻很難到達(dá)一定的精細(xì)程度，難以精確模擬滑體高速沖擊水面、初始涌浪形成及涌浪傳播過程中浪花翻滾的微觀形態(tài)，因此各測點處數(shù)值模擬的首浪高度均比水工物理模型試驗的小。隨著涌浪傳播，涌浪能量損耗之后，后期涌浪浪花高度逐漸減小，此時數(shù)值模擬與水工物理模型試驗的涌浪高度也基本一致。因此，基于完整的Navier-Stokes 方程，采用有限差分法可實現(xiàn)庫區(qū)滑坡涌浪的全過程，基本滿足精度要求；但初始涌浪形成過程的浪花形態(tài)及對岸涌浪爬升等微觀形態(tài)的模擬尚需采用更精細(xì)網(wǎng)格，以及付出更高的計算代價。

數(shù)值模擬可直觀展示滑坡涌浪傳播及其與大壩相互作用的整個過程，且與水工物理模型試驗的涌浪高度變化趨勢基本一致，數(shù)值上有一定誤差，可為庫區(qū)滑坡涌浪災(zāi)害預(yù)測和研究提供依據(jù)。

圖7 各監(jiān)測點數(shù)值模擬和模型試驗涌浪高度對比

4 壩前涌浪及其作用機理分析

4.1 壩前涌浪高度分析滑坡涌浪傳播的整個過程表明，由于大壩及庫岸的阻礙，壩前水面先呈整體緩慢抬升，首浪傳播至壩前時左岸壩肩處形成的涌浪最高。由4#測點（左岸壩肩）獲得的數(shù)據(jù)表明（圖7（d）），在t=44 s 時首浪超過壩頂；隨著壩前水體相互作用、震蕩，在t=151 s 第5 次波峰時，將形成水面最高。大壩左岸壩肩處在t=38～58 s和t=143～158 s各時刻涌浪形態(tài)變化過程分別見圖8和圖9。

圖8 在t=38～57s 各時刻左岸壩肩處水面高程

由圖8可知，在t=44 s 時壩前涌浪達(dá)到壩頂高程2460.00 m，t=52 s 時達(dá)到最大2464.10 m，t=58 s時水面下降回落至壩頂高程，表明首浪可能在t=44～58 s 時左壩肩處緩慢的漫過大壩。

由圖9可知，在第5 次波峰來臨前，壩前水面形成了巨大落差，在t=148 s 時涌浪迅速沖上壩頂，形成巨大波峰，t=151 s 水面最大高程為2468.36 m，此時超過壩頂高程約8.36 m，比首浪更高。上述分析表明由于庫區(qū)水面的反復(fù)震蕩、疊加，壩前第5 次波峰形成的涌浪最高；與柘溪塘巖光（距大壩1550 m）滑坡涌浪的第2 次波峰達(dá)到最高的情況類似。

圖9 在t=143～158s 各時刻左岸壩肩處剖面水面高程

滑坡涌浪傳播至壩前，部分水體將在兩岸壩肩處漫過壩頂，隨后壩前水面下降，與后期傳播至此的波峰相互碰撞，水面反復(fù)震蕩、疊加將形成更高的涌浪。因此，有必要在壩頂及大壩下游建立避險空間，及時預(yù)警。

4.2 滑坡涌浪對壩體水壓力作用分析隨著滑坡涌浪與大壩相互作用，滑坡涌浪對壩體作用力也不斷變化，因此分析壩面水壓力變化規(guī)律對大壩穩(wěn)定及應(yīng)力有著十分重要的意義。大壩承受的水壓力由靜水壓力P 和滑坡涌浪形成的動水壓力ΔP 組成，涌浪產(chǎn)生的動水壓力ΔP 與涌浪高度η 密切相關(guān)，故將各點的動水壓力換算為動水頭Δh=Δp γ，分析動水頭Δh 與涌浪高度η 隨時間的變化規(guī)律。通過數(shù)值模擬獲得正常水位以下202 m 深處的壩面動水頭及壩前涌浪高度隨時間的變化規(guī)律，并與水工物理模型試驗的監(jiān)測數(shù)據(jù)進行對比（圖10）。由圖10可知，數(shù)值模擬和水工物理模型試驗得出的動水頭隨涌浪高度的變化規(guī)律基本一致。動水頭和涌浪高度的波動性基本一致，隨著壩前涌浪高度增加，其動水頭也逐漸增加，波峰時涌浪產(chǎn)生的動水頭最大。在涌浪波峰、波谷處，動水頭與涌浪高度的大小關(guān)系有所差別，表現(xiàn)為波峰時Δh ＜η，而波谷時Δh ＞η。

圖10 動水頭Δh 與涌浪高度η 隨時間變化曲線

動水頭不僅與涌浪高度相關(guān)，還與涌浪的波動頻率、水深有關(guān)。定義波峰時動水頭與涌浪高度之間的折減系數(shù)為由于模型試驗的監(jiān)測點有限，在數(shù)值模擬結(jié)果中獲取在波峰時（t=45 s、85 s、158 s）折減系數(shù)沿大壩高程分布情況見圖11。

圖11 數(shù)值模擬計算折減系數(shù)隨壩高分布情況

由圖11可得出以下結(jié)論：①在t=158 s 時動水頭折減系數(shù)最大，而此時涌浪波動相對劇烈（圖10（b）），表明涌浪波動頻率越高，動水頭與涌浪高度之間的相對差值越大，折減系數(shù)越大。②折減系數(shù)沿高程的分布規(guī)律基本一致，折減系數(shù)隨著高程的降低先增后減，在0.3～0.4 倍壩高位置折減系數(shù)最大。③在t=158 s 時，由于涌浪波動頻率較高，折減系數(shù)增大，其最大值出現(xiàn)的位置也相應(yīng)抬高。上述分析表明，波峰時涌浪形成的動水頭最大，動水頭沿高程分布有一定程度折減，因此當(dāng)直接采用涌浪水面的靜水頭作為大壩水壓力計算時，結(jié)果偏于安全。

5 結(jié)論

本文以某高壩庫區(qū)3#變形體整體滑動激起的涌浪為例，建立庫區(qū)滑坡涌浪三維數(shù)值模型，并通過水工物理模型試驗對數(shù)值模擬結(jié)果進行驗證，研究了涌浪與大壩相互作用過程，得出以下結(jié)論：

（1）數(shù)值模擬得出的涌浪高度變化、水面起伏過程與水工物理模型試驗的結(jié)果基本一致，但數(shù)值上有一定誤差。隨著后期涌浪傳播及能量消耗，數(shù)值模擬與水工物理模型試驗得出的涌浪與大壩相互作用過程及高度基本一致。數(shù)值模擬能較好地反映涌浪傳播及其與大壩相互作用過程，可作為高壩大庫庫區(qū)滑坡涌浪災(zāi)害預(yù)測和預(yù)警的依據(jù)。

（2）滑坡涌浪傳播至壩前在t=44～58 s 時兩岸壩肩處漫過壩頂。在庫區(qū)水面反復(fù)震蕩、疊加作用下，在t=151 s 第5 次波峰時壩前涌浪最高，超過壩頂8.36 m，而后逐漸衰減。從涌浪災(zāi)害防治角度看，由于水體反復(fù)震蕩、疊加而形成的涌浪可能會對下游沿岸和大壩安全造成更大危害。因此，有必要在壩頂及下游采取措施，建立預(yù)警系統(tǒng)和避險空間，在極端情況下及時預(yù)警、避險，直至涌浪衰減。

（3）高壩大庫滑坡涌浪形成的動水頭與涌浪高度、波動頻率、水深有關(guān)。動水頭與涌浪高度的波動性基本一致；波動頻率越高，其與涌浪高度的差別也越大。波峰時涌浪形成的動水頭最大，但動水頭沿高程分布有一定程度的折減，折減系數(shù)隨著高程的降低先增后減。因此當(dāng)采用靜力方法計算分析壩體穩(wěn)定應(yīng)力時，其計算結(jié)果偏于安全。