管線滲漏引發地面塌陷的DEM-CFD耦合數值模擬

張小玲，王曉麗， 杜修力，許成順

(北京工業大學 城市與工程安全減災教育部重點實驗室， 北京 100124)

近幾年城市地下管線因破損滲漏誘發地面塌陷的現象層出不窮，給國家帶來了巨大的財產損失[1]。城市地陷具有突發性、隨機性、隱蔽性等特點，當這種地面塌陷發生在人群密集區時，更會造成極大的經濟損失同時也對居民的生命造成威脅，因此研究地下管線滲漏誘發地面塌陷的災變機理具有十分重要的現實意義。

地下管線滲漏誘發地面塌陷的問題，從微觀角度分析是土體在水流作用下的流失問題和地面土體逐漸失穩塌陷的機理和過程，重點是流體與土顆粒的相互耦合問題[2]。從本質上講，巖土材料都是由離散的、尺寸不一、形狀各異的顆粒或塊體組成，因此在該問題的研究上將土體視為離散單元更合理。目前針對滲流問題的研究計算效率較高的方法是土體采用離散單元法(DEM)，流體采用計算流體動力學方法(CFD)，將兩者結合起來考慮耦合作用時采用DEM-CFD計算方法。

目前，許多學者針對管線滲漏引發地面塌陷事故的外因與內因進行了研究。在外因方面，主要研究管線滲漏水壓、滲漏范圍、管線位置等因素對地面塌陷的影響。張成平等[3]通過室內試驗發現地表沉降值和沉降范圍會隨著管線滲漏水范圍的增加而增大，而管線位置與管內水壓通過影響滲漏范圍間接影響地面塌陷程度。陶高梁等[4]通過其研制的滲流破壞模型試驗裝置, 分析了泄漏孔徑和泄漏水壓力對砂土滲流破壞的影響規律，發現泄漏孔徑越小，產生空洞的臨界水壓力越大。在內因方面，主要研究土體顆粒級配、孔隙率或土顆粒與流體間的相互作用力等因素的影響。Vincens等[5]通過考慮接觸力的分布、接觸數和平均應力來分析土壤顆粒的內部穩定性。 Scholtes等[6]在恒定圍壓條件下進行了三軸試驗，發現土體發生滲透破壞后土體的剪脹性和峰值應力減小。

隨著離散元軟件的推廣，計算方法的不斷完善，應用DEM或DEM-CFD耦合計算方法研究滲流微觀問題得到推廣。Ma等[7]采用DEM-CFD計算方法分析了水力梯度與流速的關系，分析了孔隙率，粒徑組合和顆粒間滾動阻力對流動特性的影響。還有一些學者修改了經典的DEM-CFD計算方法以使流體網格更精細化或計算方法更準確。蔣明鏡等[8]基于流體的可壓縮性修改了經典DEM-CFD計算方法，通過單顆粒自由沉降速度符合理論解驗證改正方法的可行性。王胤等[9]考慮顆粒的形狀效應對DEM-CFD計算方法的影響，總結了顆粒滾動阻力對沙堆休止角與土體堆積孔隙率的影響。

目前在應用DEM-CFD耦合計算方法進行流體與土顆粒相互作用的研究主要是針對流體流速較小符合達西定律時的緩慢流動。這種情況下，由于流體壓力分布不均而作用于單位質量土體上的力(本文稱之為流體壓力梯度力) 影響很小可忽略不計，但當流體不滿足達西定律的高速滲流或者土顆粒尺寸較大、形狀不規則時，流體壓力梯度力的影響不容忽略。劉先珊等[10]通過土顆粒的運動位移間接反應了流體壓力梯度力的存在，而目前的大部分研究未考慮流體高速滲流時流體壓力梯度力對流固耦合過程的影響。

本文主要基于DEM-CFD雙向耦合計算方法，考慮高速滲流時流體壓力梯度力對流固耦合作用的影響，重點從機理角度分析流體滲流過程中土體力學特性和流體水力特性的變化規律與動態變化過程。

1 CFD-DEM 流固耦合原理

目前關于研究流固相互耦合作用的問題主要包括流體和流體之間的相互作用、土顆粒和流體的相互作用和土顆粒與土顆粒的相互作用。

1.1 流體與流體的相互作用

1.1.1 連續性方程

連續性方程指的是流體在運動過程中滿足質量守恒定律，即單位時間內通過體積元的質量凈通量等于體積元內質量的變化量。具體的表達式[11]如式(1)所示：

(1)

式中：n為土體孔隙率，v為流速矢量。

1.1.2 動量方程

Navier-Stokes方程是描述黏性不可壓縮流體動量守恒的運動方程，即單位時間內單元體動量的增加必等于單位時間凈流入單元體的動量加上單元體內流體所受合力。在流體流動過程中，不可壓縮流體動量增加來源有：黏性力產生的動量，壓差力產生的動量，體積力產生的動量和顆粒對流體作用力產生的動量。一般黏性不可壓縮流體Navier-Stokes方程表達式為：

(2)

式中：v是流速矢量;ρf為流體密度；g為重力加速度；p為流體壓力；τ為流體的黏性應力張量；fint為單位體積內顆粒與流體的作用力。

1.2 土顆粒與流體的相互作用

滲流過程中土顆粒與流體的相互作用力包括流體對顆粒的拖曳力，流體壓力梯度力，浮力，粒子運動引起的滯流阻力以及其他不穩定力等。為簡化計算，在流固耦合計算過程中只考慮拖曳力、流體壓力梯度力和浮力。

ffluid=fdrag+fgrad+ffloat

(3)

式中：ffluid是指流體對顆粒的作用力；fdrag是指流體對顆粒的拖曳力；fgrad是流體壓力梯度力；ffloat是顆粒在流體中的浮力。

以下針對滲流過程中流體對顆粒的作用力中的拖曳力、流體壓力梯度力和土顆粒所受浮力具體進行計算推導。

1.2.1 拖曳力求解

目前對于流體施加給土顆粒的拖曳力的計算,即便是形狀規則的土顆粒組成的多孔介質，也沒有明確的計算公式，主要是通過試驗數據擬合多孔介質的拖曳力，本文應用的計算公式[12]為：

(4)

式中：Cd為拖曳力系數；d為顆粒直徑；u為顆粒速度矢量；v是流速矢量；n為顆粒所在流體單元的孔隙率，經驗系數χ計算公式為：

(5)

拖曳力系數Cd計算公式為：

(6)

式中：Re為雷諾數，計算公式為：

(7)

式中：μf為流體動力黏滯系數。

1.2.2 流體壓力梯度力求解

在滲流過程中流體壓力梯度力計算公式：

(8)

當流體滿足達西定律時：

(9)

式中：uxo是流體的平均流速；k是土體的滲透系數。

當流速較大不滿足達西定律且土體孔隙率≤0.8時，流體壓力梯度▽p采用Ergun[13]公式：

(10)

當滲流過程中土體孔隙率大于0.8時，流體壓力梯度▽p采用Wen等[14]的計算公式：

(11)

1.2.3 浮力求解

滲流過程中土顆粒所受浮力的計算公式：

(12)

1.3 土顆粒與土顆粒的相互作用

基于牛頓第二定律，土顆粒所受的作用力包括土顆粒所受外力、流體對土顆粒的作用力、土顆粒自身重力，計算公式如下：

(13)

(14)

2 DEM-CFD耦合過程

2.1 耦合過程

由于CFD模塊無法模擬真實流體在土體中的流動狀態，因而基于體積平均原理[15]將流體劃分為若干個粗網格，每個粗網格代表一個流體單元，以流體單元內流體的平均特性來代表該單元的特性，該單元內的流體信息平均分配給單元內的土顆粒，相鄰流體單元通過信息傳遞來體現流體的流動性。本文在計算過程中考慮了非穩定流時流體壓力梯度力對流固耦合過程的影響，具體的DEM-CFD流固耦合計算流程見圖1。

圖1 CFD-PFC流固耦合流程圖

圖1左側為采用離散元軟件PFC3D建立的力學模型流程圖，右側為利用CFD模塊建立的流體模型流程圖，中間區域為DEM-CFD交互模塊。DEM與CFD模塊的耦合計算以兩者計算時間是否相等為判斷條件。首先DEM模塊將土體孔隙率n、顆粒的速度u、顆粒對流體的拖曳力fd等信息傳遞給CFD模塊， CFD模塊根據傳入信息求解流速v、水壓p、壓力梯度▽p等信息；若CFD模塊通過計算判斷流體N-S方程收斂，則將流體對顆粒的拖曳力、流體壓力梯度力、浮力傳遞到DEM模塊。若N-S方程不收斂，則根據流體邊界條件重新計算流速v、水壓p、壓力梯度▽p等信息再進行收斂判斷。在耦合過程中為了保證計算的穩定性，在CFD模塊設定流體100時步內為穩定，DEM模塊每20時步更新顆粒的作用力，拖曳力，孔隙率等參數。

2.2 模型尺寸驗證

目前針對管線滲漏引發地面塌陷問題的研究，多數成果都是針對二維模型進行了計算[16]，這與實際工程情況存在一定差異。本文采用離散元軟件PFC3D建立了管線滲漏引發地面塌陷的三維數值模型進行計算。考慮到計算效率與適用性，因計算模型尺寸與土顆粒尺寸不能同時滿足實際的要求，在建立三維計算模型時以土顆粒尺寸和顆粒數量為參考依據，土顆粒粒徑采用實際粒徑，生成的顆粒數量與二維計算模型顆粒生成數量接近。

為了考察三維計算模型尺寸對計算結果的影響，本文建立了不同尺寸大小的滲流計算模型，通過分析流體壓力梯度力的計算結果來確定最終的模型尺寸。圖2是在滲流深度相同(均為20 cm)的情況下，滲流入水口面積分別為15 cm×15 cm，20 cm×20 cm，30 cm×30 cm和40 cm×40 cm四種條件下流體壓力梯度力的對比結果。從圖2可以看出，當入水口面積為20 cm×20 cm，30 cm×30 cm和40 cm×40 cm時，流體壓力梯度力的峰值與最終的穩定值都比較接近；但是當入水口面積為15 cm×15 cm的數值結果與其他結果相差較大。圖3比較了當流體壓力梯度力趨于穩定時不同模型計算的總時間，由圖可以看出，所建立的模型尺寸越大，計算的顆粒數量就越多，所需的計算時間越長。因此綜合考慮各因素，本文采用20 cm×20 cm的模型尺寸來進行下一步的數值計算。

3 計算結果分析

由于離散元軟件無法反映流體在土體孔隙中的真實流動過程，因此在所建立的模型上游施加非零水壓p，下游水壓為0，上下游存在初始壓力梯度▽p，在該壓力梯度的作用下通過相鄰流體單元的信息傳遞來模擬流動過程。本文考慮的是管線破損后滲漏水對周圍土體的侵蝕作用，通過用土體上下游兩端的壓力梯度▽p代替水流從管線中滲漏而出的水壓力。圖4為滲流過程的模型示意圖。其中上游入水口為y=0的x-z平面，大小20 cm×20 cm，下游出水口為y=20 cm的x-z平面，大小20 cm×20 cm ，其他邊界為不透水邊界，水流以一定流速v從上游垂直流入。

圖2 模型尺寸影響

圖3 不同模型計算總時間對比

圖4 滲流過程模型示意圖

地表突然發生大規模塌陷，則土體內部已產生容納上方土體的空間[17]，因滲流過程產生的大尺寸裂隙或實際土體內已存在的土洞、早期人防工程等均可容納上方土體，因此本文在建模過程中考慮了地下空洞的存在，模擬自然界天然形成的空洞，空洞周圍沒有襯砌。在生成土顆粒力學模型時，考慮計算效率與顆粒數量的合理性，取計算模型大小20 cm×20 cm×20 cm，土顆粒粒徑4 mm，土顆粒數量大約11.6萬個，土體的中間位置存在半徑為5 cm的球形空洞，管線與空洞的位置與大小詳情見圖4(b)。計算模型采用的土顆粒計算參數見表1，流體與墻體計算參數見表2。流體網格劃分單元為1 cm×1 cm×1 cm，共有20×20×20個網格。根據楊斌等[18]關于多孔介質滲流規律的描述，流體從線性層流過渡到非線性層流臨界流速為0.002 3 m/s～0.007 8 m/s，從層流轉化為紊流的臨界流速為0.016 m/s～0.048 m/s。因此本文在分析流體壓力梯度力的作用時，對比了滿足達西定律的流速與不滿足達西定律的流速流體壓力梯度力的變化規律，其中滿足達西定律的流速為0.001 m/s，不滿足達西定律的流速為0.1 m/s、0.2 m/s、0.3 m/s。

表1 土顆粒計算參數

表2 流體與墻體計算參數

3.1 空洞上方土體失穩坍塌過程

圖5(a)—圖5(e)是沿滲流方向某截面的土顆粒在流體壓力梯度▽p作用下逐漸失穩直至坍塌的演化過程示意圖。

由圖5(a)可看出，在流體滲流過程中，空洞上方土顆粒最先失穩發生運動，隨著滲流過程的進行，土體失穩范圍越來越大，并逐漸向地表擴展。當土體失穩范圍延伸至地表時，地面開始出現輕微凹陷，此時空洞上方土體均產生微小位移，見圖5(b)，空洞幾何形狀發生變化。隨著土體壓力梯度的改變，空洞上方土顆粒在自重與流體滲流力等共同作用下，繼續向下運動且位移逐漸增大，塌陷范圍進一步向地表拓展，該失穩過程經過多次循環后，地表松動顆粒范圍越來越大，地表沉降凹槽深度也越來越大。由圖5對塌陷過程的模擬結果可知，土體塌陷過程中產生了傾斜的滑移面，模型表面呈現近似圓錐形的破壞形態。本文所用土顆粒粒徑為4 mm，隸屬于土的工程分類標準[19]中的粗粒土，這跟已有研究[17]中認為的粗粒土地層地陷形狀近似為圓形漏斗狀的結論基本一致。

半連續和連續模型的影響。在使用兩種模型類型的默認參數的第一個評估周期中，本文發現連續聲學模型比半連續聲學模型在開發集上的效率為36.10%。

圖5 土體坍塌過程位移變化 (單位:cm)

3.2 流體水力性質變化規律

土顆粒在滲流過程中受到的流體壓力梯度力與流壓梯度▽p有關，因此在流體速度為0.1 m/s,0.2 m/s,0.3 m/s三種工況下，分析了地表與A點土體兩端水壓差值▽p1的變化。圖6反映了該水壓差值▽p1隨滲流時間變化規律，由圖6可以看出，▽p1隨滲流時間呈現先減小后趨于穩定的變化規律，而穩定滲流的結果是在滲流過程中上下游的水頭差不隨時間變化，由此可知本文模擬的地下管線滲漏引發地面塌陷的過程是非穩定滲流過程。

圖6 土體上游與A點兩端水壓差隨時間的變化

3.3 土顆粒力學性質變化規律

本文研究管線滲漏后引發地面塌陷的機理以空洞上方土體的變化為主要研究對象，計算了空洞上方距離地表0.04 m處的一個流體單元(圖4中A點)內的土顆粒在流體作用下拖曳力fdrag、浮力ffloat、流體壓力梯度力fgrad的變化情況，由于分析流體對顆粒的作用力時流體單元較小，因此作用力數值較小，計算結果分別見圖7、圖8。

圖7 流體對土顆粒的拖曳力

圖8 土顆粒所受壓力梯度力

圖7是計算一個流體單元內所有顆粒受到的平均拖曳力fdrag隨時間變化曲線。在整個計算過程中當流速為0.1 m/s，0.2 m/s，0.3 m/s三種工況下土顆粒受到的流體拖曳力由最大值急劇減小后增大再減小趨于穩定。這是因為拖曳力的大小取決于顆粒與流體運動的相對速度，見公式(4)。在流體發生滲流過程的初始階段，由于土體剛接觸流體作用力，流體運動速度與土顆粒的運動速度差值最大，此時土顆粒受到的流體拖曳力在整個計算過程中為最大值。隨著圖6流體壓力梯度▽p1的減小，顆粒與流體的相對速度減小，土顆粒受到流體的拖曳力隨即減小。當運動的土顆粒從空洞上方的土體中流失掉落至空洞，此時流體與土顆粒的相對速度又急劇增大，因而土顆粒受到的拖曳力又急劇增大。當空洞下側土體與下游土體水壓差超過某值時，土顆粒重新發生流失過程，顆粒所受流體拖曳力開始緩慢減小。與蔣明鏡等[8]的研究中單個顆粒最終沉降速度為0.321 m/s時的拖曳力為8×10-6N相比，本文計算的拖曳力結果與蔣明鏡等[8]研究中的拖曳力結果相當，但由于本文研究的流體單元內含有多個土顆粒，因此圖7拖曳力結果數值會偏大。

當流體流速為0.1 m/s，0.2 m/s，0.3 m/s三種工況時，在滲流過程中土顆粒受到的浮力計算結果相同，具體數值大約為3.28×10-4N。

圖8是滲流過程中土顆粒受到的流體壓力梯度力隨時間變化曲線。由圖可看出，在整個滲流計算過程中，當流速為0.1 m/s，0.2 m/s，0.3 m/s時，土顆粒受到的流體壓力梯度力由最大值急劇減小，后又增大最終緩慢減小趨于某穩定值，由公式(8)—式(10)可知，流體壓力梯度力與流壓梯度▽p有關，而▽p又與顆粒與流體運動的相對速度有關，因此流體壓力梯度力變化原因與土顆粒受到的拖曳力原因類似。但當流速為0.001 m/s，流體壓力梯度力數值幾乎接近于0 N，說明流體滿足達西定律時不需要考慮流體壓力梯度力的影響；但流體高速滲流時，流體壓力梯度力與浮力數量級相同，且其數值大于拖曳力，因此在求解流體對土顆粒的作用力方程時不能忽略流體壓力梯度力的影響。

由于土顆粒的尺寸與幾何形狀影響流體壓力梯度力，因此針對土顆粒半徑為2 mm、3 mm、4 mm時的流體壓力梯度力進行了計算對比，結果如圖9所示。由圖可以看出，在計算時間內不同土顆粒半徑時的流體壓力梯度力的變化規律基本相同，但土顆粒半徑越大，流體壓力梯度力越大，最終趨于穩定的時間越快。

在滲流過程中對比流體壓力梯度力、拖曳力和浮力的計算結果可知，圖8中流速為0.001 m/s時流體滲流滿足達西定律，土顆粒所受的流體壓力梯度力數值很小，可以忽略不計；但當流速為0.1 m/s，0.2 m/s，0.3 m/s時，土顆粒所受流體壓力梯度力與浮力數量級相當，且在數值上要大于拖曳力，因此在求解力學方程中其數值不容忽視。在針對地下管線滲漏引起地面塌陷問題的研究中，很多學者都僅考慮了土顆粒所受流體拖曳力和浮力的影響而忽略了流體壓力梯度力的影響，這會導致極大的誤差。因此本文在分析高速滲流引發地面塌陷的問題時，考慮了流體壓力梯度力對土體顆粒力學性能的影響。

圖9 不同土顆粒半徑的流體壓力梯度力隨時間的變化

3.4 滲流過程中土顆粒流失規律

通過試驗手段分析土體滲流過程，只能從宏觀角度觀察到土體的侵蝕形狀，土體內部具體流失細節難以直接觀察，采用數值模擬的方法可以反映土體流失過程中的顆粒流失詳情、土體力鏈演化規律與土體內部顆粒的重組過程。在地下管線發生滲漏后，只有當土體內部流體的水力梯度高于起始水力梯度[20]，流體克服了結合水的粘滯阻力時，土體內部的流體才會發生滲流。在流體壓力梯度▽p作用下土體中的流體發生滲流的過程中，描述顆粒流失現象的參數有土顆粒流失速率、土顆粒接觸數變化量、土體孔隙率等指標。圖10、圖11、圖12分別對滲流過程中這些參數的變化規律進行了計算分析。

土顆粒流失速率是指某時間段內土顆粒流失數量占上一時間段顆粒數量的比重。圖10是施加了壓力梯度▽p后流體發生滲流的過程中，土顆粒的流失速率隨時間的變化曲線圖。圖中的點代表本文數值模擬得到的土顆粒流失速率結果，曲線為對土顆粒流失速率進行擬合得到的曲線。由圖中的結果可知，在流速為0.1 m/s，0.2 m/s和0.3 m/s三種工況下，在整個計算時間內土顆粒流失速率均大于零，且為一種波動式增長，說明在流體滲流過程中土顆粒一直處于流失狀態。在前期流失階段，顆粒流失速率成比例增加，到后期顆粒流失速率急劇增大，這是因為在滲流過程中土體內部流體對土顆粒的滲透力需要一定的時間維持。前期土顆粒排列緊密，顆粒流失過程較緩慢；隨著土顆粒的流失，土體內部形成貫通的滲流通道并不斷擴大，更多的土顆粒發生流失，流失現象更加明顯。隨著壓力梯度越大，即流體的流速越大，土顆粒的流失越快。

圖10 土顆粒流失速率隨時間的變化

圖11 土顆粒接觸數變化量隨時間變化

圖12 土體孔隙率隨時間的變化

土顆粒的接觸數是指滲流發生的過程中，某一時刻土體內部顆粒之間接觸粘結部位的總數量，接觸數變化量是指該時刻土體內部的接觸數與上一時刻接觸數的差值。接觸數變化量能夠反映短期內土體內部顆粒之間在流體壓力梯度▽p作用下的碰撞激烈程度與粘結狀態。圖11是在滲漏水不同的流速下的土顆粒接觸數變化量的變化趨勢。由圖可以看出，在流體流速v=0.1 m/s，0.2 m/s和0.3 m/s時，土顆粒接觸數的變化量均呈現先減小后增大的變化趨勢。若土顆粒的接觸數變化量大于零，說明此時間段內土體在水流作用下內部土顆粒接觸數量越來越多，顆粒粘結位置增加，土體密實度增大；若土顆粒的接觸數變化量小于零，說明土顆粒間的粘結數量減小，此時間段內土顆粒正處于流失狀態。同時，由圖可以看出，該數值模型中大約0.4 s時刻是顆粒流失量正負值的突變時刻，說明該時刻是土體內部接觸力鏈減小和顆粒間開始出現松動的臨界時刻；而0.8 s時刻是接觸數變化量變化趨勢的轉折點，說明此時刻土顆粒的流失速率達到最大。土顆粒接觸數變化量的曲線斜率后期緩慢減小到0，體現了土顆粒在水壓差的作用下土體流失速率由快到慢最終趨于穩定的變化規律。

土體孔隙率也是衡量土顆粒流失程度的重要指標之一。圖12所示為土體上下游兩端在流體壓力梯度▽p的作用下位于空洞上方土體孔隙率隨時間的變化情況。由于滲流過程中孔隙率計算區域土體位于空洞上方，該區域土顆粒流失迅速，因而計算時間相對較短。由圖可以看出，土體孔隙率前期變化較緩后迅速增加到1，這是因為前期土顆粒排列緊密，當土體上下游兩端存在流體壓力梯度時，土顆粒開始緩慢流失。隨著滲流過程的進行，土體內部孔隙逐漸形成貫穿通道，土顆粒流失速率加快，當孔隙率計算區域的土顆粒全部流失后，土體孔隙率即達到1。

4 結 論

(1) 本文對經典的CFD-DEM 耦合計算方法做進一步的改進，在高速滲流過程中引入了流體壓力梯度力對流固耦合過程的影響。流體壓力梯度力在高速滲流流速較大時，流體壓力梯度力與浮力的數量級相當其數值大于拖曳力，且土顆粒半徑越大，流體壓力梯度力越大。

(2) 在滲流過程中土顆粒流失速率在整個計算過程中緩慢增加；土顆粒接觸數變化量呈現先減小后增大的變化趨勢；土體孔隙率先緩慢增加后迅速增大。

(3) 地下管線滲漏引發地面塌陷的過程中，空洞上方的土顆粒最先失穩發生運動，松動的顆粒逐漸向地表延伸，土體內部產生傾斜滑移面，地表最終呈現圓錐形凹陷，整個塌陷過程中的塌陷模式呈圓錐形變化。