管線滲漏引發(fā)地面塌陷的DEM-CFD耦合數(shù)值模擬

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

(北京工業(yè)大學(xué) 城市與工程安全減災(zāi)教育部重點(diǎn)實(shí)驗(yàn)室， 北京 100124)

近幾年城市地下管線因破損滲漏誘發(fā)地面塌陷的現(xiàn)象層出不窮，給國家?guī)砹司薮蟮呢敭a(chǎn)損失[1]。城市地陷具有突發(fā)性、隨機(jī)性、隱蔽性等特點(diǎn)，當(dāng)這種地面塌陷發(fā)生在人群密集區(qū)時，更會造成極大的經(jīng)濟(jì)損失同時也對居民的生命造成威脅，因此研究地下管線滲漏誘發(fā)地面塌陷的災(zāi)變機(jī)理具有十分重要的現(xiàn)實(shí)意義。

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

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

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

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

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

1 CFD-DEM 流固耦合原理

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

1.1 流體與流體的相互作用

1.1.1 連續(xù)性方程

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

(1)

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

1.1.2 動量方程

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

(2)

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

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

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

ffluid=fdrag+fgrad+ffloat

(3)

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

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

1.2.1 拖曳力求解

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

(4)

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

(5)

拖曳力系數(shù)Cd計算公式為：

(6)

式中：Re為雷諾數(shù)，計算公式為：

(7)

式中：μf為流體動力黏滯系數(shù)。

1.2.2 流體壓力梯度力求解

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

(8)

當(dāng)流體滿足達(dá)西定律時：

(9)

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

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

(10)

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

(11)

1.2.3 浮力求解

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

(12)

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

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

(13)

(14)

2 DEM-CFD耦合過程

2.1 耦合過程

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

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

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

2.2 模型尺寸驗(yàn)證

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

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

3 計算結(jié)果分析

由于離散元軟件無法反映流體在土體孔隙中的真實(shí)流動過程，因此在所建立的模型上游施加非零水壓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 滲流過程模型示意圖

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

表1 土顆粒計算參數(shù)

表2 流體與墻體計算參數(shù)

3.1 空洞上方土體失穩(wěn)坍塌過程

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

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

半連續(xù)和連續(xù)模型的影響。在使用兩種模型類型的默認(rèn)參數(shù)的第一個評估周期中，本文發(fā)現(xiàn)連續(xù)聲學(xué)模型比半連續(xù)聲學(xué)模型在開發(fā)集上的效率為36.10%。

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

3.2 流體水力性質(zhì)變化規(guī)律

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

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

3.3 土顆粒力學(xué)性質(zhì)變化規(guī)律

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

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

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

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

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

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

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

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

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

3.4 滲流過程中土顆粒流失規(guī)律

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

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

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

圖11 土顆粒接觸數(shù)變化量隨時間變化

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

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

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

4 結(jié) 論

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

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

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