薛宏程 彭杏瑤 馬倩 卿正陽

摘要：

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

關(guān) 鍵 詞：

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

中圖法分類號： TV139

文獻標志碼： A

DOI：10.16232/j.cnki.1001-4179.2023.07.022

0 引 言

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

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

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

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

2.1 模型構(gòu)建

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

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

2.2 計算結(jié)果分析

2.2.1 滑坡滑動形態(tài)對比

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

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

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

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

3 討 論

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

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

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

4 結(jié) 論

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

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

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

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

參考文獻：

［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］ 鄧成進，黨發(fā)寧，陳興周.庫區(qū)滑坡涌浪傳播及其與大壩相互作用機理研究［J］.水利學報，2019，50（7）：815-823.

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

［9］ 譚海，徐青，陳勝宏，等.基于DEM-SPH耦合模型的散粒體滑坡涌浪仿真模擬［J］.巖土力學，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］ 朱圻，程謙恭，王玉峰，等.高速遠程滑坡超前沖擊氣浪三維動力學分析［J］.巖土力學，2014，35（10）：2909-2926.

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

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

［15］ 夏式偉，鄭昭煬，袁小一，等.蘆山地震湯家溝滑坡-碎屑流過程模擬［J］.山地學報，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