基于 DPM 模型的水上散粒體滑坡涌浪仿真分析

魏云鵬 閆毅志 張萬舉 李梓萌 王文雄 劉波

摘要：針對工程地質中散粒體結構的滑坡，本文以 EDEM-FLUENT 耦合方法和 DPM（Discrete Phase Modle）模型為基礎，采用 VOF 方法追蹤自由液面，模擬散粒體滑坡體沿斜坡運動所引起的涌浪產生及其傳播過程。將數值計算結果與試驗觀測數據進行對比，驗證該數值模型的有效性，并分析了 DPM 模型作用下滑坡體運動狀態對涌浪生成和傳播過程的影響。通過對數值仿真試驗結果進行分析，發現數值模擬結果與物理試驗觀測數據吻合較好，這表明 DPM 模型能夠準確模擬滑坡涌浪的產生與傳播特性。

關鍵詞：散粒體滑坡;滑坡涌浪;EDEM-FLUENT 耦合方法;DPM 模型

中圖分類號：TV39 文獻標識碼：A?? 文章編號：1006—7973（2022）03-0076-04

邊坡失穩是一種廣泛存在的地質災害，產生的滑坡災害危害性和破壞性巨大。中國是一個滑坡災害極為頻繁的國家，尤其是在我國的西南部地區，存在著大量的成因復雜、規模大和危害性強等顯著特點的大型滑坡。滑坡涌浪作為滑坡體入水后形成的次生災害，大型滑坡體高速入水激起巨大涌浪，而且在河道上下游長距離內傳播，容易造成洪水、漫壩等更大的災難，會造成人員傷亡和財產損失。

最典型的滑坡涌浪災害屬1958年07月09日由8.3級地震誘發的美國阿拉斯加 Lituya 海灣滑坡，體積約為3×107 m3的滑坡體以高速沖入水中，產生524 m 高的涌浪，引發了有記錄以來的最高涌浪，巨浪沖高對沿途植被造成了嚴重的破壞和侵蝕。1963年發生在意大利的瓦依昂水庫滑坡，高速入水的滑坡體產生的涌浪高達到175m，造成巨大的生命財產損失。我國地質災害頻發，其中滑坡災害是所有地質災害中分布最廣、危害最嚴重的一種。1985年發生在長江西陵峽新灘北岸的滑坡，大規模的滑坡體進入江中，掀起高達數十米巨大涌浪，造成大量船只被擊毀，死亡9人。2013年云南永善縣黃華鎮黃坪村發生滑坡，滑坡體積達到數十萬立方米，造成嚴重的生命和財產損失。

從以上國內外諸多的滑坡災害例子中可以看出，大規模滑坡體滑入河道，激起的涌浪對航行的船只產生嚴重威脅，同時堆積河床底部的滑坡體，會阻塞河道，形成堰塞湖;發生在庫區的大規模滑坡會引發巨大涌浪，摧毀水工建筑物，甚至會產生漫壩危及下游人民生命財產安全。為減少人員傷亡和財產損失，對潛在滑坡進行實時監控和預警，有必要開展滑坡涌浪以及高速遠程滑坡涌浪研究，預測和評估災害范圍，為滑坡涌浪災害產生的實時預報和災后緊要防范提供科學依據。

隨著計算機技術的迅猛發展，數值模擬方法成為研究滑坡涌浪問題的重要手段，該方法具有全面、準確分析問題的優點，以及處理滑坡涌浪研究中遇到的復雜耦合問題的優勢。目前針對滑坡涌浪各個過程的數值模擬研究已經取得了較大的進展。現有的滑坡涌浪模型多將滑體簡化為剛性塊體從滑面滑入水中，宋新遠基于 Navier-Stokes 方程，結合 RNG k-ε湍流模型并采用 VOF 法跟蹤非線性自由表面流場模擬了涌浪的產生過程，模擬結果與試驗結果和實例觀測結果吻合較好;但是其假設滑坡體為固定形狀的單一塊體，忽略了滑坡體在滑動過程中的變形對涌浪產生的影響。徐文杰以室內模型試驗為基礎采用 CEL 算法建立數值分析模型，并與模型試驗結果和 SPH 方法模擬結果進行了對比分析，之后進一步研究了滑坡體形狀、體積、滑面摩擦角等與涌浪的關系，其研究也主要是針對單個塊體滑坡，并沒有考慮散體滑坡。而自然界中的滑坡通常變形顯著，如果僅僅將滑坡體假定為剛性塊體過于簡單，對涌浪的產生也有很大影響，因此有必要在數值模型中考慮滑坡體的變形效應。

為了減小將滑體簡化為剛性塊體帶來的誤差，更好地描述滑坡體運動規律及其產生涌浪運動特征。本文針對水上散粒體滑坡，采用 EDEM-FLUENT 耦合方法，結合 DPM 模型對散粒體滑坡產生的滑坡涌浪特征進行研究。主要考慮離散體和水體之間相互作用的情況下，通過物理實驗驗證耦合算法的有效性，并對涌浪的產生和傳播特性進行了數值研究。

1控制方程

本文數值模擬計算中，采用基于歐拉-拉格朗日框架下的 DPM 方法對水上散粒體滑坡涌浪進行研究。在 DPM 模型中，水和空氣為連續相，采用 Navier-Stokes 方程描述，并采用 VOF 方法對自由液面進行追蹤。曳力是固體顆粒在流場中運動時受到的主要作用力，本文選用 Morsi 等提出的 Spherical 曳力模型。

1.1連續相控制方程

根據質量守恒定律推導出流動控制方程，即連續方程：

（1）

式中：為拉普拉斯算子。

根據牛頓第二定律推導出流體動量方程，控制方程中的流體速度、壓力和密度，為每個體積單元內的平均狀態變量。動量方程表達式：

（2）式中：為流體單元的平均速度; 為單元控制

體內流體與顆粒之間的相互作用力。

由于滑坡體下滑引起河道、水庫等水的運動是湍流運動，本文選取雷諾平均 NS 模型（RANS）方法中的 RNG k-ε模型。

1.2顆粒離散元控制方程

本文在顆粒接觸力和體積力的基礎上，還引入了

流體和顆粒相互作用力。顆粒的運動控制方程如下：

（3）

（4）

式中：為顆粒的平移線速度; 為作用在顆粒上的接觸力; 為顆粒所受的流體作用力;為粒受到顆粒（或墻體）的轉動力矩; 為顆粒總的接觸數。對于顆粒間接觸作用，本文采用 Hertz-Mindlin 非線性接觸模型。

1.3 EDEM-FLUENT 耦合方法

流體相和顆粒相間的相互作用力包括拖曳力、壓差力、浮力和虛質量力等。根據研究需要，本文僅考慮了浮力和拖曳力：

（5）

（6）

式中：為顆粒直徑;為曳力系數。

將計算離散相的 EDEM 軟件和計算流體相的 FLUENT 軟件耦合起來。在求解耦合系統時，在一些初始條件和增量荷載下， EDEM 軟件首先求解粒子系統，然后計算合力。然后將相互作用力傳遞到流體系統中，通過 FLUENT 求解流體系統。當流體的所有狀態變量都得到求解后，這些相互作用力就會更新并傳回到 EDEM 中，以求解下一時間步的粒子系統。

2算例驗證

2.1模型設置

本文主要對水上滑坡涌浪進行驗證，采用Viroulet 等的水上散粒體滑坡涌浪物理試驗來進行數值仿真研究。初始模型如圖1所示，圖中水深為0.15 m，滑坡斜面與水平面夾角為45°，實驗中通過布置的4個波高儀來記錄產生的涌浪。滑坡體的截面為圖中0.14 m×0.14 m 的三角形，初始 t =0時刻滑坡體剛好位于自由水面之上，滑坡體的前緣速度為零，啟動后滑坡體從靜止開始沿斜面往下滑。物理試驗中初始時刻通過擋板維持滑坡體穩定，抽去擋板后滑坡體沿著傾角45°的斜面下滑產生涌浪（本文所采用的滑坡體顆粒參數源于譚海）。

2.2邊界條件

模型上方上部為壓力出口（pressout）邊界條件，模型底部和兩側壁以及滑體四周設置為固壁邊界條件，固壁邊界為無滑移條件，近壁區處理采用標準壁面函數法。

2.3試驗結果分析

規定 t=0s 時滑坡體開始運動且前緣到達自由水面，

圖2中（a）、（b）、（c）情況中上圖為 Viroulet 等的試驗結果，中間圖為 Si 等的試驗結果，下圖為本文DPM 模型數值仿真結果，散粒體滑坡在下滑的過程中伴隨著顯著變形，在水體的拖曳力作用下，滑坡體前緣逐漸變厚，呈球形鼓泡狀。并且滑坡體入水后以較大的速度沿斜坡運動到水槽底部，并繼續向前運動，最終靜止堆積在水槽底部。進一步對比滑體的水平運動距離可知，本文的仿真試驗結果中滑坡體向前滑移的距離與前兩者試驗結果基本吻合。

圖3所示為涌浪產生區域流體相體積分數云圖、滑坡體的速度云圖以及周圍水體的速度矢量圖。從流體相體積分數云圖可知，滑坡體下滑時周圍的水體被排開，導致水面抬升形成涌浪，并且隨著時間推移近場形成的涌浪將逐漸向外傳播。從滑坡體速度云圖可知，滑坡體沿斜面下滑時周圍水體獲得較大的速度，使得在滑坡體周圍產生順時針渦流，這是因為在相互作用的過程中能量從滑坡體傳遞至了庫水中，引起周圍水體的速度變化;隨著滑坡體停止運動后，周圍水體的速度開始逐漸變小，表明滑坡體運動對其周圍水體的速度有很大影響。

3結論

（1）通過 DPM 模型可直觀地看到滑坡體的變形和顆粒的運動軌跡，對比剛性滑坡體在下滑的過程中不發生變形，本文數值試驗中滑坡體在下滑過程中運動形態與原物理試驗保持一致。

（2）與剛體滑塊滑坡體運動到水槽底部轉角處停止運動不同，可變形散粒體滑坡在運動至轉角處后還會沿著水槽底部繼續向前運動，但由于受到水體拖曳力等阻力作用滑坡體開始減速變慢，最終靜止堆積在斜面和水槽底部的轉角處。

（3）可變散粒體形滑坡體下滑過程中由于內部顆粒的摩擦和碰撞，一部分能量耗散在滑坡體內部的變形中，導致只有較少的能量從可變形滑坡體傳遞至庫水中。因此，可變形滑坡涌浪自由水面的變形沒有相應的剛體滑坡涌浪的劇烈。

參考文獻：

[1]黃潤秋.20世紀以來中國的大型滑坡及其發生機制[J].巖石力學與工程學報 ，2007（03）：433-454.

[2] Fritz H M ，? Hager W H ，? Minor H E . Lituya Bay case： Rockslide impact and wave run-up[J]. ence of Tsunami Hazards， 2001， 19（1）：3-19.

[3]鐘立勛.意大利扎依昂水庫滑坡事件的啟示[J].中國地質災害與防治學報，1993，5（2）：77.84.

[4]劉世凱.長江西陵峽新灘滑坡涌浪高度衰減因素初探[J].水利水電技術 ，1987（09）：11-14.

[5]蔣權 ， 陳希良 ， 肖江劍 ， 邱向東 ， 李攀峰.云南黃坪庫區滑坡運動及其失穩模式的離散元模擬[J].中國地質災害與防治學報 ，2018，29（03）：53-59.

[6]宋新遠.大型滑坡災害數值模擬研究[D].上海交通大學 ，2009.

[7]徐文杰.基于 CEL 算法的滑坡涌浪研究[J].工程地質學報 ，2012，20（03）：350-354.

[8]徐文杰.滑坡涌浪影響因素研究[J].工程地質學報 ，2012，20（04）：491-507.

[9] Morsi S A ，? Alexander A J . An investigation of particle trajectories? in? two-phase? flow? systems[J]. Journal? of Fluid Mechanics， 1972， 55.

[10] H.P. Zhu，Z.Y . Zhou，R.Y . Y ang，A.B. Yu. Discrete particle simulation of particulate systems： Theoretical developments[J]. Chemical Engineering Science，2006，62（13）.

[11] S. Viroulet，A.Sauret，O. Kimmoun，C. Kharif. Granularcollapse into water： toward tsunami landslides[J]. Journal? of Visualization，2013，16（3）：

[12]譚海.基于無網格方法的滑坡及滑坡涌浪研究[D].武漢大學 ，2019.

[13] PengfeiSi，HuabinShi，Xiping Yu. A general numericalmodel for surface waves generated by granular material intruding into a water body[J]. Coastal Engineering，2018，142：