考慮阻力作用的碎石土滑坡涌浪數(shù)值模擬方法研究

薛宏程 彭杏瑤 馬倩 卿正陽

摘要：

滑坡涌浪是庫區(qū)水災(zāi)害的風(fēng)險(xiǎn)來源之一，開展涌浪預(yù)測研究對(duì)提升航道安全，減小近岸洪水風(fēng)險(xiǎn)具有重要的工程應(yīng)用價(jià)值，但滑坡體入水形成涌浪的動(dòng)力過程十分復(fù)雜，目前對(duì)其進(jìn)行準(zhǔn)確模擬仍是棘手的問題?；诘刃Я黧w假設(shè)，將非均質(zhì)的碎石土滑坡視為連續(xù)的特殊流體，通過自定義函數(shù)編寫滑坡體沿滑動(dòng)面運(yùn)動(dòng)的摩擦阻力程序，利用Realizable k-ε湍流模型與VOF法建立了碎石土滑坡涌浪數(shù)值模擬模型。研究結(jié)果表明：該方法可較好地模擬碎石土滑坡體水下運(yùn)動(dòng)與堆積形態(tài)，且獲得的涌浪波幅精度較高；阻力模型中摩擦系數(shù)與內(nèi)摩擦系數(shù)同滑坡體的變形速率分別成反比關(guān)系，而湍流阻力系數(shù)增大僅能減緩滑坡后緣的變形速率。研究成果可為定量評(píng)估碎石土滑坡涌浪災(zāi)害風(fēng)險(xiǎn)、制定庫區(qū)涌浪災(zāi)害緊急避險(xiǎn)預(yù)案提供參考。

關(guān) 鍵 詞：

滑坡涌浪； 數(shù)值模擬； 等效流體； k-ε湍流模型； VOF模型

中圖法分類號(hào)： TV139

文獻(xiàn)標(biāo)志碼： A

DOI：10.16232/j.cnki.1001-4179.2023.07.022

0 引 言

碎石土滑坡是山區(qū)水庫滑坡的主要類型之一，在庫區(qū)沿岸分布十分廣泛，主要由碎石與土的混合物組成，結(jié)構(gòu)松散，在強(qiáng)降雨作用后極易失穩(wěn)，一旦高速滑入水庫、湖泊或河道后，形成的涌浪會(huì)使災(zāi)害影響區(qū)域急劇擴(kuò)大，嚴(yán)重威脅航行船舶、水工建筑物及沿岸城鎮(zhèn)的安全。2003年7月13日，三峽庫區(qū)支流千將坪滑坡及其形成的涌浪造成13人死亡，經(jīng)濟(jì)損失高達(dá)8 000 萬元［1］；2008年11月23日，三峽庫區(qū)龔家方發(fā)生崩滑，對(duì)岸涌浪高度達(dá)13 m，涌浪傳播至4.2 km外的巫山新縣城時(shí)浪高仍超過1 m，使?？看a頭的多只船舶受損［2］；2015年6月24日，三峽庫區(qū)紅巖子滑坡形成的涌浪打翻船舶13艘，造成2人死亡，4人受傷［3］。為了對(duì)此類次生災(zāi)害進(jìn)行預(yù)測和早期防范，準(zhǔn)確獲得涌浪傳播特征與影響范圍是工程安全運(yùn)行中亟需解決的問題。

數(shù)值模擬方法具有可重復(fù)性強(qiáng)、數(shù)據(jù)監(jiān)測靈活等優(yōu)勢，已成為涌浪災(zāi)害預(yù)警的重要研究手段之一［4-5］。在模擬巖質(zhì)滑坡沖擊水體形成的涌浪時(shí)，常將滑坡體簡化為剛體，再通過經(jīng)驗(yàn)公式計(jì)算給出剛體下滑速度［6-7］。但是碎石土滑坡在運(yùn)動(dòng)過程中不斷變形，滑坡下滑速度無法準(zhǔn)確模擬，且滑坡入水后與水體發(fā)生強(qiáng)烈的摻混作用，極大增加了涌浪模擬難度。一些學(xué)者嘗試采用離散元法（Discrete Element Method，DEM）將土石滑坡劃分為若干個(gè)離散元單元陣來模擬其運(yùn)動(dòng)過程，該方法可以較真實(shí)地還原滑坡體中的非線性大變形特征。徐文杰［8］通過構(gòu)建SPH-DEM流固耦合算法，再現(xiàn)了意大利瓦依昂（Vajont）庫區(qū)滑坡涌浪災(zāi)害鏈的全過程；譚海等［9］引入浮力和拖曳力實(shí)現(xiàn)了DEM模塊與SPH模塊的耦合，發(fā)現(xiàn)此模型更適合模擬水上散粒體滑坡形成的涌浪特征。但是離散元法在提高精度的同時(shí)將面臨控制參數(shù)較多、粒子分散度不能太大、計(jì)算量巨大等難點(diǎn)。Hungr等［10-11］基于連續(xù)介質(zhì)理論，將碎石土滑坡視為一種特殊流體，認(rèn)為滑坡運(yùn)動(dòng)形態(tài)與其內(nèi)部黏性和底部摩擦阻力有關(guān)。在此基礎(chǔ)上，朱圻等［12］、袁小一等［13］、高楊等［14］、夏式偉等［15］基于上述方法分別對(duì)牛圈溝、文家溝、雞尾山、湯家溝發(fā)生的滑坡事件進(jìn)行了模擬驗(yàn)證，計(jì)算結(jié)果與實(shí)際情況吻合較好。雖然上述方法能夠較好地模擬碎石土滑坡的運(yùn)動(dòng)過程與堆積形態(tài)，但是滑坡涌浪問題的復(fù)雜性在于滑坡運(yùn)動(dòng)過程中除受氣相影響外，還要充分考慮滑坡與水體的相互作用。針對(duì)這類典型的三相流問題如何使計(jì)算快速收斂，且能夠準(zhǔn)確模擬涌浪的形態(tài)特征，還有待進(jìn)一步研究和探索。

本文基于Hungr的理論，通過自定義函數(shù)編寫阻力模型程序控制滑坡運(yùn)動(dòng)過程，結(jié)合Realizable k-ε湍流模型與流體體積法（Volume of Fluid，VOF），對(duì)滑坡、空氣與水三相共同作用下涌浪形成過程進(jìn)行模擬，同時(shí)還討論了滑坡阻力模型中摩擦系數(shù)、內(nèi)摩擦系數(shù)、湍流阻力系數(shù)3個(gè)參數(shù)對(duì)滑坡形態(tài)特征的影響。

2 模型驗(yàn)證與結(jié)果分析

2.1 模型構(gòu)建

驗(yàn)證算例源于Fritz在蘇黎世聯(lián)邦理工學(xué)院開展的散粒體滑坡涌浪模型試驗(yàn)［16］。該試驗(yàn)在長11 m、深1 m、寬0.5 m的水槽中進(jìn)行，即涌浪在寬度方向上的傳播可忽略不計(jì)。水槽的前端設(shè)置了長3 m且坡度固定為45°的滑坡體滑動(dòng)面，滑坡散粒體通過氣動(dòng)裝置啟動(dòng)，該裝置可以給滑坡提供初始速度。如圖1所示，試驗(yàn)以水面與滑動(dòng)面交界點(diǎn)為原點(diǎn)O，在滑動(dòng)面1號(hào)（距原點(diǎn)0.4 m）與2號(hào)（距原點(diǎn)0.07 m）斷面位置安裝了激光測距傳感器，以獲得入水前滑坡厚度隨時(shí)間變化的情況，同時(shí)在①與②區(qū)域用工業(yè)相機(jī)捕捉涌浪形態(tài)。

根據(jù)上述試驗(yàn)建立二維幾何模型如圖2（a）所示，模型采用四邊形網(wǎng)格，并對(duì)整個(gè)計(jì)算區(qū)域各相交界面附近的網(wǎng)格進(jìn)行了加密處理（見圖2（b）），網(wǎng)格總數(shù)約3萬個(gè)。值得說明的是，本文計(jì)算前對(duì)計(jì)算區(qū)域內(nèi)的網(wǎng)格開展了敏感性分析，發(fā)現(xiàn)進(jìn)一步加密網(wǎng)格對(duì)計(jì)算精度影響較小，說明目前的網(wǎng)格密度足以滿足計(jì)算要求。模型邊界條件包括速度進(jìn)口、壓力進(jìn)口、自由出流與固壁邊界?，F(xiàn)選取Fritz試驗(yàn)中的兩組工況進(jìn)行模擬驗(yàn)證，其中滑坡體密度ρs均為1 620 kg/m3，水體密度ρw為998.2 kg/m3；壁面摩擦系數(shù)f′根據(jù)Fritz的滑坡材料（PP-BaSO4）取0.445，其余參數(shù)設(shè)置見表1。

2.2 計(jì)算結(jié)果分析

2.2.1 滑坡滑動(dòng)形態(tài)對(duì)比

圖3（a）和圖3（b）分別給出了1號(hào)與2號(hào)斷面處激光測距傳感器在兩組試驗(yàn)中捕獲的滑坡厚度s隨時(shí)間t的變化情況，其中滑坡前緣與水面接觸的時(shí)刻記為t=0 s。從圖3中可以看出，工況一中數(shù)模計(jì)算得到的滑坡厚度隨時(shí)間的變化情況與試驗(yàn)結(jié)果基本吻合，較好地反映了散粒體滑動(dòng)的真實(shí)狀態(tài)；工況二中滑坡體較厚且滑動(dòng)速度快，滑坡體表層的變形主要受剪切應(yīng)力影響，使得計(jì)算得到的滑坡變形速率略大于試驗(yàn)結(jié)果。不過總體來看，滑坡厚度的微小差別對(duì)涌浪形態(tài)的影響可以忽略不計(jì)。

2.2.2 近場涌浪形態(tài)對(duì)比

圖4為工況一數(shù)模計(jì)算得到的涌浪形成過程與Fritz的試驗(yàn)結(jié)果對(duì)比（上方為實(shí)際形態(tài)，下方為模擬形態(tài)，下同），數(shù)模結(jié)果中深藍(lán)色代表空氣相，青色代表水相，紅色代表滑坡相。從圖4（a）中可以看出在相對(duì)時(shí)間T=0.93時(shí)涌浪形態(tài)存在一定差異，主要是由于滑坡在水下的滑動(dòng)速度比試驗(yàn)略快；從圖4（b）中可以看出計(jì)算得到的滑坡堆積形態(tài)與試驗(yàn)基本一致，但滑坡沖入水中使兩相摻混區(qū)域的動(dòng)力黏度會(huì)迅速增大，導(dǎo)致空氣很難摻混進(jìn)水相中，故無法較好地模擬出試驗(yàn)中白色的水氣摻混現(xiàn)象；圖4（c）中，速度梯度減小促使滑坡在水下逐漸停止運(yùn)動(dòng)，滑坡與水體兩相邊界附近的動(dòng)力黏度則逐漸增大，且遠(yuǎn)遠(yuǎn)大于水體的動(dòng)力黏度值，導(dǎo)致部分卷入水中且靠近滑坡邊界的空氣無法逸出，出現(xiàn)了一層很薄的空氣滯留在滑坡表面。

由于工況二滑坡厚度與滑動(dòng)速度較工況一顯著增加，其形成的涌浪高度也更高，滑坡體與水體接觸后呈現(xiàn)出明顯的分離狀態(tài)［17］。圖5（b）與圖5（c）中涌浪自由面與試驗(yàn)結(jié)果相比略有差異，主要體現(xiàn)為摻氣效果不佳，出現(xiàn)該問題的原因依然是水體表面黏度過大導(dǎo)致空氣難以摻入，無法模擬出水花飛濺的形態(tài)。但總體來看，摻氣效果不佳和空氣滯留的影響僅局限于滑坡相的邊界附近，不會(huì)對(duì)涌浪的總體形態(tài)和傳播衰減產(chǎn)生較大影響。值得說明的是，對(duì)滑坡體邊界層區(qū)域的網(wǎng)格進(jìn)行適當(dāng)加密，可以使滯留氣體的范圍減小。

3 討 論

改變阻力模型中的摩擦系數(shù)f′、湍流阻力系數(shù)ζ，及內(nèi)摩擦系數(shù)K均會(huì)對(duì)滑坡運(yùn)動(dòng)特征及涌浪形態(tài)產(chǎn)生影響，故基于上述算例對(duì)3個(gè)參數(shù)進(jìn)行單因素影響分析。以工況一為例，當(dāng)湍流阻力系數(shù)ζ為定值2 000 m/s2，內(nèi)摩擦系數(shù)K為定值0.5，摩擦系數(shù)f′分別取值0.1，0.3，0.5時(shí)，可得到2號(hào)斷面處滑坡厚度隨時(shí)間變化的情況。從圖6可以看出，摩擦系數(shù)f′越大，對(duì)滑坡阻力作用越明顯，減緩了滑坡的變形速率。

當(dāng)摩擦系數(shù)f′取定值0.5，內(nèi)摩擦系數(shù)K取定值0.1，湍流阻力系數(shù)ζ分別取100，500，1 500 m/s2時(shí)，滑坡體在2號(hào)斷面處隨時(shí)間變化的情況如圖7所示，從圖7中可以看出改變湍流阻力系數(shù)ζ對(duì)滑坡厚度變化的影響較小，但湍流阻力系數(shù)ζ增大會(huì)減緩滑坡后緣的變形速率。

當(dāng)摩擦系數(shù)f′取定值0.5，湍流阻力系數(shù)ζ取定值100 m/s2，內(nèi)摩擦系數(shù)K分別取0.1，0.5，0.9，從圖8可以看出內(nèi)摩擦系數(shù)K對(duì)滑坡變形速率影響較大。由于增大內(nèi)摩擦系數(shù)K值會(huì)使滑坡內(nèi)部流速梯度減小，故內(nèi)摩擦系數(shù)K越大，滑坡體的變形速率就越慢。

4 結(jié) 論

（1） 本文的數(shù)值模擬方法能夠較好地模擬碎石土滑坡沖擊水體形成涌浪的全過程，通過獲得涌浪波高、流場等信息，可以為山區(qū)水庫滑坡涌浪災(zāi)害預(yù)測提供依據(jù)。

（2） 在提高涌浪模擬精度的同時(shí)，該模擬方法相比離散元法具有計(jì)算量適中、輸入控制參數(shù)較少等優(yōu)勢。

（3） 模型中的摩擦系數(shù)f′、湍流阻力系數(shù)ζ及內(nèi)摩擦系數(shù)K均會(huì)對(duì)滑坡運(yùn)動(dòng)特征產(chǎn)生影響，其中摩擦系數(shù)與內(nèi)摩擦系數(shù)同滑坡體的變形速率成反比關(guān)系，而湍流阻力系數(shù)ζ的增大僅能減緩滑坡后緣的變形速率。

（4） 值得特別說明的是，阻力模型會(huì)使滑坡與水、氣交界區(qū)域的動(dòng)力黏度迅速增大，導(dǎo)致在模擬水花飛濺等強(qiáng)摻氣現(xiàn)象時(shí)效果不佳，同時(shí)會(huì)使小部分氣體滯留在滑坡與水的交界面處難以逸出水面，針對(duì)上述問題還需進(jìn)一步優(yōu)化阻力模型。此外，在實(shí)際工程中，庫岸碎石土滑坡體特征千變?nèi)f化，還應(yīng)對(duì)碎石土滑坡的流變特性開展進(jìn)一步研究，為模擬計(jì)算提供更加準(zhǔn)確的控制參數(shù)取值范圍。

參考文獻(xiàn)：

［1］ YIN Y P，HUANG B L，CHEN X T，et al.Numerical analysis on wave generated by the Qianjiangping landslide in Three Gorges Reservoir，China［J］.Landslides，2015，12（2）：355-364.

［2］ 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.

［3］ XIAO L L，WANG J J，WARD S N，et al.Numerical modeling of the June 24，2015，Hongyanzi landslide generated impulse waves in Three Gorges Reservoir，China［J］.Landslides，2018，15（12）：2385-2398.

［4］ ALEXANDRE P，PHILIPPE H，STPHANE A.Landslide tsunamis：Comparison between depth-averaged and Navier-Stokes models［J］.Coastal Engineering，2021，170：104022.

［5］ ZHANG G，CHEN J，QI Y，et al.A WCSPH two-phase mixture model for tsunami waves generated by granular landslides［J］.Computers and Geotechnics，2022，144：104657

［6］ MAGDALENA I，F(xiàn)IRDAUS K，JAYADI D.Analytical and numerical studies for wave generated by submarine landslide［J］.Alexandria Engineering Journal，2022，61：303-7313

［7］ 鄧成進(jìn)，黨發(fā)寧，陳興周.庫區(qū)滑坡涌浪傳播及其與大壩相互作用機(jī)理研究［J］.水利學(xué)報(bào)，2019，50（7）：815-823.

［8］ 徐文杰.滑坡涌浪流-固耦合分析方法與應(yīng)用［J］.巖石力學(xué)與工程學(xué)報(bào)，2020，39（7）：1420-1433.

［9］ 譚海，徐青，陳勝宏，等.基于DEM-SPH耦合模型的散粒體滑坡涌浪仿真模擬［J］.巖土力學(xué)，2020（增2）：1-11.

［10］ HUNGR O，EVANS S G，BOVIS M J，et al.A review of the classification of landslides of the flow type［J］.Environmental & Engineering Geoscience，2001，7（3）：221-238.

［11］ HUNGR O，MCDOUGALL S.Two numerical models for landslide dynamic analysis［J］.Computers & Geosciences，2009，35（5）：978-992.

［12］ 朱圻，程謙恭，王玉峰，等.高速遠(yuǎn)程滑坡超前沖擊氣浪三維動(dòng)力學(xué)分析［J］.巖土力學(xué)，2014，35（10）：2909-2926.

［13］ 袁小一，許強(qiáng)，程謙恭，等.高速遠(yuǎn)程滑坡-碎屑流超前沖擊氣浪分析［J］.水文地質(zhì)工程地質(zhì)，2016，43（6）：113-119.

［14］ 高楊，李濱，王國章.雞尾山高速遠(yuǎn)程滑坡運(yùn)動(dòng)特征及數(shù)值模擬分析［J］.工程地質(zhì)學(xué)報(bào)，2016，24（3）：425-434.

［15］ 夏式偉，鄭昭煬，袁小一，等.蘆山地震湯家溝滑坡-碎屑流過程模擬［J］.山地學(xué)報(bào)，2017（4）：102-109.

［16］ FRITZ H M.Initial phase of landslide generated impulse waves［D］.Zurich：ETH Zürich，2002.

［17］ XUE H C，MA Q，DIAO M J，et al.Propagation characteristics of subaerial landslide-generated impulse waves［J］.Environmental Fluid Mechanics，2019，19（1）：203-230.

（編輯：胡旭東）

Numerical simulation method of gravel soil landslide surge considering resistance effect

XUE Hongcheng1，2，PENG Xingyao1，MA Qian3，QING Zhengyang1

（1.College of River and Ocean Engineering，Chongqing Jiaotong University，Chongqing 400074，China； 2.National Engineering Research Center for Inland Waterway Regulation，Chongqing Jiaotong University，Chongqing 400074，China； 3.Southwest Waterway Engineering Institute，Chongqing Jiaotong University，Chongqing 400016，China）

Abstract：

Landslide surge is one of the risk sources of water disasters in a reservoir area.It is of great engineering application value to carry out surge prediction research to improve channel safety and reduce near-shore flood risk.However，the dynamic process of landslide body entering water to form surge is very complicated.At present，it is still a difficult problem to accurately simulate it.Based on the assumption of equivalent fluid，the heterogeneous gravel soil landslide was regarded as a continuous special fluid，and a friction resistance program of the landslide body moving along the sliding surface was written by a custom function.The Realizable k-ε turbulence model and VOF method were used to establish the numerical simulation model of gravel soil landslide surge.The results showed that this method can better simulate the underwater movement and accumulation form of gravel soil landslide，and the accuracy of surge amplitude was high.In the resistance model，the friction coefficient and the internal friction coefficient were inversely proportional to the deformation rate of the landslide body，while the increasing of the turbulent flow resistance coefficient could only slow down the deformation rate of the trailing edge of the landslide.The research results can provide a reference for quantitatively evaluating the risk of surge disaster of gravel soil landslide and formulating the emergency plan of surge disaster in a reservoir area.

Key words：

landslide surge；numerical simulation；equivalent fluid；k-ε turbulence model；VOF model