王梅力，祖福興，徐績清，王平義，韓林峰

(1.重慶交通大學(xué) 國家內(nèi)河航道整治工程技術(shù)研究中心，重慶 400074；2.重慶交通大學(xué) 建筑與城市規(guī)劃學(xué)院，重慶 400074；3.重慶市交通規(guī)劃勘察設(shè)計(jì)院，重慶 400074)

滑坡涌浪是滑坡體滑入江河湖海后，將所攜帶能量在水體中進(jìn)行充分轉(zhuǎn)換而激起的巨大波浪。滑坡涌浪是滑坡的次生災(zāi)害，但對涌浪生成區(qū)域和傳播區(qū)域所造成的災(zāi)害非常巨大，有時(shí)甚至超過了滑坡本身[1-4]。河道型水庫滑坡涌浪的一個(gè)顯著特點(diǎn)就是河道邊界對涌浪傳播的作用。當(dāng)滑坡涌浪傳播到達(dá)河道邊界時(shí)，部分涌浪被反射而形成反射波，而反射波與相繼而來的入射波會(huì)發(fā)生疊加，疊加后的波高會(huì)很大，往往對河道中的船舶及岸邊的建筑造成二次災(zāi)害。而河道的邊界往往是非常復(fù)雜的，河道岸邊的坡度各不相同，引起涌浪的反射率不同，形成的反射周期也會(huì)不同；直線河岸與曲線河岸會(huì)影響波浪的反射角度，反射波的方向會(huì)不同，不同的平面曲率和圓弧角度都會(huì)影響反射波周期。

目前國內(nèi)外有關(guān)河道邊界對滑坡涌浪周期的研究較少，Law 等[5]定性地對滑坡涌浪的傳播規(guī)律進(jìn)行了研究，結(jié)合試驗(yàn)數(shù)據(jù)分析推導(dǎo)出滑坡涌浪的傳播距離和衰減程度成反比的規(guī)律；Risio 等[6-7]建立了圓形岸坡的庫區(qū)滑坡涌浪模型，以此為基礎(chǔ)研究了圓形岸坡涌浪的爬高規(guī)律，對滑坡涌浪的三維流場進(jìn)行了深入研究；王育林等[8]對峽谷型河道的危巖崩滑所產(chǎn)生的靜態(tài)和動(dòng)態(tài)破壞進(jìn)行了分析，對滑坡給水工建筑物、航道及航行船舶帶來的危害進(jìn)行了論述；胡小衛(wèi)[9]考慮到庫岸滑坡的多發(fā)性與危害性，通過研究水庫滑坡成因和產(chǎn)生涌浪的各方面影響因素，以三峽水庫為依托，建立了水槽概化模型，通過控制富裕水深、有效接觸面積和滑坡坡度3 個(gè)因素，計(jì)算滑坡涌浪的各項(xiàng)特性；楊飛鵬[10]采用數(shù)值模擬方法，研究了在兩種截面形狀以及含有彎道情況下滑坡涌浪的傳播與衰減規(guī)律；謝海清等[11]利用FLOW-3D 計(jì)算仿真工具建立了數(shù)值模型，用于對狹窄型庫區(qū)河道的滑坡涌浪的產(chǎn)生、傳播進(jìn)行模擬研究，對滑坡產(chǎn)生的初始浪高和傳播過程中的浪高衰減規(guī)律進(jìn)行了公式擬合；黃筱云等[12]利用FLOW-3D 計(jì)算流體力學(xué)軟件，建立了數(shù)值模型模擬滑坡體沿著斜坡下滑后產(chǎn)生涌浪以及傳播的過程，用以研究V 型河道中產(chǎn)生的滑坡涌浪的爬高及傳播特點(diǎn)。本文基于FLOW-3D 計(jì)算流體力學(xué)軟件建立三維滑坡涌浪數(shù)值模型，模擬分析滑坡體入水點(diǎn)位于河道直線段、彎道凸岸頂點(diǎn)、彎道凹岸頂點(diǎn)情形下，復(fù)雜河道邊界對滑坡涌浪周期的影響規(guī)律，對深入揭示滑坡涌浪傳播特性及其災(zāi)害防控具有重要的意義。

1 數(shù)值模擬研究方法

本文采用的水動(dòng)力學(xué)控制方程是連續(xù)性方程和不可壓縮黏性流體運(yùn)動(dòng)的Navier-Stokes 方程，將流體視為不可壓縮黏性流體，即：

流體的動(dòng)量守恒定律由運(yùn)動(dòng)方程（N-S 方程）來描述，即：

式中：ρ 為流體密度（kg/m3）；Ax，Ay，Az分別代表流體在x，y，z 方向單元面內(nèi)流體可流過區(qū)域的面積分?jǐn)?shù)；u，v，w 分別為x，y，z 方向上的速度分量（m/s）；VF為單元內(nèi)流體可流動(dòng)區(qū)域的體積分?jǐn)?shù)；Gx，Gy，Gz分別為流體在x，y，z 方向上的質(zhì)量力加速度（m/s2）；fx，fy，fz分別為流體x，y，z 方向上的黏滯力加速度（m/s2）。

FLOW-3D 中的紊流模型包含普朗特混合長度模型、k-ε 方程模型、LES 模型、ε 方程模型、RNG k-ε 模型5 種紊流模型來模擬紊流。本文采用RNG k-ε 紊流模型[13]計(jì)算涌浪的產(chǎn)生和傳播。采用有限差分法對控制方程進(jìn)行離散。

VOF 法是一種處理復(fù)雜自由表面的有效方法，其基本思想是：在任何一個(gè)計(jì)算區(qū)域內(nèi)，水和氣體的體積分?jǐn)?shù)之和為1，定義F 為計(jì)算區(qū)域內(nèi)水占區(qū)域的流體分?jǐn)?shù)，在FLOW-3D 中，對于某個(gè)計(jì)算網(wǎng)格單元而言：F=1表 示計(jì)算網(wǎng)格單元內(nèi)充滿水體， F=0表 示計(jì)算網(wǎng)格單元內(nèi)充滿氣體， 0

2 模型的建立與驗(yàn)證

2.1 模型的建立

選取與長江三峽庫區(qū)萬州某河段原型一致的物理模型[14]。以試驗(yàn)中的1 個(gè)工況為例，物理模型尺寸為560 m（河寬）×51.8 m（水深），對岸岸坡為20°。河道凹岸彎曲半徑為1 120 m，河道凸岸彎曲半徑為560 m，河道彎曲角度為90°，滑坡面傾角選取40°。滑坡體尺寸選取70 m×35 m×14 m。構(gòu)建的河道和滑坡體三維模型如圖1(a)所示。模型采用六面體結(jié)構(gòu)化網(wǎng)格，在滑塊落水處采用漸變網(wǎng)格加密，x 方向網(wǎng)格最小尺寸0.1 m，最大0.2 m；y 方向網(wǎng)格尺寸均為0.1 m；z 方向尺寸均為0.1 m，如圖1(b)所示。

圖1 滑坡涌浪數(shù)值模擬三維模型及網(wǎng)格劃分Fig.1 3D model and mesh of landslide surge numerical simulation

網(wǎng)格塊包含的區(qū)域?yàn)橛?jì)算區(qū)域，由于該網(wǎng)格塊為長方體，因此6 個(gè)面都需要給定邊界條件，該區(qū)域上表面為自由表面，采用壓力邊界條件，給定水面壓力p=1.013×105Pa；由于物理模型試驗(yàn)是在靜水中進(jìn)行，因此其余表面邊界條件采用無滑移壁面WALL 條件，水工建筑物視為混凝土材質(zhì)，糙率為0.014。

2.2 模型驗(yàn)證

以物理模型試驗(yàn)結(jié)果[14]對數(shù)值模擬結(jié)果進(jìn)行驗(yàn)證。選擇在順直河道進(jìn)行模型試驗(yàn)，河道實(shí)際水位取145 m，對應(yīng)河道實(shí)際水深51.8 m，滑坡面傾角選擇40°，滑坡體尺寸選取70 m×35 m×14 m，波高測點(diǎn)共布置16 個(gè)，如圖2 所示。各測點(diǎn)物理試驗(yàn)換算成原型的結(jié)果與數(shù)值模擬試驗(yàn)結(jié)果對比見圖3。

圖2 凸岸滑坡涌浪測點(diǎn)平面布置（單位：m）Fig.2 Layout plan of swell measuring points of convex landslide (unit: m)

從圖3 可以看出：數(shù)模計(jì)算結(jié)果與物模實(shí)測結(jié)果的最大偏差為4.94%，平均偏差為4.06%。說明FLOW-3D對滑坡涌浪的產(chǎn)生及傳播的模擬具有較高精度，可用于復(fù)雜邊界對滑坡涌浪周期的影響研究。

3 河道邊界對滑坡涌浪周期的影響

3.1 岸坡坡度對滑坡涌浪周期的影響規(guī)律

為了模擬分析岸坡坡度對滑坡涌浪周期的影響規(guī)律，選取滑坡體處河對岸岸坡坡度為20°和40°進(jìn)行對比。以圖2 滑坡入水點(diǎn)向?qū)Π蹲龃咕€，在垂線上靠近滑坡入水點(diǎn)處設(shè)置A1 測點(diǎn)（對應(yīng)圖2 中1#測點(diǎn)），在垂線上靠近對岸坡腳處設(shè)置A2 測點(diǎn)（對應(yīng)圖2 中5#測點(diǎn)）。A1 測點(diǎn)靠近滑坡入水點(diǎn)，此處的涌浪周期主要受入射波控制，反射波起次要作用。A2 測點(diǎn)靠近河對岸岸坡坡腳處，此處的涌浪周期主要受反射波控制，入射波起次要作用。試驗(yàn)完成后，提取A1 測點(diǎn)和A2 測點(diǎn)隨時(shí)間的水位變化值，利用傅里葉自相關(guān)離散型光譜頻率函數(shù)[15]，通過Matlab 編程，繪制A1 測點(diǎn)和A2 測點(diǎn)頻譜圖（圖4）。

從圖4 可以看出，入射波起主導(dǎo)作用的A1 測點(diǎn)，包括兩個(gè)控制波周期，即主控制波周期T11（T11=1/0.09=11.11 s）和輔控制波周期T12（T12=1/0.18=5.56 s）。反射波起主導(dǎo)作用的A2 測點(diǎn)，包括3 個(gè)控制波周期，即T21（T21=1/0.04=25.00 s），T22（T22=1/0.11=9.09 s）和T23（T23=1/0.18=5.56 s），3 個(gè)控制波周期中T22的波譜強(qiáng)度最大。A2 測點(diǎn)的T23周期和A1 測點(diǎn)的T12周期相同，涌浪波譜強(qiáng)度前者比后者低，說明A2 測點(diǎn)的T23周期的涌浪是由A1 測點(diǎn)的T12周期的涌浪傳播衰減而來。而A2 測點(diǎn)的T21和T22周期的涌浪是由A1 測點(diǎn)的T11周期演變而來。A1 測點(diǎn)的T11周期的涌浪到達(dá)對岸岸坡后分成了兩部分：一部分爬高后反射回來，經(jīng)爬高反射后造成這部分周期延長，就成為了A2 測點(diǎn)的周期為T21的波浪；另一部分破碎后反射回來，經(jīng)破碎反射后造成這部分周期縮短，就成了A2 測點(diǎn)的周期為T22的波浪。對比20°和40°兩種不同的岸坡坡度，前者的破碎后反射波浪較后者的大，表明岸坡坡度越小，涌浪破碎數(shù)量越多。

圖4 不同岸坡坡度作用下A1 及A2 測點(diǎn)頻譜Fig.4 Spectra of point A1 and point A2 under the action of different bank slope angles

3.2 彎曲河道凹岸對滑坡涌浪周期的影響規(guī)律

3.2.1 凹岸平面曲率對滑坡涌浪周期的影響規(guī)律

選取凹岸彎曲半徑分別為1 330 和1 540 m，即平面曲率1/1 330 和1/1 540 進(jìn)行對比試驗(yàn)?；氯胨c(diǎn)選取彎道凸岸中間點(diǎn)臨水處，以滑坡入水點(diǎn)向?qū)Π蹲龃咕€，在垂線上靠近滑坡入水點(diǎn)處設(shè)置B1 測點(diǎn)，在垂線上靠近對岸坡腳處設(shè)置B2 測點(diǎn)。利用前述方法繪制1#和2#測點(diǎn)頻譜圖（圖5）。

從圖5 可見，入射波起主導(dǎo)作用的B1 測點(diǎn)，包括兩個(gè)控制波周期，即主控制波周期T11（T11=1/0.08=12.50 s）和輔控制波周期T12（T12=1/0.18=5.56 s）。反射波起主導(dǎo)作用的B2 測點(diǎn)，包括3 個(gè)控制波周期，即T21（T21=1/0.04=25.00 s），T22（T22=1/0.11=9.09 s）和T23（T23=1/0.18=5.56 s），3 個(gè)控制波周期中T22的波譜強(qiáng)度最大。B2 測點(diǎn)的T23周期和B1 測點(diǎn)的T12周期相同，涌浪波譜強(qiáng)度前者比后者低，說明B2 測點(diǎn)的T23周期的涌浪是由B1 測點(diǎn)的T12周期的涌浪傳播衰減而來。而B2 測點(diǎn)的T21周期和T22周期的涌浪是由B1 測點(diǎn)的T11周期演變而來。B1 測點(diǎn)的T11周期的涌浪到達(dá)對岸岸坡后分成了兩部分：一部分爬高后反射回來，經(jīng)爬高反射后造成這部分周期延長，就成為了B2 測點(diǎn)的周期為T21的波浪；另一部分破碎后反射回來，經(jīng)破碎反射后造成這部分周期縮短，就成了B2 測點(diǎn)的周期為T22的波浪。對比平面曲率1/1 330 和1/1 540 的凹岸河道邊界，兩者的最大主控制波周期均是由T11（T11=1/0.08=12.50 s）變化到T22（T22=1/0.11=9.09 s），表明不同平面曲率的凹岸河道邊界，對涌浪周期的影響是相同的，均是整體性略微縮短。

3.2.2 凹岸圓弧角度對滑坡涌浪周期的影響規(guī)律 選取凹岸圓弧角度分別為30°和90°進(jìn)行對比試驗(yàn)。B1 測點(diǎn)和B2 測點(diǎn)的位置同3.2.1，利用前述方法繪制B1 測點(diǎn)和B2 測點(diǎn)頻譜圖（圖6）。

圖5 不同凹岸平面曲率作用下B1 及B2 測點(diǎn)頻譜Fig.5 Spectra of point B1 and point B2 under the action of plane curvature of different concave banks

從圖6 可以看出，入射波起主導(dǎo)作用的B1 測點(diǎn)，包括兩個(gè)控制波周期，即主控制波周期T11（T11=1/0.09=11.10 s）和輔控制波周期T12（T12=1/0.18=5.56 s）。反射波起主導(dǎo)作用的B2 測點(diǎn)，包括3 個(gè)控制波周期，即T21（T21=1/0.04=25.00 s），T22（T22=1/0.1=10.00 s）和T23（T23=1/0.18=5.56 s），3 個(gè)控制波周期中T22的波譜強(qiáng)度最大。B2 測點(diǎn)的T23周期和B1 測點(diǎn)的T12周期相同，涌浪波譜強(qiáng)度前者比后者低，說明B2 測點(diǎn)的T23周期的涌浪是由B1 測點(diǎn)的T12周期的涌浪傳播衰減而來。而B2 測點(diǎn)的T21周期和T22周期的涌浪是由B1 測點(diǎn)的T11周期演變而來。B1 測點(diǎn)的T11周期的涌浪到達(dá)對岸岸坡后分成了兩部分：一部分爬高后反射回來，經(jīng)爬高反射后造成這部分周期延長，就成為了B2 測點(diǎn)的周期為T21的波浪；另一部分破碎后反射回來，經(jīng)破碎反射后造成這部分周期縮短，就成了B2 測點(diǎn)的周期為T22的波浪。對比圓弧角度為30°和90°的凹岸河道邊界，兩者最大的主控制波周期均是由T11（T11=1/0.09=11.10 s）變化到T22（T22=1/0.1=10.00 s），表明不同圓弧角度的凹岸河道邊界，對涌浪周期的影響是相同的，均是整體性略微縮短。

圖6 不同凹岸圓弧角度作用下B1 及B2 測點(diǎn)頻譜Fig.6 Spectra of point B1 and point B2 under action of different concave bank arc angles

3.3 彎曲河道凸岸對滑坡涌浪周期的影響規(guī)律

選取凸岸圓弧角度分別為60°和90°進(jìn)行對比試驗(yàn)。將滑坡入水點(diǎn)由凸岸中間坡腳臨水處移至凹岸中間坡腳臨水處，以滑坡入水點(diǎn)向?qū)Π蹲龃咕€，在垂線上靠近滑坡入水點(diǎn)處設(shè)置C1 測點(diǎn)，在垂線上靠近對岸坡腳處設(shè)置C2 測點(diǎn)。同前述方法，繪制C1 測點(diǎn)和C2 測點(diǎn)頻譜圖（圖7）。

從圖7 可以看出，入射波起主導(dǎo)作用的C1 測點(diǎn)，包括兩個(gè)控制波周期，即第一控制波周期T11（T11=1/0.08=12.50 s）和第二控制波周期T12（T12=1/0.2=5.00 s）。反射波起主導(dǎo)作用的C2 測點(diǎn)，包括3 個(gè)控制波周期，即T21（T21=1/0.04=25.00 s），T22（T22=1/0.11=9.09 s）和T23（T23=1/0.18=5.56 s）。涌浪到達(dá)對岸岸坡后，由于凸岸邊界的特殊性，受到的影響異常激烈。但總的來說分成了兩部分：一部分爬高后反射回來，經(jīng)爬高反射后造成這部分周期延長；另一部分破碎后反射回來，經(jīng)破碎反射后造成這部分周期縮短。

圖7 凸岸作用下C1 及C2 測點(diǎn)頻譜Fig.7 Spectra of point C1 and point C2 under convex bank action

4 結(jié) 語

（1）選取與長江三峽庫區(qū)萬州某河段為原型一致的物理模型，基于計(jì)算流體力學(xué)FLOW-3D 軟件，建立了滑坡涌浪三維數(shù)值模型，經(jīng)與物理模型試驗(yàn)驗(yàn)證符合性較好。按照滑坡入水點(diǎn)位于河道直線段、彎道凸岸頂點(diǎn)、彎道凹岸頂點(diǎn)3 種工況，計(jì)算模擬分析了復(fù)雜河道邊界對滑坡涌浪周期的影響規(guī)律。

（2）靠近滑坡入水點(diǎn)的涌浪周期主要受入射波控制，反射波起次要作用；靠近滑坡體對岸岸坡坡腳處的涌浪周期主要受反射波控制，入射波起次要作用。對比20°和40°兩種不同的岸坡坡度，前者的破碎后反射波浪較后者的大，表明岸坡坡度越小，涌浪破碎數(shù)量越多。

（3）不同平面曲率的凹岸河道邊界，對涌浪周期的影響相同，均是整體性略微縮短。不同圓弧角度的凹岸河道邊界，對涌浪周期的影響也均是整體性略微縮短。

（4）由于凸岸邊界的特殊性，涌浪到達(dá)對岸岸坡后，總體上分成了兩部分：一部分爬高后反射回來，延長了波周期；另一部分破碎后反射回來，這部分的波周期縮短。