顆粒離散相模型的改進及其在滑坡涌浪數值模擬應用驗證

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

關鍵詞：滑坡涌浪；流固耦合；CFD-DEM耦合方法；數值模擬

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

郝悅彤，肖鴻.顆粒離散相模型的改進及其在滑坡涌浪數值模擬應用驗證［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

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

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

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

1流固耦合數值方法

1.1計算流體運動學方法

1.1.1控制方程

對于流體相，采用以下基本方程：連續性方程見（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）方法捕捉自由面，對涌浪的形成及演變過程中多相流體交界面進行追蹤。通過針對多相的體積分數的連續性方程求解來完成對相之間的界面跟蹤，對于第q相流體，

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

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

式中：腳標對應不同種類的流體。

1.2離散元方法

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

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

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

單個顆粒的旋轉運動方程見式（7）：

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

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

1.3虛擬多孔介質球孔隙率模型

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

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

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

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

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

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

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

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

其中：

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

Rep=（12）

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

fp=f（ε）fp0（13）

其中fp0見式（14）：

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

Cp0=（0.63+）2

f（ε）為冪函數形式，見式（16）：

f（ε）=ε2-Χ

其中，見式（17）：

Χ=3.7-0.65exp-

1.5耦合計算流程

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

2單顆粒沉降模擬

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

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

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

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

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

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

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

流體計算域尺寸為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方向上數目分別是20×20×30，共12 000個（顆粒按初始時刻z方向的位置每10層分為一組，紅色、綠色、藍色分別表示上、中、下3組顆粒），顆粒粒徑為0.002 m，顆粒間球心間距略大于顆粒粒徑為0.002 01 m，顆粒密度為2 500 kg/m3。

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

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

4顆粒堆積體坍塌-涌浪試驗模擬

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

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

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

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

5結論

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

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

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

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

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

參考文獻：

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

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

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

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

［5］WANG T，ZHANG F S，FURTNEY 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耦合計算中的體積分數算法［J］.清華大學學報（自然科學版），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模型最優網格尺寸的理論推證［J］.泥沙研究，2023，48（5）：1-7，34.

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

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

［13］彭愷然，劉紅帥，平新雨，等.CFD-DEM耦合模擬中拖曳力模型精度［J］.吉林大學學報（地球科學版），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］肖華波，王澤皓，石偉明，等.猴子巖水庫色玉滑坡涌浪災害鏈CFD-DEM耦合數值模擬［J］.工程科學與技術，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.

（責任編輯：程茜）