摘要：滑坡涌浪鏈生災(zāi)害的演化過程涉及滑坡體-河流之間復(fù)雜的相互作用，流固耦合方法在顆粒尺度下的應(yīng)用正日益凸顯。為了深入研究這一問題，研究采用顆粒尺度的CFD（Computational Fluid Dynamics）-DEM（Discrete Element Method）流固耦合方法。通過編寫UDF，以改進計算流體力學(xué)軟件FLUENT中離散相模型只能進行動量交換而無法考慮排水效應(yīng)的缺點。這不僅改善了滑坡體與流體間的體積置換問題，而且克服了網(wǎng)格顆粒臨界尺寸比的限制。在顆粒堆積體坍塌-涌浪試驗?zāi)M中，數(shù)值模擬計算得到的最大涌浪高度為8.67 cm，與物理試驗結(jié)果較為一致，證明了改進后的CFD-DEM流固耦合模型能有效還原顆粒運動狀態(tài)及涌浪的波動變化，能夠滿足滑坡涌浪過程模擬中滑坡體粗顆粒材料建模和高分辨率真實地形網(wǎng)格的需求。

關(guān)鍵詞：滑坡涌浪；流固耦合；CFD-DEM耦合方法；數(shù)值模擬

中圖分類號：TV642.4文獻標識碼：A文章編號：1001-9235（2024）07-0111-07

郝悅彤，肖鴻.顆粒離散相模型的改進及其在滑坡涌浪數(shù)值模擬應(yīng)用驗證［J］.人民珠江，2024，45（7）：111-117.

Improvement of Discrete Phase Model for Particles and Its Validation for Application in Numerical Simulation of Landslide Surges

HAO Yuetong，XIAO Hong*

（State Key Laboratory of Hydraulics and Mountain River Engineering，Sichuan University，Chengdu 610065，China）

Abstract：The evolution of landslide surge chain-generated hazards involves complex landslide-stream interactions，and the application of fluid-solid coupling methods at the granular scale is becoming more and more prominent.In order to study this problem in depth，this study adopts the CFD（Computational Fluid Dynamics）-DEM（Discrete Element Method）fluid-solid coupling method at the granular scale.The UDF is written to improve the disadvantage that the discrete phase model in the Computational Fluid Dynamics software FLUENT can only perform momentum exchange but cannot consider the drainage effect.This not only improves the volume replacement problem between landslides and fluids，but also overcomes the limitation of the critical size ratio of mesh particles.In the simulation of particle accumulation body collapse-surge test，the maximum surge height obtained by numerical simulation is 8.67 cm，which is more consistent with the physical test results，proving that the improved CFD-DEM fluid-solid coupling model proposed in this paper can effectively restore the particle motion state and the fluctuation change of the surge，and it can meet the needs of the coarse-grained material modeling of landslides in the simulation of the landslide surge process and the high-resolution real topography grid in the simulation of landslide surge process.

Keywords：landslide surge;fluid-solid coupling;CFD-DEM coupling method;numerical simulation

庫岸邊坡由于庫水位波動、暴雨、地震以及人為施工活動等影響容易發(fā)生滑坡災(zāi)害［1］，滑坡涌浪作為滑坡體入水后的次生災(zāi)害，所產(chǎn)生的破壞性和影響范圍都超越了滑坡災(zāi)害本身［2］。大規(guī)?；麦w高速入水激起大量水體，水體回落形成涌浪在河道、庫區(qū)等水域內(nèi)傳播。一方面，涌浪的沖擊會對水工建筑物造成破壞甚至產(chǎn)生漫壩或潰壩的風(fēng)險，對下游人民生命財產(chǎn)安全形成威脅；另一方面，滑坡體入水最終堆積在河道形成堰塞體會進一步增加潰壩漫壩的風(fēng)險［3］。

隨著計算機硬件水平的不斷提高，數(shù)值分析方法在研究庫岸災(zāi)害方面的應(yīng)用變得愈發(fā)重要［4］。河床沖刷、泥石流等諸多巖土工程問題往往涉及復(fù)雜的流體-顆粒耦合作用，在對其機理進行的研究中，顆粒尺度的流固耦合方法正發(fā)揮著越來越重要的作用［5-6］。利用離散元方法（Discrete Element Method，DEM）構(gòu)建巖土體，結(jié)合常用于流體流動建模的計算流體動力學(xué)（Computational Fluid Dynamics，CFD）方法，Tsuji等［7］提出了一種顆粒尺度的流固耦合方法，即CFD-DEM耦合方法，這種方法成功用于模擬流化床問題。在CFD-DEM耦合方法中，流體運動由Euler框架下的連續(xù)方法描述，固體顆粒在考慮粒間碰撞以及流體-顆粒作用的基礎(chǔ)上由Lagrange框架下的離散方法描述其微觀特性，在空間和時間尺度上都遠小于連續(xù)方法描述的宏觀特性，需要建立宏觀尺度（連續(xù)介質(zhì)）和微觀尺度（離散介質(zhì)）的橋梁，也就是流體體積分數(shù)［8］。目前，CFD軟件FLUENT中的離散相模型對顆粒處理尚不完善，僅考慮了動量交換而未考慮顆粒與流體之間的體積置換，對于滑坡涌浪問題的模擬計算會產(chǎn)生較大誤差。此外，傳統(tǒng)的CFD-DEM耦合方法的使用受限于網(wǎng)格顆粒臨界尺寸比［9-10］，難以滿足滑坡涌浪過程模擬中滑坡體粗顆粒材料建模和高分辨率真實地形網(wǎng)格的需求［11］。

本文旨在優(yōu)化FLUENT中CFD-DEM耦合數(shù)值模擬方法，并驗證其進行滑坡涌浪數(shù)值模擬的可靠性。首先引入虛擬多孔介質(zhì)球孔隙率計算模型，不僅考慮到各相間體積置換，而且可以解決實際應(yīng)用過程中顆粒尺寸大于網(wǎng)格尺寸時產(chǎn)生計算不穩(wěn)定的缺點；其次通過單顆粒沉降模擬、顆粒團入水及出水模擬和顆粒堆積體坍塌-涌浪試驗?zāi)M，驗證了模型的準確性和有效性。研究成果可以為相關(guān)地區(qū)重大水電工程開發(fā)運行的安全保障提供參考。

1流固耦合數(shù)值方法

1.1計算流體運動學(xué)方法

1.1.1控制方程

對于流體相，采用以下基本方程：連續(xù)性方程見（1）：

+++=0（1）

動量方程見（2）：

式中：u、v、w分別為x、y、z 3個方向上的速度分量；ρ為密度；t為時間；μ為動力黏度；p為壓強；fx、fy、fz分別為流體單元在3個方向上所受的體積力，一般考慮為重力；fpfx、fpfy、fpfz分別為流體所受顆粒-流體相互作用在3個方向上的大小。

1.1.2 VOF方法

使用VOF（Volume of Fluid）方法捕捉自由面，對涌浪的形成及演變過程中多相流體交界面進行追蹤。通過針對多相的體積分數(shù)的連續(xù)性方程求解來完成對相之間的界面跟蹤，對于第q相流體，

式中：αq為第q相流體在某個計算單元的體積分數(shù)值；pq為從第p相到第q相的質(zhì)量輸移；qp為從第q相到第p相的質(zhì)量輸移，默認右側(cè)源項sα為0。

當αq=0時，代表該單元格沒有該種流體；當αq=1，代表該單元格充滿該種流體；當0lt;αq＜1，代表該單元格為該流體和其他流體的交界面。根據(jù)單元格上的αq值，分配計算域中每個單元格的屬性及其存儲的變量，見式（4）、（5）。

式中：腳標對應(yīng)不同種類的流體。

1.2離散元方法

對于顆粒相，根據(jù)牛頓第二定律，每個粒子切向的運動方程為式（6）：

mi=（Fn，k，i+Ft，k，i）+mig+Fi（6）

式中：mi為顆粒i的質(zhì)量；Fn，k，i為相鄰顆粒k產(chǎn)生的法相接觸力；Ft，k，i為相鄰顆粒k產(chǎn)生的切向接觸力；Fi為流體作用在顆粒i上的力。

單個顆粒的旋轉(zhuǎn)運動方程見式（7）：

Ii=aεijlniFt，k，l（7）

式中：Ii為顆粒i的轉(zhuǎn)動慣量；ωi為顆粒i的轉(zhuǎn)動角速度；a為顆粒的半徑；ni為顆粒i表面的單位法向量。

1.3虛擬多孔介質(zhì)球孔隙率模型

計算網(wǎng)格孔隙率是計算顆粒相和流體相之間相互作用力的關(guān)鍵［12］。傳統(tǒng)CFD-DEM使用的方法有中心孔隙率計算方法（Centre Void Fraction Method）和劃分孔隙分數(shù)法（Divided Void Fraction Method）。前者誤差較大且顆粒體積大于網(wǎng)格體積時會出現(xiàn)流場不連續(xù)的情況，導(dǎo)致整個CFD-DEM模擬計算的數(shù)值不穩(wěn)定；后者雖能得到相對光滑的體積分數(shù)場，但在非均勻網(wǎng)格中較難實現(xiàn)，同時也無法解決顆粒體積大于網(wǎng)格體積時的數(shù)值不穩(wěn)定問題。在應(yīng)用過程中需要滿足網(wǎng)格尺寸與顆粒粒徑尺寸比大于4的限制［8］。

本文引入虛擬多孔介質(zhì)球孔隙率模型，將顆?？紤]為3～5倍大的多孔介質(zhì)球，將顆粒在原本位置體積占有的影響平均到比原有位置處顆粒體積大幾倍的網(wǎng)格上，最后再將所有顆粒對其周圍網(wǎng)格的影響疊加總和，從而完成了顆粒拉格朗日法到流體計算域歐拉場的映射轉(zhuǎn)換。虛擬多孔介質(zhì)球孔隙率模型對于實際顆粒粒徑尺寸接近于網(wǎng)格尺寸時產(chǎn)生的數(shù)值誤差較小，對于多孔介質(zhì)球粒徑尺寸為實際顆粒粒徑3～5倍情況下的計算穩(wěn)定性更好。

每個顆粒影響區(qū)域下覆蓋著很多個流體計算域網(wǎng)格，每個顆粒的影響域都會在其覆蓋流體計算網(wǎng)格上產(chǎn)生一個體積占有值，代表著每個網(wǎng)格被顆粒i所占有的體積比值（圖1）。對所有顆粒執(zhí)行相同的操作，最終每個網(wǎng)格上都有被不同顆粒所影響的體積比值，將這些值在每個流體網(wǎng)格上疊加，就能得到流場中的顆粒所占有體積的信息，見式（8）：

式中：i為第i個顆粒；di，j為顆粒與周圍流體網(wǎng)格中心的距離；Vp，i為第i個顆粒的體積；Vcell，j為第j個流體網(wǎng)格的體積；a為多孔介質(zhì)球與實際顆粒的半徑尺寸比，這里取值為4。

對于每個網(wǎng)格，被固體顆粒占據(jù)的體積與孔隙率之和為1，孔隙率見式（9）：

1.4流體與顆粒之間的相互作用

流體與顆粒之間的相互作用是流-固耦合模型連接流場信息與顆粒信息的主要內(nèi)容，流體與顆粒之間相互作用力計算的是否準確，是流-固耦合成敗的關(guān)鍵［13］。模擬設(shè)置的顆粒密度遠大于液相密度，流體重力、虛擬質(zhì)量力、Basset力、Saffman力和Magnus力在所有受力中所占量級較小，研究中將其忽略，主要考慮相間拖曳力的作用。本文考慮目前常用的2種拖曳力模型并進行驗證。

a）Ergun、Wen和Yu模型。結(jié)合計算密集顆粒堆積體的Ergun公式以及計算稀疏顆粒堆積體的Wen和Yu公式，得到Ergun、Wen和Yu模型，見式（10）、（11）［14］：

其中：

式中：fp為流體對顆粒產(chǎn)生的拖曳力；ε為網(wǎng)格孔隙率；dp為顆粒直徑；ufi為顆粒球心處的流體速度，通常采用顆粒球心所在流體網(wǎng)格中的流速；up為顆粒速度；CP為拖曳力系數(shù)；Rep為顆粒雷諾數(shù)，見式（12）：

Rep=（12）

b）Di Felice模型。Di Felice［15］通過對顆粒堆積體內(nèi)滲流實驗數(shù)據(jù)進行整理和擬合，提出了Di Felice模型公式，此模型認為顆粒堆積體內(nèi)流體作用于顆粒的拖曳力可以根據(jù)單顆粒在流體內(nèi)平穩(wěn)沉降時受到的拖曳力fp0乘一個關(guān)于堆積體孔隙率ε的函數(shù)f（ε）來獲得，見式（13）：

fp=f（ε）fp0（13）

其中fp0見式（14）：

其中，Cp0為拖曳力系數(shù)，見式（15）：

Cp0=（0.63+）2

f（ε）為冪函數(shù)形式，見式（16）：

f（ε）=ε2-Χ

其中，見式（17）：

Χ=3.7-0.65exp-

1.5耦合計算流程

流體相和顆粒相相互作用的耦合需要得到流體對顆粒的作用和顆粒對流體的作用。因此耦合方案由2個部分構(gòu)成：一是獲得顆粒周邊流體的性質(zhì)，為顆粒模型提供邊界條件；二是獲得并更新相間信息，如體積分數(shù)。本文的CFD-DEM耦合計算步驟依賴于ANSYS FLUENT的DPM離散相模型的應(yīng)用，其計算流程見圖2。

2單顆粒沉降模擬

球形顆粒的基本沉降是流體-顆粒相互作用中的重要問題，選擇此類顆粒沉降問題，以驗證流體和顆粒的相互作用力是否計算準確。選取的計算域尺寸為4 mm×4 mm×8 mm，初始水面高4 mm，球形顆粒從水面上方1 mm的位置處由靜止狀態(tài)自由落體。顆粒直徑0.1 mm，顆粒密度2 500 kg/m3?？諝庀嗝芏? kg/m3，運動黏度1×10-5 Pa·s；水相密度1 000 kg/m3，運動黏度0.002 Pa·s。網(wǎng)格尺寸設(shè)置為0.1 mm，網(wǎng)格尺寸和顆粒粒徑尺寸比為1。

對于低雷諾數(shù)，顆粒的沉降速度規(guī)律符合Stokes定律使用范圍下的理論解，顆粒入水后的速度隨時間的變化關(guān)系見式（18）：

對2種拖曳力模型進行對比驗證，結(jié)果見圖3。

由圖3可知，虛擬多孔介質(zhì)球孔隙率模型可以應(yīng)用于計算網(wǎng)格顆粒尺寸比為1時的情況，有效消除了傳統(tǒng)模型中網(wǎng)格尺寸與顆粒粒徑比大于4的限制。同時，在該情形下，采用Di Felice拖曳力模型的結(jié)果與Stokes理論解具有較好的一致性，最終沉降速度為0.004 87 m/s。因此，本文選取Di Felice拖曳力模型。

在0.015～0.020 s模型與理論解相差較大，這是由于Stokes理論解的計算分為水面上自由落體和水面下顆粒沉降，水面為定值；而數(shù)值模擬時，水面會被顆粒擾動，水面處液相和顆粒相之間有相互作用。在拖曳力的驗證中，主要依據(jù)最終沉降速度來驗證拖曳力計算的準確性。本文中顆粒最終達到了理想的穩(wěn)定沉降速度，可以驗證模型的可靠性。

3顆粒入水及出水模擬驗證

顆粒團入水及出水涉及復(fù)雜的水面運動情況以及流體-顆粒之間相互作用。顆粒團進入水中，會排開相應(yīng)體積的水，水面上升；顆粒浮出水面時，顆粒原本占據(jù)的體積會被液體填充，自由面高度會下降。通過模擬顆粒團入水及出水的過程，驗證改進后的模型對于固體顆粒、水和空氣三相體積處理的守恒性。

流體計算域尺寸為0.05 m×0.05 m×0.15 m，初始時刻自由水面位置為z=0.05 m處，上方空氣相密度1 kg/m3，運動黏度1×10-5 Pa·s；下方水相密度1 000 kg/m3，運動黏度1×10-3 Pa·s。初始時刻顆粒像柵格一樣在水面上方排列，其x、y、z方向上數(shù)目分別是20×20×30，共12 000個（顆粒按初始時刻z方向的位置每10層分為一組，紅色、綠色、藍色分別表示上、中、下3組顆粒），顆粒粒徑為0.002 m，顆粒間球心間距略大于顆粒粒徑為0.002 01 m，顆粒密度為2 500 kg/m3。

DEM時間步長設(shè)置為2.0×10-5 s，CFD時間步長設(shè)置為1.0×10-4 s。計算2 s后顆粒沉降結(jié)束自由面恢復(fù)平靜，此時更改水體密度為7 500 kg/m3，黏度不變，再計算2 s顆粒浮出水面并達到穩(wěn)定。

圖4顯示，未改進的模型在T=2 s時，顆粒完全入水，水面仍為初始液面高度，未考慮顆粒相的體積排水作用。改進模型后，T=2 s時顆粒完全入水，水面上升的高度為顆粒團全部固體顆粒的體積，水面達到理論值0.070 1 m，T=4 s時顆粒團上浮運動結(jié)束，水面上升的高度為被淹沒顆粒所排開水的體積，水面為理論值0.058 0 m。證明本文模型有效改善了FLUENT中離散相模型無法考慮顆粒排水效應(yīng)的缺點。

4顆粒堆積體坍塌-涌浪試驗?zāi)M

顆粒堆積體坍塌-涌浪試驗可以驗證CFD-DEM流固耦合模型模擬滑坡-涌浪災(zāi)害的有效性［16］。Robbe-Saule等［17］所做的室內(nèi)水槽試驗包含了顆粒堆積體失穩(wěn)、涌浪形成和傳播等類似滑坡涌浪災(zāi)害鏈的過程，被廣泛用作標準檢驗?zāi)Ｐ汀?/p>

試驗裝置左側(cè)為顆粒堆積區(qū)，高0.41 m，寬0.165 m，總體積0.102 m3；右側(cè)為一個長2.00 m、高0.30 m、寬0.15 m的玻璃水槽，水槽中充滿深度為0.05 m的水。數(shù)值模擬計算參數(shù)見表1，DEM碰撞參數(shù)見表2。

圖5為顆粒堆積體坍塌過程不同時刻水槽試驗結(jié)果（每組狀態(tài)左圖）與數(shù)值模擬結(jié)果（每組狀態(tài)右圖）的對比。通過對比圖中顆粒堆積體坍塌過程及水面涌浪運動狀態(tài)可以看出，在整個演進過程中，數(shù)值模擬結(jié)果與水槽試驗結(jié)果相差較小。當顆粒與水面相撞時會形成初始涌浪，隨著更多的顆粒侵入水體，涌浪高度逐漸增加；隨后，涌浪在傳播過程中逐漸變陡，波峰在重力影響下濺落到水體表面。對比水槽試驗結(jié)果和數(shù)值模擬結(jié)果可知，CFD-DEM流固耦合模型能夠較為準確地捕捉這一現(xiàn)象，初步驗證了該模型解決滑坡涌浪問題的可靠性。

圖6顯示，未改進的模型由于未考慮顆粒的排水效應(yīng)會低估涌浪高度，改進模型后可縮減涌浪高度模擬的誤差。水槽試驗中測量的最大涌浪高度為8.30 cm，數(shù)值模擬計算得到的最大涌浪高度為8.67 cm，結(jié)果較為一致。數(shù)值模擬形成涌浪達到最大涌浪高度的時間早于水槽試驗結(jié)果約0.1 s，可能是由于物理試驗提拉擋板的時間誤差及其對顆粒的擾動造成的。同時數(shù)值模擬將不規(guī)則顆粒簡化為單一粒徑球形顆粒，低估了顆粒間以及顆粒與流體間的作用力，造成結(jié)果偏大。該研究結(jié)果與肖華波等［16］數(shù)值模擬結(jié)果類似。通過對顆粒坍塌運動過程和涌浪高度演化過程的對比，發(fā)現(xiàn)顆粒運動及涌浪傳播模擬結(jié)果與物理試驗結(jié)果較為吻合，證明了本文所提出改進后的CFD-DEM流固耦合模型能較好還原顆粒運動狀態(tài)及涌浪的波動變化。

5結(jié)論

本文對FLUENT中CFD-DEM流固耦合模型進行改進，并進行多角度驗證，得到以下主要結(jié)論。

a）通過引入虛擬多孔介質(zhì)球孔隙率模型，改進了FLUENT離散相模型無法做到考慮顆粒體積分數(shù)的缺點，克服了傳統(tǒng)CFD-DEM模型網(wǎng)格顆粒臨界尺寸限制在實際應(yīng)用方面的困難。

b）通過單顆粒沉降模擬，驗證了網(wǎng)格顆粒尺寸比為1時，采用Di Felice拖曳力模型的結(jié)果與理論解具有較好的一致性，最終沉降速度為0.004 87 m/s。

c）通過多顆粒入水及出水模擬，顆粒完全入水后，水面上升高度達到理論值0.070 1 m，驗證了虛擬多孔介質(zhì)球孔隙率模型可以有效改進FLUENT離散相模型，考慮顆粒體積的排水效應(yīng)，從而實現(xiàn)強耦合。

d）通過顆粒堆積體坍塌-涌浪試驗?zāi)M，驗證了改進后的CFD-DEM流固耦合模型能有效還原顆粒運動狀態(tài)及涌浪的波動變化，數(shù)值模擬計算得到的最大涌浪高度為8.67 cm，與物理試驗結(jié)果較為一致，驗證了該模型解決滑坡涌浪問題的可靠性。

參考文獻：

［1］周家文，陳明亮，瞿靖昆，等.水庫滑坡災(zāi)害致災(zāi)機理及防控技術(shù)研究與展望［J］.工程科學(xué)與技術(shù)，2023，55（1）：110-128.

［2］魏云鵬.基于CFD-DEM方法的散粒體滑坡涌浪研究［D］.昆明：昆明理工大學(xué)，2022.

［3］王雷，解明禮，黃會寶，等.不同規(guī)?；氯胨T發(fā)涌浪災(zāi)害特征差異性分析［J］.人民珠江，2024，45（2）：18-28.

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

［5］WANG T，ZHANG F S，F(xiàn)URTNEY J，et al.A review of methods，applications and limitations for incorporating fluid flow in the discrete element method［J］.Journal of Rock Mechanics and Geotechnical Engineering，2022，14：1005-1024.

［6］ZHANG P Y，MU L L，HUANG M S.A coupled CFD-DEM investigation into hydro-mechanical behaviour of gap-graded soil experiencing seepage erosion considering cyclic hydraulic loading［J］.Journal of Hydrology，2023，624.DOI：10.1016/j.jhydrol.2023.129908.

［7］TSUJI Y，KAWAGUCHI T，TANAKA T.Discrete particle simulation of two-dimensional fluidized bed［J］.Powder Technology，1993，77（1）：79-87.

［8］劉德天，傅旭東，王光謙.CFD-DEM耦合計算中的體積分數(shù)算法［J］.清華大學(xué)學(xué)報（自然科學(xué)版），2017，57（7）：720-727.

［9］ZHAO T，HOULSBY G T，UTILI S.Investigation of granular batch sedimentation via DEM-CFD coupling［J］.Granular Matter，2014，16：921-932.

［10］姚鵬，王遠，冒劉鵬，等.CFD-DEM模型最優(yōu)網(wǎng)格尺寸的理論推證［J］.泥沙研究，2023，48（5）：1-7，34.

［11］李東陽，年廷凱，吳昊，等.滑坡-堵江-涌浪災(zāi)害鏈模擬的DEM-CFD耦合分析方法及其應(yīng)用［J］.工程科學(xué)與技術(shù)，2023，55（1）：141-149.

［12］吳亮.含固體顆粒的兩相流界面變化的數(shù)值研究：DEM-VOF方法的實現(xiàn)［D］.天津：天津大學(xué)，2018.

［13］彭愷然，劉紅帥，平新雨，等.CFD-DEM耦合模擬中拖曳力模型精度［J］.吉林大學(xué)學(xué)報（地球科學(xué)版），2021，51（5）：1400-1407.

［14］ERGUN S.Fluid Flow Through Packed Columns［J］.Chemical Engineering Progress，1952，48：89-94.

［15］DI FELICE R.The Voidage Function for Fluid-Particle Interaction Systems［J］.International Journal of Multiphase Flow，1994，20（1）：153-159.

［16］肖華波，王澤皓，石偉明，等.猴子巖水庫色玉滑坡涌浪災(zāi)害鏈CFD-DEM耦合數(shù)值模擬［J］.工程科學(xué)與技術(shù)，2023，55（1）：161-170.

［17］ROBBE-SAULE M，MORIZE C，HENAFF R，et al.Experimental in-vestigation of tsunami waves generated by granular col-lapse into water［J］.Journal of Fluid Mechanics，2021，907.DOI：/10.1017/jfm.2020.807.

（責(zé)任編輯：程茜）