城市地區(qū)管道滲漏引發(fā)地面塌陷成災機理的CFD-DEM聯(lián)合計算模擬分析

王越林，陸燁

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

近年來，有關(guān)城市地面塌陷的事故頻繁發(fā)生，其中，一些事故造成了巨大的經(jīng)濟損失甚至人員傷亡。地面塌陷是由于地質(zhì)條件(如巖溶形成)和人類活動(如采礦、建筑、管道滲漏)造成的[1-2]。在這些影響因素中，管道滲漏引起的事故大約占到55%。特別是城市地區(qū)，管道通常埋深較淺，每天承受各種荷載(如交通、施工、挖掘、打樁等)，管道在各種荷載作用下會由于自然原因或強制磨損而開裂甚至斷裂。某些城市地下水位較高，在水力梯度作用下，土體被侵蝕并沖入到無水流或重力流(供水或排水)管道。隨著滲漏侵蝕的擴散，擾動區(qū)從管道附近向上發(fā)展，最終導致地表塌陷。

在實際應用中，往往很難檢測到管道的滲漏，并且無法觀察滲漏侵蝕的過程。針對這個問題，一些學者進行了相關(guān)研究。Guo等[3]開展了隧道滲流侵蝕的室內(nèi)模型試驗，研究重力流管道破裂引起土體侵蝕的形成過程。王帥超[4]通過室內(nèi)模型試驗方法對管道破損滲漏導致的地面塌陷問題進行研究，發(fā)現(xiàn)地下管道裂縫不同形狀與尺寸對地下空洞發(fā)展規(guī)律的影響。除室內(nèi)模型試驗外，部分學者還通過數(shù)值模擬進行滲漏侵蝕研究。滲漏侵蝕問題涉及到兩種介質(zhì)：土體(固體)和水(流體)，兩者在地表塌陷的形成過程中都發(fā)揮著重要作用。因此，在相關(guān)研究中，采用離散單元法(DEM)和計算流體動力學(CFD)分別對土體和水進行模擬并進行聯(lián)合計算的數(shù)值模擬，得到越來越廣泛的應用。Shamy等[5]首次將CFD-DEM流固耦合方法引入到巖土工程問題研究中，對土坡滲流和飽和土體震動液化問題進行了分析。Tao等[6]利用CFD-DEM方法研究了堤壩中管道沖蝕的詳細過程。鄭剛等[7]利用PFC-CFD方法研究災害演化過程中的土顆粒應力鏈和縫隙兩側(cè)結(jié)構(gòu)受力的發(fā)展變化情況，認為漏水、漏砂過程是一個應力拱不斷建立、破壞的連續(xù)變化過程。周健等[8]采用離散元模擬開展了土體滲流及液化問題的研究。蔣明鏡等[9]將描述流體體變-壓力非線性關(guān)系的Tait狀態(tài)方程(EOS)和計算流體動力學(CFD)相結(jié)合，建立模擬弱可壓縮流體的CFD-DEM耦合計算模塊。通過單顆粒水下自由沉降和一維固結(jié)試驗算例來驗證CFD-DEM耦合的可行性。雖然，以上研究實現(xiàn)了二維到三維CFD-DEM流固耦合模擬方法的轉(zhuǎn)變，但仍無法用于模擬大型的、具有復雜邊界條件的實際工程。

鑒于此，筆者基于上海地區(qū)地面塌陷的工程背景，提出一種Fluent和PFC聯(lián)合計算方法，主要利用ANSYS中的Fluent模塊與顆粒流程序(PFC)進行聯(lián)合計算，并建立實際塌陷模型。在建模基礎(chǔ)上，研究上海北部地區(qū)地層(上部黏性土層、下部砂性土層)中管道滲漏引起滲流侵蝕的過程，以及不同地下管道裂縫尺寸和不同地下管道破裂位置作用下砂性與黏性土中地層變形和地面塌陷規(guī)律。

1 CFD-DEM聯(lián)合計算理論

離散元(DEM)法將模擬對象劃分為離散單元，并假設其為剛性圓盤或球體，單元間按實際情況選擇合理的接觸本構(gòu)關(guān)系且滿足運動方程，用時步迭代方法求解各單元運動方程，繼而求得不連續(xù)體的整體運動形態(tài)。離散單元法的特點是允許單元間的相對運動，單元間無須滿足連續(xù)介質(zhì)力學的位移連續(xù)與變形協(xié)調(diào)條件，其運動方程的時步迭代通常采用顯式的中心差分格式，無須建立大型剛度矩陣，計算速度快，所需儲存空間小，尤其適用于求解巖土體大變形及非線性問題。

計算流體動力學(CFD)是在計算機上數(shù)值求解流體動力學基本方程的方法，通過數(shù)值計算和圖像顯示進行流場分析、流場計算和流場預測，在比較短的時間內(nèi)，能預測模擬對象性能，并通過改變各種參數(shù)，以更加深刻地理解問題產(chǎn)生的機理，為實驗提供指導，節(jié)省實驗所需的人力、物力和時間。

1.1 流體對土顆粒的作用力

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

fb=-

(1)

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

(2)

1.2 計算流程

為模擬水對土體的滲透流蝕作用，提出一種Fluent和PFC聯(lián)合計算方法。采用DEM軟件PFC求解固體顆粒運動方程，用CFD程序Fluent求解流體動力方程。將Fluent求得的三維流場導入到PFC程序中，將顆粒與流體相互作用力作為耦合紐帶，實現(xiàn)相應的流固聯(lián)合計算，計算流程如圖1所示。網(wǎng)格劃分采用Fluent中的前處理器Gambit實現(xiàn)，由于結(jié)構(gòu)網(wǎng)格面對復雜幾何外形時生成困難，以及耗費大量人工，自動化程度不高等缺點，因此，復雜模型可采用非結(jié)構(gòu)化網(wǎng)格進行劃分。同時，李曉蛟等[10]對PFC-Fluent聯(lián)合計算的有效性進行了驗證。

圖1 PFC-Fluent聯(lián)合計算流程Fig.1 PFC-Fluent joint computing

2 土體與管道模型的建立

上海地區(qū)管道埋深通常在地面以下3 m以內(nèi)，地下水位在0.5～1.0 m，選用上海北部地區(qū)典型地層情況，即上部1.5 m厚黏性土層、下部2.5 m厚砂性土層進行數(shù)值建模，結(jié)合不同地下管道裂縫尺寸和不同地下管道破裂位置的影響，進行CFD-DEM聯(lián)合計算，并探究砂性土與黏性土地層中地面塌陷的產(chǎn)生和發(fā)展規(guī)律。

模型中模型槽尺寸為10 m×6 m×4 m(長×寬×高)，在模型槽內(nèi)自下而上按線性接觸模型分別生成高度為2.5 m的砂性土層和1.5 m的黏性土層，其孔隙率均為0.41，并對所有土顆粒施加重力，在重力作用下顆粒沉積、固結(jié)，并達到平衡狀態(tài)。最后在黏性土層上部覆蓋0.3 m厚的混凝土，代表低等級路面。在土顆粒的生成中，如果按照實際土顆粒尺寸在模型槽中生成，其數(shù)目將達到億級，這對計算機來說是不現(xiàn)實的，因此，對10 m×6 m×4 m的模型進行初步研究，結(jié)果表明，管道頂部滲漏引起的顆粒位移場是軸對稱的。在滿足極限計算能力的前提下，將模型槽尺寸更改為10 m×1 m×4 m，同時，對土顆粒尺寸進行擴大，在盡可能生成較多數(shù)目顆粒的前提下，生成顆粒的最大半徑為0.061 4 m，最小半徑為0.018 6 m，平均半徑為0.04 m，總數(shù)目約6萬顆。砂性土和黏性土的參數(shù)取值皆參考已有文獻[11-12]的取值，具體宏觀和細觀參數(shù)見表1、表2。

表1 土體模型宏觀參數(shù)Table 1 Macroscopic parameters of soil model

表2 土體模型細觀參數(shù)Table 2 Meso parameters of soil model

主要針對城市地下排水管道內(nèi)滲引發(fā)的地面塌陷問題進行研究。城市地下排水管道多為混凝土管道，內(nèi)徑大于500 mm時，采用鋼筋混凝土管，且地面塌陷多發(fā)生在大孔徑管道破裂處，因此，模擬的管道直徑為1 000 mm，位于地面以下3 m處。在PFC中無法對管道裂縫進行模擬，借助Rhino5.0進行管道模擬，并將其作為墻單元導入PFC中。其中，管道中心位置位于地面以下3 m。管道上裂縫假設為圓孔型，共分為5種：管道頂部開孔，孔徑分別為3倍、5倍和8倍平均粒徑；管道右側(cè)腰部開孔，孔徑為5倍平均粒徑；管道底部開孔，孔徑為5倍平均粒徑，見圖2。

圖2 土體與管道模型

流場采用CFD程序Fluent進行計算，對于地下水滲流這種多孔介質(zhì)流，一般采用Fluent程序中的多孔介質(zhì)模型進行模擬。Fluent程序中的多孔介質(zhì)模型是在標準流動方程中增加了一個動量源項，包括慣性損失項和黏性損失項。對于地下水滲流，由于流速較小，一般可忽略它的慣性損失項，只考慮黏性損失項。黏性損失項的大小表現(xiàn)為黏性阻力系數(shù)，土體的黏性阻力系數(shù)定義為

(3)

式中：Dp為平均粒徑；n為孔隙率；1/α為黏性阻力系數(shù)。

流場區(qū)域為地下水位以下，不包括管道內(nèi)部，取上海地區(qū)地面以下0.5 m水位為標準流場模型，地下水位表面為壓力進口邊界，進口壓力為0，管道裂縫為壓力出口邊界，出口壓力為0，其余邊界均為不透水邊界。將計算出的黏性阻力系數(shù)輸入Fluent程序中，得出的流場見圖3。由圖3可知，在水力梯度的影響下，地下水流入管道頂部裂縫，且上部黏性土層中流場流速較小，下部砂性土層中流場流速較大，均符合流體力學性質(zhì)。

圖3 流場矢量圖(單位：m/s)

3 CFD-DEM聯(lián)合計算結(jié)果分析

將Fluent計算出的流場導入PFC3D土體模型后，進行CFD-DEM聯(lián)合計算。管道頂部的滲漏模擬中，裂縫尺寸為5倍平均粒徑，地下水位0.5 m。圖4為不同計算周期下的水土流失過程。在滲漏初期，土體擾動區(qū)位于滲漏附近。隨著時間的推移(計算次數(shù)的增加)，擾動區(qū)不斷往橫向、縱向拓展，最終延伸到地面，形成沉降土槽。在此階段，如果混凝土路面受到荷載(一般為交通荷載)反復作用，則會導致路面塌陷，造成經(jīng)濟損失，甚至付出生命代價。

圖5(a)為沉降位移值達到0.02 m時采樣下限的擾動區(qū)包絡線圖。結(jié)果表明，當計算步數(shù)為20萬時，擾動區(qū)跨徑(擾動區(qū)到達地表時，引發(fā)的地面塌陷區(qū)域直徑大小)被限制在約1 m寬(管徑大小)范圍內(nèi)。隨著水土流失的加劇，擾動區(qū)跨徑擴大到管徑的5倍左右。為了描述滲漏引起的地表沉降，繪制了不同計算步數(shù)下的地表沉降圖，見圖5(b)。結(jié)果表明，隨著計算步數(shù)的增加，地表沉降值也逐漸增大，其最大值發(fā)生在滲漏點(管道頂部)。當計算步達20萬步時，地表沉降槽范圍在4 m左右，最大地表沉降值近0.3 m，這時已經(jīng)形成較為明顯的地面塌陷；隨著計算步數(shù)的增加，其造成的地面破壞程度和經(jīng)濟損失也隨之加大，當計算步達到80萬步時，地表沉降槽范圍已擴大到6 m，最大地表沉降值也增加至1 m左右。

圖4 地層中不同管道裂縫尺寸土顆粒位移云圖(單位：m)Fig.4 Cloud map of soil particle displacement in different pipeline crack sizes in

圖5 管道頂部5倍粒徑開孔時擾動區(qū)域包絡線與地表沉降曲線Fig.5 Enveloping line and surface subsidence curve of disturbance area when opening holes at the top of pipeline

3.1 不同管道裂縫尺寸的影響

通過將管道頂部圓形裂縫直徑設置為3倍平均粒徑、5倍平均粒徑和8倍平均粒徑，以研究裂縫開口尺寸的影響。圖6、圖7顯示了在計算步數(shù)為80萬步時，不同管道裂縫尺寸引起的土顆粒位移規(guī)律、擾動區(qū)情況和地表沉降情況，得出結(jié)論如下：

1)隨著裂縫尺寸的增加，土體內(nèi)部土顆粒的流失越來越嚴重，當5倍平均粒徑與8倍平均粒徑開孔的工況已經(jīng)形成較大的空洞區(qū)時，3倍平均粒徑開孔工況還未產(chǎn)生明顯的空洞。

2)隨著裂縫尺寸的增加，擾動區(qū)跨徑增加速度加快，當計算步數(shù)小于20萬時，由于3倍平均粒徑開孔工況未產(chǎn)生明顯的擾動區(qū)，所以，其地表擾動區(qū)跨徑接近于0。同時，其8倍平均粒徑開孔造成的擾動區(qū)跨徑增加速度明顯大于5倍平均粒徑開孔的。當計算步數(shù)大于20萬步時，擾動區(qū)跨徑增加速度隨著開孔增大而增大。

3)隨著裂縫尺寸的增加，地表沉降值與地表沉降速度隨之加大，當計算步數(shù)為80萬步時，8倍平均粒徑、5倍平均粒徑開孔造成的地表沉降速度分別為0.015、0.013 m/萬步，而3倍平均粒徑開孔造成的地表沉降速度只有0.006 m/萬步。當3倍平均粒徑開孔產(chǎn)生地表沉降最大值達0.5 m時，8倍平均粒徑開孔產(chǎn)生地表沉降最大值接近1.2 m。管道通常以漸進的方式破裂，這意味著土體侵蝕不會在管道開裂早期造成嚴重的沉降問題，但會隨著管道滲漏的不斷擴大而加速。如果在早期可以檢測到管道破裂問題，則可以避免地面塌陷和經(jīng)濟損失。

圖6 不同管道裂縫尺寸下土顆粒位移云圖(單位：m)Fig.6 Cloud map of soil particle displacement under

圖7 不同管道裂縫尺寸下擾動區(qū)域包絡線與地表沉降曲線Fig.7 Envelope and surface subsidence curve of disturbed area under different crack sizes of

3.2 不同管道裂縫位置的影響

將尺寸為5倍平均粒徑的裂縫分別設置在管道的頂部、右側(cè)腰部和底部，可以得到不同裂縫位置下的影響。圖8、圖9顯示了在計算步數(shù)為80萬步時，管道上不同位置發(fā)生滲漏時的土體位移和地面沉降情況，得出結(jié)論如下：

1)對于土顆粒流失量而言，當管道右側(cè)腰部發(fā)生滲漏時，土體擾動區(qū)和地表沉降槽均向右移動，但擾動區(qū)范圍和地面沉降槽寬度仍與頂部發(fā)生滲漏情況相似。當滲漏發(fā)生在管道底部時，由于土壓力和水的重力作用，土顆粒流失較少，土體擾動主要發(fā)生在管道下面，并沒有延伸到地面，所以，地表沉降幾乎沒有發(fā)生。

2)對于擾動區(qū)跨徑增加速度而言，當計算步數(shù)小于40萬步時，管道頂部開孔造成的擾動區(qū)跨徑增加速度明顯大于腰部開孔的，但當計算步數(shù)大于40萬步時，腰部開孔造成的擾動區(qū)跨徑增加速度是最大的。

圖8 地層中不同管道裂縫位置下土顆粒位移云圖(單位：m)Fig.8 Cloud map of soil particle displacement under different pipeline crack positions in

3)對于地表沉降值與地表沉降速度而言，管道頂部開孔造成的地表沉降最大值大于腰部開孔，大于底部開孔，分別為1.0、0.6、0.3 m，其地表沉降速度亦是如此。由此可見，管道底部產(chǎn)生的裂縫引發(fā)的擾動區(qū)域和地表沉降值較小，而頂部和腰部產(chǎn)生裂縫會導致更大的地表沉降和擾動區(qū)域跨度。

通過采用Fluent和PFC聯(lián)合計算方法對復合地層進行了數(shù)值模擬，結(jié)合不同地下管道裂縫尺寸和不同地下管道破裂位置的影響，得出地層變形和地面塌陷的規(guī)律與王帥超[4]和張成平等[13]文章中的模型試驗結(jié)果相似。

4 結(jié)論

基于計算流體動力學(CFD)和離散元(DEM)理論，提出一種Fluent和PFC聯(lián)合計算數(shù)值模擬方法模擬上海地區(qū)混合地層中管道滲漏引起的工程問題。采用PFC程序?qū)ν翆舆M行建模，用FLUENT程序?qū)Φ叵滤鲌鲞M行計算，然后導入PFC中進行聯(lián)合計算。在建立數(shù)值模型的基礎(chǔ)上，研究了管道不同裂縫尺寸和不同裂縫位置對地面沉降的影響，得到以下結(jié)論：

1)采用的計算方法能夠合理模擬地層中地下水問題，進一步證實了土體滲透侵蝕破壞過程。這種方法可以處理大型的、具有復雜邊界和幾何斷面的工程情況。

2)當管道發(fā)生滲漏時，在水力梯度作用下，土顆粒被沖進管道，土層在裂縫周圍首先產(chǎn)生擾動區(qū)，并逐漸向上延伸，最后形成路面塌陷。隨著擾動區(qū)跨徑的增大，路面塌陷程度也隨之增大。

3)隨著裂縫尺寸的增加，土體流失增加，擾動區(qū)跨徑擴大速度、地表沉降值與地表沉降速度都增加，尤其是當裂縫尺寸為8倍平均粒徑開孔時，其沉降槽范圍是裂縫尺寸為3倍平均粒徑時的2.3倍左右，地表沉降最大值約是裂縫尺寸為3倍平均粒徑時的2.5倍。

4)管道頂部開裂造成的地表沉降最大值和地表沉降速度均大于腰部開孔和底部開孔的地表沉降最大值和地表沉降速度。在計算初期(小于40萬步)3種裂縫位置對擾動區(qū)跨徑影響程度依次為：頂部開孔>腰部開孔>底部開孔；隨計算時間增長(大于40萬步)，裂縫腰部開孔的擾動區(qū)影響范圍將超過頂部開孔造成的擾動區(qū)影響范圍。而底部開孔的裂縫，土體擾動主要發(fā)生在管道下面，并沒有延伸到地面，地表沉降幾乎沒有發(fā)生。因此，管道底部產(chǎn)生的裂縫引發(fā)的擾動區(qū)域和地表沉降值較小，而頂部和腰部產(chǎn)生裂縫會導致更大的地表沉降和擾動區(qū)域跨度。