王越林，陸燁

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

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

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

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

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

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

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

1.1 流體對(duì)土顆粒的作用力

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

fb=-

(1)

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

(2)

1.2 計(jì)算流程

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

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

2 土體與管道模型的建立

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

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

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

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

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

圖2 土體與管道模型

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

(3)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

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

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

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