李東陽,年廷凱*,吳 昊,張彥君
(1.大連理工大學(xué) 土木工程學(xué)院,遼寧 大連 116024;2.中國科學(xué)院 水利部成都山地災(zāi)害與環(huán)境研究所,四川 成都 610299;3.中國地質(zhì)調(diào)查局武漢地質(zhì)調(diào)查中心,湖北 武漢 430205)
中國西南地區(qū),尤其是青藏高原地區(qū),溝壑林立且河流發(fā)育,一直是滑坡類災(zāi)害的高發(fā)區(qū)[1–2]?;麦w阻塞河流所形成的堰塞壩可誘發(fā)如上游洪澇、潰壩洪水及泥石流等一系列次生災(zāi)害[3–4]。同時(shí),滑坡體侵入河流造成的涌浪會(huì)導(dǎo)致災(zāi)害影響范圍顯著擴(kuò)大[5–6]?;篓C堵江–涌浪災(zāi)害鏈具有致災(zāi)機(jī)理復(fù)雜、致災(zāi)范圍廣及破壞性強(qiáng)等顯著特點(diǎn),可對(duì)沿岸居民生命財(cái)產(chǎn)安全及流域內(nèi)水利設(shè)施安全造成嚴(yán)重威脅。
滑坡–堵江–涌浪災(zāi)害鏈演化過程涉及復(fù)雜的滑坡–河流多相動(dòng)力演化及耦合作用,從數(shù)值角度對(duì)其演化過程進(jìn)行模擬是一個(gè)極大的挑戰(zhàn)[7]。一些研究基于連續(xù)力學(xué)理論的雙層深度平均方法將滑坡體及河流用兩相流體建模[7–9]。然而,上述研究假定滑坡體為連續(xù)介質(zhì)流體,對(duì)固–液耦合作用的描述過于簡化,且忽略了滑坡體具有孔隙的非連續(xù)特性[10]。相較而言,基于非連續(xù)介質(zhì)理論的離散單元法(discrete element method,DEM)通過顆粒對(duì)土體建模,在處理土體大變形問題中具有突出優(yōu)勢(shì)[11]。DEM與計(jì)算流體動(dòng)力學(xué)(computational fluid dynamics,CFD)的耦合方法能夠有效處理土–水耦合作用問題,已成功應(yīng)用于水下滑坡模擬[12–13]。然而,傳統(tǒng)的固–液兩相DEM–CFD耦合方法并不足以描述滑坡涌浪的形成及傳播過程。此外,局部平均的DEM–CFD耦合方法受限于網(wǎng)格顆粒臨界尺寸比,難以滿足滑坡堵江及涌浪過程模擬中滑坡體粗顆粒材料的建模及高分辨率真實(shí)地形網(wǎng)格的要求,這些都限制了DEM–CFD耦合數(shù)值方法在滑坡–堵江–涌浪災(zāi)害鏈模擬中的應(yīng)用。
本文旨在通過數(shù)值模擬方法研究滑坡–堵江–涌浪災(zāi)害鏈演化過程。首先,通過在局部平均DEM–CFD耦合方法中引入流體體積法(volume of fluid,VOF)建立了一套能夠描述河流自由水面演化的擴(kuò)展DEM–CFD耦合數(shù)值框架,實(shí)現(xiàn)了對(duì)滑坡涌浪形成及傳播過程的追蹤。然后,提出了一種新的局部孔隙率計(jì)算方法,即虛擬球模型,來克服網(wǎng)格顆粒臨界尺寸比的限制。在對(duì)該方法的準(zhǔn)確性和有效性進(jìn)行驗(yàn)證后,在模型尺度上模擬了滑坡堵江及涌浪演化過程,并分析其滑坡–河流耦合作用機(jī)制。最后,將其應(yīng)用至2018年10月11日金沙江白格滑坡堵江事件的數(shù)值反演。
本研究所提出的DEM–CFD耦合數(shù)值方法中,滑坡體通過DEM建模,并通過求解牛頓第二定律來模擬滑坡體的運(yùn)移過程[14]:
式(1)~(2)中,mi、vi和ωi分別為顆粒質(zhì)量、速度及角速度,g為重力加速度,F(xiàn)ij和Mij分別為顆粒的接觸作用力及扭矩,F(xiàn)f為流體對(duì)顆粒的耦合作用力。
對(duì)于流體相,通過CFD求解顆粒相體積分?jǐn)?shù)修正后的局部平均Navier–Stokes方程[15]:
式(3)~(4)中, ε為局部孔隙率,u、 ρf、 μf分別為流體速度、密度及黏度,p為流體壓力,fpf為表征顆粒對(duì)流體阻力的動(dòng)量源項(xiàng),式(4)右側(cè)其余3項(xiàng)則依次為壓力梯度項(xiàng)、應(yīng)力項(xiàng)及重力項(xiàng)。
在本研究中,將VOF法引入傳統(tǒng)的局部平均DEM–CFD耦合方法中,從而允許對(duì)多相流體介質(zhì)建模,涌浪的形成及演化過程通過對(duì)多相流體交界面的追蹤來實(shí)現(xiàn)。本研究的擴(kuò)展DEM–CFD耦合方法中,每個(gè)流體單元被分為水相和空氣相,并用 αw和 αa分別表示相應(yīng)相的體積分?jǐn)?shù)。其中,水相的體積分?jǐn)?shù)滿足0 ≤αw≤1,則空氣相的體積分?jǐn)?shù)可表示為[16]:
因此, 0 <αw<1被定義為自由液面的存在。式(4)中的流體密度 ρf及黏度 μf,參考兩相流體性質(zhì),通過線性插值計(jì)算。通過求解基于流體速度的輸運(yùn)方程,自由液面運(yùn)動(dòng)可計(jì)算如下[17]:
式中,ur為兩相介質(zhì)相對(duì)速度,ur=u1?u2。值得注意的是,方程中已經(jīng)考慮顆粒相的存在。左側(cè)第1項(xiàng)中ε的時(shí)間導(dǎo)數(shù)考慮到了體積分?jǐn)?shù)變化對(duì)液相體積分?jǐn)?shù)的影響。式(6)的前兩項(xiàng)與式(4)具有相同的平流形式,第3項(xiàng)為人工項(xiàng),用于減少數(shù)值耗散,以達(dá)到使自由液面更為清晰的目的。
本研究中,流體與顆粒之間的耦合作用力Ff包括浮力Fb和拖曳力Fd,其中,浮力Fb基于壓力梯度求解,可表示為:
式中,Vp為顆粒體積。拖曳力是由流體域顆粒間的相對(duì)運(yùn)動(dòng)引起的,作用于顆粒質(zhì)心,其方向與相對(duì)運(yùn)動(dòng)方向相反。通過Di Felice模型來計(jì)算拖曳力[18]:
局部孔隙率 ε是顆粒相在流體控制方程中的重要表征,對(duì)于計(jì)算流體對(duì)顆粒的拖曳力至關(guān)重要。然而,在局部平均DEM–CFD耦合數(shù)值方法中,顆粒尺寸接近或大于最小流體網(wǎng)格尺寸時(shí),會(huì)導(dǎo)致顆粒相占據(jù)的體積分?jǐn)?shù)接近甚至大于網(wǎng)格體積,從而導(dǎo)致流體域孔隙度場(chǎng)不連續(xù)并使數(shù)值收斂困難[19]。因此,對(duì)于最小網(wǎng)格尺寸與最大顆粒粒徑的比值,建議保持在4倍以上[20]。但對(duì)于真實(shí)案例模擬而言,復(fù)雜地形建模顯然要求高分辨率的地形網(wǎng)格,因此限制了傳統(tǒng)DEM–CFD耦合方法在滑坡–堵江–涌浪災(zāi)害鏈演化模擬中的應(yīng)用。為此,本研究提出了一種基于虛擬球模型的局部孔隙率計(jì)算的改進(jìn)方法,如圖1所示。
圖1 虛擬球模型示意圖Fig. 1 Schematic of the virtual sphere model
在擴(kuò)展DEM–CFD耦合系統(tǒng)中,孔隙率計(jì)算的網(wǎng)格范圍通過虛擬球判定,其質(zhì)心與真實(shí)顆粒相同,直徑為4倍真實(shí)顆粒直徑。真實(shí)顆粒對(duì)孔隙率計(jì)算的貢獻(xiàn)體積被均勻劃分至虛擬球確定的相關(guān)網(wǎng)格。此外,同一網(wǎng)格被允許關(guān)聯(lián)至不同虛擬球。因此,某一網(wǎng)格單元的局部孔隙率可通過式(1)計(jì)算得出[21]:
本研究中,擴(kuò)展的DEM–CFD耦合方法基于離散元軟件EDEM及CFD軟件Fluent實(shí)現(xiàn),其中耦合模塊基于C語言開發(fā)。其計(jì)算流程如圖2所示。
圖2 DEM–CFD耦合計(jì)算流程Fig. 2 Numerical solution strategy for the extended coupled DEM–CFD model
通過DEM對(duì)滑坡體建模,并求解其運(yùn)移過程。而在CFD模塊中,首先通過引入的VOF模型追蹤河流自由液面,其次通過壓力隱式算子分裂(pressure-implicit with splitting of operators,PISO)算法迭代求解流體相控制方程,進(jìn)而求解計(jì)算域內(nèi)流體演化。DEM及CFD的時(shí)間步長分別滿足Rayleigh時(shí)間及Courant–Friedrichs–Lewy(CFL)條件的約束。在耦合模塊中,基于虛擬球模型計(jì)算局部孔隙率,并進(jìn)一步計(jì)算包括浮力及拖曳力在內(nèi)的耦合作用力。隨后,將其傳輸至DEM模塊進(jìn)行接觸力計(jì)算,而顆粒對(duì)流體的作用力則通過顯式及隱式動(dòng)量源項(xiàng)代入流體控制方程,從而實(shí)現(xiàn)雙向耦合計(jì)算。
為了驗(yàn)證本研究中擴(kuò)展DEM–CFD耦合方法的有效性及準(zhǔn)確性,采用單顆粒從空氣入水沉降模型,并將顆粒沉降速度的數(shù)值解與理論解進(jìn)行對(duì)比。單顆粒沉降驗(yàn)證案例如圖3所示。從圖3(a)可以看出:計(jì)算域?yàn)?.1 m×0.1 m×0.2 m的立方體,并采用4 mm尺寸的正六面體網(wǎng)格均勻離散。其中,上半部分為空氣,下半部分為水。直徑為 1 mm 的顆粒在距水面0.05 m處生成并從靜止?fàn)顟B(tài)開始下落,模擬所需具體參數(shù)見表1。
圖3 單顆粒沉降驗(yàn)證Fig. 3 Verification of single particle sedmentation
表1 單顆粒沉降模擬參數(shù)Tab. 1 Parameters for single particle sedimentation simulation
根據(jù)式(12),采用顯式正向差分法分別計(jì)算不同時(shí)刻顆粒的沉降速度理論解。顆粒沉降速度–時(shí)間曲線如圖3(b)所示,由圖3(b)可見,顆粒沉降過程的理論解及數(shù)值解表現(xiàn)出良好的一致性,并都達(dá)到了0.134 m/s的穩(wěn)定沉降速度。
保持顆粒直徑不變,將網(wǎng)格尺寸設(shè)定為1~5 mm以改變網(wǎng)格尺寸與顆粒直徑比,分別采用傳統(tǒng)的中心點(diǎn)模型[22]和虛擬球模型計(jì)算孔隙率,并進(jìn)行顆粒沉降過程模擬,結(jié)果如圖4所示。
圖4 中心點(diǎn)模型和虛擬球模型中尺寸效應(yīng)對(duì)比Fig. 4 Comparison of size ratio in central model and virtual sphere model
由圖4可知:在中心點(diǎn)模型中,隨著尺寸比的不斷減小,顆粒最終沉降速度的誤差不斷增大,這一結(jié)果反映了局部平均DEM–CFD方法的臨界尺寸效應(yīng)。與之相反,虛擬球模型在不同網(wǎng)格尺寸下的模擬結(jié)果相對(duì)穩(wěn)定,其結(jié)果與理論解的相對(duì)誤差在1%以內(nèi),且計(jì)算的顆粒最終沉降速度的精度最大能提高 40%。結(jié)果表明,虛擬球模型有效地消除了局部孔隙率計(jì)算中網(wǎng)格尺寸與顆粒直徑比的限制。
滑坡體延坡面快速運(yùn)移,隨后侵入并阻塞河流,同時(shí)激發(fā)涌浪延河道及對(duì)岸邊坡傳播。在堵江過程中,碎屑物質(zhì)進(jìn)入河谷,隨后到達(dá)對(duì)岸并沉積下來形成堰塞壩?;篓C堵江–涌浪演化過程模擬模型如圖5所示。由圖5可知,滑坡–堵江–涌浪災(zāi)害鏈?zhǔn)且环N典型的斜坡運(yùn)動(dòng)–河流系統(tǒng)的耦合作用現(xiàn)象,從模型尺度上建立災(zāi)害鏈演化的簡化模擬模型,有助于理解滑坡–河流相互作用過程。
圖5 滑坡–堵江–涌浪演化過程模擬模型Fig. 5 Numerical modelling of the landslide river blocking and impulse wave evolution
野外調(diào)查及統(tǒng)計(jì)數(shù)據(jù)表明,坡度在30°~45°的V型峽谷為滑坡堵江災(zāi)害的高發(fā)區(qū)域[23]。因此,模型設(shè)置如圖5所示,研究區(qū)域?yàn)? m×4 m的V型峽谷,其兩側(cè)岸坡均設(shè)置為35°,滑坡運(yùn)動(dòng)方向設(shè)定為與河道垂直。共計(jì)520 kg的滑坡體,初始位置設(shè)置在模型中段頂部,其初始速度為7 m/s。河流水深為0.3 m,流速為0.5 m/s。在 DEM 和 CFD 模塊中,岸坡表面均施加無滑移邊界條件,CFD模塊中頂部邊界設(shè)為開放邊界。計(jì)算所需參數(shù)見表2,其中DEM參數(shù)已在之前的研究中經(jīng)過標(biāo)定和驗(yàn)證[24]。
表2 滑坡壩形成過程模擬參數(shù)Tab. 2 Simulation parameters for landslide dam formation process
圖6為滑坡–堵江–涌浪演化的模擬結(jié)果。由圖6可知:滑坡體侵入河流后,與河水、河谷之間的劇烈動(dòng)能交換導(dǎo)致涌浪被激發(fā)。此外,碰撞下的能量耗散導(dǎo)致滑坡體前緣速度驟降,隨后對(duì)中后部滑坡體產(chǎn)生阻礙作用致使其開始沿河道方向擴(kuò)散。與此同時(shí),滑坡體前緣受到中后部滑坡體的推擠作用,在河谷地形控制下轉(zhuǎn)向,并混合涌浪一起沿對(duì)岸山谷爬升。當(dāng)滑坡體–涌浪的混合體沿對(duì)側(cè)山谷坡面爬升達(dá)到最大高度后便開始回落,滑坡體在河谷中穩(wěn)定并沉積形成滑坡壩。在此過程中,河流被完全阻塞。隨著上游來水持續(xù)被阻,壩前水位下降,而壩后水位不斷上升并形成堰塞湖。
圖6 滑坡–堵江–涌浪演化模擬結(jié)果Fig. 6 Simulation results of the river blockage and impulse wave evolution
圖7顯示了滑坡–堵江–涌浪災(zāi)害鏈過程模擬中的耦合作用力演化過程。由圖7可以觀察到:在滑坡體侵入河流的瞬間,顆粒–流體耦合作用力出現(xiàn)顯著的短暫下降,這是由于河流對(duì)滑坡體的制動(dòng)作用導(dǎo)致前緣速度降低。隨著滑坡體持續(xù)侵入河流,耦合作用力繼續(xù)上升,直到滑坡體到達(dá)谷底?;麦w與谷底之間強(qiáng)烈碰撞產(chǎn)生的能量耗散致使耦合作用力不斷減小。此后,隨著堵江的完成及蓄水的增加,耦合作用力緩慢增長。
圖7 耦合作用力演化過程Fig. 7 Evolution process of the coupling forces
圖8為滑坡壩幾何形態(tài)及涌浪在對(duì)岸最大爬升高度的演化過程。由圖8可知:值得注意的是,涌浪在對(duì)岸的傳播要略微晚于壩體,體現(xiàn)了滑坡體對(duì)涌浪傳播的驅(qū)動(dòng)作用。此外,從壩體沉積長度的演化過程可以發(fā)現(xiàn),壩體上游的沉積長度顯著小于壩體下游方向,上游一側(cè)壩體的延伸速度也明顯低于下游一側(cè)??梢姾恿鲗?duì)滑坡體的運(yùn)移沉積有明顯的干擾作用,即河水能顯著增強(qiáng)滑坡體順流方向的運(yùn)移,同時(shí)抑制逆流方向的運(yùn)移。
圖8 壩體形態(tài)及涌浪高度演化過程Fig. 8 Evolution process of the dam geometry and the wave height
2018年10月11日凌晨,西藏江達(dá)縣與四川省白玉縣交界處的白格村金沙江右岸邊坡(東經(jīng)98°42′17.98″,北緯31°4′56.41″)突發(fā)大規(guī)?;瑒?dòng)[25]。超過2.3×107m3的失穩(wěn)高位滑坡體高速滑動(dòng)并侵入金沙江,激發(fā)了最大高度約160 m的涌浪沿對(duì)岸快速攀爬侵蝕,在對(duì)岸留下了最大標(biāo)高達(dá)3 040 m的涌浪侵蝕區(qū)域,明顯擴(kuò)大了災(zāi)害影響范圍[26–27]。隨后,滑坡體阻塞金沙江并形成堰塞湖,形成了順河長度達(dá)1 500~2 000 m,橫河寬度約510 m的堰塞壩,并在潰決前阻水達(dá)2.9×108m3[26]。
本研究以四川省測(cè)繪地理信息局提供的2017年實(shí)測(cè)數(shù)據(jù)的10 m分辨率災(zāi)前數(shù)字高程模型和2018年10月12日高性能無人機(jī)航空攝影測(cè)量的1 m分辨率災(zāi)后數(shù)字高程模型為基礎(chǔ)建立模擬模型?!?0·11”白格滑坡堵江模擬模型如圖9所示,滑坡體被離散為2~5 m的顆粒,并通過DEM建模。CFD計(jì)算域內(nèi)最小網(wǎng)格尺寸為2.5 m,考慮到災(zāi)害發(fā)生在金沙江汛期,參考調(diào)查信息設(shè)置金沙江初始水位高程為2 890 m,上部由空氣填充,頂部為開放邊界。根據(jù)水文數(shù)據(jù),確定入流量為1 680 m3/s[28]。在上游設(shè)置質(zhì)量流入入口,并提前進(jìn)行充分時(shí)間的模擬,以在耦合計(jì)算前達(dá)到穩(wěn)定的初始河流條件。DEM參數(shù)則通過經(jīng)驗(yàn)公式、文獻(xiàn)[27,29–30]數(shù)據(jù)及參數(shù)標(biāo)定獲得,相關(guān)計(jì)算參數(shù)見表3。其中,摩擦系數(shù)可通過滑坡的高程落差及最大水平運(yùn)動(dòng)距離之比來確定,對(duì)于白格滑坡,其值可設(shè)定為0.45。此外,滑坡體的密度、泊松比及剪切模量等材料參數(shù)可參考前人研究的結(jié)果或結(jié)論。進(jìn)一步地,基于壩體沉積形態(tài)開展參數(shù)反演,從而確定模擬參數(shù)取值。
圖9 “10·11”白格滑坡堵江模擬模型Fig. 9 Numerical modelling for the “10·11” Baige landslide simulation
表3 白格滑坡案例模擬參數(shù)Tab. 3 Simulation parameters for Baige landslide event
圖10為“10·11”白格滑坡–堵江–涌浪演化過程。由圖10可見:在11 s左右,滑坡體前緣已抵達(dá)金沙江面,并有少量侵入。隨后,滑坡體高速?zèng)_擊河道,在滑坡體與河床的劇烈碰撞下,滑坡體前緣速度驟然下降,并沿對(duì)岸向上攀升。河流在滑坡體的劇烈擾動(dòng)下激發(fā)涌浪,涌浪隨后沿對(duì)岸岸坡向上傳播。后部滑坡體持續(xù)侵入河道,一方面推動(dòng)滑坡體前緣繼續(xù)向上攀升,另一方面在狹窄山谷地形及前部滑坡體的制約下,沿著河道方向向上下游擴(kuò)展。在33 s時(shí),可以觀察到河道中的滑坡體與坡面上滑坡體之間存在明顯的速度差異。此時(shí),涌浪波已經(jīng)到達(dá)對(duì)岸,并保持在滑坡體的前方。值得注意的是,在上游方向沿山谷的初始延伸明顯大于下游方向,上游側(cè)山谷夾角小于下游側(cè),因此可知,地形約束加劇了滑坡體延河谷方向的運(yùn)動(dòng)。在44 s左右,涌浪波到達(dá)對(duì)岸最大高程,并開始回落到河道中。同時(shí),可以觀察到涌浪波也沿河道向上游和下游方向傳播(圖10(e))。隨后,大部分滑坡體沉積在山谷中,形成了滑坡壩,而兩側(cè)的壩體前緣繼續(xù)沿山谷方向緩慢延伸(圖10(f))。
圖10 “10·11”白格滑坡–堵江–涌浪演化過程Fig. 10 River blockage and impulse wave evolution of the “10·11” Baige landslide
“10·11”白格滑坡災(zāi)害的滑坡速度演化過程及最大涌浪高度演化過程如圖11所示。由圖11可以發(fā)現(xiàn):其速度演化包括啟動(dòng)、加速、減速及沉積4個(gè)階段,最大平均速度為47 m/s。此外,其涌浪高度演化曲線峰值平緩,可見受地形影響,涌浪爬高在對(duì)側(cè)岸坡上不同位置先后達(dá)到峰值。最大涌浪爬高為海拔3 033 m,超過水面高程153 m,十分接近調(diào)查所得的160 m。
圖11 滑坡速度及涌浪高度演化過程Fig. 11 Evolution of the landslide velocity and the wave height
“10·11”白格滑坡事件模擬結(jié)果如圖12所示。通過觀察圖12中的累積滑坡運(yùn)移路徑及涌浪侵蝕面積可以發(fā)現(xiàn),其模擬得到的滑坡最終沉積形態(tài)整體呈弓形,順河方向沉積長度為1 759.01 m,橫河方向沉積長度為544.14 m,均與現(xiàn)場(chǎng)調(diào)查值一致。模擬得到的最大壩頂高程為海拔2 990 m,與海拔2 985 m的實(shí)際最大壩頂高程十分接近。此外,通過與實(shí)際災(zāi)害影響邊界對(duì)比可進(jìn)一步發(fā)現(xiàn),整體模擬災(zāi)害影響范圍與實(shí)地調(diào)查獲得的災(zāi)害邊界高度吻合,進(jìn)一步說明了模擬結(jié)果的可靠性。
圖12 “10·11”白格滑坡事件模型結(jié)果Fig. 12 Simulation result of the “10·11” Baige landslide event
本文提出了一種擴(kuò)展DEM–CFD耦合數(shù)值方法,并將其應(yīng)用至滑坡堵江及涌浪過程模擬和演化機(jī)制研究,進(jìn)一步模擬了“10·11”金沙江白格滑坡堵江事件。具體結(jié)論如下:
1)通過引入VOF模型及虛擬球模型,所提出的擴(kuò)展DEM–CFD耦合方法實(shí)現(xiàn)了對(duì)河流自由液面演化的追蹤,并克服了網(wǎng)格顆粒臨界尺寸的限制,實(shí)現(xiàn)了耦合作用力的精確求解,能夠滿足滑坡體粗顆粒材料的建模及高分辨率地形網(wǎng)格的要求。
2)模擬結(jié)果表明,滑坡壩形成過程涉及強(qiáng)烈的滑坡–河流相互作用現(xiàn)象?;麦w驅(qū)動(dòng)涌浪在對(duì)岸爬升及河道內(nèi)傳播;而河流顯著增加了滑坡體的動(dòng)能耗散,且在上下游方向上對(duì)滑坡體的沉積產(chǎn)生相反的影響。
3)所提出的擴(kuò)展DEM–CFD耦合數(shù)值方法準(zhǔn)確再現(xiàn)了“10·11”金沙江白格滑坡堵江及涌浪演化過程,模擬得到的滑坡運(yùn)移路徑、沉積形態(tài)及涌浪侵蝕面積與實(shí)地調(diào)查結(jié)果保持了較高的一致性,從而進(jìn)一步證明了方法的有效性。