一種CFD-DEM 流固耦合方法在滲流導致城市地面沉降問題中的應用

李曉蛟， 陸 燁， 武亞軍

(上海大學土木工程系， 上海200444)

城市局部區域地面沉降會對周邊構建筑物及使用人群造成影響， 帶來巨大經濟損失， 而地下水滲流是造成地面沉降的主要原因. 對于城市局部區域地面沉降這一問題， 現場監測有較大難度， 且現場外部環境隨時間變化， 影響測量結果的因素較多. 模型試驗一方面難以揭示地面沉降的細觀機理， 另一方面也不能得到土體全場位移規律. 因此， 數值模擬方法可作為現場監測和模型試驗的補充， 對沉降過程及機理進行深入研究.

在以往的研究中， 對于城市地面沉降的數值模擬主要采用有限元方法(finite element method， FEM)[1-2]或離散元方法(discrete element method， DEM). 有限元方法能夠較好地模擬宏觀尺度下連續介質的變形， 但是在模擬細觀尺度下離散介質的運動規律時有所不足， 即有限元方法要滿足變形協調條件， 其內部細觀結構的尺寸要遠小于模擬對象的尺寸[3]， 這對于研究滲流作用下顆粒的遷移規律顯然不再適用. 而Cundall等[4]提出的離散元方法則非常適合用于模擬細觀尺度下土顆粒的遷移運動， 并應用于多種巖土問題[5-8]. 但是， 當模擬對象中涉及滲流或孔壓時， 如何處理流體是需要解決的問題. Li等[8]曾提出一個簡化的流體-顆粒耦合方法， 該方法中顆粒之間的孔隙被認為是通過顆粒接觸之間的管子連接的流域， 通過建立網格來表征流域的體積以及管子的長度和直徑. 這種方法對少量顆粒系統比較簡便， 但在三維情況下， 孔隙體積計算將變得繁瑣.

計算流體動力學(computational fluid dynamics， CFD)技術是目前廣泛使用的模擬流體的計算方法， 不僅可以進行三維流場的計算， 還可以模擬復雜形狀的流場[9]. Zeghal等[10]首次在巖土工程問題的研究中引入了CFD-DEM 流固耦合方法， 分析了土坡滲流和飽和土體震動液化問題， 取得了較好的結果. 周健等[11]和羅勇等[12]也運用CFD-DEM 方法開展了土體滲流等問題的研究. 蔣明鏡等[13]利用PFC2D， 借助Fish 函數進行二次開發研究流固耦合問題， 通過單顆粒水下自由沉降和一維固結試驗算例來驗證CFD-DEM 耦合的可行性， 取得了較好的結果. 王胤等[14]提出了一種考慮粒間滾動阻力的CFD-DEM 模擬方法， 能夠較好地反映顆粒的形狀效應， 建立了流體與固體顆粒的動量傳遞模型， 能夠很好地描述巖土問題中水-土結合的力學性質， 但該方法的應用范圍具有一定局限性， 例如DEM 模型不能識別具有復雜形狀的CFD 模型. 從以上研究可以看出， 已有的CFD-DEM 流固耦合模擬方法雖然實現了二維到三維的跨越， 在巖土工程中的應用也越來越多， 但大部分研究工作基于結構化網格， 所模擬的流場形狀十分規則單調(如方形或圓柱形)， 對于大型的、具有復雜形狀的實際工程， 模擬起來仍有一定難度.

鑒于此， 本工作采用一種 CFD-DEM 方法， 即顆粒流程序 (particle flow code， PFC)與計算流體力學程序(Fluent)流固耦合方法， 通過單顆粒在水中自由沉降算例驗證本方法的正確性；依據本方法建立具有復雜流場形狀的基坑臨空面滲漏與河道邊坡滲流兩個算例， 在設置不同地下水位的前提下， 通過Fluent 求解滲流場， 然后將滲流場導入PFC 土顆粒模型中， 實現耦合計算， 最終模擬出不同水力梯度作用下滲流導致地面沉降的整個動態過程. 為便于區分，本工作提出的CFD-DEM 方法在下文中統一稱為PFC-Fluent 方法.

1 PFC-Fluent 計算理論基礎

本工作中的流固耦合程序采用PFC3D 模擬固相顆粒的細觀行為， 采用均化流體CFD 技術模擬液相流體的運動. 對于固相顆粒， 采用牛頓第二定律求解顆粒的平移和轉動方程; 對于液相流體， 通過求解平均的Navier-Stokes 方程模擬孔隙流體的運動; 流體對顆粒的作用力由Ergun 方程表達[15]. 每個流體單元包含一定數量的土顆粒， 流體流動產生的滲透力以體力的形式作用于土顆粒.

PFC-Fluent 耦合程序是單向耦合， 對于地下水滲流導致的地面沉降問題， 土體在地下水滲流力作用下被沖刷破壞， 而地下水位基本保持不變， 流體對土顆粒的影響遠遠大于土顆粒對流體的影響， 因此在這種情況下， 采用PFC-Fluent 耦合程序是可行的.

1.1 固相運動方程

固相顆粒的運動由顆粒離散元來模擬. 作用在單個顆粒上的主要有體力(顆粒的自重)、其他與之接觸的顆粒通過接觸點施加的接觸力， 以及流體作用在該顆粒上的力. 在這些力的作用下， 顆粒的運動滿足牛頓第二定律， 平動和轉動方程為

式(1)和(2)中：u 為土顆粒速度矢量; t 為時間; fmech為顆粒重力與接觸力的合力; ffluid為流體對單個土顆粒的合力; m 為土顆粒質量; g 為重力加速度; ω 為角速度矢量; M 為力矩; I 為慣性矩.

1.2 液相運動方程

對于密度不變不可壓縮流體的運動， 可以由平均的Navier-Stokes 連續方程和動量方程來表達， 即

式(3)和(4)中：n 為孔隙率; ?為哈密頓算子; v 為流體速度矢量; ρf為流體密度; ?p 為壓力梯度; τ 為黏性應力張量; fd為單位體積內流體與顆粒之間的作用力.

流體網格內孔隙率n 的計算方法有兩種：質心法和多面體法. 本工作采用質心法， 即當顆粒球心位于網格內時即識別為整個顆粒都處于網格內. 質心法的優點是計算效率高， 缺點是隨著顆粒的運動， 孔隙率會產生跳躍， 因此在模擬中設定的網格尺寸相對顆粒來說足夠粗， 以削弱這種影響， 相對于整個流場， 網格又是足夠細的， 這就在保征準確性的前提下提高了計算效率.

1.3 流體對土顆粒的作用力

在固-液兩相飽和介質中， 流體對顆粒的作用力從大類上可分為靜水力和動水力. 靜水力即為顆粒靜置于流體中的浮力作用， 流體對顆粒的浮力大小與該顆粒周圍的流體壓力梯度有關， 表達式為

式中：fb表示單個顆粒受到的浮力; r 表示顆粒半徑.

動水力包括拖拽力、上舉力和虛擬質量力， 后兩種力的大小和拖拽力相比可以忽略. 拖拽力的表達式為

式中：fdrag為土顆粒受到的拖拽力; Cd為拖拽力系數， 適用于整個區域的近似公式

其中Re 為雷諾數，

μ 為黏滯系數.

2 PFC-Fluent 耦合流程與程序實現

采用PFC 軟件求解固體顆粒運動方程， 用Fluent 軟件求解流體動力方程. 將Fluent 求得的三維流場導入PFC 程序中， 將顆粒與流體相互作用力作為耦合紐帶， 實現相應的流固耦合計算. PFC-Fluent 耦合計算流程如圖1 所示.

圖1 PFC-Fluent 耦合方法流程圖Fig.1 Flow chart of PFC-Fluent coupled method

3 數值模擬驗證

首先， 通過3 個不同尺寸顆粒在水下自由沉降的算例， 驗證PFC-Fluent 耦合程序的準確性. 然后， 通過不同水位差的兩組算例來驗證PFC-Fluent 耦合計算的可行性. 兩組數值算例如下：基坑開挖時， 圍護結構出現滲流導致坑外地層沉降; 河道邊坡在滲流作用下發生破壞.

3.1 算例一 3 個不同尺寸顆粒水下自由沉降

3.1.1 自由沉降速度理論解

斯托克定律表明， 小球在水中重力作用下的沉降速度是一定的. 采用PFC-Fluent 耦合程序模擬3 個不同尺寸顆粒在水中重力作用下的自由落體運動， 以驗證本方法的正確性. 在黏性流體中， 球形顆粒在重力作用下下落的運動方程如下：

當下落的球體達到一個穩定的速度， 加速度為0， 式(9)變為

3.1.2 PFC-Fluent 數值模擬

將半徑為1.0， 1.5， 2.0 mm 的顆粒置于黏滯系數為0.001 的流體中， 在重力作用下自由下落. 顆粒與液體的密度分別為2 650， 1 000 kg/m3， 重力加速度為9.8 m/s2. 將上述參數代入式 (7)， (8)， (10)， 并用 Matlab 數值求解得到自由沉降速度理論解分別為 0.245， 0.324，0.391 m/s.

耦合程序中顆粒速度變化過程如圖2 所示， 最終速度達到穩定值， 與理論解一致， 表明本工作采用的耦合程序可以實現對顆粒施加流體作用.

圖2 不同時刻顆粒下落速度Fig.2 Drop velocity of particle at different time

3.2 算例二 基坑臨空面滲漏

本算例考慮基坑開挖完成后還未澆筑底板時， 由于降水坑內外水位差較大， 同時基坑圍護結構在臨近坑底處出現裂縫， 水土在水力梯度作用下向基坑臨空面滲漏， 導致坑外土體流失從而產生地面沉降.3.2.1 數值建模

基于顆粒流軟件PFC3D 平臺建立基坑邊坡土體模型的基本流程如下.

第 1 步， 采用 wall 命令生成一個尺寸為 6.0 m×1.0 m×9.0 m(長×寬×高)的模型槽， 在模型槽內按線性接觸模型生成高度為6.5 m， 孔隙率為0.4 的砂土層. 土體通過一系列球形顆粒來模擬， 對所有土顆粒施加重力， 在重力作用下顆粒沉積、固結， 并達到平衡狀態. 在砂土層上部按平行黏結點接觸模型生成高度為1.5 m， 孔隙率為0.4 的黏土層， 施加重力， 同樣達到平衡狀態. 基坑模型中顆?？倲禐?7 721 個， 顆粒最大半徑為0.06 m， 最小半徑為0.02 m， 平均半徑為0.04 m， 顆粒與墻體間的摩擦系數取值為0， 以減小邊界效應的影響. 模型中參數取值參考類似地質條件[16]， 如表1 所示.

表1 PFC3D 模型細觀參數Table 1 Mesoscopic parameters of PFC3D model

第2 步， 由于本工作主要考慮臨空面滲漏對地面沉降的影響， 不考慮基坑圍護結構變形的影響， 故采用wall 命令在距模型槽右側邊界2.5 m 處生成一面無厚度墻作為基坑的圍護結構，圍護結構總長7 m， 刪除圍護結構右側高度為3.5 m 區域內的土顆粒， 圍護結構插入比為1. 針對實際工程中基坑圍護結構產生裂縫的情況， 在圍護結構離坑底0.5 m 處開一圓孔. 若圓孔半徑太小， 土顆粒會堵塞裂縫， 故經過不斷試算， 得到圓孔直徑為5 倍土粒平均粒徑可引起滲流，因本工作重點不對圓孔尺寸進行對比， 故只取這一種尺寸進行計算. 基坑土體模型如圖3 所示.

第 3 步， 打開 PFC 流體模塊， 將 Fluent 軟件計算出的流場 (3.2.2 節詳述)導入 PFC 軟件中， 即可進行流固耦合計算.

3.2.2 流場計算

對于多孔介質的孔隙率及相關阻力系數， 孔隙率取值與土體孔隙率一致; 對于土中水滲流這種流速較低的層流狀態， 阻力系數只考慮黏性阻力系數， 根據孔隙率與土顆粒粒徑計算出黏性阻力系數值[13].

本算例模擬了3 種不同地下水位， 其中坑底地下水位保持不變， 均為坑底地面以下0.5 m， 坑外地下水位依次取地面以下1.5， 0.5 和0 m. 流場進口采用壓力進口邊界， 進口壓力為0 kPa， 出口采用壓力出口邊界， 出口壓力為0 kPa， 其余邊界設置為墻邊界， 水在重力作用下流動. 計算出的流場如圖4 所示， 將其導入PFC 建立的土體模型中進行耦合計算.

圖4 基坑開挖流場流速矢量圖Fig.4 Velocity vector diagram of flow field for excavation

3.2.3 模擬結果

在算例二中， 造成地面沉降的原因有兩方面， 一方面是土顆粒在地下水滲流作用下從基坑圍護結構破損處流出(流出的土顆粒被程序自動刪除)， 繼而上部土體往下運動， 最終造成坑外地面沉降; 另一方面是土顆粒在水力梯度作用下從高水位向低水位運動， 導致坑外地面沉降.本算例分別模擬了3 種不同地下水位， 通過不同時步土顆粒位置坐標變化得出地表沉降位移如圖5 所示. 從圖中可以看出， 隨計算時步的增加， 地表沉降值越來越大; 基坑邊緣處地表沉降值最大， 沉降可達0.35 m; 坑外水位越高， 地表沉降值越大. 水位高度在地面以下0.5 m 與在地面以下0 m 的地表沉降比較接近， 原因可能是二者水力梯度比較接近.

圖5 基坑臨空面滲漏引起的地表沉降Fig.5 Ground settlements caused by leaking on retaining structure of excavation

在3 種不同水位高度下， 土顆粒位移云圖比較一致， 故僅選取水位高度在地面以下0 m 這一種情況下的不同時步的土顆粒位移云圖(見圖6)， 圖中各顏色表示土顆粒位移大小，藍色表示位移為0 m， 紅色表示位移為1.0 m 及以上， 中間每隔0.1 m 用不同顏色表示. 模擬結果表明， 隨著時步不斷增加， 基坑圍護結構開孔處土顆粒不斷流失， 上部土體不斷往下移動， 最終形成地表沉降. 同時， 基坑底部在滲流作用下有輕微隆起.

圖6 水位高度在地面以下0 m 時的顆粒位移Fig.6 Soil particle displacements for underground water level of 0 m

3.3 算例三 河道邊坡滲流

本算例考慮河道邊坡在水力梯度作用下發生滲流破壞誘發邊坡外側的地面沉降.

3.3.1 數值建模

基于顆粒流軟件PFC3D 平臺建立河道邊坡土體模型的基本流程如下.

第1 步， 采用wall 命令生成一個尺寸為12.5 m×1.0 m×9.0 m(長×寬×高)的模型槽， 在模型槽內按線性接觸模型生成高度為6.5 m， 孔隙率為0.4 的砂土層. 土體通過一系列球形顆粒來模擬， 對所有土顆粒施加重力， 在重力作用下顆粒沉積、固結， 并達到穩定狀態. 在砂土層上部按平行黏結點接觸模型生成高度為1.5 m， 孔隙率為0.4 的黏土層， 施加重力， 達到穩定狀態.模型槽中顆粒總數為131 308 個， 土顆粒尺寸與算例二一致. 顆粒與墻體間的摩擦系數取值為0， 以減小邊界效應的影響. 模型中參數取值如表1 所示.

第 2 步， 河道邊坡是一個傾斜角為 45°， 長 6.4 m 的護坡墻， 由 2 000 個半徑為 0.1 m 的圓球按平行黏結接觸模型生成. 借助Fish 語言刪除邊坡右側高度為3.5 m 內的土顆粒(見圖7).

第 3 步， 打開 PFC 流體模塊， 將 Fluent 軟件計算出的流場導入 PFC 軟件中， 即可進行流固耦合計算.

3.3.2 流場計算

河道邊坡模型的流場同樣用Fluent 進行計算， 由于算例三與算例二的土體顆粒級配、孔隙率一樣， 故本算例中流體網格尺寸和黏性阻力系數與算例二一致. 本算例模擬了河道水位的3 種高度， 其中地下水位保持不變， 為地面以下0.5 m， 河道水位高度依次取為2， 1 和0 m. 流場進口采用壓力進口邊界， 進口壓力為0 kPa， 出口采用壓力出口邊界， 根據水位深度， 出口壓力依次取為0， 9.81， 19.62 kPa， 其余邊界設置為墻邊界. 計算出的流場如圖8 所示， 將其導入PFC 建立的土體模型中進行耦合計算.

圖 7 河道PFC 模型Fig.7 PFC model for river embankment

圖8 河道滲流流場流速矢量圖Fig.8 Velocity vector diagram of flow field for river embankment

3.3.3 模擬結果

算例三的地表沉降位移如圖9 所示. 當河道水位高度為2 m 時， 運行80 萬步后地表沉降位移最大值為0.34 m; 當河道水位高度為1 m 時， 運行80 萬步后地表沉降位移最大值為0.58 m; 當河道水位高度為0 m 時， 運行80 萬步后地表沉降位移最大值為1.14 m. 可見， 邊坡兩側水位高度相差越大， 地表沉降值也越大; 距河道坡頂越近， 地表沉降值越大.

圖9 河道邊坡滲流引起的地表沉降Fig.9 Ground settlements caused by seepage under river embankment

在3 種不同河道水位下， 運行80 萬步時不同水位土顆粒位移云圖如圖10 所示(圖中顏色表示意義與圖6 一致). 模擬結果表明， 河道邊坡兩側水位高度相差越大， 土顆粒運動越迅速，地表沉降值也越大.

圖10 河道邊坡滲流引起的土顆粒位移Fig.10 Soil particle displacements caused by seepage under river embankment

河道邊坡滲流算例主要用來觀察在初始流場作用下滲流對土顆粒位移的影響， 最終模擬結果符合河道滲流規律. 河道邊坡兩側水位高度差和滲透系數是影響流場分布的重要因素， 在本算例中二者無明顯變化. 土顆粒雖然發生了較大位移， 但仍然處于現有流場范圍內， 且流體網格尺寸相對顆粒尺寸是足夠粗的， 當有土顆粒流出時， 周圍其他土顆粒也會流入， 因此孔隙率變化并不是特別明顯. 本算例主要研究的是固體顆粒的位移規律， 流體對顆粒的影響遠大于顆粒對流體的影響， 因此采用PFC-Fluent 是有一定合理性的.

4 結 論

(1) 在Navier-Stokes 方程與Ergun 方程基礎上， 通過計算流體力學程序Fluent 求解滲流場， 將求得的滲流場按節點、單元的形式導入顆粒流軟件PFC 中， 實現了PFC 與Fluent 的固液耦合.

(2) 通過3 個不同尺寸顆粒在水下自由沉降算例， 驗證了PFC-Fluent 方法的準確性; 通過基坑臨空面滲漏與河道邊坡滲流兩個算例， 驗證了PFC-Fluent 方法的可行性， 并得出了在不同水位高度下地面沉降規律和土體全場位移云圖， 對地面沉降災害的預防與治理具有一定的指導意義.

(3) PFC-Fluent 方法的優勢在于可從細觀角度展現滲流引起地面沉降的整個動態過程，并能實現對各種復雜形狀流場的模擬求解計算. 這也表明本方法可應用于大型、復雜、實際工程中地面沉降的流固耦合計算， 并從細觀層面揭示其機理.

(4) PFC-Fluent 方法是單向耦合的， 當土顆粒對流體的作用極大地影響流場分布時需進行改進. 接下來將對本方法作進一步的改進， 流場的更替與雙向耦合將是后續研究的重點.