基于流固耦合方法的滑坡涌浪數(shù)值模擬

司鵬飛，田利勇

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

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

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

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

1 數(shù)值模型

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

(1)

(2)

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

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

(3)

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

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

(4)

(5)

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

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

(6)

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

2 模型驗(yàn)證

2.1 模型設(shè)置

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

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

圖2 涌浪爬高、波高測(cè)點(diǎn)布置(單位：m)

2.2 邊界條件

2.3 結(jié)果分析

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

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

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

圖3 不同測(cè)點(diǎn)涌浪波高計(jì)算值與試驗(yàn)值對(duì)比

圖4 不同測(cè)點(diǎn)涌浪爬高計(jì)算值與試驗(yàn)值對(duì)比

3 不同滑坡體運(yùn)動(dòng)狀態(tài)下的涌浪特性

3.1 勻加速運(yùn)動(dòng)

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

mgisinβ=mai

(7)

vi=ait

(8)

式中：β為斜坡傾角。

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

圖5 滑坡體勻加速運(yùn)動(dòng)與流固耦合作用下涌浪波高對(duì)比

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

3.2 流固耦合+勻(加)速運(yùn)動(dòng)

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

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

圖6 滑坡體分段運(yùn)動(dòng)狀態(tài)下的涌浪波高對(duì)比

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

4 結(jié) 論

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

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

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