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

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

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

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

中圖分類號(hào)：TV39 文獻(xiàn)標(biāo)識(shí)碼：A?? 文章編號(hào)：1006—7973（2022）03-0076-04

邊坡失穩(wěn)是一種廣泛存在的地質(zhì)災(zāi)害，產(chǎn)生的滑坡災(zāi)害危害性和破壞性巨大。中國是一個(gè)滑坡災(zāi)害極為頻繁的國家，尤其是在我國的西南部地區(qū)，存在著大量的成因復(fù)雜、規(guī)模大和危害性強(qiáng)等顯著特點(diǎn)的大型滑坡?；掠坷俗鳛榛麦w入水后形成的次生災(zāi)害，大型滑坡體高速入水激起巨大涌浪，而且在河道上下游長距離內(nèi)傳播，容易造成洪水、漫壩等更大的災(zāi)難，會(huì)造成人員傷亡和財(cái)產(chǎn)損失。

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

從以上國內(nèi)外諸多的滑坡災(zāi)害例子中可以看出，大規(guī)?；麦w滑入河道，激起的涌浪對(duì)航行的船只產(chǎn)生嚴(yán)重威脅，同時(shí)堆積河床底部的滑坡體，會(huì)阻塞河道，形成堰塞湖;發(fā)生在庫區(qū)的大規(guī)?；聲?huì)引發(fā)巨大涌浪，摧毀水工建筑物，甚至?xí)a(chǎn)生漫壩危及下游人民生命財(cái)產(chǎn)安全。為減少人員傷亡和財(cái)產(chǎn)損失，對(duì)潛在滑坡進(jìn)行實(shí)時(shí)監(jiān)控和預(yù)警，有必要開展滑坡涌浪以及高速遠(yuǎn)程滑坡涌浪研究，預(yù)測和評(píng)估災(zāi)害范圍，為滑坡涌浪災(zāi)害產(chǎn)生的實(shí)時(shí)預(yù)報(bào)和災(zāi)后緊要防范提供科學(xué)依據(jù)。

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

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

1控制方程

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

1.1連續(xù)相控制方程

根據(jù)質(zhì)量守恒定律推導(dǎo)出流動(dòng)控制方程，即連續(xù)方程：

（1）

式中：為拉普拉斯算子。

根據(jù)牛頓第二定律推導(dǎo)出流體動(dòng)量方程，控制方程中的流體速度、壓力和密度，為每個(gè)體積單元內(nèi)的平均狀態(tài)變量。動(dòng)量方程表達(dá)式：

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

體內(nèi)流體與顆粒之間的相互作用力。

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

1.2顆粒離散元控制方程

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

流體和顆粒相互作用力。顆粒的運(yùn)動(dòng)控制方程如下：

（3）

（4）

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

1.3 EDEM-FLUENT 耦合方法

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

（5）

（6）

式中：為顆粒直徑;為曳力系數(shù)。

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

2算例驗(yàn)證

2.1模型設(shè)置

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

2.2邊界條件

模型上方上部為壓力出口（pressout）邊界條件，模型底部和兩側(cè)壁以及滑體四周設(shè)置為固壁邊界條件，固壁邊界為無滑移條件，近壁區(qū)處理采用標(biāo)準(zhǔn)壁面函數(shù)法。

2.3試驗(yàn)結(jié)果分析

規(guī)定 t=0s 時(shí)滑坡體開始運(yùn)動(dòng)且前緣到達(dá)自由水面，

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

圖3所示為涌浪產(chǎn)生區(qū)域流體相體積分?jǐn)?shù)云圖、滑坡體的速度云圖以及周圍水體的速度矢量圖。從流體相體積分?jǐn)?shù)云圖可知，滑坡體下滑時(shí)周圍的水體被排開，導(dǎo)致水面抬升形成涌浪，并且隨著時(shí)間推移近場形成的涌浪將逐漸向外傳播。從滑坡體速度云圖可知，滑坡體沿斜面下滑時(shí)周圍水體獲得較大的速度，使得在滑坡體周圍產(chǎn)生順時(shí)針渦流，這是因?yàn)樵谙嗷プ饔玫倪^程中能量從滑坡體傳遞至了庫水中，引起周圍水體的速度變化;隨著滑坡體停止運(yùn)動(dòng)后，周圍水體的速度開始逐漸變小，表明滑坡體運(yùn)動(dòng)對(duì)其周圍水體的速度有很大影響。

3結(jié)論

（1）通過 DPM 模型可直觀地看到滑坡體的變形和顆粒的運(yùn)動(dòng)軌跡，對(duì)比剛性滑坡體在下滑的過程中不發(fā)生變形，本文數(shù)值試驗(yàn)中滑坡體在下滑過程中運(yùn)動(dòng)形態(tài)與原物理試驗(yàn)保持一致。

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

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

參考文獻(xiàn)：

[1]黃潤秋.20世紀(jì)以來中國的大型滑坡及其發(fā)生機(jī)制[J].巖石力學(xué)與工程學(xué)報(bào) ，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].中國地質(zhì)災(zāi)害與防治學(xué)報(bào)，1993，5（2）：77.84.

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

[5]蔣權(quán) ， 陳希良 ， 肖江劍 ， 邱向東 ， 李攀峰.云南黃坪庫區(qū)滑坡運(yùn)動(dòng)及其失穩(wěn)模式的離散元模擬[J].中國地質(zhì)災(zāi)害與防治學(xué)報(bào) ，2018，29（03）：53-59.

[6]宋新遠(yuǎn).大型滑坡災(zāi)害數(shù)值模擬研究[D].上海交通大學(xué) ，2009.

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

[8]徐文杰.滑坡涌浪影響因素研究[J].工程地質(zhì)學(xué)報(bào) ，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]譚海.基于無網(wǎng)格方法的滑坡及滑坡涌浪研究[D].武漢大學(xué) ，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：