基于PFC-CFD的土體管涌數(shù)值模擬研究

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2 南京長(zhǎng)江大堤管涌案例

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

圖1 南京長(zhǎng)江大堤管涌案例

3 管涌計(jì)算模擬

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

圖2 管涌計(jì)算模型

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

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

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

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

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

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

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

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

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

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

4 結(jié) 論

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

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

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

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