岳 霞，戴福初，朱 杰，朱雨軒

(北京工業(yè)大學(xué) 建筑工程學(xué)院，北京 100124)

1 研究背景

滑坡是最具危害性的地質(zhì)災(zāi)害之一，而滑坡涌浪是滑坡災(zāi)害由點(diǎn)到面的拓展，受災(zāi)范圍更廣、造成的損失更大[1]。隨著國內(nèi)外水利事業(yè)的大規(guī)模發(fā)展，越來越多的水電站修建在高山之中，庫水位的頻繁變動(dòng)導(dǎo)致滑坡活動(dòng)加劇，一旦發(fā)生滑坡涌浪將嚴(yán)重威脅庫區(qū)人民的生命財(cái)產(chǎn)安全。例如，1963年意大利北部的Vajont水庫滑坡，庫區(qū)左岸27億m3的巖土體以25～30 m/s的速度沖入水庫，造成高達(dá)125 m的沖擊波，2 500萬m3的洪水順流而下，造成下游隆加羅市2 000余人喪生[2-3]。目前庫區(qū)沿岸存在大量潛在滑坡，因此開展滑坡涌浪研究，對潛在滑坡涌浪災(zāi)害進(jìn)行預(yù)測具有重要的意義。

滑坡運(yùn)動(dòng)預(yù)測是建立在對滑坡運(yùn)動(dòng)機(jī)理理解的基礎(chǔ)上，Hungr[4]和Zhao等[5]將現(xiàn)有的滑坡運(yùn)動(dòng)模型分為2類，將滑坡理想化為單點(diǎn)的集中質(zhì)量模型和連續(xù)介質(zhì)模型。第一個(gè)集中質(zhì)量模型是由Heim[6]提出的“雪橇模型”，模型假定滑坡運(yùn)動(dòng)過程中的能量損失都是由摩擦造成的。Sassa等[7]基于塊體理論和滑坡運(yùn)動(dòng)過程中材料的連續(xù)性原理提出了滑坡運(yùn)動(dòng)模型，模型能對流動(dòng)性土質(zhì)滑坡運(yùn)動(dòng)進(jìn)行準(zhǔn)確反演。汪發(fā)武等[8-9]、郭劍等[10-12]、沈偉等[13-16]分別對Sassa模型中的表觀摩擦角、地轉(zhuǎn)偏向力對滑體運(yùn)動(dòng)的影響、側(cè)壓力的選取、滑床的侵蝕裹挾等問題進(jìn)行了研究，使模型的輸入?yún)?shù)更容易獲取、模擬效果更加準(zhǔn)確。本文采用的類流體滑坡模型是以Saasa的滑坡運(yùn)動(dòng)模型為基礎(chǔ)，并考慮了滑坡-涌浪之間的強(qiáng)烈耦合作用力。

現(xiàn)有的流體模型已經(jīng)比較完善，但滑坡涌浪運(yùn)動(dòng)交界面處滑坡與流體之間復(fù)雜的耦合作用力仍是滑坡涌浪模型計(jì)算的難點(diǎn)，一種理想的簡化方法是預(yù)先定義滑坡的運(yùn)動(dòng)，但這種方法也只能考慮單相耦合。Quecedo等[17-18]針對這個(gè)問題提出采用多相流模型來反映流體與滑坡的相互作用，計(jì)算中所有相均被視為流體。Mao等[19]建立了由一般本構(gòu)方程計(jì)算的多相流模型模擬滑坡、流體、空氣，進(jìn)一步深化了對該多相流模型的認(rèn)識(shí)。Pastor等[20-22]在水深積分方程和SPH計(jì)算方法的基礎(chǔ)上建立滑坡涌浪計(jì)算模型，模型中滑坡體的孔隙壓力會(huì)隨著滑坡的運(yùn)動(dòng)而改變，能更好地反映滑體的巖土特性。Ward等[23]針對類流體滑坡提出了與SPH方法類似的“海嘯球”模型，其特點(diǎn)是將水球柱作為基本單元而非單個(gè)小球，得到的結(jié)果更加準(zhǔn)確，之后肖麗莉等[24-26]、王佳佳等[27-28]又提出了可重復(fù)分割的“Tsunami Squares”模型，模型使計(jì)算精度不受塊體大小的限制。

本文根據(jù)Sassa滑坡運(yùn)動(dòng)模型、COMCOT海嘯模型、流體力學(xué)建立了流動(dòng)性滑坡涌浪動(dòng)力學(xué)模型，該模型以多相流理論和淺水方程為基礎(chǔ)，將滑坡、水體兩相視為流體，然后由水深積分方程計(jì)算波高的變化。郭劍等[10]也曾做過此類研究，但尚未考慮到滑體與水體之間的強(qiáng)烈耦合作用力。本文在Sassa土質(zhì)滑坡運(yùn)動(dòng)模型和郭劍流動(dòng)性滑坡涌浪動(dòng)力學(xué)模型的基礎(chǔ)上推導(dǎo)了涉水滑坡與水體之間的強(qiáng)烈耦合作用力，然后根據(jù)開源COMCOT程序進(jìn)行程序的編寫，并將COMCOT開源程序的開放邊界、移動(dòng)邊界、波浪破碎應(yīng)用到模型中，使模型能夠考慮滑坡涌浪強(qiáng)烈耦合作用力、波浪破碎等因素對涌浪高度的影響。

2 數(shù)學(xué)模型

2.1 控制方程

2.1.1 滑坡模型

類流體滑坡模型的控制方程是由動(dòng)量方程和連續(xù)方程組成，由流體連續(xù)性原理可知，控制單元體的厚度變化由兩側(cè)的流量差決定，因此模型的連續(xù)方程為

(1)

式中：h表示滑體厚度；Q1x、Q1y分別表示滑坡在x、y方向的體積通量。

本文的滑坡模型與Sassa模型的主要不同點(diǎn)在于涉水滑坡運(yùn)動(dòng)需要考慮滑坡與水體之間的耦合作用力，下面根據(jù)x方向滑坡的受力推導(dǎo)出滑坡x方向的動(dòng)量方程，單元體的受力情況如圖1所示。

圖1 滑坡控制單元體x方向的受力分析示意圖Fig.1 Forces acting on the slice of landslide in x direction

由圖1可知，滑坡單元體主要受到重力、滑床的支持力、滑體的側(cè)壓力、滑床的摩擦力、水體的拖曳力，其中滑坡重力G不會(huì)在xOy平面產(chǎn)生分力，滑床的支持力是滑體所受的作用力在垂直于滑面方向的投影，側(cè)壓力、摩擦力、拖曳力的作用方向均平行于滑面，在垂直于滑面方向的投影為0，支持力主要由豎直向下的重力產(chǎn)生，因此，滑床的支持力在x方向的分力FNx為

(2)

其中，

q=tan2α+tan2β。

式中：α、β分別為滑體底面與xOz、yOz平面的夾角；m為滑坡的質(zhì)量；g為重力加速度；ρs為滑坡的密度。

滑坡側(cè)壓力主要由滑坡周圍的控制單元體產(chǎn)生，本文利用靜止土壓力法計(jì)算滑坡單元體受到的側(cè)壓力P，即

(3)

式中：dh為滑體厚度在x方向的增量；k為側(cè)壓力系數(shù)，k=1-sinθ，θ為滑坡體的內(nèi)摩擦角。

滑坡運(yùn)動(dòng)時(shí)滑坡底部受到的摩擦力由摩爾-庫倫準(zhǔn)則計(jì)算，因此基底摩擦力R和其在x方向的分力Rx分別為：

(4)

Rx=DxR。

(5)

式中：R為滑坡運(yùn)動(dòng)時(shí)的摩擦力；Ns為重力沿滑面方向的分量；φa為滑坡的有效摩擦角；c為滑坡的黏聚力；Dx為速度的方向向量，負(fù)號表示阻力方向與速度方向相反，其中，v1z=v1xtanα+v1ytanβ，v1x、v1y、v1z分別為滑坡在x、y、z方向的平均速度。

根據(jù)流體力學(xué)[29]可知，水下滑坡運(yùn)動(dòng)時(shí)，拖曳力的大小與滑體-水體之間相對運(yùn)動(dòng)速度的平方成正比，方向與滑坡運(yùn)動(dòng)的方向相反，水下滑坡運(yùn)動(dòng)受到拖曳力Fw為

Fw=CDρw|v1x-v2x>|(v1x-v2x)S=

(6)

式中：CD為水體的黏滯阻力系數(shù)；ρw為水體的密度；v2x、v2y分別為水體x、y方向的速度；S為滑坡與水體的接觸面積；q0=tan2α0+tan2β0，α0、β0分別表示滑坡上表面與xOz平面、yOz平面的夾角。

由牛頓第二定律可知Fx=max，由流體動(dòng)力學(xué)可知類流體滑坡的加速度ax為

(7)

聯(lián)立式(2)—式(7)可得x方向滑坡單元體的動(dòng)量方程為

(8)

同理，可得y方向滑坡單元體的動(dòng)量方程為

(9)

2.1.2 涌浪模型

涌浪模型的控制方程也是由連續(xù)方程(式(10))和動(dòng)量方程(式(11)、式(12))組成，式(11)、式(12)右邊的第3項(xiàng)是在COMCOT海嘯模型控制方程的基礎(chǔ)上根據(jù)流體動(dòng)力學(xué)與滑坡模型中滑坡涌浪的作用力與反作用力推導(dǎo)得到的，其它項(xiàng)均與COMCOT海嘯模型相同。具體方程如下所示。

(10)

(11)

(12)

式中：H=h1+η，η表示水面位移，h1表示水底到靜止水面的距離；n為曼寧系數(shù)；Q2x、Q2y分別表示涌浪x、y方向的體積通量。

2.2 離散形式

滑坡、涌浪模型的控制方程都是通過有限差分形式離散，對流項(xiàng)采用三階迎風(fēng)差分法離散，動(dòng)量方程右側(cè)的其他項(xiàng)采用中心差分法離散。圖2顯示了一個(gè)離散單元的計(jì)算方式，滑坡厚度或涌浪波高h(yuǎn)位于單元中心處，滑坡或涌浪x、y方向的通量位于單元網(wǎng)格線上，滑坡厚度或涌浪波高h(yuǎn)的值由x、y方向的通量計(jì)算。

圖2 離散網(wǎng)格Fig.2 Discrete grids

2.3 邊界條件

滑坡、涌浪體積通量只在濕網(wǎng)格(總水深>0，滑坡厚度>0)計(jì)算更新，然后由連續(xù)方程結(jié)合移動(dòng)邊界條件確定計(jì)算區(qū)域下一個(gè)時(shí)間層的水面位移、滑坡運(yùn)動(dòng)、干網(wǎng)格(總水深等于0，滑坡厚度等于0)。移動(dòng)邊界的具體形式如圖3所示，若i點(diǎn)滑坡厚度(水深)>0，i+1點(diǎn)滑坡厚度(水深)等于0，且i點(diǎn)原始地形加上滑坡厚度(水深)后的高程大于i+1點(diǎn)的高程，則邊界移動(dòng)到i+1點(diǎn)(圖3(b)移動(dòng)邊界)，否則邊界保持在i點(diǎn)(圖3(a)未移動(dòng)邊界)。

圖3 移動(dòng)邊界方案Fig.3 Scheme of moving boundary

3 模型驗(yàn)證

3.1 龔家方滑坡

2008年11月23日重慶市巫山縣發(fā)生的龔家方滑坡，位于三峽大壩上游120 km、巫山縣城下游4 km處。已有研究[30-31]對龔家方滑坡涌浪的現(xiàn)場調(diào)查表明，龔家方滑坡發(fā)生時(shí)，水庫蓄水位為173 m，激起的初始涌浪高度約10 m，對岸涌浪爬坡浪高達(dá)12 m，涌浪的具體沿岸爬坡高度如圖4[26]、表1[26, 30]所示。已有研究[10，26]表明，龔家方滑坡在平面上呈簸箕形，前、后緣高程分別為120 m、400 m，滑體長290 m，寬40～60 m，平均厚度20～30 m，入水體積約38×104m3?；禄瑒?dòng)后，殘留滑體分布于高程30～220 m之間，且堆積體的最大厚度可達(dá)90 m?；乱灾酗L(fēng)化和強(qiáng)風(fēng)化的三疊系灰?guī)r、泥灰?guī)r為主，粒徑多<30 cm，D50=25 cm，本文的類流體滑坡涌浪動(dòng)力學(xué)模型可以應(yīng)用于龔家方滑坡涌浪的模擬。

表1 龔家方滑坡涌浪爬坡高度一覽表[26, 30]Table 1 Observed climbing height of Gongjiafanglandslide surge[26, 30]

3.2 參數(shù)選取

龔家方滑坡發(fā)生后的數(shù)字高程模型(DEM)數(shù)據(jù)來源于地理空間數(shù)據(jù)云平臺(tái)(http:∥www.gscloud.cn)，滑坡發(fā)生前的地形參考滑坡兩側(cè)高程對其進(jìn)行復(fù)原，高程<173 m的點(diǎn)設(shè)為水體覆蓋區(qū)域，水深為地形高程與173 m的差值。模擬的計(jì)算時(shí)間步長為0.1 s，總計(jì)算時(shí)間為3 000 s，差分網(wǎng)格的實(shí)際大小為30 m×30 m，網(wǎng)格個(gè)數(shù)為 394×145，滑坡、水體的初始狀態(tài)均為靜止。根據(jù)已有研究資料[10, 32-34]和反復(fù)試算得到模擬需要的參數(shù)，其中滑體與滑帶的黏聚力會(huì)隨著下滑速度的增大而逐漸喪失[26]，結(jié)合Sassa模型在實(shí)際案例中的應(yīng)用[8]和反復(fù)試算認(rèn)為龔家方滑坡運(yùn)動(dòng)過程中的黏聚力c值為0。本文考慮的是滑坡涌浪運(yùn)動(dòng)過程，表2給出了滑坡涌浪運(yùn)動(dòng)狀態(tài)下的具體參數(shù)。

表2 龔家方滑坡涌浪模擬參數(shù)Table 2 Parameters for modeling of Gongjiafanglandslide surge

圖5 不同時(shí)刻滑坡體厚度分布Fig.5 Variation of landslide thickness with time

3.3 模擬結(jié)果

圖5顯示了不同時(shí)間滑坡體厚度的分布情況。由圖5(a)可知滑坡的初始最大厚度為34 m；由圖5(b)可知滑坡滑動(dòng)初期，滑坡厚度大的地方受力大，最先開始運(yùn)動(dòng)，但滑坡整體的面積未發(fā)生較大改變；由圖5(d)可知在20 s時(shí)，滑坡運(yùn)動(dòng)到平緩地形時(shí)，滑體開始堆積，最大滑坡厚度位于高程140 m地形較平緩處；圖5(f)給出了滑坡最終堆積體形態(tài)，堆積體后緣高程在230 m左右，與220 m的野外觀測值接近，堆積體的最大厚度為55 m左右。由此可知模型能較好地反映滑坡運(yùn)動(dòng)特征和滑體的堆積形態(tài)。

圖6顯示了滑坡發(fā)生區(qū)不同時(shí)刻的河道剖面，在滑坡運(yùn)動(dòng)停止之前時(shí)，滑坡發(fā)生區(qū)河底的高程一直在上升，水體的高程也會(huì)隨之上升，此時(shí)近岸的波幅都比較高。由圖6(b)可知只有當(dāng)滑坡的入水體積和滑速達(dá)到一定條件時(shí)河道的涌浪才會(huì)向外傳播，在7 s時(shí)近岸產(chǎn)生的波高還未滿足向外傳播的條件；由圖6(c)可知在17 s時(shí)，滑坡產(chǎn)生的涌浪已經(jīng)傳播到河道中部，此時(shí)近岸的最大涌浪爬坡高度為24.3 m；由圖6(e)可知在35 s時(shí)，滑坡基本停止運(yùn)動(dòng)，此時(shí)的涌浪波已經(jīng)到達(dá)了對岸，并在對岸產(chǎn)生了一定的爬坡高度；由圖6(f)可知在300 s時(shí)，滑坡發(fā)生區(qū)河段的涌浪高度已經(jīng)降低到2 m以下，河道水面已無較大的災(zāi)害涌浪。

圖7顯示了不同時(shí)刻的河道涌浪波高，由圖7(a)可知在40 s時(shí)，產(chǎn)生的第一個(gè)涌浪波以半圓形的形式向河道傳播，并且已經(jīng)到達(dá)了河道對岸，在40s以后受地形影響涌浪開始沿著河道向上下游進(jìn)行傳播；由圖7(b)可知在240 s時(shí)，涌浪已經(jīng)到達(dá)巫山城區(qū)，涌浪在巫山城區(qū)的最大爬坡高為1.6 m；由圖7(c)可知在1 000 s時(shí)，涌浪已經(jīng)傳播到整個(gè)計(jì)算區(qū)域，計(jì)算域內(nèi)的涌浪已經(jīng)降低到1.5 m以下，而滑坡發(fā)生區(qū)河段的水面波動(dòng)仍然比較劇烈；由圖7(d)可知在計(jì)算結(jié)束3 000 s時(shí)，計(jì)算域內(nèi)的最大波高在0.7 m左右，但達(dá)到完全穩(wěn)定還需要一定的時(shí)間。同時(shí)由于岸邊波能的集聚，爬坡涌浪要比河道內(nèi)的涌浪波更高，模擬得到的涌浪傳播情況與已有研究[24]比較吻合。

圖6 滑坡發(fā)生區(qū)河道橫剖面變化過程Fig.6 Variations of landslide (gray), water (green), and amplified wave (red) at cross-section ofriver channel in landslide area

圖7 不同時(shí)刻的河道涌浪波高Fig.7 Plan map showing variation of wave height of the surge at different instances

圖8 涌浪沿岸爬坡高度Fig.8 Climbing height of the surge along the reservoir bank

圖8顯示了涌浪沿岸爬坡高度模擬值與現(xiàn)場觀測值的對比。為了使結(jié)果能更準(zhǔn)確地反映爬坡值的具體變化，涌浪爬坡高度模擬值是對每一個(gè)岸邊點(diǎn)都進(jìn)行了取值，涌浪爬坡高度的現(xiàn)場觀測值是結(jié)合圖4、表1中的數(shù)據(jù)得到的。

由圖8(a)可知近岸的最大涌浪爬坡高度為24.3 m，由圖8(b)可知滑坡正對岸的涌浪爬坡高度在10～18.9 m之間，模擬結(jié)果與現(xiàn)場觀測吻合較好。由沿程涌浪爬坡高度曲線可知，隨著沿程距離的增加，涌浪爬坡高度不斷減小，但也存在波動(dòng)現(xiàn)象，主要與觀測點(diǎn)的具體地形條件有關(guān)。

4 結(jié) 論

本文基于Sassa準(zhǔn)三維滑坡運(yùn)動(dòng)模型和開源COMCOT海嘯模型提出了將滑坡-涌浪運(yùn)動(dòng)視為一個(gè)連續(xù)過程、可以反映滑坡與水體之間強(qiáng)烈耦合作用力的流動(dòng)性滑坡涌浪動(dòng)力學(xué)模型，并利用三峽庫區(qū)龔家方滑坡涌浪數(shù)據(jù)進(jìn)行了模型驗(yàn)證，取得了以下認(rèn)識(shí)：

(1)通過對滑體和水體的受力分析，推導(dǎo)了模型的動(dòng)量方程，由物質(zhì)的連續(xù)性原理得到了連續(xù)方程；采用有限差分法對模型的控制方程進(jìn)行了離散，并利用Fortran語言編寫了差分方程的計(jì)算程序，從而建立了流動(dòng)性滑坡涌浪動(dòng)力學(xué)模型。

(2)為了驗(yàn)證模型的有效性，將模型應(yīng)用到龔家方滑坡涌浪案例進(jìn)行模擬，得到滑坡堆積體后緣高程在230 m左右，近岸涌浪的最大爬坡為24.3 m，正對岸的涌浪爬坡高度在10～18.9 m之間，兩岸的涌浪爬坡高度與野外觀測值情況比較吻合。

(3)通過模型輸出的滑坡厚度分布情況、庫區(qū)各點(diǎn)涌浪的瞬時(shí)高度、庫岸邊坡的最大涌浪爬坡高度、河道剖面滑坡涌浪運(yùn)動(dòng)情況等，可以全面了解滑坡涌浪運(yùn)動(dòng)的全過程，下一步將利用模型對其他庫區(qū)的潛在滑坡涌浪進(jìn)行模擬，為滑坡涌浪災(zāi)害預(yù)警提供參考依據(jù)。