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

王越林　陸燁

摘 要：近年來，城市地區(qū)地面沉降和塌方事故頻繁發(fā)生，管道滲漏是導(dǎo)致此類事故的主要因素之一。當(dāng)管道破裂時，管道內(nèi)流體為重力流，管道周圍土體會被沖蝕進(jìn)管道，最終導(dǎo)致沉降或地面塌陷。為探討管道滲漏時土體在水流作用下的遷移過程，基于計算流體動力學(xué)（CFD）和離散元（DEM）理論，提出一種Fluent和PFC聯(lián)合計算數(shù)值模擬方法。采用顆粒流程序（PFC3D）對地層進(jìn)行建模，用ANASYS中的FLUENT模塊計算地下水流場并導(dǎo)入PFC3D，進(jìn)行聯(lián)合計算。采用此流固耦合方法，針對上海北部地區(qū)地層情況，即上部1.5 m厚黏性土層、下部2.5 m厚砂性土層進(jìn)行數(shù)值建模，研究地下管道裂縫尺寸和地下管道裂縫位置對管道滲漏引發(fā)的地下空洞區(qū)發(fā)展及地面塌陷的影響。數(shù)值模擬結(jié)果表明：裂縫尺寸越大，土體流失量越大，砂性地層內(nèi)土體擾動區(qū)域范圍越廣，由于上部黏性土層的存在，地表沉降值不大;土體擾動區(qū)域范圍以裂縫位置為起點(diǎn)不斷擴(kuò)展，越往上部延伸，擾動區(qū)域范圍越大，且擾動區(qū)域范圍遠(yuǎn)遠(yuǎn)大于空洞區(qū)范圍。

關(guān)鍵詞：離散元;地面塌陷;地面沉降;管道滲漏;數(shù)值模擬

中圖分類號：TU434? ?文獻(xiàn)標(biāo)志碼：A ??文章編號：2096-6717（2021）02-0060-08

Abstract：In recent years， land subsidence and collapse accidents occur frequently in urban areas， and the leakage of pipelines is one of the main factors leading to such accidents. When the pipeline breaks and the fluid in the pipeline is of gravity flow， the soil around the pipeline will be eroded into the pipeline， which will eventually lead to settlement or ground collapse. To explore the migration process of soil under action of water flow during pipeline leakage， the authors propose a Fluent and PFC combined computational numerical simulation method based on computational fluid dynamics （CFD） and discrete element （DEM） theory. This method mainly uses particle flow program （PFC3D） to model the formation. The FLUENT module in the ANASYS is used to calculate the groundwater flow field， which will be imported into the PFC3D to conduct the joint computation. Using this flow-solid coupling method， this paper conducts a numerical modeling of the stratigraphic situation in the northern part of Shanghai， and this domain is the upper 1.5 m thick clay layer and the lower 2.5 m thick sandy soil layer. Then the influences of the crack size and the crack position of underground pipeline on the development of underground cavity area caused by pipeline leakage and the ground collapse are explored. The numerical simulation results show that the larger the crack size is， the greater the soil loss is， the wider the soil disturbance area is in sandy strata， and the surface settlement will not demonstrate large value due to the existence of the upper clay soil layer. The soil disturbance area starts from the crack location and continues to expand. The higher it extends， the larger the disturbance area is， and the disturbance area is much larger than the cavity area.

Keywords：discrete element; ground collapse; land subsidence; pipe leakage; numerical simulation

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

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

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

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

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

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

1.1 流體對土顆粒的作用力

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

1.2 計算流程

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

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)合計算，并探究砂性土與黏性土地層中地面塌陷的產(chǎn)生和發(fā)展規(guī)律。

模型中模型槽尺寸為10 m×6 m×4 m（長×寬×高），在模型槽內(nèi)自下而上按線性接觸模型分別生成高度為2.5 m的砂性土層和1.5 m的黏性土層，其孔隙率均為0.41，并對所有土顆粒施加重力，在重力作用下顆粒沉積、固結(jié)，并達(dá)到平衡狀態(tài)。最后在黏性土層上部覆蓋0.3 m厚的混凝土，代表低等級路面。在土顆粒的生成中，如果按照實(shí)際土顆粒尺寸在模型槽中生成，其數(shù)目將達(dá)到億級，這對計算機(jī)來說是不現(xiàn)實(shí)的，因此，對10 m×6 m×4 m的模型進(jìn)行初步研究，結(jié)果表明，管道頂部滲漏引起的顆粒位移場是軸對稱的。在滿足極限計算能力的前提下，將模型槽尺寸更改為10 m×1 m×4 m，同時，對土顆粒尺寸進(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。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

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

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

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

參考文獻(xiàn)：

[1] 尚彥軍， 龐文利， 王開洋， 等. 太古代片麻狀花崗巖淺埋隧道變形破壞及地面塌陷分析： 以集寧隧道為例[J]. 工程地質(zhì)學(xué)報， 2018， 26（6）： 1574-1584.

SHANG Y J， PANG W L， WANG K Y， et al. Deformation damage of shallow tunnels and caveinduced ground collapse in Archeozoic gneissic granites： the Jining tunnel as case example [J]. Journal of Engineering Geology， 2018， 26（6）： 1574-1584. （in Chinese）

[2] 宋許根， 劉秀敏， 陳從新， 等. 程潮鐵礦西區(qū)采空區(qū)地表塌陷機(jī)制與變形規(guī)律初探[J]. 巖石力學(xué)與工程學(xué)報， 2018， 37（Sup2）： 4262-4273.

SONG X G， LIU X M， CHEN C X， et al. Study on the mechanism of surface collapse and ground deformation in western mined area of Chengchao iron mine [J]. Chinese Journal of Rock Mechanics and Engineering， 2018， 37（Sup2）： 4262-4273. （in Chinese）

[3] GUO S， SHAO Y， ZHANG T Q， et al. Physical modeling on sand erosion around defective sewer pipes under the influence of groundwater [J]. Journal of Hydraulic Engineering， 2013， 139（12）： 1247-1257.

[4] 王帥超. 城市地下管道滲漏引起的路面塌陷機(jī)理分析與研究[D]. 鄭州： 鄭州大學(xué)， 2017.

WANG S C. Analysis and study on subsidence mechanisms of road caused by leakage of urban underground pipeline [D]. Zhengzhou： Zhengzhou University， 2017.（in Chinese）

[5] EL SHAMY U， ZEGHAL M. Coupled continuum-discrete model for saturated granular soils [J]. Journal of Engineering Mechanics， 2005， 131（4）： 413-426.

[6] TAO H， TAO J L. Quantitative analysis of piping erosion micro-mechanisms with coupled CFD and DEM method [J]. Acta Geotechnica， 2017， 12（3）： 573-592.

[7] 鄭剛， 戴軒， 張曉雙. 地下工程漏水漏砂災(zāi)害發(fā)展過程的試驗(yàn)研究及數(shù)值模擬[J]. 巖石力學(xué)與工程學(xué)報， 2014， 33（12）： 2458-2471.

ZHENG G， DAI X， ZHANG X S. Experimental study and numerical simulation of leaking process of sand and water in underground engineering [J]. Chinese Journal of Rock Mechanics and Engineering， 2014， 33（12）： 2458-2471. （in Chinese）

[8] 周健， 周凱敏， 姚志雄， 等. 砂土管涌-濾層防治的離散元數(shù)值模擬[J]. 水利學(xué)報， 2010， 41（1）： 17-24.

ZHOU J， ZHOU K M， YAO Z X， et al. Numerical simulation of piping-filter prevention in sandy soil by discrete element method [J]. Journal of Hydraulic Engineering， 2010， 41（1）： 17-24. （in Chinese）

[9] 蔣明鏡， 張望城. 一種考慮流體狀態(tài)方程的土體CFD-DEM耦合數(shù)值方法[J]. 巖土工程學(xué)報， 2014， 36（5）： 793-801.

JIANG M J， ZHANG W C. Coupled CFD-DEM method for soils incorporating equation of state for liquid [J]. Chinese Journal of Geotechnical Engineering， 2014， 36（5）： 793-801. （in Chinese）

[10] 李曉蛟， 陸燁， 武亞軍. 一種CFD-DEM流固耦合方法在滲流導(dǎo)致 城市地面沉降問題中的應(yīng)用[J]. 上海大學(xué)學(xué)報（自然科學(xué)版）， 2020， 26（5）： 842-852.

LI X J， LU Y， WU Y J. The application of CFD-DEM coupling method in seepage-induced ground subsidence [J]. Journal of Shanghai University， 2020， 26（5）： 842-852. （in Chinese）

[11] 羅勇， 龔曉南， 吳瑞潛. 顆粒流模擬和流體與顆粒相互作用分析[J]. 浙江大學(xué)學(xué)報（工學(xué)版）， 2007， 41（11）： 1932-1936.

LUO Y， GONG X N， WU R Q. Analysis and simulation of fluid-particles interaction with particle flow code [J]. Journal of Zhejiang University（Engineering Science）， 2007， 41（11）： 1932-1936. （in Chinese）

[12] 陳建峰， 李輝利， 周健. 黏性土宏細(xì)觀參數(shù)相關(guān)性研究[J]. 力學(xué)季刊， 2010， 31（2）： 304-309.

CHEN J F， LI H L， ZHOU J. Study on the relevance of macro-micro parameters for clays [J]. Chinese Quarterly of Mechanics， 2010， 31（2）： 304-309. （in Chinese）

[13] 張成平， 岳躍敬， 蔡義. 管線滲漏水范圍對淺埋隧道圍巖變形和破壞的影響[J]. 巖石力學(xué)與工程學(xué)報， 2015， 34（2）： 392-400.

ZHANG C P， YUE Y J， CAI Y. Influence of pipeline leakage range on ground deformation and failure during shallow tunnelling [J]. Chinese Journal of Rock Mechanics and Engineering， 2015， 34（2）： 392-400. （in Chinese）

[14] 張?zhí)悾?周愛國， 馮小銘， 等. 南京市地面塌陷發(fā)育特征及防治對策[J]. 中國安全科學(xué)學(xué)報， 2011， 21（3）： 3-8.

ZHANG T L， ZHOU A G， FENG X M， et al. Study on distribution characteristics of ground collapse in Nanjing City and its countermeasures [J]. China Safety Science Journal， 2011， 21（3）： 3-8. （in Chinese）

[15] 周健， 姚志雄， 張剛. 砂土滲流過程的細(xì)觀數(shù)值模擬[J]. 巖土工程學(xué)報， 2007， 29（7）： 977-981.

ZHOU J， YAO Z X， ZHANG G.Mesomechanical simulation of seepage flow in sandy soil [J]. Chinese Journal of Geotechnical Engineering， 2007， 29（7）： 977-981. （in Chinese）

[16] CUNDALL P A. A computer model for simulating progressive large-scale movements in blocky rock systems [C]// ISRM Symposium， 1971： 129-136.

[17] ZHANG D M， MA L X， ZHANG J， et al. Ground and tunnel responses induced by partial leakage in saturated clay with anisotropic permeability [J]. Engineering Geology， 2015， 189： 104-115.

[18] MAIR R J.Tunnelling and geotechnics： new horizons [J]. Géotechnique， 2008， 58（9）： 695-736.

[19] RAVICHANDRAN N， MACHMER B， KRISHNAPILLAI H， et al. Micro-scale modeling of saturated sandy soil behavior subjected to cyclic loading [J]. Soil Dynamics and Earthquake Engineering， 2010， 30（11）： 1212-1225.

[20] CUNDALL P A， STRACK O D L. A discrete numerical model for granular assemblies [J].Géotechnique， 1979， 29（1）： 47-65.

[21] 周健， 池永. 土的工程力學(xué)性質(zhì)的顆粒流模擬[J]. 固體力學(xué)學(xué)報， 2004， 25（4）： 377-382.

ZHOU J， CHI Y. Simulating soil properties by particle flow code [J]. Acta Mechanica Solida Sinica， 2004， 25（4）： 377-382. （in Chinese）

[22] 倪小東， 趙帥龍， 王媛， 等. 巖體水力劈裂的細(xì)觀PFC-CFD聯(lián)合分析[J]. 巖石力學(xué)與工程學(xué)報， 2015， 34（Sup2）： 3862-3870.

NI X D， ZHAO S L， WANG Y， et al. Numerical analysis of hydraulic fracturing of rock mass on mesoscopic level by coupled PFC-CFD method [J]. Chinese Journal of Rock Mechanics and Engineering， 2015， 34（Sup2）： 3862-3870. （in Chinese）

[23] 羅勇. 土工問題的顆粒流數(shù)值模擬及應(yīng)用研究[D]. 杭州： 浙江大學(xué)， 2007.

LUO Y. Simulation of soil mechanical behaviors using discrete element method based on particle flow code and it application [D]. Hangzhou： Zhejiang University， 2007.（in Chinese）

（編輯 王秀玲）