崩塌落石產生涌浪的流固耦合運動分析

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

1.武漢地質調查中心，武漢 430205

2.中國地質大學工程學院，武漢 430074

3.中國地質環境監測院，北京 100081

0 引言

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

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

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

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

1 剪刀峰斜坡概況

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

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

圖1 2007年剪刀峰發生崩塌后照片Fig.1 Jiandao Peak photo after rockfalling in 2007

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

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

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

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

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

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

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

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

2.2 流固耦合運動結果分析

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

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

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

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

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

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

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

3 討論

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

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

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

圖3 巖塊運動瞬時圖片Fig.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）由于概率法能更真實準確地預測滾石的運動特性，在硬件容許和算法優化的條件下可考慮將概率方法引入耦合數值模型中。

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

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

4 結論及建議

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

1）三維流體力學與固體力學耦合分析涌浪可獲得從微觀到宏觀的直觀結果，這一方向將是涌浪研究的發展趨勢。

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

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

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

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]秦志英，陸啟韶.基于恢復系數的碰撞過程模型分析[J].動力學與控制學報，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]何思明，吳永，李新坡.滾石沖擊碰撞恢復系數研究[J].巖土力學，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].安徽農業科學，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].巖土力學，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，Freitas 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]章廣成，向欣，唐輝明.落石碰撞恢復系數的現場試驗與數值計算[J].巖石力學與工程學報，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，Ferná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.