游碧波，周翠英

(中山大學(xué)工學(xué)院∥巖土工程與信息技術(shù)研究中心，廣東 廣州 510275)

管涌問(wèn)題是水利工程界普遍關(guān)注的重要研究課題。根據(jù)國(guó)內(nèi)外的研究現(xiàn)狀，管涌現(xiàn)象的機(jī)理研究主要分為理論解析模型研究、室內(nèi)外實(shí)驗(yàn)研究以及數(shù)值模擬等方面。由于管涌問(wèn)題涉及到水與土顆粒的相互作用，其機(jī)理十分復(fù)雜，傳統(tǒng)的理論模型和實(shí)驗(yàn)研究均存在某些局限[1]。隨著計(jì)算機(jī)技術(shù)的不斷發(fā)展和各種數(shù)值模擬方法的應(yīng)用，以數(shù)值分析技術(shù)模擬滲透變形的發(fā)生及其發(fā)展過(guò)程，成為研究滲透變形問(wèn)題的一個(gè)主要方向。根據(jù)現(xiàn)有數(shù)值模擬成果，以有限元為代表的連續(xù)介質(zhì)方法在管涌研究中得到了廣泛的應(yīng)用[2～5]。采用連續(xù)介質(zhì)方法模擬管涌時(shí)，基本是以滲流場(chǎng)的模擬為基礎(chǔ)，將離散的顆?？闯墒沁B續(xù)的流體加以研究。連續(xù)介質(zhì)方法依賴(lài)高度簡(jiǎn)化的、規(guī)定性質(zhì)的本構(gòu)方程而對(duì)顆粒物性參數(shù)、粒徑、形狀及其分布等對(duì)管涌演化的影響則難以考慮；并且管涌區(qū)域的參數(shù)在不斷變化，連續(xù)介質(zhì)方法對(duì)管涌通道及管涌通道過(guò)渡區(qū)的網(wǎng)格劃分問(wèn)題難以處理，因此只能模擬穩(wěn)定滲流場(chǎng)的狀況。管涌的實(shí)質(zhì)是細(xì)顆粒從粗顆粒組成的骨架孔隙中流失，土體顆粒和水流相互作用的力學(xué)機(jī)制十分復(fù)雜，連續(xù)介質(zhì)方法對(duì)反映堤基管涌這類(lèi)顆粒運(yùn)移引起的模型和參數(shù)不斷調(diào)整的復(fù)雜問(wèn)題，以及實(shí)現(xiàn)對(duì)管涌發(fā)展演化過(guò)程的模擬等均存在一定的局限性。

離散單元法與連續(xù)介質(zhì)方法相比，在分析具有離散性質(zhì)的物質(zhì)方面表現(xiàn)出了優(yōu)越性。其基本思想是將不連續(xù)體分離成若干離散單元的集合，單元之間無(wú)須滿足位移連續(xù)與變形協(xié)調(diào)條件，運(yùn)動(dòng)方程的迭代采用顯式的中心差分格式，速度快，所需空間小，適用于求解大變形和非連續(xù)問(wèn)題[6]。顆粒元方法作為離散元的一種，是專(zhuān)門(mén)用來(lái)求解固體力學(xué)大變形和顆粒流動(dòng)問(wèn)題的數(shù)值方法，它通過(guò)接觸方程和運(yùn)動(dòng)方程將顆粒之間的復(fù)雜的相互作用由計(jì)算機(jī)來(lái)完成，因而可以對(duì)土顆粒相互之間的細(xì)觀作用形象化的表現(xiàn)出來(lái)。目前，顆粒元在巖土、礦冶、農(nóng)業(yè)、食品、化工、制藥和環(huán)境等領(lǐng)域均已廣泛地應(yīng)用[7]。對(duì)于水-土相互作用問(wèn)題，目前已有學(xué)者開(kāi)始采用顆粒元方法對(duì)離散顆粒與水的相互作用進(jìn)行分析研究[8]，也有學(xué)者開(kāi)始采用顆粒流程序針對(duì)管涌試驗(yàn)進(jìn)行數(shù)值模擬分析，并得出了一些有價(jià)值的研究成果[9]?？梢钥闯?，以顆粒元為代表的離散元方法，為在細(xì)觀尺度上解釋和分析管涌現(xiàn)象提供了一條途徑。

管涌發(fā)生時(shí)，管涌出口附近地層的顆粒流失將引起管涌的迅速發(fā)展，導(dǎo)致管涌通道的形成，甚至發(fā)生堤基破壞。對(duì)管涌口附近地層的顆粒流失規(guī)律、流速分布的分析研究，弄清管涌口附近顆粒流失的規(guī)律，對(duì)于分析險(xiǎn)情、指導(dǎo)搶險(xiǎn)，具有重要價(jià)值，對(duì)于進(jìn)一步研究管涌演化規(guī)律具有重要意義。但在管涌口出現(xiàn)時(shí)，第一要?jiǎng)?wù)是搶險(xiǎn)堵漏，因而難以對(duì)管涌口附近地層的相關(guān)細(xì)觀參數(shù)進(jìn)行測(cè)量和現(xiàn)場(chǎng)分析，此方面目前尚未見(jiàn)相關(guān)研究。本文擬通過(guò)顆粒流程序平臺(tái)PFC2D(Particle Flow Code in Two Dimension)，采用流-固耦合的分析模型對(duì)管涌口附近的顆粒流失現(xiàn)象進(jìn)行模擬，并通過(guò)FISH語(yǔ)言開(kāi)發(fā)程序，對(duì)顆粒流失量、孔隙率、流速等計(jì)算結(jié)果進(jìn)行分析。

1 顆粒流基本理論

顆粒流是基于離散元的原理提出的專(zhuān)門(mén)用于模擬圓形顆粒介質(zhì)的運(yùn)動(dòng)及其相互作用的一種數(shù)值方法[10]。它采用圓盤(pán)模型模擬顆粒，通過(guò)接觸方程和運(yùn)動(dòng)方程將顆粒之間的復(fù)雜的相互作用由計(jì)算機(jī)來(lái)完成，它除了細(xì)觀模擬上的優(yōu)勢(shì)之外，隨著力的平衡過(guò)程，土體的宏觀運(yùn)動(dòng)也可以通過(guò)土顆粒的微觀運(yùn)動(dòng)表示出來(lái)。

1.1 基本假設(shè)

根據(jù)文獻(xiàn)[10]中對(duì)顆粒流方法的相關(guān)闡述，采用顆粒流進(jìn)行模擬時(shí)具有如下基本假設(shè)：①顆粒單元為剛性體；②接觸發(fā)生在很小的范圍內(nèi)，即點(diǎn)接觸；③接觸特性為柔性接觸，接觸處允許有一定的“重疊”量；④“重疊”量的大小與接觸力有關(guān)，與顆粒大小相比，“重疊”量很??；⑤接觸處有特殊的粘結(jié)強(qiáng)度；⑥顆粒單元為圓盤(pán)形(或球形)。

1.2 計(jì)算循環(huán)過(guò)程

顆粒流方法以牛頓第二定律與力—位移定律為基礎(chǔ)，對(duì)模型進(jìn)行循環(huán)計(jì)算。圖1給出了計(jì)算過(guò)程循環(huán)圖。

圖1 顆粒流計(jì)算過(guò)程循環(huán)圖(改自文[10])

PFC2D的計(jì)算過(guò)程可分為以下4步：①接觸檢索；②利用接觸本構(gòu)關(guān)系計(jì)算作用于顆粒上的不平衡力與不平衡力矩；③顆粒滿足運(yùn)動(dòng)方程(即牛頓第二運(yùn)動(dòng)定律)；④更新顆粒位置，進(jìn)行下一步迭代。

1.3 接觸模型

對(duì)材料本構(gòu)特性的模擬，PFC2D是通過(guò)接觸本構(gòu)模型來(lái)實(shí)現(xiàn)。每一顆粒的接觸本構(gòu)有：①接觸剛度模型；②滑動(dòng)模型；③粘結(jié)模型。接觸剛度模型提供了接觸力和相對(duì)位移的彈性關(guān)系；滑動(dòng)模型則強(qiáng)調(diào)切向和法向接觸力使得接觸顆?？梢园l(fā)生相對(duì)移動(dòng)；粘結(jié)模型是限制總的切向和法向力使得顆粒在粘結(jié)強(qiáng)度范圍內(nèi)發(fā)生接觸，包括接觸連接和平行連接[11]。本文的模型在模擬砂性土顆粒之間接觸作用時(shí)采用接觸剛度模型，黏性土顆粒之間采用粘結(jié)模型中的平行粘結(jié)模型。粘結(jié)模型能同時(shí)傳遞力和力矩，發(fā)生在接觸顆粒間圓形或方形有限范圍內(nèi)，可以描述顆粒之間有限范圍內(nèi)有粘接作用特性，由法向剛度、切向剛度、法向強(qiáng)度、切向強(qiáng)度和連接半徑等參數(shù)進(jìn)行定義，相對(duì)應(yīng)的作用力以總接觸力與力矩表示。

2 管涌口數(shù)值模型的建立

2.1 生成模型

北江大堤位于廣東省北部，是防御北江和西江洪水侵襲廣州和珠江三角洲地區(qū)的重要屏障，是廣東省最重要的一條堤防，也是全國(guó)確保的七大堤防之一。北江大堤石角段歷史上曾多次決堤出險(xiǎn)，是北江大堤的主要險(xiǎn)段，吸引了一些學(xué)者對(duì)此段堤基的管涌問(wèn)題開(kāi)展相關(guān)研究[12-13]。石角堤段的遙堤在0+650斷面樂(lè)排河閘附近未采取任何防滲措施，堤身、堤基滲透性較大，一旦該堤擋水將存在發(fā)生滲透變形的危險(xiǎn)[14](其地層斷面見(jiàn)圖2)。本文基于此管涌危險(xiǎn)堤段的地層條件，建立模型，模擬發(fā)生管涌時(shí)逸出口處的顆粒流失現(xiàn)象。

首先選取模擬對(duì)象的范圍，考慮計(jì)算規(guī)模和效率的關(guān)系，集中選取該堤基管涌口周?chē)欢ǚ秶鷥?nèi)的顆粒進(jìn)行建模。由圖2可以看出，堤內(nèi)地層結(jié)構(gòu)為典型二元堤基，覆蓋層為粉質(zhì)黏土，下面是中粗砂層。據(jù)此建立數(shù)值模型尺寸為10 m×6 m，其中覆蓋粉質(zhì)黏土層厚度為1 m，下臥透水中粗砂層厚度為5 m，模型在左、右、底側(cè)分別設(shè)置PFC2D透水邊界墻作為邊界條件，以模擬實(shí)際的水流邊界。

在該尺寸的模型中，如果按實(shí)際土顆粒的大小模擬，得到的顆粒數(shù)量遠(yuǎn)遠(yuǎn)超過(guò)百萬(wàn)甚至千萬(wàn)，目前的計(jì)算機(jī)和軟件還無(wú)法滿足要求。這里，我們關(guān)心的是管涌發(fā)展的演化機(jī)制，而不是為了模擬堤基土的顆粒數(shù)量，所以在模擬過(guò)程中對(duì)土顆粒的尺寸進(jìn)行了同比放大，同時(shí)將堤基簡(jiǎn)化為雙層計(jì)算，這樣即克服了計(jì)算機(jī)容量和速度的限制，又能滿足計(jì)算精度的要求，同時(shí)，對(duì)模擬該處管涌的規(guī)律也有一定的參考意義。顆粒試樣尺寸分別為：粉質(zhì)黏土采用粒徑為0.015 m的顆粒，中粗采用粒徑為0.05 m的顆粒，編制FISH函數(shù)由PFC2D隨機(jī)生成器生成。顆粒生成的具體過(guò)程如下。

首先在墻體所在范圍內(nèi)充滿流體，在墻體范圍內(nèi)采用自由下落法生成顆粒，定義顆粒的相關(guān)參數(shù)和重力加速度，首先生成下層的砂土顆粒，讓顆粒在重力作用下達(dá)到平衡狀態(tài)，這樣可以消除內(nèi)部的不平衡力；待顆粒靜止后，刪除多余顆粒，再用同樣方法生成上層黏土顆粒，此時(shí)土體狀態(tài)即為模型初始狀態(tài)。這種顆粒生成方法的優(yōu)點(diǎn)是可以保證試樣生成過(guò)程中顆粒通過(guò)相互作用不斷調(diào)整，較易達(dá)到自然受力平衡；最后在黏土層的中間部位，刪除寬度為0.1 m的顆粒帶，以模擬管涌口。最終生成的模型顆粒總數(shù)約為12 500個(gè)，模型如圖3所示。

圖2 石角遙堤(樂(lè)排河閘)0+650斷面圖

圖3 管涌口數(shù)值模型

2.2 參數(shù)的選取

砂土顆粒采用接觸剛度模型，黏土顆粒采用平行粘結(jié)接觸模型。本文的數(shù)值模擬強(qiáng)調(diào)的是顆粒的運(yùn)移規(guī)律與實(shí)際情況相近。

因此，在建立實(shí)際工程模型之前，先構(gòu)造并運(yùn)行簡(jiǎn)化的測(cè)試模型，并設(shè)置流體，在左側(cè)和下側(cè)施加水頭，模擬在實(shí)際的水力邊界條件，顆粒在水壓力下發(fā)生起動(dòng)，反復(fù)調(diào)整顆粒相關(guān)參數(shù)，使得數(shù)值模型與實(shí)際管涌具有大致相同的顆粒起動(dòng)條件。粉質(zhì)黏土覆蓋層參數(shù)的確定參考前人試驗(yàn)中通過(guò)建立模型進(jìn)行粘接力試驗(yàn)確定，這樣，可保證其粘結(jié)剛度與實(shí)際情況接近。顆粒相關(guān)參數(shù)如表1所示。

表1 顆粒流模型細(xì)觀參數(shù)表

需要指出的是數(shù)值模型采用的各項(xiàng)參數(shù)，如剛度、摩擦系數(shù)、平行粘結(jié)的各項(xiàng)指標(biāo)并非真實(shí)堤基土的物理指標(biāo)，而僅僅是顆粒流數(shù)值模擬的自定義指標(biāo)，與真實(shí)土層物理指標(biāo)沒(méi)有直接的關(guān)系，而且迄今還沒(méi)有能從細(xì)觀力學(xué)指標(biāo)反應(yīng)宏觀力學(xué)指標(biāo)的關(guān)系公式。本文參數(shù)取值遵循顆粒流數(shù)值模擬的基本途徑：首先構(gòu)建并運(yùn)行簡(jiǎn)化模型，這樣可對(duì)力學(xué)細(xì)觀的概念有深入的了解，通過(guò)這種前期簡(jiǎn)化模型的運(yùn)行，可得計(jì)算模型的幾何參數(shù)及初始條件，分析簡(jiǎn)化模型的結(jié)果，將其與實(shí)測(cè)結(jié)果或宏觀現(xiàn)象進(jìn)行對(duì)比，大致確定出各參數(shù)的取值范圍，一般還需加以修改，多次試算即可得出計(jì)算參數(shù)的取值。

3 計(jì)算結(jié)果分析

在模擬墻體上施加壓力差，形成向右向上的水流壓力，在程序中設(shè)置測(cè)量圈和內(nèi)置函數(shù)來(lái)求取顆粒的孔隙率和滲透率指標(biāo)，測(cè)量圈的設(shè)置如圖4所示。并編寫(xiě)了試樣各個(gè)部分的流體cell孔隙率的FISH函數(shù)，用來(lái)記錄管涌口附近孔隙率及水流流速的變化規(guī)律，以下對(duì)模擬的結(jié)果進(jìn)行詳細(xì)討論和分析。

圖4 測(cè)量圈的設(shè)置

3.1 孔隙率的變化

圖5記錄了管涌演化過(guò)程中各個(gè)測(cè)量圈內(nèi)的孔隙率變化情況。

圖5 各測(cè)量圈內(nèi)顆?？紫堵首兓€圖

對(duì)圖5中的各條曲線進(jìn)行分析可知，孔隙率總體上均呈逐漸增大的趨勢(shì)，并且隨著顆粒流失結(jié)束而趨于穩(wěn)定，該結(jié)果與和相關(guān)研究的結(jié)論是一致的[14]。管涌口兩側(cè)靠近覆蓋層的測(cè)量圈1和3的孔隙率變化基本呈加速上升的態(tài)勢(shì)，這是由于在管涌發(fā)展過(guò)程中，此區(qū)域形成了滲流集中的通道，導(dǎo)致顆粒加速流失。曲線1比曲線3上升得更快說(shuō)明左側(cè)通道發(fā)展比右側(cè)更為迅速，這顯然是受到左側(cè)水壓力大的影響。

靠近管涌口附近的測(cè)量圈2的孔隙率變化非常明顯，其對(duì)應(yīng)的曲線2可劃分為兩個(gè)階段：首先是穩(wěn)定發(fā)展階段，孔隙率和流速呈平緩上升；然后是迅速演化階段，當(dāng)管涌演化到一定階段時(shí)(12萬(wàn)步左右)，發(fā)展為加速流失階段，曲線迅速抬升，這表明此時(shí)顆粒迅速流失。

曲線4在開(kāi)始階段逐漸上升，孔隙率由0.2增加到0.3左右，到20萬(wàn)步之后，又有減小的趨勢(shì)。結(jié)合測(cè)量圈1的孔隙率變化曲線可以推測(cè)，曲線4的這種起伏變化是由于管涌顆粒流失后，在地層內(nèi)形成了細(xì)小的通道，在測(cè)量圈1內(nèi)的通道受到不斷變大的水流壓力的影響，最后演變成主要的滲流通道，因而孔隙率迅速增長(zhǎng)，而測(cè)量圈4靠近通道，其孔隙率受到通道顆粒運(yùn)移時(shí)壓力的作用，因而孔隙率減小，密實(shí)度有所增加。

曲線5在剛開(kāi)始上升較快，在15萬(wàn)步后，曲線呈水平狀態(tài)，孔隙率保持在0.27～0.28左右，在25萬(wàn)步之后，曲線又出現(xiàn)上升。分析其原因，測(cè)量圈5位于管涌口下方，離管涌口有一定的距離，剛開(kāi)始時(shí)，受到水壓力的作用，顆粒雖然不能從管涌口流失，但顆粒與流體發(fā)生復(fù)雜的相互作用，顆粒結(jié)構(gòu)發(fā)生了調(diào)整，變得較為松散，直到管涌口附近顆粒流失到一定程度時(shí)，此處的顆粒便開(kāi)始流失。

測(cè)量圈6內(nèi)的孔隙率變化過(guò)程(對(duì)應(yīng)曲線6)較為復(fù)雜，有一定的起伏，曲線最后階段才有一定的上升。這表明此位置的顆粒受到擾動(dòng)，土體變疏松后又變得密實(shí)，這可能也是由于此處復(fù)雜的水—土相互作用造成的。

3.2 孔隙流速的變化

在模擬計(jì)算過(guò)程中，在各個(gè)測(cè)量圈內(nèi)對(duì)流體網(wǎng)格內(nèi)各節(jié)點(diǎn)上的流速進(jìn)行平均，得到了管涌演化全過(guò)程中各部位的平均孔隙水流速度變化曲線，如圖6所示。

圖6 各測(cè)量圈內(nèi)孔隙水流速變化曲線圖

通過(guò)圖6中各測(cè)量圈的流速變化曲線可看出，除曲線2上升顯著外，其它曲線的變化都較為平緩；曲線2在12萬(wàn)步左右，流速一直呈波動(dòng)增長(zhǎng)的趨勢(shì)， 這正好與圖5的孔隙率圖反應(yīng)出的顆粒流失情況相對(duì)應(yīng)；曲線1、3曲線存在明顯的波動(dòng)，特別是曲線1在15萬(wàn)步之后，波動(dòng)非常明顯，并呈持續(xù)上升，這反應(yīng)了此處流體和顆粒復(fù)雜的相互作用，從下文的圖7可知，這是由于此處形成了管涌通道所造成；由于測(cè)量圈4受到測(cè)量圈1范圍內(nèi)的管涌通道的影響，在15～20萬(wàn)步時(shí)，曲線4存在較明顯波動(dòng)；曲線5、6基本呈一水平線，這說(shuō)明管涌發(fā)生后，距離管涌口較遠(yuǎn)的位置，其孔隙流速不會(huì)發(fā)生大的變化。

結(jié)合以上對(duì)于孔隙率和流速的變化曲線的分析，可以對(duì)雙層堤基管涌出口附近的演化過(guò)程進(jìn)行以下的歸納闡述：管涌口下方滲透地層的顆粒在管涌發(fā)展過(guò)程中是沿著實(shí)際水流速度最大的通道溢出管涌口，管涌的發(fā)生僅會(huì)影響管涌口附近一定范圍內(nèi)的流速重新分布，且流速分布情況會(huì)有一個(gè)從相等的平均流速狀態(tài)向局部流速集中并逐漸增大的改變過(guò)程。管涌初始階段試樣內(nèi)部各處的流速相差不大，在水壓力作用下，靠近管涌口的地層，其顆粒開(kāi)始流失，并導(dǎo)致此處流速的增大，流速的增大又促進(jìn)了此處管涌顆粒的加速流失，這種水與顆粒之間的相互作用始終貫穿于模擬的整個(gè)過(guò)程。管涌通道的形成，會(huì)導(dǎo)致通道之外的地層受到水壓力的作用，局部地層有可能變得密實(shí)，但這只是暫時(shí)現(xiàn)象，顆粒的流失最終會(huì)導(dǎo)致整體地層密實(shí)度的下降。

3.3 顆粒流失的規(guī)律

圖7記錄了管涌口顆粒逸出到管涌地層內(nèi)部通道形成的全過(guò)程。

圖7 管涌口PFC試樣的演化

從圖7可將整個(gè)模擬過(guò)程劃分為管涌初期(顆粒啟動(dòng)，開(kāi)始流失)、管涌中期(管涌口擴(kuò)大，顆粒流失加快)和管涌后期(劇烈流失，形成通道)三個(gè)階段。

0～15萬(wàn)時(shí)步為第一階段。管涌口附近的砂性土顆粒由于受到水壓力的作用開(kāi)始流出。下臥層中，中下部顆粒的位移受到上部顆粒的限制，基本未發(fā)生移動(dòng)，但由于土顆粒的不規(guī)則移動(dòng)形成了一些孔隙，這些孔隙呈分散狀分布在透水層的下半層。管涌繼續(xù)發(fā)展，覆蓋層下面靠近管涌口位置的顆粒開(kāi)始流失，出現(xiàn)了明顯的接觸沖刷。

15～30萬(wàn)時(shí)步為第二階段。由于下臥層管涌顆粒的流出以及水流沖刷的作用下，管涌口周?chē)母采w層顆粒的粘結(jié)被破壞，管涌口邊緣的黏土顆粒組強(qiáng)度喪失，拆散成單獨(dú)的顆?；蝾w粒組，并向上運(yùn)動(dòng)直到流失。宏觀上的表現(xiàn)為，隨著內(nèi)部砂土顆粒大量涌出，管涌口開(kāi)始向兩側(cè)擴(kuò)展管涌開(kāi)始加劇；顆粒的加速流失，導(dǎo)致管涌口正下方的下臥層形成了一定深度的空腔，在管涌口附近的接觸層沖刷更加明顯，覆蓋層在重力作用下產(chǎn)生局部的下陷，這也解釋了管涌口周?chē)鷷?huì)發(fā)生地陷的原因。值得關(guān)注的是，下臥砂層內(nèi)部土顆粒的不規(guī)則移動(dòng)現(xiàn)象逐漸消失，這可能是由于水壓力作用下，顆粒產(chǎn)生相互擠壓，而當(dāng)上層砂性土部分流失后，中下層顆粒之間的擠壓應(yīng)力釋放，孔隙分布逐漸變得均勻?？梢酝茰y(cè)，在此演化階段，顆粒的流失，并不會(huì)影響堤基的整體穩(wěn)定性，堤基仍然能正常的承受上部荷載。

30～40萬(wàn)時(shí)步可看作第三階段。內(nèi)部顆粒從左側(cè)(上游)大量流失，并且流失顆粒較集中的位于管涌口下部緊貼覆蓋層下側(cè)的砂性土。覆蓋層下陷更加明顯，管涌口繼續(xù)擴(kuò)大，最終透水層內(nèi)部由于顆粒的流失形成較大的空腔，并向左右擴(kuò)展，形成多條明顯的滲流通道，并且隨著通道的逐漸形成，孔隙不斷增大導(dǎo)致流速的相對(duì)集中，周?chē)念w粒迅速沿通道流失。這一現(xiàn)象也表明，隨著通道的形成，堤基的穩(wěn)定性將越來(lái)越差，隨著通道附近顆粒進(jìn)一步被掏蝕，堤基強(qiáng)度將逐漸喪失，最終導(dǎo)致堤基破壞。

北江大堤管涌的相關(guān)資料表明：剛開(kāi)始發(fā)生管涌時(shí)，管涌口水流較小，多為清水或夾雜少量細(xì)砂，流失到一定程度時(shí)，管涌口處水流會(huì)迅速增大，帶出大量的泥沙[15]。本文的模擬也反映了相同的演化規(guī)律。同時(shí)，根據(jù)北江大堤石角段在94年特大洪水時(shí)的管涌測(cè)量資料得出的結(jié)論：對(duì)于多層堤基條件下透水層發(fā)生管涌前后的對(duì)比發(fā)現(xiàn)，透水地層的密實(shí)度有所下降，結(jié)構(gòu)變得松散，砂層的有效粒徑由管涌前的0.05～0.16，變?yōu)楣苡亢蟮?.15～1.1[16]。本文采用顆粒流模擬管涌時(shí)得到的細(xì)觀參數(shù)的變化規(guī)律，與真實(shí)的管涌在宏觀和細(xì)觀上均具有一致性。

以上基于離散元原理的PFC2D程序，較好的模擬了管涌口顆粒逸出的全過(guò)程。通過(guò)模擬結(jié)果與文獻(xiàn)[14]中的實(shí)際管涌情況對(duì)比可看出：在宏觀方面，堤基管涌從發(fā)生到破壞，顆粒流出、管涌口擴(kuò)展、管涌口下方透水層顆粒流失導(dǎo)致沖刷、管涌口周?chē)耐翆铀?、管涌通道逐漸形成幾個(gè)宏觀過(guò)程得到了形象的模擬，模擬的結(jié)果與實(shí)際管涌的宏觀現(xiàn)象基本吻合；細(xì)觀方面，孔隙率和流速的變化規(guī)律也與實(shí)際情況相對(duì)應(yīng)。

4 結(jié) 語(yǔ)

本文利用離散介質(zhì)的二維顆粒流方法，通過(guò)對(duì)實(shí)際地層進(jìn)行簡(jiǎn)化，建立了基于實(shí)際工程的管涌口演化顆粒流數(shù)值模型，模擬得到了水—土相互作用作用與管涌發(fā)生的全過(guò)程，分析了管涌產(chǎn)生與發(fā)展的細(xì)觀機(jī)制。通過(guò)對(duì)管涌口附近的孔隙率、流速等細(xì)觀參數(shù)的變化規(guī)律進(jìn)行分析，觀察顆粒流失的宏觀現(xiàn)象，得到以下結(jié)論：

1)雙層堤基管涌從發(fā)生到破壞，顆粒流出、管涌口擴(kuò)展、管涌口下方透水層顆粒流失導(dǎo)致沖刷、管涌口周?chē)耐翆铀?、管涌通道逐漸形成幾個(gè)宏觀過(guò)程得到了形象的模擬，模擬的結(jié)果與實(shí)際管涌的宏觀現(xiàn)象基本吻合。

2)管涌口下方滲透地層的顆粒在管涌發(fā)展過(guò)程中是沿著實(shí)際水流速度最大的通道溢出管涌口。管涌口附近流速分布情況會(huì)有一個(gè)從相等的平均流速狀態(tài)向局部流速集中并逐漸增大的改變，并且明顯大于其他區(qū)域。距離管涌口較遠(yuǎn)的地層中，流速不會(huì)發(fā)生大的變化。滲透層的顆粒流失與孔隙流速的增加是相互影響和促進(jìn)的關(guān)系，水與顆粒的相互作用始終貫穿于整個(gè)過(guò)程。顆粒的流失最終會(huì)導(dǎo)致整體地層密實(shí)度的下降。

3)顆粒流比傳統(tǒng)連續(xù)介質(zhì)模擬方法可以更明晰地再現(xiàn)土顆粒的細(xì)觀參數(shù)變化軌跡，并且可以得到實(shí)際管涌過(guò)程中無(wú)法定量測(cè)的管涌發(fā)展過(guò)程中的土樣特性和水流特性等數(shù)據(jù)。在得到管涌演化宏觀現(xiàn)象的同時(shí)，可以通過(guò)分析土顆粒的運(yùn)動(dòng)規(guī)律，分析靜力觸探的細(xì)觀演化規(guī)律，方便地改變應(yīng)力邊界條件，調(diào)節(jié)土的力學(xué)參數(shù)進(jìn)行試驗(yàn)，因而是一種行之有效的方法。

本文僅是對(duì)管涌口的演化規(guī)律進(jìn)行了較為初步計(jì)算和分析，在模擬中，涉及到參數(shù)取值及一些細(xì)部處理的問(wèn)題時(shí)，模擬的結(jié)果在一定程度上缺乏說(shuō)服力，但對(duì)管涌口模擬這一復(fù)雜問(wèn)題而言，基于PFC2D軟件，應(yīng)用顆粒流理論及數(shù)值計(jì)算方法來(lái)模擬卻是有益的嘗試。由于目前的計(jì)算機(jī)水平限制，有些問(wèn)題尚未能開(kāi)展，比如對(duì)管涌過(guò)程中土體內(nèi)部管涌通道的測(cè)量及演化規(guī)律等，還有待開(kāi)展進(jìn)一步的研究。

參考文獻(xiàn)：

[1] 周健, 張剛.管涌現(xiàn)象研究的進(jìn)展與展望[J].地下空間,2004,24(4): 536-542.

[2] 李守德,徐紅娟,田軍. 均質(zhì)土壩管涌發(fā)展過(guò)程的滲流場(chǎng)空間性狀研究[J]. 巖土力學(xué),2005(12)：2001-2004.

[3] 李守德,張曉海,劉志祥. 基坑開(kāi)挖工程管涌發(fā)生過(guò)程的模擬[J]. 工程勘察，2003(2)：14-17,21.

[4] 胡琦,凌道盛,陳仁朋,等. 粉砂地基深基坑工程土體滲透破壞機(jī)理及其影響研究[J]. 巖土力學(xué)，2008(11)：2967-2972.

[5] 丁留謙,吳夢(mèng)喜,劉昌軍,等. 雙層堤基管涌動(dòng)態(tài)發(fā)展的有限元模擬[J]. 水利水電技術(shù)，2007(2)：36-39.

[6] 周先齊,徐衛(wèi)亞,鈕新強(qiáng),等.離散單元法研究進(jìn)展及應(yīng)用綜述[J]. 巖土力學(xué)，2007( S1)：408-416.

[7] 楊洋,唐壽高. 顆粒流的離散元法模擬及其進(jìn)展[J]. 中國(guó)粉體技術(shù)，2006(5)：38-43.

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

[9] 張剛,周健,姚志雄. 堤壩管涌的室內(nèi)實(shí)驗(yàn)與顆粒流細(xì)觀模擬研究[J]. 水文地質(zhì)工程地質(zhì)，2007(6)：83-86.

[10] Itasca Consulting Group Incorporation. PFC2D User’s Guide[R].Minneapolis, Minnesota: [s.n.], 2004.

[11] Itasca Consulting Group Incorporation.PFC2D theory and background[R].Minneapolis, Minnesota: [s.n.], 2004.

[12] 陳明東,陳國(guó)寶. 廣東省北江石角堤段堤基機(jī)械管涌原因分析[J]. 工程地質(zhì)學(xué)報(bào)，1998(4)：28-33.

[13] 曹洪,張挺,陸培炎. 北江大堤石角段強(qiáng)透水堤基滲流分析[J]. 巖石力學(xué)與工程學(xué)報(bào)，2001(6)：855-858.

[14] 周曉杰. 堤防的滲透變形及其發(fā)展的研究[D].北京：清華大學(xué)，2006.

[15] STERPI D. Effects of the erosion and transport of fine particles due to seepage flow[J].International Journal of Geomechanics，2003, 3(1): 111-122.

[16] 林叔忠.北江大堤石角段強(qiáng)透水堤基滲透變形機(jī)理初探[C]∥北江大堤石角段滲流問(wèn)題技術(shù)研討會(huì),1998.