CFD-DEM方法在砂土深基坑管涌中的應用

夏 霄，甘鵬路

（中國電建集團華東勘察設計研究院有限公司，浙江 杭州 311122）

1 概述

管涌作為一種土體內部侵蝕現象是基坑工程中的常見災害，它引起的基坑擋墻倒塌事故時有發生[1]。例如，2015 年杭州某濱江基坑開挖中發生了管涌誘發的擋墻倒塌事故，承受較高水力梯度的砂質土體從擋墻裂縫處涌出，造成坑外土體大幅沉降和擋墻的嚴重變形[2]。類似基坑擋墻薄弱處（如樁或地連墻連接處）破裂引發的管涌災害屢有發生[3-4]。

正是由于管涌細觀機制的復雜性，相應的數值分析研究也逐步展開。一類是利用有限元將土體視為連續性介質并與滲透水流進行耦合的唯象分析[5-7]，它很難捕捉管涌中土顆粒流失，重新排列、沉積等重要特征。因此，另一類基于流體力學（CFD）和散體介質（DEM）的多相耦合仿真分析方法（CFD-DEM）在管涌研究中越發常見[8]。這種方法將土體視為不同相的組合，如土骨架以及土骨架孔隙中的水流與細顆粒，它可考慮內部侵蝕或管涌過程中細顆粒從土骨架中的流失。CFD-DEM耦合仿真分析方法可捕捉基坑管涌過程中土體的主要力學特征，也可用于解決實際基坑管涌問題[9]。

本文利用CFD-DEM 耦合仿真分析方法研究細粒含量（Fc）與水力梯度（i）對間斷級配土基坑中管涌災害的影響。首先介紹了CFD-DEM 流固耦合仿真模擬方法，然后通過建立基坑管涌模型，考慮不同的Fc和i進行一系列參數分析，得到不同工況下基坑土體顆粒速度場和流場特征，驗證了該方法在基坑管涌模擬中的優越性。

2 CFD-DEM流固耦合仿真分析方法

CFD-DEM 方法由三個部分構成，分別是離散元法（DEM），計算流體力學（CFD）及CFD-DEM耦合算法。DEM主要用于對土顆粒的仿真，CFD則是仿真土中的滲透水流。CFD-DEM耦合算法主要計算流體—顆粒間相互作用力，以及實現CFD 與DEM 間信息交互。開源 DEM 軟件 LIGGGHTS 和 CFD 軟件 Open?FOAM 是本研究中具體使用的多相耦合仿真分析工具。下面簡要介紹這三個仿真分析方法的基本原理和控制方程。

2.1 DEM方法控制方程

DEM 中顆粒運動遵循牛頓第二定律,顆粒i的平動和轉動由下式計算：

式中：mi、It——顆粒i的質量和轉動慣量；

Ui、ωi——顆粒i的平動速度與轉動角速度；

本研究中顆粒間接觸模型采用Hertz 非線性接觸定律[10]。對于切向上的滑動阻力，利用Mohr-Coulomb準則計算。

2.2 CFD方法控制方程

CFD 域中，流體運動主要由連續性方程和動量方程控制，如下所示：

式中：Uf、p——流體網格的平均速度和壓強；

εf——流體網格中顆粒體的孔隙率；

fp——單位體積顆粒與流體間平均相互作用力，即

ΔV——流體網絡體積；

k——顆粒數；

ρf、μf——流體密度和粘度；

g——重力加速度，本文中為50倍的常重力速度。

2.3 流體—顆粒相互作用力計算

在管涌問題中，流體—顆粒相互作用力Ff主要包括拖曳力、粘滯力和動壓差力。本研究中也主要考慮上述三種流固相互作用力，忽略其它相互作用力對管涌特性的影響。Ff的計算如公式（3）所示：

流體對土顆粒產生的拖曳力采用Buijtenen 等（2011）[11]方法計算，如公式（4）所示，它適用于顆粒體密度較大的情況，在流體雷諾數變化較大時也成立：

式中：dp、Up——顆粒的直徑與速度；

Vp——單個顆粒的體積；

εp——流體網格中顆粒所占體積百分比，εp+εf=1；

Rep——單個顆粒的雷諾數。

作用于顆粒上的動壓差力和粘滯力計算見公式（5）和（6）[12]：

式中：τ——流體剪切力。

3 基坑管涌仿真分析工況與參數設置

3.1 仿真分析工況

本文主要通過不同的土體細粒含量Fc與水力梯度i工況來研究其對基坑管涌的影響。一般認為Fc直接決定土體組構，采用Fc=15%和35%的試樣研究兩種不同的典型土體組構對基坑管涌的影響。

為研究基坑管涌啟動前后土體力學性質演變，在仿真過程中，水力梯度由i=0.8 增至1.8（i= Δp/ρgL，Δp是基坑內處水壓差，L是模型中流體最短流動路徑，本研究中是0.5m，g是模型中采用的重力加速度值）。當i=0.8時，管涌并未啟動，而i增至1.8時，管涌開始發生并持續至仿真結束。表1總結了本研究考慮的四個管涌仿真工況。圖1是本研究中采用的土體的級配曲線，粗粒組和細粒組均為單一直徑顆粒，分別為2mm和0.5mm。考慮到CFD-DEM 耦合仿真計算效率，每一種工況的仿真時長均為1.0s。結果表明，在此較短的仿真時長內，基坑變形也可達到穩定狀態。

圖1 試樣級配曲線

表1 仿真分析工況

3.2 模型尺寸與參數設置

圖2 是基坑DEM 模型的尺寸示意圖。由于本研究主要針對管涌過程中顆粒與流體力學特性演化，因此擋墻假定為鋼性體，在仿真過程中不發生變形。為研究擋墻裂縫引起的基坑管涌災害，在擋墻10cm高度處開有2.5cm 高度的孔洞。擋墻右側和左側分別是基坑的開挖側與支護側，高度分別為50cm 和25cm，長度分別均為25cm和50cm。整個模型厚度為5cm，將實際的基坑管涌簡化成平面應變問題。結果表明這一尺寸設置可以減小邊界效應的影響。仿真模型的CFD域與DEM域范圍大小相同，保證所有顆粒均受滲透水流的作用。表2列出了模型中顆粒與流體的力學性質參數。

圖2 基坑管涌仿真DEM模型

表2 CFD與DEM參數匯總

3.3 邊界條件設置

為模擬實際中基坑擋墻兩側水壓差，在CFD 域的開挖側頂面施加一定水壓，使得基坑內外水力梯度在整個仿真過程中均維持在0.8或1.8。CFD域其余面均為不透水面。對于顆粒與DEM邊界的接觸，假定它們之間的摩擦系數等于0，法向接觸剛度是顆粒彈性模量的10倍。

3.4 仿真分析步驟

首先在基坑區域上方生成一定數量的顆粒并使其在重力作用下自由下落。在填至模型1/5高度時，將顆粒體豎向縮至一定的相對密實度，本研究中將土體相對密實度設定為Re=65%～67%。重復上述兩個步驟直至顆粒填滿基坑支護側與開挖側。在DEM模型準備完成后，立即在CFD域的支護與開挖側頂面上施加固定的水壓力差以開始基坑管涌的仿真模擬。所有顆粒及它們間的接觸信息在整個仿真過程中每隔0.01s輸出。

4 基坑管涌仿真分析主要結果

圖3是細粒含量分別為Fc=15%和35%情況下基坑管涌啟動時刻（t=0.2s，i=1.8）地層變形與顆粒速度場。在管涌啟動時刻，擋墻裂縫處和基坑開挖側坑底土體首先開始流動，具有較大的初始速度，而支護側土體相對流動速度較小。圖4 是在Fc=15%和35%時，基坑管涌仿真結束時刻地層變形與顆粒速度場。此時，支護側土體出現明顯沉降，開挖側坑底土體明顯隆起。在擋墻裂縫前方出現一條管涌通道，支護側土顆粒經此流入開挖側并逐漸趨于穩定。當土體中含較多細粒時，支護側土體沉降較大。這主要由于裂縫寬度與顆粒粒徑之比較大時，顆粒不易堵塞，更易通過裂縫，因而在相同水力條件下流動性更強，導致最終支護側顆粒流失量增大。

圖4 最終時刻(t=1.0s)基坑變形與顆粒速度場

圖5 是管涌過程基坑流場的CFD-DEM 仿真結果。基坑支護側水體在水頭差作用下通過擋墻裂縫流入開挖側，裂縫周圍水流速度達到最大。

圖5 基坑管涌過程流場仿真結果

此外，通過比較不同工況下支護側土體流失量，可以發現，對于土體細粒含量較低（Fc=15%）情況，土體流失量較多。如前分析，這主要由于裂縫寬度與顆粒粒徑之比較大時，顆粒在裂縫處不易堵塞，更易流動導致的。

進一步對比Fc=15%和35%兩種工況下的支護側地表沉降演化，可以發現，地表沉降都在靠近擋墻處最大，且在管涌啟動時刻（t=0.2s）不明顯，而在t=0.4～0.8s 時刻內發展較快，隨后沉降增量趨于減小。對于Fc=15%的情況，由于其支護側土體流失量較小，地表沉降范圍主要局限在0

然而，開挖側坑底土體隆起演化結果與支護側地表沉降演化不同。在管涌啟動時刻（t=0.2s），坑底隆起即達到最大值，隨后下降并由于支護側土體流入而最終略有上升。這一現象與圖3中開挖側土體在管涌啟動時刻即具有較大速度相吻合。這說明坑底隆起主要是管涌流對坑底土顆粒向上的拖曳力，使得土體懸浮于管涌流中引起。

5 結論

本文介紹了一種可反映流體與砂土顆粒間相互作用的CFD-DEM 方法，并利用該法對深基坑開挖過程擋墻裂縫導致的管涌問題進行了仿真分析。通過砂土基坑管涌的宏細觀仿真結果，可以發現CFD-DEM 方法可以較好重現基坑管涌的主要特征，如坑底隆起、地表沉降和管涌通道中流體流動特征。管涌過程中地表沉降逐步增大，而坑底隆起在管涌啟動時即達到最大值。管涌災害主要由滲流力對基坑開挖側土體力學性質的削弱引起，在實際工程中為了防止管涌災害的發生，除提高擋墻或止水帷幕施工質量外，仍需對開挖側土體實施良好的加固措施。