李 靜，陳健云，徐 強，孫 迅

（大連理工大學 建設工程學部，遼寧 大連 116024）

1 研究背景

汶川地震導致西南震區(qū)發(fā)生了大量的邊坡失穩(wěn)，即使當時未直接失穩(wěn)破壞，但由于強震損傷嚴重，降低了邊坡的抗滑穩(wěn)定能力，近年來在不斷發(fā)生的余震以及強暴雨等作用下，邊坡失穩(wěn)也在不斷發(fā)生。顯然，如果這一災害發(fā)生在庫區(qū)，滑坡體瞬時沖擊庫水將會引發(fā)巨大的涌浪，對擋水壩體造成沖擊壓力作用并威脅下游人民生命和財產安全。

我國西南地區(qū)工程地質條件多為高陡邊坡，滑坡涌浪問題日益突顯出來。因此，對于西南高壩大庫，特別是狹長型河谷水庫，不僅要考慮強震對壩體的損傷破壞，還要考慮強震或暴雨引發(fā)的滑坡涌浪對大壩沖擊造成的災害疊加風險。歷史上多次滑坡涌浪造成了嚴重的結構破壞，如1959年意大利龐特塞拱壩庫區(qū)大約300萬m3的滑坡引發(fā)的漫頂涌浪高達20多米。1961年拓溪水電站大壩上游右岸山體滑坡，造成約165萬m3的土石體以20 m/s的速度涌入水庫，漫過正在施工的溢流壩段壩頂?shù)挠坷烁叨葹?.6 m。1963年意大利瓦伊昂拱壩庫區(qū)內滑坡導致300萬m3的土石沖入200多米水深的狹窄河谷內，形成的涌浪在漫過壩頂時仍有100多米高，在下游距離大壩1.4 km的河口出口處立波仍然高達70 m，沖毀了大部分朗格尼亞鎮(zhèn)。

由于滑坡涌浪的巨大危害，圍繞著滑坡涌浪首浪高度與滑坡體相關參數(shù)之間的關系，基于經驗公式、物理模型試驗以及數(shù)值分析開展了大量的研究。谷建［1］采用潘家錚法分析了某水電站滑坡涌浪不同滑速下的涌浪高度。龐昌俊［2］通過試驗研究了邊坡坡度、體積和入水速度等因素對滑坡涌浪的影響。殷坤龍等［3］根據(jù)體積守恒和水下滑塊運動，將滑坡涌浪分成體積涌浪部分和沖擊涌浪部分。吳佳壕［4］對小灣水電站進行了涌浪分析。徐文杰［5］采用數(shù)值方法研究了滑塊形狀、體積、滑面摩擦系數(shù)、水面寬度等因素對滑坡涌浪及爬坡高度的影響?？娂獋惖龋?］對剛性滑塊入水涌浪進行了分析。目前的研究主要是通過經驗公式對涌浪浪高及其爬坡進行計算，使用數(shù)值模擬對涌浪過程進行自由面追蹤，而對于涌浪對擋水建筑物沖擊作用力方面的研究相對較少。黃錦林［7］通過模型試驗對壩面滑坡涌浪沖擊壓力分布進行了測試。楊艷［8］則通過水工模型試驗，研究了滑坡涌浪對架空直立碼頭的正面波壓。

本文針對滑坡體入水過程的涌浪生成及傳播過程進行數(shù)值模擬，分析不同因素對涌浪要素、傳播過程及壩面沖擊壓力特征的影響。

2 基于SPH方法的滑坡涌浪分析

2.1 SPH方法基本理論光滑粒子水動力學方法（the Smoothed Particle Hydrodynamics，SPH）將分析區(qū)域離散為光滑核函數(shù)加權的粒子（如圖1所示）［9］。將流體動量方程離散成粒子形式，則可以得到流體基于SPH方法的方程：

圖1 SPH分析模型

式中：vi為粒子i的速度；t為時間；N為粒子j支持域內粒子總數(shù)；p為壓力；ρj、mj為粒子j的密度和質量；?為梯度算子；為粒子i和粒子j之間的距離；h為光滑長度；為表征粒子j對粒子i影響的核函數(shù)［10］，滿足歸一化、緊支性及狄拉克函數(shù)等條件，本文選取三次B樣條函數(shù)；Πij為黏性力項。

對任意一個粒子i，本文采取連續(xù)性密度法，將連續(xù)方程離散成粒子形式［11］：

通過式（2）得到的粒子密度決定了粒子分配和光滑長度h的變化。

粒子壓力根據(jù)粒子密度采用流體狀態(tài)方程計算，本文采用Mie-Gruneisen狀態(tài)方程：

式中：ρ0為粒子初始密度；C及S1、S2、S3為與沖擊波波速Us與波后質點粒子速度Up相關的常數(shù)［12］；γ0為Gruneisen常數(shù)；α為γ0的一階體積修正系數(shù)；μ=ρρ0-1；E為初始內能，在常溫下取值為0。

黏性力項通過下式得到［13］:

光滑長度h的選擇關系到計算效率和計算精度，太小會降低計算效率，太大會丟失局部特性。本文通過平均密度來對光滑長度進行調節(jié)變換：

式中：h0為初始的光滑長度；ρ0為初始密度；d為空間維數(shù)。

由于在邊界上的粒子存在缺陷，并且會發(fā)生積分截斷的現(xiàn)象，本文采用耦合邊界法［14］，設置兩層固壁粒子來模擬固壁邊界。積分時間步長需滿足Courant-Friedrichs-Levy（CFL）條件。本文使用簡單經典的一階格式用于積分，時間步長由下式決定：

式中：KCFL為數(shù)值常數(shù)，一般取0.3。

庫區(qū)流體SPH粒子與滑塊有限元單元的耦合采用基于法向罰函數(shù)的點面接觸算法。

2.2 滑坡涌浪數(shù)值分析模型剖面如圖2所示，長300 m，深50 m，滑塊單元尺寸和水體粒子間距均為1 m，圖2中右邊界為固定邊界，上游為無反射邊界。滑塊體尺寸為10 m×20 m×50 m（長×寬×高），滑塊距離左端邊界2 m（如圖2所示）?；瑝K初速度和到達底部速度均為0 m/s。為了探討粒子間距的影響，本文分別對1、2和5 m粒子間距模型進行了模擬。

圖2 模型剖面

數(shù)值模型主要參數(shù)參考文獻［12］中的數(shù)值，具體如表1所示。

表1 水體Gruneisen狀態(tài)方程主要參數(shù)

圖3為涌浪產生和傳播過程。從圖3中可以看出，滑坡體入水后，水體受到滑坡體擠壓，兩側水體迅速升高，右側水體與滑坡體分離，隨后形成反向卷入波，隨滑坡體進入水體的大量空氣在坡體右下側形成空腔，進一步抬高液體表面高度。在滑塊到達水底后，由于重力作用，氣泡發(fā)生破碎，能量耗散下波浪稍有回落向右傳播。隨后回涌的水體與滑塊相撞，再次成浪向右傳播。首浪傳播至右端壩面后爬坡到最大值，然后發(fā)生明顯回落。

以上現(xiàn)象在粒子間距1 m和2 m時都可以觀察到，粒子間距5 m則過于粗糙，反向卷入波和氣體空腔現(xiàn)象不明顯，觀察不出水體的運動變化特性，自由面波形變化也由于粒子過于稀疏不能很好呈現(xiàn)。

圖3 涌浪傳播過程

圖4為滑坡體下落位移及速度隨時間的變化曲線。從圖4可以看出，滑坡體滑入水中后速度逐漸變大，然后逐漸變小。到達底部后，滑塊位移與速度均有些微的震蕩，最終趨于靜止，這個現(xiàn)象更符合實驗結果和實際情況［15］。

采用不同的粒子間距參數(shù)得到的滑塊下落位移曲線幾乎一致，都能反映滑塊的運動特性。不過粒子間距較大時，滑塊會更早的受到邊界的排斥力作用，到達水底的時間略微延后。

圖5為右側邊界的涌浪沖擊壓力包絡分布圖。由圖5可見，1 m和2 m粒子間距得到的邊界最大沖擊壓力沿高度的分布基本一致，下部略有差別。涌浪作用下，壩面最大沖擊壓力部位在水表面附近，并且遠超過靜水壓力，甚至比庫底靜水壓力還要大。5 m的粒子間距過大，不能合理反映沖擊壓力的分布規(guī)律。說明在粒子稀疏的情況下，SPH粒子的偶然性和不均勻性大大提高，不足以正確反映水體粒子的運動特性。

圖4 滑塊下落過程時間曲線

滑坡入水點不同距離的涌浪高度隨時間的變化如圖6所示。從圖6可以看出，距離滑坡入水點較近部位的初始涌浪高度較大，距離稍遠處由于受到氣泡破碎、浪體回涌的影響而小于最大浪高。由于沿程能量損失，在傳播過程中涌浪高度緩慢衰減，平穩(wěn)浪高要小于初始最大浪高。說明穩(wěn)定涌浪呈現(xiàn)近似孤立波形態(tài)傳播。

圖5 壩面涌浪沖擊壓力包絡圖

圖6 不同距離處自由面曲線

2.3 自由面波形對于經典孤立波，假定波幅遠小于水深，則水面線方程為：

式中：y為水面高程；a為波幅；L為長度標量，L=d 4d 3a，d為水深。孤立波的重力勢能為：

式中：b為滑塊寬度；l為滑塊長度；ρ為水體密度。

假定滑塊重力勢能均轉化為孤立波的能量，有：

其中：m為滑塊質量。

代入Boussinesq方程［16］，可以得到孤立波波形的解析解為：

圖7 自由面波形

圖7為應用前文所述的SPH方法計算的t=10s時自由面波形數(shù)值模擬結果與式（10）得到的解析解波形對比。從圖7可見，數(shù)值模擬結果與解析解吻合良好，說明本文數(shù)值模型分析結果是可靠的。

3 影響因素分析

3.1 滑坡體初始入水速度對涌浪及沖擊壓力的影響我國西部水電工程位于高山峽谷地區(qū)，山高坡陡，來自于不同高度的滑坡體的入水速度是不同的，本文以庫水水深為70 m為例，對滑坡體入水速度分別為0、10及20 m/s情況下涌浪波高及壩面沖擊壓力進行比較研究。滑坡入水的形式有很多，不同的入水方式會對首浪高度以及涌浪傳播產生影響。由于本文的重點在于研究涌浪對于壩體表面的沖擊作用，從水庫涌浪產生本身來講，都可以等效為邊界沖擊條件和初始條件，因此計算中采用簡化的垂直入水得到不同的入水速度和入水體積的影響。

圖8 不同入水速度影響比較

圖8分別為滑坡體不同入水速度的位移-時間曲線、沿程最大浪高幅值曲線及壩面涌浪沖擊壓力包絡分布。由圖8可見，隨著滑塊入水速度的增大，滑塊到達底部的時間大大縮短，沿程涌浪幅值增加。涌浪對壩面的沖擊壓力幅值在壩面上部和下部隨著入水速度的增加而增加，但是在壩面中部沖擊壓力隨入水速度增加呈不規(guī)則變化。

3.2 滑坡體寬度對涌浪及沖擊壓力的影響滑坡體入水形態(tài)與寬度有關，本文以40、30和20 m三種尺寸為例研究滑坡體寬度對涌浪及沖擊壓力的影響。滑塊尺寸（長×寬×高）分別為10 m×40 m×25 m、10 m×30 m×33.33 m和10 m×20 m×50 m，水深25 m。

圖9分別為不同滑坡體寬度下的位移-時間曲線、沿程最大浪高幅值曲線及壩面涌浪沖擊壓力包絡分布。從圖9可見，體積不變條件下，隨著滑坡體寬度增加，入水阻力增大，入水過程明顯變得緩慢，首浪最大波高降低，穩(wěn)定波高沿著傳播途徑的衰減速率減小。當傳播距離超過100 m后，涌浪穩(wěn)定浪高隨著滑坡體寬度的增加而增大。不同寬度下形成的涌浪對壩面的沖擊壓力分布基本上一致，幅值隨著滑坡體寬度的增加而增加。

圖9 不同滑塊寬度影響比較

3.3 庫水水深對涌浪及沖擊壓力的影響本文分別對25、50和70 m三種不同庫水深度下的涌浪及其作用進行分析，結果如圖10所示。

從圖10可見，水深不影響滑坡體在水中的運動速度。由于滑坡體體積和入水速度相同的情況下，隨著水深增加，滑坡體傳遞給水體的能量向滑坡體下部深處增加，水平向擠壓產生的首浪最大浪高減小，說明在淺水情況下，滑坡更易形成較高的首浪波高，對于距離滑坡較近的大壩沖擊是不利的。但是隨著涌浪向擋水壩面的傳播，不同水深情況下的涌浪穩(wěn)定波高逐漸趨向一致。

從壩面涌浪沖擊壓力的分布也可以看出，不同水深情況下，最大沖擊壓力在壩面中上部的分布以及最大幅值基本上一致。但是隨著水深的增加，沖擊壓力在中下部的變化呈波動性多峰值分布，說明由于滑坡沖擊在深水情況下能量傳遞的復雜性，對于壩面的沖擊壓力有不利影響。

圖10 不同庫水深度影響比較

4 涌浪對擋水重力壩的影響

4.1 不同水位下壩前涌浪形成及傳播過程經過加固后的庫區(qū)邊坡通常具有較好的抗滑穩(wěn)定能力。但是，經歷強震作用后的邊坡抗滑穩(wěn)定能力會降低。如果強震中受損的邊坡在強余震或者強暴雨發(fā)生失穩(wěn)，其次生災害主要是對強震中同樣受損的擋水壩體的涌浪沖擊作用。

在滑坡涌浪浪高較大或者壩前水位較高時，涌浪有漫過壩頂?shù)那闆r，本文基于SPH-有限元耦合分析模型，分別研究低水位和高水位兩種情況下的滑坡涌浪沖擊作用，模型如圖11所示，數(shù)值模型中采用的主要參數(shù)如表2所示。

圖11 模型剖面

表2 主要材料參數(shù)

計算壩體高度103 m，庫區(qū)水體模擬區(qū)域長300 m，水深分別為70和100 m?；麦w體積為10 m×20 m×50 m，入水速度為0 m/s。

圖12為庫區(qū)涌浪傳播過程。從圖12可以看出，高水位和低水位情況下，涌浪產生和傳播過程與前面的分析基本上類似，隨著滑坡體下落，右側庫水被推擠迅速向右運動，隨后卷入波產生、越過滑塊、向左回涌、滑塊右側下部形成氣體空腔隨滑塊向底部運動等現(xiàn)象漸次發(fā)生。然后氣體空腔減小、消失，回涌水體沖擊到左側邊界開始向上爬升，然后回落再次形成新的波浪向右傳播，在右側壩面上爬升的涌浪開始回落。

在70 m低水位下最大浪高為12.74 m，在壩前沒有越過壩頂；100 m高水位情況下，最大浪高約11.24 m，10 s時壩前涌浪已開始漫上壩頂，隨后涌浪漫過壩頂，涌入下游，沖擊大壩下游面。

圖12 庫區(qū)涌浪傳播過程

4.2 不同水位下滑坡涌浪沖擊下的壩體位移圖13為壩頂順河向位移時程。由圖13可知，兩種水位下壩體受到的流體沖擊過程類似。首先受到的是快速傳播的壓縮波作用（波速1440 m/s），這是滑坡入水對庫水的瞬間沖擊產生的，對壩面形成第一次作用，壩體在壓縮波作用下發(fā)生較大動態(tài)位移，隨后振動幅度逐漸衰減。在8.2s左右，滑坡產生的涌浪重力波（波速28～33 m/s）傳播到壩前，在涌浪沖擊下壩體振動位移再次增加，隨后振動幅度逐漸衰減。高水位下，10s左右涌浪爬坡越過壩頂沖向下游，涌浪作用下的壩體位移振動平衡位置逐漸偏移。

圖13 壩頂點位移-時間曲線

4.3 不同水位下滑坡涌浪沖擊下的壩體拉應力圖14為滑坡庫水沖擊下壩踵和下游折坡點的第一主應力變化時程曲線。從圖14可見，低水位下的壩體應力在涌浪沖擊作用下比庫水壓縮波動壓力作用下的大，而高水位下相反。這一方面是因為滑坡入水過程中，滑坡體沖擊庫水傳遞的壓縮波動能量與擠壓庫水傳遞的涌浪重力波動能量在不同的水深條件下不同，另一方面是高水位下，由于涌浪越過壩頂就降低了涌浪沖擊壓力的影響。

圖14 壩踵與下游折坡點第一主應力時程曲線

壩踵和下游折坡點應力均在壓縮波和涌浪沖擊到達時最大，隨后衰減，相比而言，高水位下由于涌浪越壩使得涌浪回涌減少，下游折坡點的應力衰減更快；下游折坡點在庫水沖擊作用下的應力增加幅度遠遠大于壩踵，說明滑坡引起的庫水沖擊壓力對壩體中上部影響更大，壩頭部位更為危險。

4.4 庫水沖擊壓力分析當前研究滑坡涌浪對擋水結構的沖擊壓力作用開展的相對較少，有限的研究也主要集中在滑坡入水引起的涌浪波特征及其作用上。實際上，滑坡對庫水的影響包括兩部分，一部分是滑坡入水對庫水邊界沖擊引起的流體體積壓縮波，波速為1440 m/s；一部分是滑坡體擠壓庫水形成的流體表面重力波，在本文中傳播波速為28～33 m/s之間；由于波速的差異，壓縮體積波首先傳播到壩面形成沖擊，隨后是涌浪沖擊作用。兩種波的幅值和壩面分布形式是不同的。圖15為大壩上游面壓力。從圖15（a）可以看出，由于壓縮波為流體的機械波動，其作用范圍在庫水靜水位以下，并且在水表面為零。而涌浪為表面重力波，其作用范圍與涌浪在壩面的爬升高度有關；從圖15（b）可以看出，其作用范圍高出靜水位線超過12 m，由于涌浪沖擊作用性質，即使在水表面也具有很大的作用幅值。

圖15 大壩上游面壓力

在縮尺模型試驗中由于庫水相似材料選擇的限制，很難同時滿足兩種波動的相似率設計要求，由于壓縮波波速很快，縮尺條件下通常只能測得涌浪的生成、傳播和作用，很難測到壓縮波的傳播特征和影響。而在實際工程中庫區(qū)尺度相對較大，大壩同時受到兩種庫水沖擊的先后作用，在實際中是不能忽視庫水壓縮波沖擊作用的影響的。實際上，由于兩者的強度、壩面分布形式均不相同，壓縮波沖擊作用對大壩也有很大的影響。特別是兩種沖擊波是具有不同周期特性的動力作用。按照壩面涌浪壓力分布采用靜力計算的方法分析壩體動力響應可能是偏于危險的。

5 結論

針對水利樞紐庫區(qū)邊坡失穩(wěn)誘發(fā)的庫水沖擊問題，分析了滑坡涌浪的產生、傳播以及對壩面沖擊壓力的作用特征，研究了滑坡入水速度、庫水水深和滑塊寬度等因素對以上特征的影響，并分別針發(fā)生漫頂和不漫頂?shù)膬煞N水位下的壩體動力響應特征進行了分析。本文得到以下結論：（1）基于SPH方法的粒子間距選擇必須滿足不同波動的波長要求，過大的粒子間距無法捕捉流體的正確波動形態(tài)。通過與解析解的對比，得到的壩面沖擊壓力作用以及對壩體位移和應力的分析是合理的，說明本文的分析模型是可靠的。（2）滑坡體參數(shù)及庫水深度等因素對于涌浪生成及傳播過程的定性特征影響不大，但是對于壩面的沖擊壓力強度和分布特征影響很大。（3）庫區(qū)滑坡入水會引發(fā)庫水機械波和重力波的沖擊作用。兩種波的傳播速度差異接近2個數(shù)量級，壩體先后受到庫水壓縮波和涌浪的沖擊作用。壓縮波沖擊作用在靜水位以下，靜水位部位沖擊壓力為零；涌浪沖擊作用范圍在涌浪自由表面以下，在靠近水位線附近強度最大。（4）壩體上部在沖擊波作用下的應力增大幅值遠遠大于中下部，不過由于靜力作用下壩踵部位為拉應力薄弱部位，因此，庫水沖擊作用下壩踵和壩頭部位都為危險薄弱部位。（5）涌浪沖擊作用起主要作用，但由于兩種沖擊作用的動力特性以及在壩面的分布形式不同，對于壩體不同部位動力響應的影響規(guī)律是不同的，在實際工程中不能忽視壓縮體積波沖擊作用對壩體動力破壞的影響。

庫水沖擊作用對壩體的影響與滑坡體入水的形式、形態(tài)、體積、速度、角度以及庫水深度、滑坡體入水點與壩體的距離、壩體與庫水的相對高度、壩體動力特性等很多因素有關，本文只研究了主要影響因素中的幾個，其他相關因素的影響還有待于進一步深入研究。

［1］谷健.某水電站庫區(qū)變形體滑坡涌浪初步分析［J］.水電站設計，2013，29（4）：56-59.

［2］龐昌俊.二維斜滑坡涌浪的試驗研究［J］.水利學報，1985（11）：54-59.

［3］殷坤龍，杜娟，汪洋.清江水布埡庫區(qū)大堰塘滑坡涌浪分析［J］.巖土力學，2008，29（12）：3266-3270.

［4］吳佳壕.山區(qū)河道型水庫庫岸坡涌浪災害研究與預測［D］.成都：成都理工大學，2011.

［5］徐文杰.滑坡涌浪影響因素研究［J］.工程地質學報，2012，20（4）：491-507.

［6］繆吉倫，陳景秋，張永祥 .基于SPH方法的立面二維涌浪數(shù)值模擬［J］.中南大學學報（自然科學版），2012，43（8）：3244-3249.

［7］黃錦林.庫岸滑坡涌浪對壩體影響研究［D］.天津：天津大學，2011.

［8］楊艷.山區(qū)河道型水庫滑坡涌浪對架空直立式碼頭結構作用影響研究［D］.重慶：重慶交通大學，2014.

［9］QIU L.Two-dimensional SPH simulations of landslide-generated water waves［J］.Journal of Hydraulic Engineer?ing，2008，134（5）：668-671.

［10］GóMEZ-GESTEIRA M，DALRYMPLE R.Using a three-dimensional smoothed particle hydrodynamics method for wave impact on a tall structure［J］.Journal of Waterway，Port，Coastal and Ocean Engineering，2004，130（2）：63-69.

［11］AMINI Y，EMDAD H，F(xiàn)ARID M.A new model to solve fluid-hypo-elastic solid interaction using the smoothed particle hydrodynamics（SPH）method［J］.European Journal of Mechanics-B/Fluids，2011，30（2）：184-194.

［12］MONAGHAN J J，GINGOLD R A.Shock simulation by the particle method SPH［J］.Journal of Computational Physics，1983，52（2）：374-389.

［13］COLAGROSSI A，LANDRINI M.Numerical simulation of interfacial flows by smoothed particle hydrodynamics［J］.Journal of Computational Physics，2003，191（2）：448-475.

［14］張姝慧.求解淺水波方程的光滑粒子流體動力學法［D］.合肥：安徽大學，2007.

［15］YIM S，YUK D，PANIZZO A，et al.Numerical simulations of wave generation by a vertical plunger using RANS and SPH models［J］.Journal of Waterway，Port，Coastal and Ocean Engineering，2008，134（3）：143-159.

［16］李孟國，王正林，蔣德才.關于波浪Boussinesq方程的研究［J］.青島海洋大學學報（自然科學版），2002，32（3）：345-354.