肖華波，王澤皓，石偉明，王東坡，歐陽朝軍

(1.中國電建集團(tuán) 成都勘測(cè)設(shè)計(jì)研究院有限公司，四川 成都 610072；2.成都理工大學(xué) 地質(zhì)災(zāi)害防治與地質(zhì)環(huán)境保護(hù)國家重點(diǎn)實(shí)驗(yàn)室，四川 成都 610059；3.中國科學(xué)院 水利部 成都山地災(zāi)害與環(huán)境研究所，四川 成都 610041)

位于河道、水庫周邊的大型滑坡、崩塌等地質(zhì)災(zāi)害一旦發(fā)生，可能造成堵江并產(chǎn)生巨大的涌浪，對(duì)人民的生命財(cái)產(chǎn)安全、船舶安全及水電站的運(yùn)行等構(gòu)成威脅[1–3]。例如：2003年7月，位于三峽庫區(qū)的千將坪地區(qū)突然發(fā)生了一起方量為1.5×107m3的滑坡，該滑坡形成的堰塞壩將青干河堵塞，造成24名人員死亡，直接經(jīng)濟(jì)損失約700萬美元[4]。2013年7月，云南省永善縣黃華鎮(zhèn)黃坪村發(fā)生山體滑坡，滑坡體沖入金沙江激起15～20 m高的涌浪，造成12人失蹤，多部車輛損毀[5]。2015年6月，重慶市巫山縣紅巖子發(fā)生滑坡，其方量約為2.3×105m3，坡體入水后產(chǎn)生5 m多高的涌浪，造成2人死亡，6人受傷，13條船只毀壞[6–7]。

數(shù)值模擬方法是國內(nèi)外開展滑坡–涌浪災(zāi)害相關(guān)研究的重要手段[8–9]，可較為快速、準(zhǔn)確、全面地分析滑坡–涌浪動(dòng)力學(xué)過程。然而，由于滑坡入水產(chǎn)生涌浪的過程伴隨滑坡體與水體之間的流固耦合作用，較為復(fù)雜，且空間尺度較大，很難用單一介質(zhì)模型進(jìn)行求解?？紤]到流體和滑坡的不同力學(xué)性質(zhì)，應(yīng)采用不同的數(shù)值模擬方法。在處理離散顆粒時(shí)，離散單元法（discrete element method，DEM）具有很好的靈活性，能夠很好地反映滑坡體滑動(dòng)、碎裂的過程，因而被廣泛應(yīng)用于滑坡災(zāi)害動(dòng)力過程的數(shù)值模擬[10–12]；計(jì)算流體力學(xué)方法（computational fluid dynamics，CFD）在水體演進(jìn)過程的計(jì)算精度方面有其獨(dú)特優(yōu)勢(shì)[13]。通過兩種方法間的耦合求解，能改善傳統(tǒng)模擬方法在研究流固耦合現(xiàn)象時(shí)的不足，考慮滑坡運(yùn)動(dòng)過程中庫水與滑坡體間的水土相互作用，還原水庫滑坡在水環(huán)境影響下的真實(shí)演進(jìn)過程，在模擬時(shí)間和空間尺度都較大的真實(shí)滑坡涌浪災(zāi)害時(shí)具有一定可靠性。

CFD–DEM流固耦合模型已被廣泛應(yīng)用于化學(xué)分析、工業(yè)制造等領(lǐng)域，目前在地質(zhì)災(zāi)害領(lǐng)域也初步開展了一些應(yīng)用。Li等[14]使用具有遠(yuǎn)端相互作用的鍵合粒子網(wǎng)絡(luò)模擬柔性防護(hù)網(wǎng)，利用計(jì)算流體力學(xué)和離散元耦合研究泥石流對(duì)柔性防護(hù)網(wǎng)沖擊的動(dòng)力響應(yīng)。Park[15]針對(duì)潰壩過程中的固液氣混合流動(dòng)，研究了顆粒密度對(duì)自由表面和顆粒運(yùn)動(dòng)的影響，并對(duì)運(yùn)動(dòng)特征進(jìn)行了分類。Shi等[16]應(yīng)用CFD–DEM的耦合方法研究滑坡壩的滲流特性，定量分析了粗顆粒與細(xì)顆粒不同的滲流破壞機(jī)理。Zhao等[17]采用CFD–DEM耦合程序?qū)ζ矫鎽?yīng)變條件下Vajont滑坡過程進(jìn)行了2維分析研究，并將模擬結(jié)果與ALE–FEM方法獲得的結(jié)果進(jìn)行了比較。

本文利用開源程序CFDEM、OpenFOAM及LIGGGHTS，基于計(jì)算流體力學(xué)方法與離散單元法構(gòu)建流固耦合模型[18–20]；將該耦合模型應(yīng)用于真實(shí)情況下3維大尺度滑坡–涌浪災(zāi)害動(dòng)力過程演進(jìn)模擬；以四川省猴子巖水庫色玉滑坡為例，對(duì)滑坡涌浪災(zāi)害的動(dòng)力學(xué)過程進(jìn)行分析，一次解決滑坡變形失穩(wěn)及涌浪形成與傳播的問題，為山區(qū)大區(qū)域尺度大型滑坡–涌浪復(fù)合災(zāi)害的分析研究提供新的技術(shù)手段。

1 CFD–DEM耦合模型

1.1 計(jì)算流體力學(xué)控制方程

CFD–DEM耦合模型中，流體包含水體與空氣兩相，采用VOF模型求解多相流的界面問題。在VOF模型中，流體界面的演化通過求解以下方程來描述：

式中，i為網(wǎng)格數(shù)，ρi、μi分別為第i個(gè)網(wǎng)格的流體密度和黏度。

在CFD–DEM耦合模型中，流體的控制方程為Navier–Stokes方程，如式（5）～（6）所示：

1.2 離散單元法控制方程

在CFD–DEM耦合模型中，顆粒的運(yùn)動(dòng)狀態(tài)由牛頓第二定律進(jìn)行描述，包括平動(dòng)及轉(zhuǎn)動(dòng)兩部分。除顆粒之間的相互作用力和力矩外，平動(dòng)及轉(zhuǎn)動(dòng)兩部分控制方程中還分別添加了流體作用于顆粒上的力和力矩，方程如下：

式（9）～（10）中，mi為 顆粒i的質(zhì)量，ui為顆粒i的平動(dòng)速度，nj為與顆粒j接觸的顆?？倲?shù)，F(xiàn)ij為顆粒i與顆粒j的相互作用力，F(xiàn)fi為流體作用在顆粒i上的合力，Ii為轉(zhuǎn)動(dòng)慣量， θi為顆粒i的旋轉(zhuǎn)角速度，Mij為顆粒j作用于顆粒i的力矩，Mif為流體作用于顆粒i的力矩，g為重力加速度。

顆粒間相互作用力的計(jì)算采用基于Hertz接觸理論的Hertz–Mindlin模型[21]。當(dāng)兩個(gè)顆粒的接觸距離小于兩者半徑之和時(shí)，將產(chǎn)生相互作用力，法向包括彈簧力和阻尼力兩項(xiàng)，切向包括剪切力和阻尼力兩項(xiàng)。該模型計(jì)算公式如式（11）所示：

式中，kn、kt分別為法向、切向彈性系數(shù)， δn為法向重疊量， δt為切向位移， γn、 γt分別為法向、切向阻尼系數(shù)，vn為法向相對(duì)速度，vt為切向相對(duì)速度。

1.3 流體–顆粒相互作用

在CFD–DEM耦合模型中，流體–顆粒間的相互作用通過在控制方程中交換相互作用力來考慮。本文主要考慮拖拽力和浮力的影響[21]，即：

式中：Cd為拖拽力系數(shù)；up為顆粒速度； χ為孔隙率的一個(gè)經(jīng)驗(yàn)修正系數(shù)，此系數(shù)使得該關(guān)系式能夠適用于更大范圍的雷諾數(shù)的情況。

1.4 CFD–DEM耦合流程

在CFD–DEM耦合框架下，應(yīng)用OpenFOAM程序求解流體運(yùn)動(dòng)方程，應(yīng)用LIGGGHTS程序計(jì)算散粒體滑坡顆粒運(yùn)動(dòng)方程；將顆粒與流體的相互作用力作為耦合紐帶，通過CFDEM耦合模塊進(jìn)行動(dòng)量交互傳遞。

CFD–DEM耦合流程為：1）建立計(jì)算模型，給定初始流場(chǎng)、初始顆粒位置、初始顆粒速度；2）利用DEM求解器計(jì)算作用于顆粒上的各個(gè)力，求解顆粒動(dòng)量方程，更新每個(gè)顆粒的位置和速度；3）完成DEM循環(huán)后，計(jì)算CFD網(wǎng)格單元中的孔隙率與流體–顆粒相互作用力，并將其傳遞到CFD循環(huán)；4）求解流體動(dòng)量方程和連續(xù)性控制方程，得到流體速度場(chǎng)與壓力場(chǎng)；5）流場(chǎng)收斂后，計(jì)算作用于顆粒上的流體作用力，代回DEM求解器；6）CFD–DEM循環(huán)計(jì)算，直至達(dá)到計(jì)算結(jié)束時(shí)間。

2 算例驗(yàn)證

2.1 顆粒堆積體坍塌–涌浪試驗(yàn)

為驗(yàn)證CFD–DEM流固耦合模型模擬滑坡–涌浪災(zāi)害的有效性，對(duì)Robbe–Saule等[25]所做的室內(nèi)水槽試驗(yàn)進(jìn)行數(shù)值模擬。該試驗(yàn)過程包含了顆粒堆積體失穩(wěn)、涌浪形成及傳播等類似滑坡入水的關(guān)鍵過程，研究了顆粒堆積體的高寬比和體積、顆粒的直徑和密度等因素對(duì)涌浪高度的影響情況。因此，被廣泛用作標(biāo)準(zhǔn)檢驗(yàn)?zāi)Ｐ汀?/p>

試驗(yàn)裝置左側(cè)為顆粒堆積區(qū)，高寬比為H0/L0；右側(cè)為一個(gè)長2.00 m、高0.30 m、寬0.15 m的玻璃水槽，其中充滿深度為0.05 m的水。通過提拉擋板釋放顆粒，通過改變顆粒堆積體的高度和寬度控制體積變量V。數(shù)值模擬使用的計(jì)算參數(shù)見表1。

表1 顆粒堆積體坍塌–涌浪試驗(yàn)計(jì)算參數(shù)Tab. 1 Calculation parameters of the particle accumulation collapse–surge experiment

2.2 對(duì)比驗(yàn)證結(jié)果

圖1為在H0/L0=2.5，V=10.2 dm3工況下，顆粒堆積體坍塌演進(jìn)過程中不同時(shí)刻的水槽試驗(yàn)結(jié)果（圖1左）與數(shù)值模擬結(jié)果（圖1右）。對(duì)比圖1中顆粒堆積體坍塌過程及水面涌浪運(yùn)動(dòng)狀態(tài)可以清晰看出，在整個(gè)演進(jìn)過程中，數(shù)值模擬結(jié)果與水槽試驗(yàn)結(jié)果基本保持一致。當(dāng)顆粒撞擊水面時(shí)會(huì)形成初始涌浪，隨著顆粒體大量侵入水體，涌浪高度隨之增加；而后，涌浪在傳播過程中逐漸變陡，波峰在重力影響下濺落到水體表面。除此之外，顆粒堆積體坍塌到水中產(chǎn)生涌浪的整個(gè)過程，涉及到顆粒、水體及空氣的復(fù)雜相互作用。由圖1試驗(yàn)結(jié)果可以觀察到，水體在短時(shí)間內(nèi)并不能完全充滿顆粒間隙，在顆粒堆積體坍塌過程中將形成一個(gè)明顯的凹形干濕邊界；對(duì)比數(shù)值模擬結(jié)果可知，CFD–DEM流固耦合模型能夠較為準(zhǔn)確地捕捉這一現(xiàn)象，初步驗(yàn)證了該模型解決流固耦合問題的可靠性。

圖1 顆粒堆積體坍塌演進(jìn)過程對(duì)比Fig. 1 Comparison of particle accumulation collapse evolution process

通過提取水槽試驗(yàn)及數(shù)值模擬中各時(shí)刻涌浪高度，并進(jìn)行對(duì)比分析，結(jié)果如圖2所示。從圖2中的涌浪峰值高度看，試驗(yàn)中測(cè)量的最大涌浪高度為8.3 cm，模擬計(jì)算得到的最大涌浪高度為8.1 cm，結(jié)果較為一致。從圖2中涌浪高度演進(jìn)過程來看，模擬結(jié)果稍微超前于試驗(yàn)結(jié)果，但計(jì)算模擬的整體趨勢(shì)及高差范圍與試驗(yàn)結(jié)果較為吻合。分析造成誤差的原因有兩點(diǎn)：一是，實(shí)際試驗(yàn)中需提拉擋板，造成時(shí)間誤差，導(dǎo)致涌浪延后產(chǎn)生，在提升過程中對(duì)顆粒堆積體也會(huì)造成擾動(dòng)，一定程度上影響涌浪產(chǎn)生過程；二是，試驗(yàn)中顆粒大多是不規(guī)則的，為簡化模擬采用粒徑單一的球形顆粒代替，會(huì)低估顆粒間的相互作用。

圖2 數(shù)值模擬與物理試驗(yàn)得到的涌浪高度對(duì)比Fig. 2 Comparison of surge height between numerical and experimental results

初始顆粒堆積體的體積對(duì)涌浪的形成起著重要作用，因此，Robbe–Saule等[25]進(jìn)行了一系列不同顆粒堆積體體積工況下的試驗(yàn)。根據(jù)對(duì)應(yīng)的試驗(yàn)參數(shù)，本文進(jìn)行顆粒堆積體體積分別為3.71、5.69、7.42、10.15 dm3的4組數(shù)值模擬，以進(jìn)一步驗(yàn)證CFD–DEM流固耦合模型在不同體積工況下的適用性。各組工況下試驗(yàn)與模擬結(jié)果的最大涌浪高度如圖3所示。由圖3可知：在一定范圍內(nèi)，最大涌浪高度隨初始顆粒堆積體體積的增大而增大；在各組顆粒堆積體體積工況下，模擬與試驗(yàn)結(jié)果的最大涌浪高度擬合度較高，誤差均在8%以內(nèi)。

圖3 不同顆粒體積工況下最大涌浪高度對(duì)比Fig. 3 Comparison of maximum surge height under different particle volume conditions

通過對(duì)顆粒坍塌運(yùn)動(dòng)過程、涌浪高度演化過程等進(jìn)行對(duì)比，發(fā)現(xiàn)顆粒運(yùn)動(dòng)及涌浪傳播模擬結(jié)果與試驗(yàn)結(jié)果吻合很好，可證明本文的流固耦合模型能有效還原顆粒運(yùn)動(dòng)狀態(tài)及涌浪的波動(dòng)變化。

3 色玉滑坡–涌浪災(zāi)害演進(jìn)過程分析

3.1 色玉滑坡–涌浪災(zāi)害概況

基于上述CFD–DEM流固耦合模型，對(duì)四川省猴子巖水庫色玉堆積體前緣滑坡進(jìn)行數(shù)值模擬研究，分析滑坡失穩(wěn)、涌浪形成和傳播的災(zāi)害演進(jìn)全過程。

猴子巖水電站位于四川省甘孜藏族自治州康定縣境內(nèi)大渡河干流上；色玉堆積體位于電站庫區(qū)右岸，距下游大壩約4.3 km。色玉堆積體順河長約650 m，橫河寬約350 m，前、后緣高程分別為1 820、2 100 m，高差約為280 m。高程約1 930～2 100 m處為較寬緩的臺(tái)狀地形，臺(tái)地前緣至坡腳地形坡度約為50～70°。堆積臺(tái)地主要由冰磧堆積物（glQ3）組成，臺(tái)地后緣覆蓋少量崩坡積堆積物（col+dlQ4）。堆積體厚度較大，一般為83～100 m；下伏基巖主要為泥盆系中上統(tǒng)河心組（D2–3h）灰、灰白色中厚層白云巖、白云質(zhì)灰?guī)r夾灰?guī)r。

2017年10月24日，色玉堆積體前緣發(fā)生滑坡，此時(shí)水庫蓄水位為1 832 m?；率Х€(wěn)后的地貌特征如圖4所示。由圖4可知，失穩(wěn)方量約80×104m3，失穩(wěn)高度最大達(dá)160 m，失穩(wěn)后在下游3.7 km處的電站進(jìn)水口形成1.7 m高的涌浪。截止目前，該堆積體又多次發(fā)生垮塌，并形成多處凹槽地形。由于該堆積體距離大壩相對(duì)較近，其失穩(wěn)涌浪對(duì)工程和周邊人類活動(dòng)安全的不利影響需重點(diǎn)關(guān)注和研究。

圖4 色玉滑坡全貌Fig. 4 Panorama of the Seyu landslide

3.2 計(jì)算模型及參數(shù)設(shè)置

通過對(duì)色玉堆積體失穩(wěn)前后地形數(shù)據(jù)進(jìn)行差值分析，確定模擬的物源。計(jì)算使用了2萬多個(gè)半徑為2 m的球形顆粒，參考Hu等[26]所做的色玉滑坡附近區(qū)域土工試驗(yàn)測(cè)試結(jié)果，設(shè)置模型中顆粒間的泊松比和楊氏模量等參數(shù)。色玉滑坡數(shù)值計(jì)算模型見圖5，計(jì)算域網(wǎng)格大小為40 m×40 m×40 m。為減小涌浪高度計(jì)算誤差，將庫水位上下20 m范圍內(nèi)的網(wǎng)格進(jìn)行加密，網(wǎng)格大小為10 m×10 m×10 m。計(jì)算參數(shù)選取見表2。

表2 色玉滑坡計(jì)算參數(shù)Tab. 2 Calculation parameters of Seyu landslide

圖5 色玉滑坡3維模型Fig. 5 3D model of Seyu landslide

3.3 模擬結(jié)果

為驗(yàn)證數(shù)值模擬的準(zhǔn)確性，將涌浪計(jì)算結(jié)果與實(shí)測(cè)數(shù)據(jù)進(jìn)行對(duì)比。圖5中，P點(diǎn)（225，–2 700）所示位置為色玉堆積體下游3.7 km的電站進(jìn)水口處，堆積體失穩(wěn)后，在該處實(shí)測(cè)到了高1.7 m的涌浪。提取P點(diǎn)的涌浪高度計(jì)算結(jié)果如圖6所示。

圖6 P點(diǎn)涌浪高度變化Fig. 6 Changes of surge height at point P

由圖6可知：在滑坡發(fā)生142 s后，首浪傳播到P點(diǎn)，高度為2.12 m，這與實(shí)測(cè)的數(shù)據(jù)較為接近。當(dāng)涌浪還沒有傳播到P點(diǎn)時(shí)，水面在±0.5 m范圍內(nèi)波動(dòng)，這是由于流體初始條件設(shè)置及部分網(wǎng)格不規(guī)則造成的誤差。通過對(duì)該處涌浪計(jì)算結(jié)果與實(shí)測(cè)數(shù)據(jù)進(jìn)行對(duì)比，可以證明色玉滑坡–涌浪災(zāi)害的數(shù)值模擬結(jié)果較為準(zhǔn)確。

色玉滑坡運(yùn)動(dòng)及涌浪傳播動(dòng)態(tài)演進(jìn)過程模擬結(jié)果如圖7所示，顯示了0～30 s時(shí)間內(nèi)滑坡及涌浪形態(tài)。由圖7可知：從滑坡體失穩(wěn)運(yùn)動(dòng)至靜止堆積的持續(xù)時(shí)間約為20 s，從滑坡入水處涌浪形成到涌浪傳播至對(duì)岸的持續(xù)時(shí)間約為10 s。t=5 s左右時(shí)，在滑坡顆粒沖擊作用下，涌浪開始逐漸形成；t=10 s，大部分滑坡顆?；?，涌浪高度進(jìn)一步增加，沿滑坡發(fā)生點(diǎn)呈圓弧狀向外擴(kuò)散；t=15 s，滑坡顆粒全部進(jìn)入水庫，近岸水體沿滑坡顆粒運(yùn)動(dòng)方向下陷，涌浪在對(duì)岸已經(jīng)開始爬升；t=20 s，滑坡顆粒不再運(yùn)動(dòng)，堆積于水庫底部，此刻對(duì)岸涌浪爬升到最大高度；t=20～30 s，涌浪開始衰退，并繼續(xù)向上下游傳播。

圖7 色玉滑坡運(yùn)動(dòng)及涌浪傳播動(dòng)態(tài)演進(jìn)過程Fig. 7 Dynamic evolution process of Seyu landslide movement and surge propagation

圖8展示了滑坡顆粒速度和涌浪速度在滑坡從失穩(wěn)運(yùn)動(dòng)到靜止堆積這段時(shí)間內(nèi)的變化情況。圖9展示了所有滑坡顆粒的最大速度及平均速度隨時(shí)間的變化曲線。

圖8 滑坡顆粒及涌浪速度Fig. 8 Velocity of the landslide particles and surge

圖9 滑坡顆粒速度隨時(shí)間變化曲線Fig. 9 Variation curves of landslide particles velocity with time

由圖8、9可知：隨著滑坡運(yùn)動(dòng)不斷發(fā)展，重力勢(shì)能轉(zhuǎn)換為動(dòng)能，滑坡顆粒速度開始處于快速增加階段，平均速度在6.0 s時(shí)達(dá)到最大值，為16.12 m/s；滑坡中部運(yùn)動(dòng)速度明顯高于兩側(cè)運(yùn)動(dòng)速度，滑坡前部運(yùn)動(dòng)速度高于滑坡后部運(yùn)動(dòng)速度。之后，由于滑坡前緣顆粒持續(xù)滑入水庫底部，受到地形影響，顆粒在水庫底部逐漸減速并開始堆積，滑坡體的平均速度也開始減小。直至滑坡發(fā)生約20 s后，平均速度減到0.04 m/s，可以認(rèn)為此時(shí)滑坡完全停止運(yùn)動(dòng)。雖然滑坡整體運(yùn)動(dòng)速度在6 s后減小，但仍有部分顆粒在激烈碰撞下達(dá)到較高的運(yùn)動(dòng)速度，單個(gè)顆粒的最大速度達(dá)到44.29 m/s，發(fā)生在10.0 s左右。

采用橫向剖面進(jìn)一步分析涌浪沿滑坡運(yùn)動(dòng)方向的傳播過程及涌浪的爬升高度，剖面線位置如圖10所示。為更加精確地顯示涌浪高度變化，在剖面上選取3個(gè)監(jiān)測(cè)點(diǎn)對(duì)各時(shí)刻的涌浪高度進(jìn)行提取分析，監(jiān)測(cè)點(diǎn)分別為滑坡入水處A點(diǎn)、河道中部B點(diǎn)及滑坡對(duì)岸處C點(diǎn)，3點(diǎn)坐標(biāo)分別為（472，1 018）、（561，946）、（658，868）。

圖10 剖面線及監(jiān)測(cè)點(diǎn)位置Fig. 10 Location of profile lines and monitoring points

圖11為剖面位置涌浪演進(jìn)過程模擬分析結(jié)果，可清楚觀察到滑坡發(fā)生后短時(shí)間內(nèi)產(chǎn)生巨大涌浪的全過程。

結(jié)合圖11可以看出，涌浪運(yùn)動(dòng)過程主要分為3個(gè)階段，具體如下：

圖11 涌浪演進(jìn)過程剖面Fig. 11 Profile of surge evolution process

1）涌浪產(chǎn)生及傳播階段：滑坡失穩(wěn)后在重力作用下沖擊水庫，受滑坡前緣顆粒擠壓，水面沿滑坡方向向下凹陷，涌浪高度快速增加并向?qū)Π秱鞑ァ?/p>

2）涌浪爬升階段：色玉滑坡方量較大，且由于河谷狹窄，涌浪在10 s后即傳播到對(duì)岸，在慣性力作用下開始沿坡面向上爬升，并在20 s左右爬升到最大高度27.32 m；在這一階段，滑坡入水處水體產(chǎn)生的凹陷開始收縮。

3）涌浪回流階段：在20 s之后，由于滑坡的運(yùn)動(dòng)趨于靜止，沖向?qū)Π兜挠坷嗽谥亓Φ淖饔孟?，開始回流。

剖面上監(jiān)測(cè)點(diǎn)的涌浪高度變化如圖12所示。

圖12 沿滑坡運(yùn)動(dòng)方向監(jiān)測(cè)點(diǎn)涌浪高度變化Fig. 12 Changes of surge height at monitoring points along the landslide movement direction

由圖12可知：A點(diǎn)顯示了滑坡發(fā)生后產(chǎn)生的初始涌浪高度，在3.5 s時(shí)達(dá)到6.84 m，隨后近岸水體在大量滑坡顆粒沖擊下迅速下陷；涌浪波峰在10.5 s時(shí)傳播至B點(diǎn)，此時(shí)河道中部的最大涌浪高度為16.62 m；C點(diǎn)為滑坡對(duì)岸附近，可顯示涌浪傳播到對(duì)岸邊時(shí)的高度，該點(diǎn)涌浪高度在15.5 s時(shí)達(dá)到最大值19.04 m?；庐a(chǎn)生的涌浪經(jīng)過A、B、C點(diǎn)之后，開始在對(duì)岸沿坡面爬升，最大爬升高度達(dá)到27.32 m。由此可見，在山區(qū)狹窄河谷地形條件下，水庫岸坡失穩(wěn)將產(chǎn)生巨大涌浪，尤其在對(duì)岸爬高最大。

為分析滑坡涌浪沿河道方向的傳播過程，在滑坡區(qū)與電站進(jìn)水口P點(diǎn)之間，等距離選取了D（525，200）、E（725，–525）、F（700，–1 272）、G（425，–1 975）、P（225，–2 700）5個(gè)監(jiān)測(cè)點(diǎn)進(jìn)行涌浪高度監(jiān)測(cè)，具體位置如圖5所示。沿河道方向各監(jiān)測(cè)點(diǎn)涌浪高度變化如圖13所示。由圖13可知：相較于滑坡涌浪近場(chǎng)區(qū)域，遠(yuǎn)場(chǎng)區(qū)域涌浪高度明顯降低。當(dāng)滑坡體高速?zèng)_擊水體時(shí)，能量主要沿滑動(dòng)方向傳遞，造成沿滑動(dòng)方向上的涌浪高度、傳播速度和能量大，沿河道方向較小。色玉滑坡失穩(wěn)入水后，涌浪傳播到各監(jiān)測(cè)點(diǎn)的時(shí)間分別為22、46、80、110、142 s。滑坡涌浪在河流彎道處EF段的傳播時(shí)間最長，在直道段傳播時(shí)間稍短。以上結(jié)果表明滑坡涌浪在沿河道傳播過程中，其傳播特征受到庫區(qū)地形條件的影響較大。各監(jiān)測(cè)點(diǎn)最大涌浪高度分別為4.26、3.75、3.15、3.09、2.12 m，均出現(xiàn)在第1個(gè)波峰，后續(xù)波峰峰值衰減較大。沿程涌浪是一個(gè)逐漸衰減過程，最大涌浪高度隨著傳播距離的增加而減小。

圖13 沿河道方向監(jiān)測(cè)點(diǎn)涌浪高度變化Fig. 13 Changes of surge height at monitoring points along the river

4 結(jié) 論

1）介紹了一種CFD–DEM耦合算法，流體的控制方程為連續(xù)性方程及動(dòng)量方程，顆粒的運(yùn)動(dòng)由牛頓第二定律控制，通過基于拖拽力模型的動(dòng)量交換項(xiàng)實(shí)現(xiàn)流體與固體之間的耦合。

2）基于顆粒堆積體坍塌試驗(yàn)，進(jìn)行了不同顆粒體積工況下的模型驗(yàn)證，模擬結(jié)果與試驗(yàn)結(jié)果吻合較好，證明了耦合模型能準(zhǔn)確模擬顆粒與流體的相互作用過程。

3）將耦合模型應(yīng)用于實(shí)際發(fā)生的大規(guī)?；氯胨C涌浪災(zāi)害鏈的工程四川省猴子巖水庫色玉滑坡–涌浪的模擬。結(jié)果表明：滑坡體從失穩(wěn)運(yùn)動(dòng)至靜止堆積的持續(xù)時(shí)間約20 s，顆粒平均速度最大達(dá)到16.12 m/s；滑坡引起的涌浪在約10 s后傳播到對(duì)岸，之后開始沿坡面向上爬升，最大爬升高度達(dá)到27.32 m。

CFD–DEM耦合模型不僅對(duì)室內(nèi)坍塌–涌浪試驗(yàn)的模擬具有高效性和準(zhǔn)確性，也可以應(yīng)用于山區(qū)河谷大規(guī)?；氯胨C涌浪災(zāi)害鏈這一復(fù)雜流固耦合問題的模擬，實(shí)現(xiàn)滑坡失穩(wěn)運(yùn)動(dòng)及涌浪形成傳播全過程分析，為山區(qū)河谷防災(zāi)減災(zāi)提供新的技術(shù)手段。