本發(fā)明涉及一種滑坡涌浪計算方法，屬于水利工程災(zāi)害預(yù)測領(lǐng)域。

背景技術(shù)：

水利水電建設(shè)的長足發(fā)展為社會發(fā)展做出了重要貢獻(xiàn)，但也帶來一系列地質(zhì)環(huán)境問題。相比于一般的邊坡失穩(wěn)問題，庫岸邊坡失穩(wěn)會對下游的人民生命財產(chǎn)安全構(gòu)成更加巨大的威脅?；麦w滑入水庫后產(chǎn)生的涌浪有可能翻越壩頂，淹沒下游地區(qū)，甚至由于涌浪的沖擊作用，大壩存在潰決的風(fēng)險。在國內(nèi)水利水電工程建設(shè)中，不乏水庫滑坡失穩(wěn)破壞并產(chǎn)生涌浪災(zāi)害的先例。如1961年3月柘溪水庫塘巖光滑坡，165萬m³的滑坡體沖入庫區(qū)后產(chǎn)生了高達(dá)21m的涌浪，造成重大損失。隨著一系列高壩大庫在我國西南地區(qū)建立，水庫庫岸邊坡穩(wěn)定問題越來越吸引人們的注意。在此背景下，數(shù)值模擬作為災(zāi)害預(yù)測的手段之一，就顯的非常必要和重要了。

目前的數(shù)值模擬基本是首先用能量守恒原理估算滑坡體入水前的滑速，而后依據(jù)經(jīng)驗(yàn)公式計算初始涌浪，最后用計算流體力學(xué)的方法求解涌浪的傳播過程，實(shí)際上只是求解了已知初始條件和邊界條件的計算流體力學(xué)問題。滑坡體入水產(chǎn)生涌浪的過程是一個復(fù)雜的水動力學(xué)問題，過程中涉及到三相介質(zhì)的相互耦合、浸入、波浪破碎等現(xiàn)象，介質(zhì)間的界面難于捕捉描述。現(xiàn)有的數(shù)值方法，包括基于淺水方程的解法，不能反映三相介質(zhì)的強(qiáng)烈相互作用。主要的界面描述方法有VOF法、Level Set法等，涌浪自由面的描述在計算流體力學(xué)中本身就是很大的難題，數(shù)值計算的準(zhǔn)確性可以得到進(jìn)一步改善。相對于涌浪產(chǎn)生來說，涌浪的傳播過程研究成果已較為豐富。但此類計算的前提是已知首浪的各種要素，如浪高、速度等。從某種意義上說，現(xiàn)有方法是用一個粗估的初始涌浪數(shù)據(jù)作為初始計算條件，即使涌浪傳播過程的模擬再精細(xì)，也很難得到合理可信的計算結(jié)果。

盡管國內(nèi)外學(xué)者對此已做了大量的研究工作，取得了豐富的研究成果，但基本上是著眼于某個方面的研究，還沒有針對從庫岸邊坡失穩(wěn)運(yùn)動到涌浪產(chǎn)生與傳播全過程系統(tǒng)的數(shù)值模擬。而實(shí)際上此類災(zāi)害的發(fā)生基本是在幾十秒到幾分鐘極短的時間內(nèi)，庫岸邊坡從失穩(wěn)、滑動到涌浪產(chǎn)生再到涌浪傳播應(yīng)該是一個完整連續(xù)的過程，每個階段的準(zhǔn)確描述都將對最后災(zāi)害預(yù)報的可信度產(chǎn)生較大影響。

技術(shù)實(shí)現(xiàn)要素：

本發(fā)明所要解決的技術(shù)問題是針對背景技術(shù)的缺點(diǎn)，將庫岸邊坡滑坡體入水到涌浪傳播各階段作為完整統(tǒng)一的過程進(jìn)行分析研究。采用原始的不可壓縮粘性流體Navier-Stokes方程描述水體、空氣、滑坡體的相互作用。采用改進(jìn)的守恒式Level Set方法，準(zhǔn)確捕捉流體自由表面，真實(shí)模擬涌浪產(chǎn)生過程，為涌浪傳播及災(zāi)害預(yù)測提供合理的前提條件。在此基礎(chǔ)上，采用沿水深積分的流體控制方程模擬涌浪傳播過程，對涌浪運(yùn)動軌跡和泛濫區(qū)域進(jìn)行數(shù)值預(yù)報，評價涌浪災(zāi)害。本發(fā)明是在已經(jīng)利用動力時程強(qiáng)度折減法進(jìn)行土質(zhì)庫岸邊坡動力穩(wěn)定分析，確定了庫岸可能的失穩(wěn)滑動范圍的前提下進(jìn)行的。

本發(fā)明為實(shí)現(xiàn)上述發(fā)明目的采用如下技術(shù)方案：

一種滑坡涌浪計算方法，包括如下步驟：

步驟1：根據(jù)已有地形、水位和滑坡體形狀位置資料，建立三維有限元模型；

步驟2：結(jié)合Navier-Stokes方程和改進(jìn)的守恒式Level Set方法，計算流體各運(yùn)動參數(shù)，捕捉自由面；

步驟3：采用沿水深積分的流體控制方程模擬涌浪傳播過程，對涌浪運(yùn)動軌跡和泛濫區(qū)域進(jìn)行數(shù)值預(yù)報。

所述步驟2中對土質(zhì)或松散介質(zhì)坡體，基于“流體化”的特點(diǎn)，以流體力學(xué)基本方程N(yùn)avier-Stokes方程配以非牛頓流體或擬流體本構(gòu)關(guān)系描述其物理特性和運(yùn)動過程，由于將滑坡體視作似流體，滑坡體、水體和空氣三相介質(zhì)可用同一本構(gòu)模型表示：

式中，τ_y0表示流體屈服應(yīng)力，φ為內(nèi)摩擦角，I_2D是應(yīng)變張量D的第二不變量，μ為黏性系數(shù)，η為表示流體特性的參量。

在利用Navier-Stokes方程計算出的流場信息的基礎(chǔ)上，采用改進(jìn)的守恒式雙重Level Set方法捕捉滑坡體、水體和空氣三者的交界面，反映三者間強(qiáng)烈的流固耦合作用，真實(shí)模擬滑坡體入水產(chǎn)生涌浪的過程，改進(jìn)的守恒式雙重Level Set方法包含兩個部分：指示函數(shù)的輸運(yùn)，和指示函數(shù)的重新初始化。

在三相流中，各點(diǎn)的物理參數(shù)用指示函數(shù)H表示，各點(diǎn)的物理參數(shù)包括密度或黏性；引入指示函數(shù)H₁、H₂，滑坡體、水體和空氣三相介質(zhì)各自區(qū)域分別記作Ω₁、Ω₂、Ω₃，滑坡體與水體、水體與空氣、滑坡體與空氣界面依次用Γ₁₂、Γ₂₃、Γ₁₃表示，P表示物理參數(shù)，則各區(qū)域參數(shù)用統(tǒng)一的方程表示：

P＝P₁+(P₂-P₁)H₁+(P₃-P₂)H₂

在Ω₁內(nèi)，H₁＝H₂＝0；在Ω₂內(nèi)，H₁＝1，H₂＝0；在Ω₃內(nèi)，H₁＝H₂＝1。

交界面的位置確定方法，引入標(biāo)記距離函數(shù)和指示函數(shù)H，或H＝0.5表示界面，重新初始化方程中包含界面的法線方向，由于指示函數(shù)H本身的局限性，不適宜用于求解遠(yuǎn)離H＝0.5處的界面法線方向，選取一個臨界值δ，當(dāng)時，利用計算界面法線方向，否則用H計算法線方向，具體過程表述如下：

(1)求解兩個指示函數(shù)的輸運(yùn)方程，得到和H；

(2)對遠(yuǎn)離界面處的標(biāo)記距離函數(shù)進(jìn)行初始化；

(3)計算界面彌散區(qū)域各點(diǎn)的法線方向；

(4)對指示函數(shù)H進(jìn)行初始化；

(5)根據(jù)H計算靠近界面處的距離標(biāo)記函數(shù)

(6)進(jìn)行物理特性插值，確定界面位置，進(jìn)行下一步。

交界面的位置確定過程中選取的臨界值δ一般可以取0.1Δ～0.2Δ，其中Δ表示網(wǎng)格尺寸，當(dāng)時表示遠(yuǎn)離界面，表示靠近界面。

以步驟2中首浪產(chǎn)生數(shù)值模結(jié)果作為前提進(jìn)行連貫的數(shù)值分析，用沿水深積分的Navier-Stokes方程即淺水方程作為控制方程求解涌浪的傳播過程，通過步驟3的計算，可以得到各時刻各測點(diǎn)的水深數(shù)值，對涌浪的運(yùn)動軌跡和泛濫區(qū)域進(jìn)行數(shù)值預(yù)報，評價涌浪災(zāi)害。

典型滑坡體截面產(chǎn)生的初始浪高計算方法如下：在滑坡體進(jìn)入水體之后的任意時刻，取指示函數(shù)H值分布圖，取兩點(diǎn)，分別為點(diǎn)A和點(diǎn)B，設(shè)H_A＜0.5，點(diǎn)A位于空氣中，H_B＞0.5，點(diǎn)B位于水體中，A點(diǎn)高程h_A，B點(diǎn)高程h_B，根據(jù)線性插值遠(yuǎn)離，容易得出，H＝0.5的點(diǎn)的高程為初始浪高h(yuǎn)＝h_O-h_Still，其中h_Still表示靜水位。

本發(fā)明采用上述技術(shù)方案具有如下有益效果：

本發(fā)明針對庫岸邊坡失穩(wěn)到涌浪產(chǎn)生與傳播全過程進(jìn)行數(shù)值模擬采用，將庫岸邊坡從失穩(wěn)、滑動到涌浪產(chǎn)生再到涌浪傳播看作是一個完整連續(xù)的過程，采用Navier-Stokes方程和改進(jìn)的守恒式Level Set方法計算出首浪的各參數(shù)，在此基礎(chǔ)上，采用淺水方程計算涌浪的傳播，對涌浪運(yùn)動軌跡和泛濫區(qū)域進(jìn)行數(shù)值預(yù)報，進(jìn)行涌浪災(zāi)害評價。

附圖說明

圖1是本發(fā)明三維有限元網(wǎng)格圖，

圖2是本發(fā)明典型斷面密度分布圖，

圖3是本發(fā)明三相流指示函數(shù)H₁、H₂示意圖，

圖4是本發(fā)明首浪高度計算示意圖，

圖5是本發(fā)明河道截面處涌浪水位隨時間變化過程(靠近滑坡體入水處)，

圖6是本發(fā)明河道截面處涌浪水位隨時間變化過程(靠近大壩)。

具體實(shí)施方式

下面結(jié)合附圖和具體實(shí)施方式，進(jìn)一步闡明本發(fā)明。應(yīng)理解下述具體實(shí)施方式僅用于說明本發(fā)明而不用于限制本發(fā)明的范圍。

本發(fā)明采用原始的不可壓縮粘性流體Navier-Stokes方程描述水體、空氣、滑坡體的相互作用。采用改進(jìn)的守恒式Level Set方法，準(zhǔn)確捕捉流體自由表面，真實(shí)模擬涌浪產(chǎn)生過程。在此基礎(chǔ)上，采用沿水深積分的流體控制方程模擬涌浪傳播過程，對涌浪運(yùn)動軌跡和泛濫區(qū)域進(jìn)行數(shù)值預(yù)報，評價涌浪災(zāi)害。

一種滑坡涌浪計算方法，包括如下步驟：

步驟1：根據(jù)已有地形、水位和滑坡體形狀位置資料，建立三維有限元模型；

步驟2：結(jié)合Navier-Stokes方程和改進(jìn)的守恒式Level Set方法，計算流體各運(yùn)動參數(shù)，捕捉自由面；

步驟3：采用沿水深積分的流體控制方程模擬涌浪傳播過程，對涌浪運(yùn)動軌跡和泛濫區(qū)域進(jìn)行數(shù)值預(yù)報。

步驟1具體如下：

根據(jù)實(shí)際地形資料，利用前處理軟件，如GID、AutoCAD、Abaqus等建立水體、空氣、滑坡體三維有限元模型。

步驟2具體如下：

步驟(2-1)對土質(zhì)或松散介質(zhì)坡體，基于“流體化”的特點(diǎn)，以流體力學(xué)基本方程N(yùn)avier-Stokes方程配以非牛頓流體或擬流體本構(gòu)關(guān)系描述其物理特性和運(yùn)動過程：

進(jìn)一步地，由于將滑坡體視作似流體，滑坡體、水體和空氣三相介質(zhì)可用同一本構(gòu)模型表示：

式中，τ_y0表示流體屈服應(yīng)力，φ為內(nèi)摩擦角，I_2D是應(yīng)變張量D的第二不變量，μ為黏性系數(shù)，η為表示流體特性的指數(shù)。

步驟(2-2)在利用Navier-Stokes方程計算出的流場信息的基礎(chǔ)上，采用改進(jìn)的守恒式雙重Level Set方法捕捉滑坡體、庫水和空氣三者的交界面，反映三者間強(qiáng)烈的流固耦合作用，真實(shí)模擬滑坡體入水產(chǎn)生涌浪的過程。改進(jìn)的守恒式雙重Level Set方法包含兩個部分：指示函數(shù)的輸運(yùn)，指示函數(shù)的重新初始化。

為了確定交界面的位置，引入標(biāo)記距離函數(shù)和指示函數(shù)H?；騂＝0.5表示界面。重新初始化方程中包含界面的法線方向，由于指示函數(shù)H本身的局限性，不適宜用于求解遠(yuǎn)離H＝0.5處的界面法線方向。選取一個臨界值δ，當(dāng)時，利用計算界面法線方向，否則用H計算法線方向。

步驟(2-2)的具體過程表述如下：

(1)在步驟(2-1)求解得到的速度場的基礎(chǔ)上，求解兩個指示函數(shù)的輸運(yùn)方程，得到和H；

(2)對遠(yuǎn)離界面處的標(biāo)記距離函數(shù)進(jìn)行初始化；

(3)計算界面彌散區(qū)域各點(diǎn)的法線方向；

(4)對指示函數(shù)H進(jìn)行初始化；

(5)根據(jù)H計算靠近界面處的距離標(biāo)記函數(shù)

(6)進(jìn)行物理特性插值，確定界面位置，進(jìn)行下一步。

進(jìn)一步地，上述過程中選取的臨界值δ一般可以取0.1Δ～0.2Δ，其中Δ表示網(wǎng)格尺寸，當(dāng)時表示遠(yuǎn)離界面，表示靠近界面。

交替計算步驟(2-1)、(2-2)，可以得到三相流運(yùn)動參數(shù)、捕捉界面位置，觀察涌浪產(chǎn)生過程。選取初始涌浪作為涌浪傳播計算的初始數(shù)值，用以計算后續(xù)的涌浪傳播。

進(jìn)一步地，在三相流中，各點(diǎn)的物理參數(shù)(如密度，黏性等)可以用指示函數(shù)H表示。引入指示函數(shù)H₁、H₂三相介質(zhì)各自區(qū)域分別記作Ω₁、Ω₂、Ω₃，界面分別用Γ₁₂、Γ₂₃、Γ₁₃表示。P表示物理參數(shù)，則各區(qū)域參數(shù)可以用統(tǒng)一的方程表示：

P＝P₁+(P₂-P₁)H₁+(P₃-P₂)H₂

在Ω₁內(nèi)，H₁＝H₂＝0；在Ω₂內(nèi)，H₁＝1,H₂＝0；在Ω₃內(nèi)，H₁＝H₂＝1。

步驟3具體如下：

以步驟2中首浪產(chǎn)生數(shù)值模結(jié)果作為前提進(jìn)行連貫的數(shù)值分析，用沿水深積分的Navier-Stokes方程即淺水方程作為控制方程求解涌浪的傳播過程。通過步驟3的計算，可以得到各時刻各測點(diǎn)的水深數(shù)值，對涌浪的運(yùn)動軌跡和泛濫區(qū)域進(jìn)行數(shù)值預(yù)報，評價涌浪災(zāi)害。

進(jìn)一步地，典型滑坡體截面產(chǎn)生的初始浪高計算方法如下：

在滑坡體進(jìn)入水體之后的任意時刻，取指示函數(shù)H值分布圖。此處以二維示意圖為例，三維情況同理。取兩點(diǎn)，分別為點(diǎn)A，點(diǎn)B，不妨設(shè)H_A＜0.5(空氣)，H_B＞0.5(水體)，A點(diǎn)高程h_A，B點(diǎn)高程h_B。根據(jù)線性插值遠(yuǎn)離，容易得出，H＝0.5的點(diǎn)的高程為初始浪高h(yuǎn)＝h_O-h_Still，其中h_Still表示靜水位。

下面舉本發(fā)明的一個具體實(shí)施例：

某滑體干密度為2100kg/m³、飽和密度為2200kg/m³、水的密度為1000kg/m³、空氣的密度為1kg/m³。

步驟1、根據(jù)已有地形、水位和滑坡體形狀位置資料，做出合理的簡化，建立三維有限元模型。

模型中地形表面坐標(biāo)是從DEM數(shù)據(jù)獲得，所以坐標(biāo)系為大地坐標(biāo)系。z軸為豎直向，向上為正。覆蓋層滑坡涌浪模型結(jié)點(diǎn)41475個，單元207081個，其中滑坡體單元17261個。三維網(wǎng)格圖如圖1。邊界條件有庫岸邊界和截取的庫水邊界。假設(shè)截取的庫水為反射邊界，即涌浪傳到此處反射向回傳播；庫岸邊界為反射邊界，即涌浪傳到岸坡上發(fā)生反射，然后繼續(xù)向?qū)Π秱鞑ゲ⑴c入射波相互迭加和影響。涌浪傳播河道視為天然河道，參考相關(guān)文獻(xiàn)，粗糙系數(shù)n取0.05。

步驟2、結(jié)合Navier-Stokes方程和改進(jìn)的守恒式Level Set方法，計算流體各運(yùn)動參數(shù)，捕捉自由面。

(2-1)對土質(zhì)或松散介質(zhì)坡體，基于“流體化”的特點(diǎn)，以流體力學(xué)基本方程N(yùn)avier-Stokes方程配以非牛頓流體或擬流體本構(gòu)關(guān)系描述其物理特性和運(yùn)動過程：

進(jìn)一步地，由于將滑坡體視作似流體，三相介質(zhì)可用同一本構(gòu)模型表示：

式中，τ_y0表示流體屈服應(yīng)力，φ為內(nèi)摩擦角，I_2D是應(yīng)變張量D的第二不變量，μ為黏性系數(shù)，η為表示流體特性的參量。

(2-2)在利用Navier-Stokes方程計算出的流場信息的基礎(chǔ)上，采用改進(jìn)的守恒式雙重Level Set方法捕捉滑坡體、庫水和空氣三者的交界面，反映三者間強(qiáng)烈的流固耦合作用，真實(shí)模擬滑坡體入水產(chǎn)生涌浪的過程，典型斷面密度分布圖如圖2。改進(jìn)的守恒式雙重Level Set方法包含兩個部分：指示函數(shù)的輸運(yùn)，指示函數(shù)的重新初始化。

為了確定交界面的位置，引入標(biāo)記距離函數(shù)和指示函數(shù)H?；騂＝0.5表示界面。重新初始化方程中包含界面的法線方向，由于指示函數(shù)H本身的局限性，不適宜用于求解遠(yuǎn)離H＝0.5處的界面法線方向。選取一個臨界值δ，當(dāng)時，利用計算界面法線方向，否則用H計算法線方向。步驟(2-2)的具體過程表述如下：

(1)在步驟(2-1)求解得到的速度場的基礎(chǔ)上，求解兩個指示函數(shù)的輸運(yùn)方程，得到和H；

(2)對遠(yuǎn)離界面處的標(biāo)記距離函數(shù)進(jìn)行初始化；

(3)計算界面彌散區(qū)域各點(diǎn)的法線方向；

(4)對指示函數(shù)H進(jìn)行初始化；

(5)根據(jù)H計算靠近界面處的距離標(biāo)記函數(shù)

(6)進(jìn)行物理特性插值，確定界面位置，進(jìn)行下一步。

進(jìn)一步地，上述過程中選取的臨界值δ一般可以取0.1Δ～0.2Δ，其中Δ表示網(wǎng)格尺寸，當(dāng)時表示遠(yuǎn)離界面，表示靠近界面。

交替計算步驟(2-1)、(2-2)，可以得到三相流運(yùn)動參數(shù)、捕捉界面位置，觀察涌浪產(chǎn)生過程，選取初始涌浪作為涌浪傳播計算的初始數(shù)值，用以計算后續(xù)的涌浪傳播。

進(jìn)一步地，在三相流中，各點(diǎn)的物理參數(shù)(如密度，黏性等)可以用指示函數(shù)H表示。引入指示函數(shù)H₁,H₂，三相介質(zhì)各自區(qū)域分別記作Ω₁,Ω₂,Ω₃，界面分別用Γ₁₂,Γ₂₃,Γ₁₃表示，如圖3。P表示物理參數(shù)，則各區(qū)域參數(shù)可以用統(tǒng)一的方程表示：

P＝P₁+(P₂-P₁)H₁+(P₃-P₂)H₂

在Ω₁內(nèi)，H₁＝H₂＝0；在Ω₂內(nèi)，H₁＝1,H₂＝0；在Ω₃內(nèi)，H₁＝H₂＝1。

步驟3、采用沿水深積分的流體控制方程模擬涌浪傳播過程。

以步驟2中首浪產(chǎn)生數(shù)值模結(jié)果作為前提進(jìn)行連貫的數(shù)值分析，用沿水深積分的Navier-Stokes方程即淺水方程作為控制方程求解涌浪的傳播過程。通過步驟3的計算，可以得到各時刻各測點(diǎn)的水深數(shù)值，對涌浪的運(yùn)動軌跡和泛濫區(qū)域進(jìn)行數(shù)值預(yù)報，評價涌浪災(zāi)害。

進(jìn)一步地，典型滑坡體截面產(chǎn)生的初始浪高計算方法如下：

在滑坡體進(jìn)入水體之后的任意時刻，取指示函數(shù)H值分布圖。此處以二維示意圖為例，三維情況同理。在靠近界面處取兩點(diǎn)，分別為點(diǎn)A，點(diǎn)B，不妨設(shè)H_A＜0.5(空氣)，H_B＞0.5(水體)，A點(diǎn)高程h_A，B點(diǎn)高程h_B。根據(jù)線性插值遠(yuǎn)離，容易得出，H＝0.5的點(diǎn)的高程為初始浪高h(yuǎn)＝h_O-h_Still，其中h_Still表示靜水位。

取t＝22.13s時刻的涌浪高度作為首浪高度，如圖4，計算結(jié)果見表1。

表1 典型斷面涌浪計算參數(shù)及結(jié)果

在靠近滑坡體入水處和靠近大壩處選取兩個截面，分別布置3個測點(diǎn)。兩截面出水深隨時間分布圖分別見圖5、圖6。T＝700s時，涌浪傳播到壩址附近并形成最大浪高為2.7m。故壩址處受涌浪影響最高水位為382.7m左右，比壩頂高程(384m)約低1.3m。計算得出的滑坡涌浪高度低于大壩高度，對大壩的沖擊速度較低，在無外界不可抗力因素的影響下，滑坡涌浪對大壩存在一定的影響，但不會產(chǎn)生危害性的破壞。

綜上所述，本發(fā)明提供了一種滑坡涌浪的計算方法，該方法將涌浪產(chǎn)生和傳播作為一個統(tǒng)一連貫的過程研究，能表現(xiàn)出滑坡體、庫水和空氣三者強(qiáng)烈耦合作用。本發(fā)明不僅對庫岸邊坡失穩(wěn)破壞及涌浪次生災(zāi)害的預(yù)測具有重大的理論意義和實(shí)用價值，而且對其它類似的滑坡、泥石流、堰塞湖等地質(zhì)災(zāi)害的數(shù)值模擬亦具有一定的借鑒意義。

本發(fā)明方案所公開的技術(shù)手段不僅限于上述技術(shù)手段所公開的技術(shù)手段，還包括由以上技術(shù)特征任意組合所組成的技術(shù)方案。

以上述依據(jù)本發(fā)明的理想實(shí)施例為啟示，通過上述的說明內(nèi)容，相關(guān)工作人員完全可以在不偏離本項(xiàng)發(fā)明技術(shù)思想的范圍內(nèi)，進(jìn)行多樣的變更以及修改。本項(xiàng)發(fā)明的技術(shù)性范圍并不局限于說明書上的內(nèi)容，必須要根據(jù)權(quán)利要求范圍來確定其技術(shù)性范圍。