司鵬飛，田利勇

(1.上海市水利工程設(shè)計研究院有限公司，上海 200061; 2.上海市水務(wù)局防汛減災(zāi)工程技術(shù)研究中心，上海 200061)

庫區(qū)滑坡是水庫蓄水運行期間普遍存在的問題?；麦w入水產(chǎn)生的涌浪破壞力巨大，由于波浪爬高及漫頂作用，會對庫區(qū)擋水建筑物、水上航行、周邊居民的生命財產(chǎn)造成嚴(yán)重威脅[1]。例如，1963年意大利的瓦依昂水庫庫岸發(fā)生滑坡，近3億m3的滑坡體以30 m/s的速度滑入水庫，掀起高達(dá)175 m的涌浪，摧毀了壩區(qū)的所有建筑物，66名水庫技術(shù)人員全部遇難，涌浪以超出壩頂約100 m的高度沖入下游河谷，沖毀5個村鎮(zhèn)，并奪去近3 000人的生命[2]。1961年我國發(fā)生了首例水庫蓄水初期誘發(fā)的大型滑坡，湖南柘溪水庫的大壩上游右岸1.5 km處的塘巖光發(fā)生了165萬m3的大型滑坡，該滑坡以近20 m/s的速度入水后，產(chǎn)生的涌浪傳播到對岸時浪高仍有21 m，并以3.6 m高的涌浪越過正在施工的壩頂，造成40余人死亡[3]。

為了研究滑坡涌浪的特性并對此類災(zāi)害進行預(yù)警及風(fēng)險評估，學(xué)者們建立了多種不同的數(shù)學(xué)模型[4-5]。針對崩塌、巨石滾落以及歷時較短、發(fā)展迅速的滑坡涌浪問題，很多學(xué)者將滑坡體簡化為剛性體處理，即不考慮滑坡體在短時間運動過程中自身的變形特性。Wu等[6]利用大渦模擬模型以及VOF水面追蹤方法，將物理模型試驗中測得的滑塊實際運動軌跡輸入到數(shù)值模型中，研究了三維塊體下滑造成的涌浪爬坡和回落現(xiàn)象；Heinrich[7]利用二維水動力NASA-VOF2D模型，計算了三角形滑塊入水產(chǎn)生的波面，在整個模擬中利用已知的試驗數(shù)據(jù)對滑塊的運動參數(shù)進行了設(shè)定；Watts等[8]設(shè)定滑坡體位移-時間曲線為對數(shù)雙曲余弦函數(shù)形式，利用Boussinesq方程模擬了滑坡體入水產(chǎn)生的涌浪傳播過程。但是以上數(shù)學(xué)模型都是將滑坡體的運動作為已知參數(shù)直接輸入到建立的水波模型中進行計算，沒有考慮滑坡體入水后與水體的相互作用。也有許多學(xué)者直接假定滑坡體入水后作勻速運動或勻加速運動[9-11]，或假設(shè)滑坡體的位移-時間曲線為正弦或余弦形式[12-13]，或者假設(shè)水體對滑坡體的作用力與滑坡體速度滿足一定的關(guān)系[14]。但是在實際情況中滑坡體的運動特性難以事先確定，特別是滑坡體入水后的速度變化規(guī)律更是難以預(yù)測，這樣的假設(shè)往往會使模擬結(jié)果與滑坡涌浪的實際特性存在一定的偏差。

本文通過建立考慮滑坡涌浪過程中滑坡體與水體相互作用的流固耦合模型，利用FLOW-3D軟件進行數(shù)值求解，研究剛性滑坡體入水生成涌浪的規(guī)模及其傳播過程，并針對滑坡體在不同運動狀態(tài)下所產(chǎn)生的涌浪特性進行對比，分析流場、涌浪爬高及其傳播過程的異同。

1 數(shù)值模型

流體的控制方程采用考慮流固耦合作用的三維N-S方程：

(1)

(2)

式中：xi、xj(i，j=1,2,3)為坐標(biāo)分量；ui為xi方向的流體速度分量；p為壓力；νt為紊動黏性系數(shù)；gi為重力加速度在xi方向的分量；VF為計算域所占網(wǎng)格單元的體積分?jǐn)?shù)，Ai為計算單元在xi方向的截面面積比例；ρ為流體的密度；t為時間。FLOW-3D采用FAVOR(fractional area/volume obstacle representation)方法定義計算單元并求得相應(yīng)單元的VF與Ai的值[15]，隨著滑坡體運動過程中位置和坐標(biāo)的更新，流體計算單元的面積和體積分?jǐn)?shù)也會被重新計算。

控制方程的空間離散采用三維矩形交錯網(wǎng)格，數(shù)值離散采用有限差分法，同時采用VOF方法[16-17]對自由水面進行追蹤處理?？坍嬜杂伤娴腣OF函數(shù)滿足方程：

(3)

式中：α為流體占計算單元空間的體積分?jǐn)?shù)，取值范圍為0～1。α=0說明計算單元內(nèi)全為空氣，α=1說明計算單元內(nèi)全為水體，當(dāng)0<α<1時，水氣混合處于自由水面。式(3)中的右端項反映了滑坡體運動對流體空間的擠占作用。

采用k-ε雙方程紊流模型模擬紊流效應(yīng)，流體的紊動黏性系數(shù)νt=Cμk2/ε，其中經(jīng)驗系數(shù)Cμ取為0.09，紊動能k和紊動能耗散率ε的守恒方程分別表示為

(4)

(5)

式中：ν為水體的運動黏性系數(shù)；經(jīng)驗常數(shù)Cε1=1.44，Cε2=1.92，σk=1.0，σε=1.0。

利用FLOW-3D中的通用剛體運動模型(general moving object, GMO)模擬滑坡體運動，可以選擇指定滑坡體的運動軌跡或是通過流固耦合計算，當(dāng)采用流固耦合模型時，滑坡體的控制方程可表示為

(6)

式中：Fi為滑坡體沿xi方向受到的合力，包括重力、流體的阻力和壁面摩擦阻力；m為滑坡體的質(zhì)量；vi為滑坡體沿xi方向的運動速度。

2 模型驗證

2.1 模型設(shè)置

采用建立的模型模擬滑坡體沿斜坡滑入水中產(chǎn)生的涌浪特性及其傳播情況，為驗證數(shù)值模型的有效性，將數(shù)值計算結(jié)果與Wu等[6]在港池中進行的滑坡涌浪試驗觀測數(shù)據(jù)進行對比。數(shù)值模型設(shè)置如圖1所示，試驗中的涌浪爬高測點和傳播過程中的涌浪波高測點布置如圖2所示，其中1號、2號用來測量涌浪沿斜坡的爬高值，3號、4號、5號、6號用來測量涌浪傳播過程中不同位置處波高的變化情況。模型區(qū)域設(shè)置為6.6 m(長)×3.7 m(寬)×3.3 m(高)，3個方向的網(wǎng)格數(shù)量分別為60個、100個和60個，共計36萬個網(wǎng)格；斜坡位于模型區(qū)域一端，長高比為2∶1?；瑝K剖面形狀為三角形，相應(yīng)尺寸為0.91 m(長)×0.61 m(寬)×0.455 m(高)，滑塊密度設(shè)為2 640 kg/m3。初始水深設(shè)為3 m，初始時刻滑塊下表面剛好與靜止水面齊平，斜坡壁面摩擦系數(shù)設(shè)為0.2，模擬時間為4 s。

圖1 模型設(shè)置示意圖 (單位：m)

圖2 涌浪爬高、波高測點布置(單位：m)

2.2 邊界條件

2.3 結(jié)果分析

利用所建立的流固耦合數(shù)學(xué)模型進行數(shù)值求解，分別將波高測點和涌浪爬高測點的數(shù)值計算結(jié)果與物理試驗結(jié)果進行對比，見圖3和圖4。

從圖3中可以看出，數(shù)值計算結(jié)果與試驗數(shù)據(jù)整體吻合較好，流固耦合模型能夠較準(zhǔn)確地模擬出滑坡入水后產(chǎn)生的首浪高度和涌浪傳播速度；4號和6號測點處第二個波峰的計算高度略大于試驗值，原因可能是由于模型中側(cè)壁設(shè)置為對稱反射邊界，而試驗中的港池側(cè)壁不一定能夠完全實現(xiàn)水波的反射，因此導(dǎo)致了波高計算值偏大。不過整體數(shù)值模擬結(jié)果還是能夠較準(zhǔn)確地反映涌浪的產(chǎn)生和傳播過程。

從圖4可以看出，模型計算所得涌浪爬高與試驗觀測值吻合較好，表明所建模型能夠準(zhǔn)確模擬滑坡入水后涌浪的傳播與爬坡過程。

圖3 不同測點涌浪波高計算值與試驗值對比

圖4 不同測點涌浪爬高計算值與試驗值對比

3 不同滑坡體運動狀態(tài)下的涌浪特性

3.1 勻加速運動

假設(shè)忽略水體對滑坡體的作用力和滑坡體所受摩擦力，即滑坡體在自身重力的作用下勻加速下滑，則滑坡體沿xi方向的加速度ai和速度vi可通過下式計算:

mgisinβ=mai

(7)

vi=ait

(8)

式中：β為斜坡傾角。

若不考慮滑坡體與水體的相互作用，將式(7)(8)得到的滑坡體運動軌跡直接作為輸入?yún)?shù)，利用所建立的數(shù)學(xué)模型進行數(shù)值求解，3號、4號波高測點的數(shù)值計算結(jié)果與2.3小節(jié)中考慮流固耦合模型時的計算結(jié)果對比如圖5所示。

圖5 滑坡體勻加速運動與流固耦合作用下涌浪波高對比

從圖5可以看出，若指定滑坡體在重力作用下作勻加速運動，則滑坡體入水產(chǎn)生的涌浪規(guī)模及其傳播速度與利用流固耦合模型計算得到的結(jié)果存在較大的差異。采用指定滑坡體作勻加速運動的計算模型所得的首浪高度要比實際情況大30%左右，而且涌浪傳播速度較流固耦合模型算出的結(jié)果也有明顯的增加，但兩種方法計算所得的整體波形基本一致，說明涌浪的產(chǎn)生和發(fā)展主要是由于滑坡體排開其運動軌跡前方一定體積的水而形成波峰，水體與滑坡體的相互作用只是改變了波峰的高度以及涌浪的傳播速度，對整體波形影響較小。通過上述研究可以得知，在滑坡涌浪模擬過程中若直接指定滑坡體的運動軌跡進行計算所得到的結(jié)果是不合理的。

3.2 流固耦合+勻(加)速運動

為了研究流固耦合作用在滑坡體不同運動階段對涌浪特性所起的作用，將滑坡涌浪過程分為兩個階段：滑坡體從開始下滑到完全浸入水中(入水階段)和滑坡體完全浸入水體后(水下階段)。入水階段采用流固耦合模型，水下階段忽略水體對滑坡體的阻力作用，采用給定滑坡體運動軌跡的方式進行計算。

首先利用流固耦合模型模擬計算滑坡體入水階段(0≤t≤1.5 s)的涌浪特性，在模擬結(jié)束后，利用所得計算結(jié)果作為忽略水體阻力的水下階段(1.5 s

圖6 滑坡體分段運動狀態(tài)下的涌浪波高對比

由圖6可知，將滑坡體完全浸入水體時作為分界點分段計算所得的結(jié)果與全過程均用流固耦合模型計算所得的結(jié)果高度吻合，并且兩種不同運動軌跡(勻速、勻加速)的計算結(jié)果也相差不大。表明當(dāng)滑坡體完全浸入水體之后，其運動狀態(tài)對涌浪的特性及傳播過程幾乎沒有影響，說明滑坡體與水體的相互作用對涌浪特性的影響主要發(fā)生在滑坡體接觸水面到完全浸沒的過程中，因此在此階段采用流固耦合模型來考慮滑坡體與水體相互作用對涌浪產(chǎn)生的影響是非常必要的。

4 結(jié) 論

a. 建立的滑坡涌浪三維流固耦合模型具有較高的精度，能夠準(zhǔn)確模擬出涌浪的波高、傳播速度、沿岸爬高等特征參數(shù)。

b. 若不考慮滑坡體與水體的相互作用，直接給定滑坡體運動軌跡進行數(shù)值模擬所得的涌浪波高及傳播速度要遠(yuǎn)遠(yuǎn)大于實際情況。

c. 滑坡體與水體相互作用對涌浪特性的影響主要體現(xiàn)在滑坡體入水階段，在此階段運用流固耦合模型進行數(shù)值求解，而在水下階段，根據(jù)前一階段流固耦合模型計算所得結(jié)果指定滑坡體入水后的運動軌跡所計算得到的結(jié)果與實際觀測結(jié)果非常吻合。表明本文提出的滑坡體入水階段采用流固耦合模型，水下階段通過直接給定滑坡體運動軌跡的方法，既能較好地模擬出涌浪的產(chǎn)生與傳播過程，又可以簡化計算，對實際大型滑坡涌浪問題的模擬研究具有借鑒意義。