陳聞瀟 石 崇 單治鋼 陳 亮 倪小東 李汪洋

(1.河海大學(xué) 巖石力學(xué)與堤壩工程教育部重點實驗室,南京 210098;2.河海大學(xué) 巖土工程科學(xué)研究所,南京 210098;3.中國電建集團(tuán) 華東勘測設(shè)計研究院有限公司,杭州 310014)

管涌是導(dǎo)致河道堤防、壩基等失穩(wěn)破壞的重要因素[1-2].在滲透水流的作用下,填充在土體骨架中的細(xì)顆粒逐漸流失,管涌通道的直徑將不斷擴(kuò)大,導(dǎo)致堤壩發(fā)生不均勻沉陷,最終導(dǎo)致崩塌[3-5].

近年來,越來越多的學(xué)者對管涌現(xiàn)象展開了研究.理論計算方面:施倩等[6]建立了管涌破壞的數(shù)學(xué)模型,介玉新等[7]研究了管涌破環(huán)的時間過程模擬.室內(nèi)試驗方面:曹洪等[8]研究了管涌通道邊壁松散層對管涌發(fā)展的影響,Mohamed Elkholy等[9]研究了土壤組成對土堤管涌的影響.數(shù)值模擬方面:常利營等[10]研究了砂土滲流規(guī)律的適用性,Y.Guo等[11]研究了不同顆粒形狀對管涌的影響.離散元方法如顆粒流(PFC)[12],因為能夠較好地模擬大變形破壞,而逐漸被較多的學(xué)者采用.周建等[13]通過PFC模擬了砂土管涌的基料-濾層系統(tǒng),從細(xì)觀角度揭示管涌發(fā)展過程中顆粒的運動特性和濾層防治機(jī)理.倪小東等[14]基于PFC對室內(nèi)砂槽試驗進(jìn)行了模擬.張剛等[15]通過PFC從細(xì)觀角度研究了管涌的形成和發(fā)展機(jī)制.但學(xué)者們的模擬大多基于室內(nèi)試驗或小尺度模型,針對現(xiàn)場管涌案例建立模型的離散元分析還有待進(jìn)一步研究.

本文基于PFC模擬管涌現(xiàn)象的發(fā)展,滲流流體的CFD計算通過python中的fipy偏微分方程求解器求解達(dá)西公式來實現(xiàn),通過將PFC-CFD結(jié)合來實現(xiàn)土體-流體的相互作用計算.參考南京長江大堤管涌案例,建立管涌計算模型,分析隨著管涌發(fā)展,滲透流速、水壓力和土體孔隙率的變化規(guī)律和動力學(xué)特性,并研究堤壩位移的變化趨勢.

1 基于PFC-CFD的多孔介質(zhì)滲流計算方法

利用PFC-CFD模塊計算三維條件下的多孔介質(zhì)流動,多孔介質(zhì)中的低雷諾數(shù)流動通?？梢酝ㄟ^達(dá)西定律描述:

式中:v為流體的速度矢量;K為滲透矩陣;μ為流體黏度;ε為孔隙率矩陣;p為流體壓力.通常假定流體的壓縮性很小,可以忽略不計,即認(rèn)為流體不可壓縮:

穩(wěn)態(tài)不可壓縮滲流方程可以通過對方程(1)兩邊同時取散度導(dǎo)出:

將公式(3)代入(2),得泊松方程:

式(4)邊界條件為入口處有

其中:vin為入口速度;出口處有p=0.這個方程通過隱式求解可以很容易得出流體的壓力場.求解方案基于穩(wěn)態(tài)流,即流入量與流出量相等,一旦壓力已知,流體速度可以由方程(1)直接獲得.

流體方程(1)和(2)是在粗流體網(wǎng)格單元[16]上求解.通過計算PFC3D顆粒與流體單元之間的重疊量確定孔隙率.考慮流動受顆粒運動的影響,滲透系數(shù)由PFC3D模型的孔隙率計算.

當(dāng)re為PFC顆粒平均半徑時,第i個粗網(wǎng)格的滲透系數(shù)Ki可以通過其孔隙率εi和Kozeny-Carman關(guān)系[17]來計算:

2 南京長江大堤管涌案例

2016年7月3 日晚,南京長江揚子江段江堤背水坡坡腳出現(xiàn)滲水現(xiàn)象,本文通過建立PFC-CFD計算模型,實現(xiàn)管涌發(fā)展的動態(tài)分析.如圖1所示,出險段長江大堤堤頂高程12.6 m,堤頂寬8.0 m,背水坡堤腳高程約8.5 m,長江下關(guān)最高潮位為9.67 m.如果滲水不斷發(fā)展,滲透水流逐步掏走通道中的土體,通道不斷擴(kuò)大,大堤就會坍塌.因此南京市防辦及相關(guān)部門及時采取應(yīng)急處置,滲水情況得到明顯控制,大堤安全總體可控.本文基于該案例進(jìn)行PFC-CFD方法的管涌過程模擬分析.

圖1 南京長江大堤管涌案例

3 管涌計算模擬

如圖2所示建立邊坡模型,在此基礎(chǔ)上考慮管涌現(xiàn)象的發(fā)展及其對堤壩的穩(wěn)定性影響.

圖2 管涌計算模型

由于耦合方法是基于三維模型,為兼顧問題的三維性,采用假三維模型來進(jìn)行平面問題的模擬,假三維厚度為3 m.模型尺寸如圖2(a)所示,整體模型長50.9 m,高12.4 m,土體的孔隙率在0.36左右.由于計算能力受限,將顆粒粒徑擴(kuò)大10倍,管涌影響區(qū)域內(nèi)細(xì)顆粒和粗顆粒半徑分別為0.04～0.08 m和0.12～0.24 m.管涌影響區(qū)域之外均為0.12～0.24 m半徑的顆粒,為提高計算效率,共生成顆粒63644個.考慮顆粒的尺寸效應(yīng)[18],需要對流體的黏滯系數(shù)擴(kuò)大10倍以滿足相似性原則.

如圖2(b)所示,基于python調(diào)用gmsh軟件建立非均勻六面體網(wǎng)格,共生成3612個網(wǎng)格,并將網(wǎng)格信息導(dǎo)入到python中的fipy偏微分方程求解器中,通過fipy實現(xiàn)公式(1)～(6)的計算.通過PFC中的CFD模塊為顆粒施加流體-顆粒相互作用力,在每一個時間步更新孔隙率和滲透系數(shù),并將數(shù)據(jù)傳遞回python中的fipy求解器中,隨之進(jìn)行下一個時間步的計算.模型采用接觸粘結(jié)模型,該模型能夠較好地模擬土體的力學(xué)響應(yīng)和破壞特性[19].堤防土體細(xì)觀力學(xué)參數(shù)見表1.

表1 計算模型細(xì)觀力學(xué)參數(shù)

管涌計算結(jié)果如圖3所示,如果從細(xì)顆粒逐漸啟動流失到管涌逐漸形成開始模擬,則需要模擬的時間過長難以接受.故從長江潮水位9.67 m開始模擬,此時已經(jīng)具有一定的初始滲透流速.從圖3(a)、3(b)可見,土體中各流體單元孔隙率的不同導(dǎo)致滲透系數(shù)不同,引起了各測點流場中壓力及流速的不均勻.隨著管涌的發(fā)展,各個測點的滲透性也在不斷地增強(qiáng),滲透流速不斷增大.

圖3 管涌模型計算結(jié)果

從圖3(b)可見,管涌由背水側(cè)的細(xì)小顆粒流失而不斷發(fā)展,滲透流速不斷變大,進(jìn)而導(dǎo)致整個滲流通道內(nèi)的土體開始流失.近水側(cè)監(jiān)測點1在15萬～20萬時間步時,孔隙率急劇上升,表明這部分的土體已經(jīng)大部分流失而接近破壞.背水側(cè)土體雖有流失,但同時水流也將近水側(cè)流失顆粒攜帶而來,故孔隙率只在較小幅度內(nèi)變化.從圖3(c)可見前10萬時步時,孔隙率的變化幅度在0.1左右,10萬～40萬時步,孔隙率的變化幅度達(dá)到了0.2～0.4.

從圖3(d)可見,在滲流初期,水壓力隨著滲透路徑呈現(xiàn)近似地均勻分布,隨著管涌發(fā)展,靠近滲流入口處的土體逐漸流失,水壓力上升明顯.從圖3(e)可見,初始各個位置的水壓力梯度是近似均勻的,在2萬～20萬步內(nèi),沿滲透路徑0～10 m處的水壓力梯度逐步下降,是因為該處的土體不斷流失(圖3(c)),而10～20 m處壓力梯度逐步提高.計算至40萬步時,0～15 m處的壓力梯度陡然下降,這是由于此處已達(dá)到土體的臨界啟動速度而發(fā)生土體流失(圖3(c)).此時的壓力梯度集中在了20～30 m處,隨著管涌的進(jìn)一步發(fā)展,這一部分土體也將會達(dá)到臨界啟動速度并流失.

模型對土體的流失量進(jìn)行了測量,圖3(f)為模擬過程中流失量隨時間的變化曲線,20萬時間步時流失量已經(jīng)達(dá)到了2.5 m3,在40萬時間步流失量已經(jīng)迅速增加至9 m3,巨大的流失量表明土體中可能已經(jīng)產(chǎn)生了空洞,這對于整個堤壩的安全是十分不利的.隨著管涌的進(jìn)一步發(fā)展,流失量將會繼續(xù)增大,這可能會導(dǎo)致堤壩發(fā)生破壞.

圖4 模型位移云圖(單位:m)

由圖4位移云圖可以看出,滲流初始階段堤壩位移主要發(fā)生在堤壩內(nèi)部滲流通道附近,在2萬、10萬、20萬、40萬時步下的最大位移分別為0.042、0.068、0.093 6和0.148 m.可見位移是逐漸擴(kuò)大的,主要發(fā)生在堤壩內(nèi)部和堤內(nèi)坡腳處,這對于堤壩土體是十分不利的,如不及時治理將會導(dǎo)致更進(jìn)一步的破壞.

4 結(jié) 論

本文采用PFC-CFD數(shù)值模擬方法,基于南京長江大堤管涌案例,建立了土體管涌的數(shù)值模擬模型,研究了管涌滲透流體變化規(guī)律,分析了土體孔隙度及位移隨管涌發(fā)展的變化趨勢.主要得出以下結(jié)論:

1)采用PFC-CFD的模擬方法可以較好地用于土體管涌現(xiàn)象的模擬,可以較好地呈現(xiàn)滲透水流與土體的相互作用.

2)隨著管涌發(fā)展,靠近滲流入口處的土體因為管涌滲透作用而率先流失,沿滲流路徑上滲透流速隨著管涌發(fā)展而不斷加大,水壓力和壓力梯度出現(xiàn)分布變化,土體流失量迅速增加.

3)堤壩位移主要發(fā)生于堤壩內(nèi)部滲流通道附近,隨著土體不斷流失,堤內(nèi)坡腳處的位移不斷增大,將導(dǎo)致堤壩破壞.