崩塌落石產(chǎn)生涌浪的流固耦合運動分析

黃波林，王世昌，殷躍平，劉廣寧，陳小婷

1.武漢地質(zhì)調(diào)查中心，武漢 430205

2.中國地質(zhì)大學(xué)工程學(xué)院，武漢 430074

3.中國地質(zhì)環(huán)境監(jiān)測院，北京 100081

0 引言

庫區(qū)崩塌落石不僅對滾石運動路徑上的物體存在威脅，而且入水后會產(chǎn)生涌浪，威脅航道安全。例如在瑞士的Lake Uri，由于1.6萬 m3的危巖崩塌后會威脅下方公路和形成涌浪，1992年有關(guān)人員對其進行了爆破，蘇黎世大學(xué)的水力學(xué)實驗室（VAW laboratory）對涌浪進行了監(jiān)測[1]。

落石運動過程十分復(fù)雜，受眾多因素影響。如斜坡的微地貌、坡面的物質(zhì)組成、落石的形態(tài)等條件都會影響落石的運動過程[2-4]。控制滾石動力學(xué)特征的主要參數(shù)是碰撞恢復(fù)系數(shù)[5]。秦志英等[6]、何思明等[7]、沈均等[8]、楊海清等[9]在理論上對碰撞恢復(fù)系數(shù)進行了分析。C.Mangwandi 等[10]、Rocscience[11]、A.Azzoni等[12]、R.W.Day[13]、章廣成等[14]通過現(xiàn)場試驗給出了碰撞恢復(fù)系數(shù)的建議值。目前有許多模擬落石碰撞過程的專業(yè)軟件，如 CRSP[15]，Rocfall[11]，STONE[16]，CADMA[17]等。但這些軟件只能計算滾石的運動特征，不能與流體耦合進行涌浪計算。

滾石以一定速度入水后會形成涌浪，滾石的幾何形態(tài)與入水姿態(tài)都會對涌浪產(chǎn)生影響，具有強烈的三維特性。進行三維流體力學(xué)與固體力學(xué)耦合分析可獲得從微觀到宏觀的直觀結(jié)果，是目前涌浪研究的熱點，也是未來研究的發(fā)展趨勢。但由于軟硬件制約，國內(nèi)鮮有文獻對此進行報導(dǎo)。在國外，B.H.Choi等[18]采用RANS系統(tǒng)進行了三維孤立波爬高分析；Silvia等[19]利用移動墻采用三維淺水波模型模擬了瓦伊昂滑坡涌浪事件；F.Montagna等[20]利用Flow-3D進行了三維海島滑坡涌浪研究，并與物理試驗進行對比分析；M.Pastor等[21]耦合三維SPH法和FEM法進行了幾個災(zāi)難性滑坡涌浪的研究；Stéphane Abadie等[22]采用兩相流的 N-S方程模擬研究了三維滑塊入水產(chǎn)生的涌浪。國外的研究多集中在概化性三維滑塊和典型滑坡涌浪實例上進行研究，尚無研究者對崩塌落石涌浪全過程進行流固耦合分析。

長江三峽庫區(qū)風(fēng)景秀麗、山高坡陡，是著名的旅游景點，也是崩塌落石的多發(fā)區(qū)。為給三峽庫區(qū)巫峽剪刀峰一帶危巖防治風(fēng)險管理提供科學(xué)依據(jù)，筆者以三峽庫區(qū)巫峽剪刀峰為原型，嘗試在計算流體力學(xué)軟件（computing fluid dynamic，CFD）中采用6自由度的運動碰撞模型（general moving objects collision model，GMO）耦合計算崩塌落石形成的涌浪效應(yīng)，為類似崩塌落石涌浪研究提供借鑒。

1 剪刀峰斜坡概況

三峽庫區(qū)巫峽剪刀峰位于重慶市巫山縣長江北岸，巫山新縣城下游16km。2006年至今的野外調(diào)查、巡查發(fā)現(xiàn)，剪刀峰一帶區(qū)域崩塌落石現(xiàn)象時有發(fā)生。該區(qū)域斜坡位于神女溪-官渡口向斜北翼，由三疊系嘉陵江組三段（T1j3）灰?guī)r構(gòu)成同傾順向坡，平均坡度超過50°。斜坡巖體剖面X節(jié)理發(fā)育，坡體被X節(jié)理和層面裂隙切割成菱形結(jié)構(gòu)體。塊體沿X節(jié)理和層面發(fā)生破壞，脫離母巖后向長江或兩側(cè)沖溝方向崩落。塊體的體積較小，但是落石高差往往超過400m，勢能巨大，對過往船只仍然構(gòu)成較大威脅。由于沿X節(jié)理破壞，該區(qū)域形成了三角面和V型溝的典型山體地貌。

剪刀峰一帶分別于2007年7月3日和7月23日發(fā)生過兩次崩塌，兩次崩塌的危巖源自同一地點。崩塌體呈板狀，東西長20m，厚0.35m，高80m，體積560m3（圖1）。7月23日第二次崩塌落石入江時，其下游50m處有一艘載有50余名乘客的客船正貼近懸崖往上游方向行使，所幸落石沒有擊中船只[23]。

圖1 2007年剪刀峰發(fā)生崩塌后照片F(xiàn)ig.1 Jiandao Peak photo after rockfalling in 2007

該區(qū)域長約2.5km，坡高險峻，因為是著名的旅游景點，不適合進行大規(guī)模危巖清除或攔網(wǎng)防護等防治工作。該斜坡段無人居住，其主要危害對象為航道過往船只，危害的主要形式為落石打擊和涌浪；因此有必要對落石及其產(chǎn)生的涌浪進行研究，為該區(qū)域的地質(zhì)災(zāi)害防治及風(fēng)險管理提供依據(jù)。

2 崩塌涌浪的流固耦合運動分析

庫岸崩塌落石產(chǎn)生涌浪是一個很復(fù)雜的過程。首先，石塊在陸地斜坡上進行滑動、彈跳；然后落入水體，石塊在水體和水下斜坡中運動，與水體相互作用產(chǎn)生涌浪，出現(xiàn)涌浪傳播。為了解決這一復(fù)雜的流固耦合運動問題，筆者在國內(nèi)首次引入Flow-3D軟件進行建模、分析。

2.1 Flow-3D介紹及耦合模型建立

Flow-3D是一款由Flow Sciences公司開發(fā)的通用流體動力學(xué)計算軟件，始于1980年的Los Alamos National Laboratory。在物質(zhì)守恒、動量守恒、能量守恒等歐拉方程框架內(nèi)，F(xiàn)low-3D采用了有限體積差分法逼近離散化計算域進行求解。該軟件有大量的模型用于模擬相變、非牛頓流體、孔隙介質(zhì)流、表面張力效應(yīng)、兩相流等。Flow-3D采用FAVOR （fractional area／volume obstacle representation）和 VOF（volume-of-fluid）技術(shù)來求解三維瞬時Navier-Stokes方程，能夠提供極為真實且詳盡的自由液面（free surface）流場信息。FAVOR和VOF技術(shù)使得在歐拉網(wǎng)格內(nèi)能夠定義固體邊界，能夠在計算流體響應(yīng)固體邊界時追蹤流體邊界。采用這一方法，固體物質(zhì)獨立生成網(wǎng)格，能夠高效率且精確地定義幾何外型。Flow-3D的FAVOR和VOF技術(shù)，使得它在描述自由液面流動方面具有獨特的準(zhǔn)確性和真實性[24]。

Flow-3D有許多不同的湍流模型用來模擬湍流，包括普朗特混合長度模型（Prandtl mixing length model）、k-ε 方 程 （k-ε model）、RNG 方 程（RNG scheme）和 LES模型（large eddy simulation model）。同時，在Flow-3D中有一個特殊的GMO碰撞計算模型[24]，能夠提供使用者預(yù)測移動對象在流體內(nèi)運動的狀況。GMO模擬剛體運動，可以是指定運動方式或與流動耦合計算。指定運動方式時流動受物體運動影響，而物體運動不受流體影響。與流動耦合時物體運動和流動是動態(tài)耦合的（兩者互相影響）。兩種方式中運動物體都可以有6個自由度，計算時可以有多種類型的運動物體，且可以相互碰撞。碰撞分析可以選擇采用彈性碰撞、部分塑性碰撞、完全塑性碰撞3種。彈性碰撞是指運動過程中運動物體間碰撞沒有能量損失。完全塑性碰撞是指運動物體間碰撞后，能量完全損失掉，這一碰撞分析采用總體摩擦系數(shù)和總體碰撞恢復(fù)系數(shù)來控制。碰撞恢復(fù)系數(shù)處于0和1之間，0代表完全塑性 ，1代表完全彈性。筆者采用了k-ε方程和GMO模型的耦合模型進行分析計算。

圖2 耦合模型Fig.2 Coupling model map

以剪刀峰斜坡河谷為例（圖2）建立了一個X方向1.6km長，Y方向1.5km寬，Z方向1km高的斜坡模型。根據(jù)調(diào)查分析，剪刀峰發(fā)生1萬m3的崩塌落石概率較小。而涌浪與塊體入水體積有密切的正相關(guān)性，假定剪刀峰發(fā)生1萬m3落石入水，以計算可能的最大涌浪高度和風(fēng)險。因此，建立了一個長80m，寬30m，厚5m的板狀GMO巖體，其密度為2600kg／m3。這一假設(shè)的塊體體積為1.2萬m3，初始重心高度652m，初始狀態(tài)為靜止。根據(jù)章廣成等[14]對灰?guī)r區(qū)落石碰撞恢復(fù)系數(shù)的研究，計算采用的碰撞恢復(fù)系數(shù)為0.72，總體摩擦系數(shù)為0.57。采用5m的計算網(wǎng)格進行離散，共有1827萬個網(wǎng)格單元。k-ε湍流模型邊界條件為：X方向（河流動方向）為靜水壓力邊界（定水頭邊界），Y方向為流出邊界，Zmax方向為自由表面邊界（水壓力為0，空氣界面），Zmin方向為不透水 Wall邊界（隔水邊界）。k-ε湍流模型初始條件為靜止水面，水位為175m。GMO模型的邊界條件為：Z方向為重力加速度，三維斜坡體為封閉邊界（不可進入的區(qū)域），其他區(qū)域為開區(qū)域（可以進入的區(qū)域）。GMO模型的初始條件為靜止?fàn)顟B(tài)。計算模擬板狀巖體從剪刀峰處崩塌-滾動-彈跳-入水-涌浪的全過程，模擬時間設(shè)置為120s。

2.2 流固耦合運動結(jié)果分析

該耦合計算模型在LENOVO THINK工作站上計算約9h，形成了152Gb的結(jié)果數(shù)據(jù)文件。GMO運動過程模擬顯示（圖3），塊體在陸地斜坡上經(jīng)歷了滑動、翻轉(zhuǎn)、彈跳等運動。巖塊從靜止?fàn)顩r開始下落后，有3次較大的彈跳發(fā)生在陸地斜坡上，有1次彈跳是入水之后發(fā)生的。由于陸地斜坡陡峻，在陸地斜坡上巖塊只有一小段進行了滑動或滾動，大部分是碰撞后彈跳飛過的。其運動軌跡顯現(xiàn)三維性，并不是在一個平面上運動。由于未考慮空氣阻力問題，運動過程中，巖塊的能量損失主要來源于碰撞和滾動時的摩擦。

雖然每次碰撞彈跳都會造成總能量損失和速度的暫時降低（圖4），但總體來看，隨著巖塊的下降，勢能轉(zhuǎn)化為動能，動能不斷增加，速度不斷增加。20 s左右速度達到最大，約80m／s，能量也達到最大，約9.30×1011J。20s后巖塊入水，巖塊入水是該過程的一個轉(zhuǎn)折點。入水后，GMO的總能量開始持續(xù)下降（圖4a），能量一部分傳遞給水體，一部分碰撞水下斜坡?lián)p失了。

巖塊進入水體后，與水體相互作用，兩者的運動相互影響。對巖塊而言，入水后其運動方向與軌跡發(fā)生較大的變化，由開始南南東的斜向下運動方向為主，慢慢轉(zhuǎn)化為東西方向擺動向下運動。瞬時過程圖（圖4b）顯示為：在波浪和水的作用下，巖塊在水中發(fā)生左右飄動下沉。運動速度圖上顯示為XY方向速度出現(xiàn)方向性擺動。巖塊在陸地477m高差的斜坡上運動了20s，在水下140m高差斜坡運動了15～20s。這是由于流固耦合作用，流體讓巖塊不易停止沉底，巖塊的水下運動時間明顯延長。運動物體完全停止后，流體獲得的總動能為6.08×1010J，與入水前巖塊的總能量相比較，能量傳遞率約為6.54%。

流體和運動物體的相互作用，極大地消耗了運動物體的能量，改變了運動物體的原有運動韻律，延緩了運動固體的下沉?xí)r間，延長了流-固能量交換過程。這一流固耦合過程，由于能量交換時間過長而不利于涌浪形成。

巖塊入水后，水體沿石塊四周壁面向上涌起飛濺，激起了較大的水花（水舌），最大高度達到了15.60m，而后散落入水中。以入水點為圓點，開始形成一個新月半圓錐狀的孤立涌浪波，在傳播過程中形成最大涌浪波約10.0m，距離岸邊約150.0 m。在涌浪的開始階段，涌浪波形成一圈一圈的環(huán)狀水波進行傳播（圖5a），衰減非常快。波的運動方向也為環(huán)狀，直至環(huán)狀波傳播至對岸，波的運動方向才開始轉(zhuǎn)為放射狀沿河道上下游進行傳播，同時反射波開始與涌浪波疊加（圖5b）。波浪傳播至河道中線（弘深線）時，波高衰減為2.40m。傳播至對岸時，波高已衰減至0.89m左右，最大爬高為1.90 m。波浪沿河道傳播500.0m后的最大浪高為1.50 m左右，1km外降低至1.0m以下。

借鑒國家海洋局2009年11月發(fā)布的《風(fēng)暴潮、海浪、海嘯和海冰災(zāi)害應(yīng)急預(yù)案》對計算區(qū)域航道進行涌浪風(fēng)險預(yù)警分區(qū)：當(dāng)波浪大于3m時為航道紅色預(yù)警區(qū)；當(dāng)波浪為2～3m時為航道橙色預(yù)警區(qū)；當(dāng)波浪為1～2m時為航道黃色預(yù)警區(qū)；當(dāng)波浪小于1m時為藍色預(yù)警區(qū)。根據(jù)這一標(biāo)準(zhǔn)，在入水區(qū)附近（離岸200m內(nèi)）為航道紅色預(yù)警區(qū)，在離剪刀峰岸線200～350m時為航道橙色預(yù)警區(qū)，在離剪刀峰岸線350m后的航道為黃色預(yù)警區(qū)。因此，河道中心線以南（上行航道）的涌浪較小，風(fēng)險較低。

另外，筆者采用的崩塌落石為1.2萬m3，實際崩塌如此大方量的可能性較低，產(chǎn)生的涌浪高度也將比本文所得結(jié)果低，預(yù)警等級相應(yīng)會降低。因此，如果將目前全河面航道的北航道往南移350m左右，將紅色-橙色預(yù)警區(qū)河面作為避讓區(qū)，則航道內(nèi)的涌浪預(yù)警將從紅色預(yù)警區(qū)降為黃色或藍色預(yù)警區(qū)，涌浪風(fēng)險大大降低，航道安全度將得到提高。

3 討論

通過對崩塌落石涌浪案例的流固耦合分析研究，發(fā)現(xiàn)有以下問題值得探討和思考：

1）就目前PC機硬件而言，小的GMO體意味著更小的離散網(wǎng)格和更多的網(wǎng)格單元。上千萬或上億的網(wǎng)格單元常常會導(dǎo)致一些PC機硬件無法滿足計算要求；即便是可以滿足，其所需的計算時間也是很漫長的。這也是筆者沒有采用更小尺寸落石塊體的重要原因之一。同樣的原因，也導(dǎo)致在網(wǎng)格的大小與結(jié)果的質(zhì)量之間需要取舍平衡。

2）對斜坡為多物質(zhì)組成時，應(yīng)分別建立多組件的斜坡模型，以設(shè)立不同的碰撞恢復(fù)系數(shù)。崩塌落石是一個很復(fù)雜的過程，斜坡與滾石特性和初始狀態(tài)的稍微改變都會強烈改變落石的動力特征。因此對未知滾石運動過程最好的評價方法可能是先采用概率的辦法進行軌跡預(yù)測，然后對大概率事件或高風(fēng)險事件進行確定性耦合分析。對已知滾石停止點的落石碰撞過程模擬，需要通過工程地質(zhì)判斷調(diào)整輸入?yún)?shù)，以達到預(yù)期的結(jié)果[11]。

圖3 巖塊運動瞬時圖片F(xiàn)ig.3 Moving instantaneous picture of rock block

圖4 巖塊運動過程線Fig.4 Process line of rock block

圖5 涌浪傳播及波速矢量圖Fig.5 Impulsive wave propagation and wave velocity vector map

3）由于概率法能更真實準(zhǔn)確地預(yù)測滾石的運動特性，在硬件容許和算法優(yōu)化的條件下可考慮將概率方法引入耦合數(shù)值模型中。

4）筆者選用了工程中應(yīng)用最廣泛的k-ε湍流模型，其基本思想是用低階關(guān)聯(lián)量和平均流體性質(zhì)來模擬未知的高階關(guān)聯(lián)量。今后可考慮應(yīng)用LES大渦模型對實際涌浪中大量存在的高Re數(shù)湍流激流現(xiàn)象進行有效求解。

5）一般來說，一個從上方滾下的大塊石經(jīng)過彈跳、翻滾、滑動可能分解成許多小塊石。碰撞分解的過程在本文中并未涉及。塊石彈跳、碰撞、分解然后入水產(chǎn)生涌浪，全過程數(shù)值模擬涉及到斷裂力學(xué)、碰撞固體力學(xué)、運動學(xué)和流體力學(xué)，是目前崩滑體涌浪研究領(lǐng)域研究的最前沿。

4 結(jié)論及建議

筆者首次嘗試使用耦合GMO碰撞模型的CFD軟件預(yù)測了崩塌落石產(chǎn)生涌浪的全過程，得到了以下結(jié)論和建議。

1）三維流體力學(xué)與固體力學(xué)耦合分析涌浪可獲得從微觀到宏觀的直觀結(jié)果，這一方向?qū)⑹怯坷搜芯康陌l(fā)展趨勢。

2）以剪刀峰河道為原型，建立了一個k-ε湍流模型和GMO碰撞模型的流固耦合模型，網(wǎng)格單元為1827萬個，假定的GMO方量為1.2萬m3。

3）通過流固耦合計算得到，巖塊經(jīng)過3次彈跳后入水，入水時運動物體的速度達到最大，約為80 m／s，動能也達到最大，約9.3×1011J。入水后，巖塊動能持續(xù)降低，能量開始傳遞給水體。在水體作用下，巖塊的運動方向與軌跡發(fā)生較大變化，下沉?xí)r間明顯延長。運動物體完全停止后，流體獲得的總動能為6.08×1010J，能量傳遞率約為6.54%。

4）根據(jù)航道內(nèi)涌浪大小分析，在入水區(qū)附近（離岸200m內(nèi)）航道為紅色預(yù)警區(qū)，在離剪刀峰岸線200～350m內(nèi)航道為橙色預(yù)警區(qū)，離剪刀峰岸線350m后的航道為黃色預(yù)警區(qū)。建議將本區(qū)段北航道線往南移350m左右，航道內(nèi)涌浪風(fēng)險大大降低，航道安全度得到提高。

5）流固耦合分析涌浪是一個新方法、新手段，需要與其他方法進行對比分析，方法自身也尚還有很多方面需要完善，有待深入研究。

（References）：

[1]Valentin Heller.Landslide Generated Impulse Waves：Prediction of Near Field Characteristics[D].Zürich：Eidgenossische Technische Hochschule Zürich，2007.

[2]Bozzolo D，Pamini B.Simulation of Rock-Falls Down a Valley Side[J].Acta Mechanica，1986，63（1）：113-130.

[3]Dourrier F，Dorren L，Nicot F，et al.Toward Objective Rockfall Trajectory Simulation Using a Stochastic Impact Model[J].Geomorphology，2009，110（3／4）：68-79.

[4]Giani G，Giacomini A，Migliazza，et al.Experimental and Theoretical Studies to Improve Rock Fall Analysis and Protection Work Design[J].Rock Mechanics and Rock Engineering，2004，37（5）：369-389.

[5]Chau K T，Wond R H C，Lee C F.Rockfall Problems in Hong Kong and Some New Experimental Results for Coefficient of Restitution[J].Int J Rock Mech Min Sci，1996，35：662-663.

[6]秦志英，陸啟韶.基于恢復(fù)系數(shù)的碰撞過程模型分析[J].動力學(xué)與控制學(xué)報，2006，4（4）：294-298.Qin Zhiying，Lu Qishao.Analysis of Impact Process Model Based on Restitution Coefficient[J].Journal of Dynamics and Control，2006，4（4）：294-298.

[7]何思明，吳永，李新坡.滾石沖擊碰撞恢復(fù)系數(shù)研究[J].巖土力學(xué)，2009，30（3）：623-627.He Siming，Wu Yong，Li Xinpo.Research on Restitution Coefficient of Rock Fall[J].Rock and Soil Mechanics，2009，30（3）：623-627.

[8]沈均，何思明，吳永.滾石對墊層材料的沖擊特性研究[J].安徽農(nóng)業(yè)科學(xué)，2009，37（17）：8286-8288.Shen Jun，He Siming，Wu Yong.Study on the Impact Properties of Rock-Fall on Cushion Material[J].Journal of Anhui Agricultural Sciences，2009，37（17）：8286-8288.

[9]楊海清，周小平.邊坡落石運動軌跡計算新方法[J].巖土力學(xué)，2009，30（11）：3411-3416.Yang Haiqing，Zhou Xiaoping.A New Approach to Calculate Trajectory of Rockfall[J].Rock and Soil Mechanics，2009，30（11）：3411-3416.

[10]Mangwandi C，Cheong Y，Adams M，et al.The Coefficient of Restitution of Different Representative Types of Granules [J].Chemical Engineering Science，2007，62（1／2）：437-450.

[11]Rocscience.Rocfall User Manual：Statistical Analysis of Rockfalls[M].[S.l.]：Rocscience Inc，2002：51-59.

[12]Azzoni A，F(xiàn)reitas M.Experimentally Gained Parameters，Decisive for Rock Fall Analysis[J].Rock Mechanics and Rock Engineering，1995，28（2）：111-124.

[13]Day R W.Case Studies of Rockfall in Soft Versus Hard Rock[J].Environmental and Engineering Geoscience，1997，3（1）：133-140.

[14]章廣成，向欣，唐輝明.落石碰撞恢復(fù)系數(shù)的現(xiàn)場試驗與數(shù)值計算[J].巖石力學(xué)與工程學(xué)報，2011，30（6）：1266-1273.Zhang Guangcheng，Xiang Xin，Tang Huiming.Field Test and Numerical Calculation of Restitution Coefficient of Rockfall Collision[J].Chinese Journal of Rock Mechanics and Engineering，2011，30（6）：1266-1273.

[15]Jones C L，Higgins J D，Andrew R D.Colorado Rockfall Simulation Program，Version 4.0（for Windows）[R].Colorado：Colorado Department of Transportation，Colorado school of Mines，Colorado Geological Survey，2000.

[16]Guzzetti F，Crosta G，Detti R.STONE：A Computer Program for the Three-Dimensional Simulation of Rock-Falls[J].Computers and Geosciences，2002，28（9）：1079-1093.

[17]Azzoni A，Barbera G，Zaninetti A.Analysis and Prediction of Rockfalls Using a Mathematical Model[J].International Journal of Rock Mechanics and MiningSciences and Geomechanics Abstracts，1996，33（4）：178.

[18]Choi B H，Pelinovsky E，Kim D C，et al.Two-and Three-Dimensional Computation of Solitary Wave Runup on Non-Plane Beach[J].Nonlin Processes Geophys，2008，15：489-502.

[19]Silvia B，Marco P.Shallow Water Numerical Model of the Wave Generated by the Vajont Landslide[J].Environmental Modelling ＆ Software，2011（26）：406-418.

[20]Montagna F，Bellotti G，Risio M D.3DNumerical Modeling of Landslide-Generated Tsunamis Around a Conical Island[J].Nat Hazards，2011（58）：591-608.

[21]Pastor M，Herreros I，F(xiàn)ernández M J A，et al.Modelling of Fast Catastrophic Landslides and Impulse Waves Induced by Them in Fjords，Lakes and Reservoirs[J].Engineering Geology，2009，109：124-134.

[22]Stéphane Abadie，Denis Morichon，Stéphan Grilli，et al.Numerical Simulation of Waves Generated by Landslides Using a Multiple-Fluid Navier-Stokes Model[J].Coastal Engineering，2010（57）：779-794.

[23]Huang B L，Chen L D，Peng X M，et，al.Assessment of the Risk of Rockfalls in Wu Gorge，Three Gorges，China[J].Landslides，2010，7（1）：1-11.

[24]Flow Science Inc.Flow-3DUser Manual[M].New Mexico：Flow Science Inc，2012.