基于SPH-DEM并行計(jì)算的滑坡涌浪模擬

0 引言

滑坡災(zāi)害若發(fā)生于水庫、湖泊、河道等水體附近，會激起巨大涌浪并導(dǎo)致次生災(zāi)害。涌浪不僅直接損害鄰近的碼頭、橋梁、公路等重要基礎(chǔ)設(shè)施，其傳播過程中引發(fā)的水面驟變還會將岸邊船只、行人卷入水中，造成的后果往往遠(yuǎn)超過滑坡本身[1]。尤其當(dāng)滑坡體接近大壩時(shí)，涌浪可能翻越大壩并對項(xiàng)體結(jié)構(gòu)造成破壞。因此，研究水庫滑坡及其誘發(fā)的涌浪對于預(yù)測和評估災(zāi)害范圍具有重要意義。

對滑坡涌浪的常見研究方法有模型實(shí)驗(yàn)法、經(jīng)驗(yàn)公式法和數(shù)值模擬法[2]。模型實(shí)驗(yàn)的結(jié)論與規(guī)律一般適用于特定對象或案例，如Harbitz等[3]采用三維縮尺模型實(shí)驗(yàn)對挪威西部多個(gè)不穩(wěn)定巖質(zhì)滑坡涌浪的產(chǎn)生、傳播和爬坡過程進(jìn)行了研究;黃波林等[4]通過白鶴灘庫區(qū)王家山滑坡涌浪縮尺物理模型，對低 Fr 數(shù)涉水滑坡涌浪的產(chǎn)生、傳播特征進(jìn)行了試驗(yàn)分析；黃錦林等5使用縮尺模型實(shí)驗(yàn)對樂昌峽水庫的滑坡涌浪風(fēng)險(xiǎn)進(jìn)行了評估。

經(jīng)驗(yàn)公式法方面，近些年相關(guān)學(xué)者針對不同工程展開分析，也總結(jié)了眾多經(jīng)驗(yàn)公式。如任坤杰[通過量綱分析與回歸分析給出了整體滑坡體的首浪高度經(jīng)驗(yàn)公式；殷坤龍等推導(dǎo)出適合于三峽庫區(qū)滑坡的涌浪計(jì)算系列公式；馬倩等8根據(jù)試驗(yàn)結(jié)果，結(jié)合經(jīng)典波浪爬坡經(jīng)驗(yàn)公式推導(dǎo)出了新滑坡涌浪爬高公式；Heller等[在對前人實(shí)驗(yàn)數(shù)據(jù)綜合分析的基礎(chǔ)上提出了波浪產(chǎn)生參數(shù) P 以綜合表示涌浪源特征等。數(shù)值模擬法一般采用N-S、SWE等作為控制方程進(jìn)行求解，以獲得滑坡涌浪物理場在時(shí)間尺度上的變化規(guī)律。由于數(shù)值方法和流體動力理論基本完備，計(jì)算結(jié)果準(zhǔn)確且可視化效果好，因此數(shù)值模擬法成為目前主流的研究方法。目前報(bào)道的工程案例多采用成熟的商業(yè)軟件，如鄧成進(jìn)等[1°]通過FLOW3D軟件研究了網(wǎng)格尺寸對庫區(qū)滑坡涌浪傳播計(jì)算精度的影響;張萬舉等[11]基于FLUENT軟件模擬矩形滑塊在“V”“U”形庫區(qū)斷面下的滑坡涌浪傳播過程；彭輝等[12]基于FLOW3D模擬了彎曲河道型庫岸滑坡涌浪傳播過程及其與下游水利設(shè)施大壩相互作用。雖然這些軟件在某些方面帶來便利和效率，但在關(guān)鍵技術(shù)上可能存在“卡脖子”的風(fēng)險(xiǎn)，不利于國產(chǎn)軟件的創(chuàng)新和發(fā)展。

鑒于此，本文采用自主研發(fā)的分布式并行SPH-DEM程序構(gòu)建計(jì)算模型，以白日埡滑坡為案例，對滑坡涌浪災(zāi)害的動力學(xué)過程進(jìn)行分析，旨在解決滑坡變形失穩(wěn)及涌浪形成與傳播問題的同時(shí)，系統(tǒng)分析涌浪的產(chǎn)生與傳播規(guī)律，以期為水庫滑坡－涌浪復(fù)合災(zāi)害的分析研究提供技術(shù)參考。

1數(shù)值模型理論方法

1.1 流體動力學(xué)控制方程

dρ/dt+ρ??u=0

式中： 為流體粒子速度矢量， ρ 為流體顆粒的密度， 為壓力， 為外力， u 是運(yùn)動黏滯系數(shù)。本文采用投影法對控制方程進(jìn)行離散。因?yàn)樗俣葓隹梢苑纸鉃闊o散分量和無旋分量的梯度，而壓力場為非旋標(biāo)量場，因此速度場被分解為兩部分： ① 由流動黏度和外力共同作用產(chǎn)生的中間速度; ② 由壓力梯度驅(qū)動的壓力速度。將兩部分速度分別計(jì)算并疊加后，即獲得粒子的速度場分布。

1. 2 光滑粒子流體動力學(xué)（SPH）

SPH是一種近似方法，其本質(zhì)是將場函數(shù)及其導(dǎo)數(shù)的積分表示為支持域中所有相鄰粒子相應(yīng)數(shù)值的求和，即函數(shù) f（x） 可以表示為

式中： W 是平滑核函數(shù)， h 是平滑長度， x_j∣，α_D=21/（16πh³） 。 x 代表粒子位置矢量， 表示粒子的微元體積。對公式（3）分別取散度和梯度，可以得到SPH的表達(dá)式為

式中： m_i 為粒子 i 的質(zhì)量， W_ij=W（x_i-x_j，h） 。

將公式（5）代人式（6）中，可得拉普拉斯算子 Δ= abla?（?f） 的表達(dá)形式為

可以證明，公式（7）能保證向量場在投影中的散度為零，即為一種精確的投影格式。

1.3 離散單元模型（DEM）

假定固體粒子都是質(zhì)量均勻的標(biāo)準(zhǔn)球體，DEM通過計(jì)算接觸顆粒之間的力和接觸力矩，更新顆粒位置和速度，其控制方程表示為

式中： r_i 是粒子 i 的半徑， I_i 是粒子 i 的慣性矩； 是粒子 χ_i 的角速度， 是與 i 粒子作用邊界的外法線單位向量方向， 是從 i 粒子的位置指向 j 粒子位置的單位向量， 是粒子的位置矢量， 是粒子 i 的速度， F_：j 是粒子 j 對粒子 χ_i 的合力， 是邊界 對粒子 i 的合力， 是邊界給粒子 i 的切向力， 是粒子 i 和 j 之間的切向力。關(guān)于切向力具體計(jì)算方法可參考各種摩擦模型[13-14]，這里不再贅述。

1.4 SPH-DEM耦合算法

流固耦合問題的核心在于計(jì)算流體與固體間的相互作用力。固體粒子所受的浮力、拖拽力、掏蝕力等，本質(zhì)上源于流體的壓力和黏性力。因此，可將DEM粒子視作流體粒子，納入與其相鄰流體顆粒的N-S方程中。通過SPH求和計(jì)算，對作用于固體顆粒表面的作用力進(jìn)行積分，以求得流體對固體的作用力；再疊加固體粒子自身的重力，即可獲取其受力狀態(tài)和運(yùn)動軌跡。根據(jù)牛頓第三運(yùn)動定律，流體粒子會受到固體粒子的反作用力，以保持系統(tǒng)動量守恒，這一守恒性在N-S方程的計(jì)算中已得到體現(xiàn)。利用方程（6）計(jì)算動量方程（2）中的壓力 ，即固體粒子 i 所受周圍流體粒子 f 的作用力 F_if

另外，固體粒子 i 受到周圍流體粒子 f 對其產(chǎn)生的黏性力 可采用方程（12）計(jì)算：

式中： ，由此，DEM運(yùn)動微分方程（8）改寫為

1.5 分布式MPI并行算法

負(fù)載平衡是并行系統(tǒng)充分發(fā)揮效能的關(guān)鍵。分布式MPI并行的關(guān)鍵是合理地將計(jì)算負(fù)荷分配到每個(gè)獨(dú)立的處理器或者進(jìn)程中。理想情況下，當(dāng)進(jìn)程數(shù)增加時(shí)，并行計(jì)算時(shí)間將線性減少。對于SPH和DEM算法而言，負(fù)載平衡的關(guān)鍵在于各進(jìn)程的粒子數(shù)大致相同，因此本次研究采取的策略是將計(jì)算域分解為幾個(gè)子區(qū)域，每個(gè)子區(qū)域都對應(yīng)一個(gè)進(jìn)程，每個(gè)進(jìn)程中包含的粒子數(shù)基本一致。

無網(wǎng)格粒子法并行計(jì)算的一個(gè)難點(diǎn)是粒子在計(jì)算域內(nèi)隨時(shí)間發(fā)生自由遷移，導(dǎo)致進(jìn)程間的粒子數(shù)量發(fā)生變化，因此要在計(jì)算中動態(tài)調(diào)整計(jì)算域，以保證進(jìn)程間計(jì)算負(fù)荷的平衡。本文采用的動態(tài)區(qū)域調(diào)整算法如圖1所示。為了將所有粒子均勻地分配到它們對應(yīng)的進(jìn)程，首先通過二分法將粒子分為兩部分；然后反復(fù)劃分每一部分，直到所有進(jìn)程接收到相同數(shù)量的粒子；最后整個(gè)計(jì)算域被均勻地劃分為 1～n 的子區(qū)域。每個(gè)子區(qū)域都對應(yīng)一個(gè)粒子塊，這些粒子塊動態(tài)存儲在相應(yīng)進(jìn)程中，可以方便增加或者刪除粒子，為后續(xù)計(jì)算做準(zhǔn)備。

粒子相互作用是SPH或DEM計(jì)算中最耗時(shí)的模塊。此模塊通過插值鄰近粒子來計(jì)算每個(gè)粒子的受力。由于計(jì)算復(fù)雜性，此模塊占總執(zhí)行時(shí)間的 90% 以上。為了實(shí)現(xiàn)有效的加速效果，需要解決粒子相互作用模塊的并行可擴(kuò)展性問題。基于此，本文提出一種重疊計(jì)算和通信的策略，通過隱藏通信延遲，提高粒子相互作用模塊的并行效率（圖2）。

由于SPH只計(jì)算平滑半徑內(nèi)的粒子，同時(shí)DEM相鄰粒子之間距離最大不超過 h ，因此邊界距離為 2h ，將粒子劃分為邊界幽靈粒子和中心粒子，可以保證中心粒子及其相鄰粒子完全落入某個(gè)獨(dú)立進(jìn)程中。并行計(jì)算可以大致分為兩個(gè)階段。以進(jìn)程2為例，由于中心粒子不需要與相鄰進(jìn)程通信，因此第一階段首先完成中心粒子相互作用力的計(jì)算。為了最小化通信開銷，在計(jì)算的同時(shí)，通過異步傳輸與其他進(jìn)程聯(lián)系，即由MPI_Irecv和MPI_Isend接收和發(fā)送幽靈粒子，當(dāng)中心粒子計(jì)算完成時(shí)，邊界粒子也被成功發(fā)送和接收。第二階段主要是更新計(jì)算落入邊界部分的幽靈粒子。這種策略充分利用了幽靈粒子傳輸導(dǎo)致的空閑等待計(jì)算資源，有效提高了程序的并行效率。

2 算例驗(yàn)證

Heinrich[15]曾在實(shí)驗(yàn)室（圖3）對滑坡涌浪開展實(shí)驗(yàn)，滑坡過程用一個(gè)剛性楔形體從 45^° 的斜面上自由下滑來模擬，楔形體的密度為 2000kg/m³ ，上表面初始位置與水面齊平

將此算例的 SPH-DEM計(jì)算結(jié)果與SPH-SPH，SPH－DDA以及物理實(shí)驗(yàn)結(jié)果進(jìn)行比較（圖4）。其中SPH-SPH是基于開源軟件DualSPHysics計(jì)算的結(jié)果，而SPH-DDA是Wang等[1的研究結(jié)果。

可以看出，本文采用的SPH-DEM模型計(jì)算精度最高，流體受到楔形體擠壓后產(chǎn)生的液面變形非常接近物理實(shí)驗(yàn)結(jié)果。特別是在 t=1.0 s時(shí)波谷和波峰的形成，規(guī)律和物理實(shí)驗(yàn)完全一致，這從側(cè)面也證實(shí)了液體受到楔形體的反作用力大小是合理的，即利用本文SPH-DEM方法描述流體在重力作用下產(chǎn)生的波動現(xiàn)象及其相互影響具有良好精度，可以應(yīng)用于滑坡涌浪模擬以及分析滑坡入水后涌浪傳播特性。

3白日埡滑坡涌浪數(shù)值模型驗(yàn)證

白日埡崩塌體所在庫區(qū)位置為深“V”形河谷，平均水深達(dá) 90m ，兩岸距離約 700m 。滑坡產(chǎn)生的涌浪高度受滑坡體積、形狀、滑速以及水深等約束明顯。涌浪在爬高過程中也會受到河道地形、涌浪形態(tài)及岸坡條件等影響。涌浪的產(chǎn)生、爬升方向與涌浪在河道中的傳播方向幾乎垂直，能量傳播和耗散隨著傳播距離和深度的變化而變化。這些特征使得白日埡滑坡涌浪呈現(xiàn)典型的三維特征，使用三維模擬能比較準(zhǔn)確地捕捉滑坡涌浪的全過程。

3.1 計(jì)算模型

白日埡崩塌體滑坡建模采用SRTM地形數(shù)據(jù)，具有較高的垂直精度，分辨率為1arc－second（約30m）。計(jì)算區(qū)域如圖5所示，庫區(qū)總長度約 7km ，其中上游 3km ，下游 4km 。水庫正常蓄水位 470.00m ，按照SPH影響域長度 h=5m 來填充水粒子，共產(chǎn)生約400萬顆粒子。水體湍流模型采用的是較為常用的 k 一 ε 模型，可以描述涌浪復(fù)雜運(yùn)動中的能量耗散。

根據(jù)監(jiān)測資料分析，白日埡堆積體前緣受庫水浸泡軟化作用和地下水作用，產(chǎn)生指向河床方向且偏向上游的滑動變形，在正常蓄水位附近安全系數(shù)最小（約為1.05）。因此，將滑面設(shè)定為與河床基本垂直，略微偏向上游的方向。局部滑動的滑面設(shè)定為橢球面，潛在滑體體積約175萬 m³ 。使用半徑為 2.5m 的球形離散單元顆粒，共2.7萬顆。固體顆粒單元的密度和泊松比分別取 2 650kg/m³ 和0.32，其余參數(shù)如剪切模量、滾動摩擦因數(shù)和阻尼系數(shù)等，通過不斷重復(fù)的迭代試算來確定。

3.2 三維滑坡涌浪過程分析

圖6顯示了在 0～120s 內(nèi)，白日埡滑坡堆積體滑坡和涌浪過程。滑坡運(yùn)動大概持續(xù)了 100s ，由此產(chǎn)生的涌浪約在150s傳播至大壩。在此期間，大約形成了3次涌浪過程：第一次最大涌浪高度約 30m ，第二次 20m ，第三次 15m ，每次間隔約 40s 。

在5s左右時(shí)，滑坡前緣抵達(dá)水邊，此時(shí)碎屑體的速度約為 10m/s ，滑坡碎石剛剛接觸到水面，因此還未形成明顯的擠壓效果，水面形態(tài)也沒有顯著變化。在10s時(shí)，在滑坡顆粒沖擊作用下，動能持續(xù)傳給接觸到的水體，由于未與滑體接觸的水體速度很小，水體前后的速度差導(dǎo)致水面產(chǎn)生擠壓，浪頭開始逐漸形成，首浪最高處約 35m ，入水處形成了明顯“水坑”。

在20s時(shí)，隨著大量的滑坡碎石快速入水，沖擊波導(dǎo)致水面呈現(xiàn)典型的圓弧狀向外擴(kuò)散，攜帶大量動能向?qū)Π兑约八闹軅鞑ィ胨帯八印鄙疃戎饾u增大。在35s時(shí)，涌浪第一次到達(dá)對岸并開始爬升，在此過程中，一部分動能耗散在與岸坡的摩擦及相互作用中，另一部分則轉(zhuǎn)換為重力勢能。當(dāng)涌浪的浪尖到達(dá)海拔約 500m 處，即最大爬高 30m 時(shí)，涌浪速度幾乎降為零，隨即重力勢能轉(zhuǎn)換為動能，水體開始向反方向運(yùn)動。

在40s時(shí)，涌浪從對岸高處加速滑落后反向傳播，并向上下游傳播。在50s時(shí)，反向波峰到達(dá)滑坡近岸處，與河道中正在傳播的波浪相互沖擊，兩個(gè)波的峰值相互重疊，導(dǎo)致局部浪高的顯著增加，形成強(qiáng)烈的沖擊，產(chǎn)生所謂的“瘋狗浪”，對周圍船只和岸坡的破壞非常大。

在 60s 時(shí)，涌浪在滑坡點(diǎn)附近岸坡爬高 20m ，達(dá)到高程 490.00m 處；70s時(shí)，第二波波峰到達(dá)對岸，涌浪爬高約 20m;80s 時(shí)，第二波涌浪開始消退，再次反向傳播；90s時(shí)，涌浪再次到達(dá)滑坡附近岸坡并開始爬高，此刻由于大部分滑體已經(jīng)入水，涌浪產(chǎn)生的疊加效益已不明顯。

3.3 涌浪對下游河道及大壩的影響

在滑坡區(qū)與壩前，沿下游河道中泓線選取了5個(gè)代表性監(jiān)測點(diǎn)進(jìn)行涌浪高度監(jiān)測，具體位置如圖7所示。其中A點(diǎn)和C點(diǎn)在河道轉(zhuǎn)彎處，B點(diǎn)和D點(diǎn)為直線河道的中點(diǎn)，E點(diǎn)為壩前監(jiān)測點(diǎn)。

各點(diǎn)水位隨時(shí)間變化關(guān)系如圖8所示。可以看出，離滑坡體較近的監(jiān)測點(diǎn)A其最大涌浪高度只有4m，遠(yuǎn)低于初始浪高 35m ，說明滑坡產(chǎn)生的能量主要沿滑動方向傳遞，沿河道方向較少。另外，涌浪波在傳播初期包含了較寬頻率范圍的分量，其中高頻分量具有較短的波長和較高的頻率。例如A點(diǎn)在不到5s時(shí)有一個(gè)小的波峰，即為高頻波形成。隨著傳播距離的增加，高頻分量由于衰減較快，其能量比例逐漸降低，但其存在使得涌浪波的波面變得復(fù)雜和不規(guī)則。涌浪波低頻分量的波動能量衰減得相對較慢，能夠保持較長時(shí)間的穩(wěn)定傳播，逐漸占據(jù)主導(dǎo)地位。這種能量分布的變化使得涌浪波在遠(yuǎn)距離傳播時(shí)，更多地表現(xiàn)出低頻波的特性。正如圖中各監(jiān)測點(diǎn)出現(xiàn)的典型波峰波谷，主要為低頻波。

白日埡滑坡失穩(wěn)入水后，涌浪傳播到各監(jiān)測點(diǎn)的時(shí)間分別為第40，55，80，105秒和第130秒。由此可以確定AB段傳播時(shí)間15s，波浪傳播速度 42m/s ;BC段傳播時(shí)間 25s ，速度 25m/s ;CD段傳播時(shí)間 25s ，速度 40m/s ;DE段傳播時(shí)間 25s ，速度 40m/s 。可以看出，滑坡涌浪在河流彎道處C的傳播速度較慢，在直道段傳播速度較快。這是因?yàn)樵谄街焙佣危坷瞬ǖ牟ㄩL較長，周期也較長，長周期的涌浪儲存更多能量，移動速度更快。

另外，江坪河水電站大壩壩頂高程 476.00m ，在正常蓄水位 470.00m 下，此案例計(jì)算的滑坡涌浪在壩前E點(diǎn)最大壅高水位為 472.80m ，不會產(chǎn)生涌浪翻壩或漫頂?shù)那樾巍?/p>

為分析壩址處地形受到涌浪的影響，在江坪河大壩左岸、壩中和右岸分別選取3個(gè)水位觀測點(diǎn)展開分析。由于左岸的坡度要比右岸大，因此左岸受到涌浪影響更加明顯，產(chǎn)生的壅高會更大，3處的水位波動過程如圖9所示。由于水體重力作用，靠近山體兩側(cè)的水位有波動的現(xiàn)象，特別是左岸水位，周期性比較明顯。大壩中點(diǎn)處的水位由于受到兩側(cè)山體邊界的影響較小，大部分時(shí)間水位比較一致。可以看出，3處水位在130s的時(shí)候有比較明顯的抬升，表明潰壩涌浪已經(jīng)傳至壩前，由于左岸比右岸要陡峭，因此水位波動更加明顯，最高處達(dá)到 474.00m 。但3處水面壅高都未超過江坪河大壩的壩頂高程，不會產(chǎn)生翻壩涌浪，大壩整體安全。

4并行效率

并行可擴(kuò)展性通常用于評估大規(guī)模程序的并行性能。在總體計(jì)算規(guī)模保持不變的條件下，具有優(yōu)越并行效率的程序其計(jì)算時(shí)間將隨著計(jì)算節(jié)點(diǎn)或進(jìn)程的增加而線性減少，這一規(guī)則也被稱為程序的強(qiáng)可擴(kuò)展性。本文評估了在128進(jìn)程上模擬400萬粒子的并行效率，如圖10所示。

可以看出，當(dāng)進(jìn)程數(shù)量小于20時(shí)，可以獲得更好的線性加速比，此時(shí)并行效率接近 80% 。但是，當(dāng)進(jìn)程數(shù)量增加到128時(shí)，程序的并行效率下降到 70% 。出現(xiàn)這種現(xiàn)象的原因是因?yàn)檫吔缌Ｗ优c中心粒子的比例過大。假設(shè)該比例為 τ ，隨著進(jìn)程數(shù)量的增加，由于中心粒子的數(shù)量隨著進(jìn)程數(shù)量的增加而減少， τ 變得越來越大。在理想情況下，存在一個(gè)閾值 τ₀ ，它表示當(dāng)邊界粒子完成傳輸時(shí)，中心粒子剛好完成計(jì)算的臨界比例。顯然，為了獲得線性加速比，根據(jù)通信隱藏策略，在粒子相互作用和系統(tǒng)更新階段，邊界幽靈粒子的傳輸總是需要在中心粒子的計(jì)算完成之前結(jié)束。這樣CPU才能始終連續(xù)執(zhí)行粒子的計(jì)算任務(wù)，這是確保并行程序獲得線性加速比的關(guān)鍵。

5結(jié)語

針對滑坡涌浪問題，本文首先介紹了一種基于SPH-DEM的耦合模型，通過將固體粒子視為流體粒子來計(jì)算流體與固體之間的作用力，然后根據(jù)粒子法的計(jì)算特性開發(fā)了分布式并行程序，通過采用高效的區(qū)域分解算法和重疊粒子計(jì)算與信息交互策略，提高了程序并行效率。將該程序應(yīng)用于白日埡崩塌體誘發(fā)的涌浪模擬分析，計(jì)算結(jié)果表明，在垂直河道方向，滑坡涌浪最大爬高約 30m ，此后由于能量耗散，爬高逐次下降。在順河道方向，涌浪傳播過程受河道地形的影響較大，浪高衰減很快，在壩前最大壅高水位為474.00m ，不會產(chǎn)生涌浪翻壩或漫頂風(fēng)險(xiǎn)。

參考文獻(xiàn)：

[1］殷坤龍，張宇，汪洋.水庫滑坡涌浪風(fēng)險(xiǎn)研究現(xiàn)狀和災(zāi)害鏈風(fēng)險(xiǎn)管控實(shí)踐[J].地質(zhì)科技通報(bào)，2022，41（2）：1-12.

[2] 陳世壯，徐衛(wèi)亞，石安池，等.高壩大庫滑坡涌浪災(zāi)害鏈研究綜述[J].水利水電科技進(jìn)展，2023，43（3）：83-93.

[3] HARBITZCB，GLIMSDALS，LσVHOLTF，etal.Rockslidetsunamisincomplex fjords：from anunstable rock slope at Akernesetto tsunamiriskinwesternNorway[J].CoastalEngineering，2014，88：101-122.

[4] 黃波林，胡劉洋，李仁江，等.水庫區(qū)低Froude數(shù)的典型涉水滑坡涌浪縮尺物理模型試驗(yàn)研究：以王家山滑坡為例[J」.巖石力學(xué)與工程學(xué)報(bào)，2023，42（8）：1989-1909.

[5]黃錦林，練繼建，張婷.滑坡涌浪作用下樂昌峽大壩安全評估［J].水利水電技術(shù)，2013，44（11）：93-97.

[6] 任坤杰，韓繼斌，陸虹.滑坡涌浪首浪高度試驗(yàn)研究[J].人民長江，2012，43（2）：43-45.

[7]殷坤龍，劉藝梁，汪洋，等.三峽水庫庫岸滑坡涌浪物理模型試驗(yàn)[J].地球科學(xué)：中國地質(zhì)大學(xué)學(xué)報(bào)，2012，37（5）：1067-1074.

[8] 馬倩，薛宏程，任昱，等.水庫滑坡涌浪爬坡特征試驗(yàn)研究[J].人民長江，2020，51（2）：195-199.

[9] HELLERV，HAGERWH.Wave types of landslide generated impulsewaves[J].OceanEngineering，2011，38（4）：630-640.

[10］鄧成進(jìn)，黨發(fā)寧，陳興周，等.庫區(qū)滑坡涌浪三維數(shù)值模擬分析[J].水利水運(yùn)工程學(xué)報(bào)，2020（6）：64-71.

[11]張萬舉，閆毅志，李梓萌，等.V、U形庫區(qū)斷面下滑坡涌浪傳播的數(shù)值模擬研究[J].水電能源科學(xué)，2022，40（8）：101-104.

[12］彭輝，黃亞杰.彎曲河道型庫岸滑坡涌浪傳播及其與大壩的相互作用[J].長江科學(xué)院院報(bào)，2022，39（10）：66-71.

[13]唐岳灝，姜清輝.SPH-DEM流固耦合方法模擬金沙江白格滑坡應(yīng)用研究[J].長江工程職業(yè)技術(shù)學(xué)院學(xué)報(bào)，2023，40（2）：1-8.

[14]唐岳灝，姜清輝.金沙江白格滑坡殘留體穩(wěn)定性及堵江風(fēng)險(xiǎn)分析[J].水利水電快報(bào)，2023，44（5）：38-44.

[15]HEINRICH P.Nonlinear water waves generated by submarine and aer-iallandslides[J].Journal ofWaterwayPortCoastal amp; OceanEngi-neering，1992，118（3）：249-266.

[16]WANGW，CHEN GQ，ZHANG H，et al. Analysis of landslide-gen-erated impulsive waves using a coupled DDA -SPH method[J].Engi-neeringAnalysiswithBoundaryElements，2016，64：267-279.

（編輯：胡旭東）

Landslide surge simulation based on SPH -DEM parallel computing

TANG Yuehao1 ，JIANG Qinghui2 （1.ScholofWaterEnginering，ChangjangstituteofTecnlogyWuhan4012，ina；.olofCilEngieing，Wu han University，Wuhan 430072，China）

Abstract：To solve the strong coupling between landslide body and water，and thecomplexfluid-solidcoupling problem causedbytheevolutionofnonlinear waterwaves，thesmothed particleflowmethod（SPH）wasused tosimulatethehydrodynamic characteristicsof water，andthediscreteelement method（DEM）wasused tosimulatethelargedeformationanddiscontinuity of landslides.Through constructing SPH-DEMfluid-solidcouplingalgorithm，anddevelopingaparalelcomputingprogram basedon distributed MPI，the large-scale high-performancecalculationof the SPH-DEMcoupling processwasrealized.This program wasapliedtothesimulation of Bairiya landslidesurge in JiangpingheReservoirof Hubei Province.The disasterchain process of landslidedevelopment，surge formation，propagationalongriverchannelandimpactonthedam wasreproducedTheresultsshowedthatthe maximumsurgeheight infrontofthedam was474.OOm，whichwaslowerthanthedam elevationof476.00 m.This program achieved more than 70% parallel acceleration efficiency when the processes number was 128.The research resultscan efectivelyanalyzethedynamic processoflandslidesurgedisaster，andthedeveloped paralel programcanprovidereference for other domestic similar software.

Key words：landslide surge；fluid-solid coupling；discrete element method（DEM）；smoothed particleflow method（SPH）：Bairiya landslide； Jiangpinghe Reservoir