蔡袁強(qiáng)，張志祥，曹志剛，嚴(yán)舒豪

?

不均勻級配砂土滲蝕過程的細(xì)觀數(shù)值模擬

蔡袁強(qiáng)1, 2，張志祥1，曹志剛1，嚴(yán)舒豪1

(1. 浙江大學(xué) 建筑工程學(xué)院，浙江 杭州，310000；2. 浙江工業(yè)大學(xué) 建筑工程學(xué)院，浙江 杭州，310000)

為了從細(xì)觀角度研究不均勻級配砂土發(fā)生滲蝕的機(jī)理，采用離散元(DEM)與計(jì)算流體動力學(xué)(CFD)耦合模擬的方法，依據(jù)Kézdi判別準(zhǔn)則，對9組不同顆粒級配的土樣滲蝕過程進(jìn)行數(shù)值模擬，分析內(nèi)部穩(wěn)定土體與易滲蝕土體在滲蝕過程中細(xì)觀參數(shù)的變化，并探討初始流速對滲蝕發(fā)展過程的影響。研究結(jié)果表明：土體的滲蝕過程可以分為起始、發(fā)展、過渡、延續(xù)4個階段；土體的顆粒級配直接影響滲蝕的發(fā)展，間隙比越大，土體內(nèi)部越不穩(wěn)定，容易被侵蝕，當(dāng)間隙比大于4時，土體滲蝕起始階段時間短，發(fā)展階段滲蝕速率大，延續(xù)階段滲蝕率仍有穩(wěn)定增長趨勢；在相同間隙比的情況下，細(xì)粒質(zhì)量分?jǐn)?shù)不影響內(nèi)部穩(wěn)定土體滲蝕曲線的起始階段和發(fā)展階段，但細(xì)粒越多延續(xù)階段的滲蝕率越大；流體邊界條件是重要的外部因素，對于易滲蝕土體，超過臨界流速，延續(xù)階段的滲蝕率仍保持增長趨勢，而對于內(nèi)部穩(wěn)定土體，流速對延續(xù)階段滲蝕率逐漸穩(wěn)定的趨勢無影響。

滲蝕；顆粒級配；間隙比；細(xì)粒質(zhì)量分?jǐn)?shù)；流速

長期滲流導(dǎo)致的土體內(nèi)部侵蝕已成為土石壩和土質(zhì)邊坡破壞的最常見原因之一[1]。在滲流力的作用下，細(xì)顆?？朔ｉg摩擦力并進(jìn)一步被水流拖動形成一條條孔隙通道，長此以往，土骨架間被部分掏空，由此引起土體局部沉降以及應(yīng)力重新分布，給工程帶來巨大的隱患[2]。此前，人們通過單元體模型研究滲蝕機(jī)理。滲蝕發(fā)生的原因可以歸結(jié)為內(nèi)部因素和外部因 素[3]。內(nèi)部因素是土體發(fā)生滲蝕的根本原因，包括顆粒級配分布、孔隙分布、顆粒形狀和孔隙形狀。人們已提出多種方法來評估土體的內(nèi)部穩(wěn)定性[4?7]。MORACI等[8?9]通過向上或向下的一維滲流試驗(yàn)分析顆粒級配分布的影響；MAROT等[10]驗(yàn)證了粗顆粒的形狀對土體滲蝕過程的影響。同時，外部因素對于滲蝕的觸發(fā)也很關(guān)鍵。TERZAGHI等[11]提出了內(nèi)部穩(wěn)定土體臨界水力梯度的理論方程；SKEMPTON等[12]認(rèn)為，內(nèi)部不穩(wěn)定土體的臨界水力梯度約為理論值的1/5~1/3；MOFFAT等[13]證實(shí)了土體的有效應(yīng)力是滲蝕發(fā)展的重要因素，并提出了應(yīng)力梯度空間中的水力?力學(xué)路徑。LIANG等[14]設(shè)計(jì)了一種新的水力?力學(xué)耦合滲蝕試驗(yàn)，研究各向同性和各向異性應(yīng)力條件下的滲蝕過程。此外，水力加載歷史對無黏性土的滲蝕有重要影響[15]。由于土體滲蝕發(fā)生在微觀孔隙或顆粒層面，從宏觀的角度來看，試驗(yàn)方法存在一定局限性，主要以經(jīng)驗(yàn)分析為主。隨著計(jì)算機(jī)技術(shù)的發(fā)展，流固耦合數(shù)值模擬方法層出不窮，并不斷被應(yīng)用到滲蝕研究中，人們得以從微觀角度研究無黏性土的滲蝕發(fā)生機(jī)理。ZOU等[16]將計(jì)算流體動力學(xué)結(jié)合顆粒流方法(CFDPFM)用于模擬細(xì)顆粒組被水流侵蝕到粗顆粒組的瞬態(tài)過程；CUI等[17]介紹了一種二維離散元與格子玻爾茲曼法耦合的計(jì)算方法(DEMLBM)，研究了地下管道局部泄漏導(dǎo)致的土體滲蝕；WANG等[18]進(jìn)一步提出了三維黏結(jié)顆粒耦合格子玻爾茲曼法(BPLBM)，研究了土壩中土體滲蝕的過程；常利營等[19]采用三維顆粒流程序(CFDEM)模擬了不同土體之前的接觸沖刷現(xiàn)象，從細(xì)觀角度分析了接觸沖刷的發(fā)生和發(fā)展過程。上述研究表明：離散元算法(DEM)結(jié)合計(jì)算流體動力學(xué)程序(CFD)可用于模擬土體的滲流現(xiàn)象，并且對于水流作用下細(xì)顆粒在粗顆?？紫吨幸苿拥那闆r模擬效果較好。目前，從細(xì)微觀層面對不同級配土體的滲蝕發(fā)生機(jī)理的研究較少，因此，本文采用PFC3D軟件結(jié)合多孔介質(zhì)下的達(dá)西流方程實(shí)現(xiàn)流固耦合計(jì)算，從細(xì)觀角度模擬級配不均勻土體的滲流侵蝕現(xiàn)象，分析不同土體內(nèi)部結(jié)構(gòu)下滲蝕的發(fā)生機(jī)制與發(fā)展規(guī)律。

1 數(shù)值計(jì)算原理

本文采用TSUJI等[20]提出的平均體積粗網(wǎng)格法計(jì)算流體和顆粒的相互作用。該方法中，流體區(qū)域被劃分為多個結(jié)構(gòu)化網(wǎng)格，根據(jù)流體邊界條件，每個流體網(wǎng)格賦予包含在該網(wǎng)格內(nèi)的顆粒以體積力。孔隙率和流體拖曳力通過每個流體單元中顆粒屬性的平均值計(jì)算。在三維顆粒流程序(PFC3D)中所作的假設(shè)如下：1) 流體單元邊長大于PFC顆粒直徑，理想的計(jì)算條件是流體單元不小于顆粒平均直徑的3倍；2) 流體特性在流體單元上是分段呈線性的；3) 流體單元不移動。對一定范圍內(nèi)的孔隙率和雷諾數(shù)，流體固體之間相互作用力的計(jì)算方程是準(zhǔn)確且連續(xù)的，同時紊流的影響也包含在流固耦合相中[21]。

1.1 流固耦合方程

在流固耦合問題中，顆粒的運(yùn)動仍然滿足牛頓第二運(yùn)動定律，考慮到流體對固體的作用，顆粒除了受到機(jī)械外力、顆粒間相互作用力，還受到額外的流體作用力，此時顆粒的運(yùn)動方程為

式中：為顆粒速度矢量；為顆粒質(zhì)量；mech為外力和顆粒間接觸力的合力；fluid為流體作用在顆粒上的合力，包括拖曳力以及由流體壓力梯度引起的力；為重力加速度；為顆粒角速度矢量；為顆粒的轉(zhuǎn)動慣量；為作用在顆粒上的合力矩。

根據(jù)包含顆粒的流體單元邊界條件可知，流體施加在顆粒上的拖曳力drag是針對每個顆粒單獨(dú)定 義的：

式中：0為單一顆粒上的拖曳力；為顆粒所在流體單元的孔隙率，為考慮局部孔隙率的經(jīng)驗(yàn)系數(shù)，此修正項(xiàng)可使式(3)適用于高孔隙率系統(tǒng)和低孔隙率系統(tǒng)以及大范圍的雷諾數(shù)[22?23]。單一顆粒上的拖曳力0的表達(dá)式為

式中：d為拖曳系數(shù)；f為流體密度；為顆粒半徑；為流體速度矢量；其中拖曳系數(shù)d表達(dá)式為

式中：p為顆粒的雷諾數(shù)。

經(jīng)驗(yàn)系數(shù)根據(jù)雷諾數(shù)p計(jì)算得到：

顆粒的雷諾數(shù)可由下式計(jì)算：

式中：f為流體的動力黏度。

作用在流體上每單位體積的體積力b可以表示為

最終流體作用在顆粒上的合力fluid可以表示為

1.2 多孔介質(zhì)中的達(dá)西流

經(jīng)典流體力學(xué)中，低雷諾數(shù)的多孔介質(zhì)流動由達(dá)西定律描述[21]，其方程為

式中：為土體模型的滲透系數(shù)。

一般假設(shè)流體是不可壓縮的，則有

對式(10)取兩側(cè)的散度可得：

將式(12)代入式(11)可得

土體模型的滲透系數(shù)與孔隙率之間的關(guān)系通過Kozeny-Carman方程描述：

CFD模塊自動將流體?顆粒相互作用力作用于PFC顆粒，通過不斷更新PFC模型中的孔隙率和滲透率信息以及CFD模塊中的流體速度場來完成雙向耦合，重新計(jì)算得到具有新的孔隙率和滲透率信息的穩(wěn)態(tài)流體流場[24]。PFC?CFD流固耦合計(jì)算流程見圖1。

圖1 PFC?CFD流固耦合計(jì)算流程圖

2 滲蝕數(shù)值模型

圖2 不均勻顆粒滲蝕數(shù)值模型

在實(shí)際情況下，發(fā)生滲蝕的土體的細(xì)顆粒較小。在模擬三維顆粒模型的條件下，為保證計(jì)算機(jī)的計(jì)算效率，并且能在合理的時間內(nèi)獲得有效的模擬效果。在數(shù)值模型建立的過程中，對實(shí)際的砂土顆粒直徑放大倍，那么，在滲蝕模型尺寸相同的條件下，要保證模型包含足夠多的顆粒來消除顆粒的尺寸效應(yīng)，當(dāng)顆粒尺寸擴(kuò)大倍后，在相同的流體流動條件下，與顆粒尺寸相關(guān)的雷諾數(shù)就會發(fā)生改變。由于雷諾數(shù)發(fā)生改變，研究區(qū)域就會出現(xiàn)紊流，而實(shí)際的滲流試驗(yàn)以達(dá)西流為主。當(dāng)土顆粒和流體發(fā)生相對運(yùn)動且為達(dá)西流時，重力場下的力學(xué)特性和1.0重力場下的力學(xué)特性相似，基于此，當(dāng)顆粒尺寸發(fā)生改變時，考慮與離心試驗(yàn)類似的相似原理[26]，則式(7)的原型方程為

式(7)的模型方程為

代入原型方程，得：

考慮原型與模型的雷諾數(shù)以及流體與顆粒之間的相對速度相等，則有

將式(14)代入式(13)得

對式(16)，其原型方程為

式(16)的模型方程為

考慮原型與模型的孔隙率以及流體壓力梯度相等，則有

因此，當(dāng)顆粒直徑放大倍時，流體的密度相應(yīng)放大倍，動力黏度放大2倍。

為了研究顆粒級配曲線對土體內(nèi)部穩(wěn)定性的影響，在數(shù)值模擬中生成9組不同間隙比(粗顆粒最小直徑與細(xì)顆粒最大直徑之比)和不同細(xì)粒質(zhì)量分?jǐn)?shù)的不均勻顆粒模型，分別編號為G2-10，G2-20，G2-30，G3-10，G3-20，G3-30，G4-10，G5-10和G6-10，顆粒分布曲線如圖3所示。滲蝕模型級配參數(shù)見表1。土體初始孔隙率均為0.5。依據(jù)Kézdi判別準(zhǔn)則，G2和G3組為內(nèi)部穩(wěn)定土體，G4組為臨界穩(wěn)定土體，G5和G6組為易滲蝕土體。本文采用的流體加載方式為一維加載，即在模型左側(cè)施加流體速度inlet，右側(cè)邊界出口壓力設(shè)為0 Pa。為了能讓細(xì)顆粒流出，墻體右側(cè)采用類似過濾網(wǎng)形式墻體，其孔徑略小于粗顆粒的最小直徑。表2所示為滲蝕數(shù)值模型各個計(jì)算單元的物理參數(shù)。由于Kézdi判別準(zhǔn)則評估土體內(nèi)部穩(wěn)定性只考慮了土體的顆粒級配曲線，未考慮顆粒的形狀以及材料的各向異性，因此，本文的滲蝕模擬沒有賦予顆粒重力場，構(gòu)建各向同性應(yīng)力狀態(tài)的模型[25]。

1—G2-10；2—G2-20；3—G2-30；4—G3-10；5—G3-20；6—G3-30；7—G4-10；8—G5-10；9—G6-10。

3 數(shù)值模擬結(jié)果

3.1 滲蝕過程中細(xì)顆粒流失質(zhì)量隨時間的變化

在模型左側(cè)邊界施加固定流體速度0.5 mm/s，計(jì)算時間為120 s，若細(xì)顆粒被沖出右側(cè)邊界，則被自動刪除，將流失的細(xì)顆粒數(shù)與原有細(xì)顆粒數(shù)之比定義為滲蝕率er。圖4(a)所示為G4-10組土樣在滲蝕發(fā)展過程中滲蝕率隨時間的變化。由圖4(a)可知：該滲蝕曲線可以分為4個部分，對應(yīng)滲蝕發(fā)展的4個階段即起始階段(A)、發(fā)展階段(B)、過渡階段(C)和延續(xù)階段(D)，其中，，和為關(guān)鍵時間點(diǎn)，這與文獻(xiàn)[3]中滲蝕模型試驗(yàn)結(jié)果較吻合(見圖4(b))。

表1 滲蝕數(shù)值模型級配參數(shù)

圖5和圖6所示分別為5組不同間隙比土體在滲蝕發(fā)展過程中滲蝕率和細(xì)顆粒流失速率隨時間的變化。從圖5和圖6可以看出：土體的間隙比越大，滲蝕起始階段的時間越短，在加速發(fā)展的過程中，滲蝕速率越大，細(xì)顆粒流失質(zhì)量也越多；G6-10和G5-10這2組易滲蝕土體達(dá)到最大滲蝕速率的時間點(diǎn)最早，且其在發(fā)展階段的滲蝕速率比其他組的明顯更大；內(nèi)部穩(wěn)定土體與易滲蝕土體的滲蝕曲線在延續(xù)階段體現(xiàn)出顯著差異，G6-10和G5-10組的滲蝕率仍有較大的發(fā)展趨勢，而G2-10和G3-10組的滲蝕率逐漸趨于穩(wěn)定。

圖7所示為在間隙比相同而細(xì)粒質(zhì)量分?jǐn)?shù)不同情況下的內(nèi)部穩(wěn)定土體滲蝕曲線。從圖7可以看出：間隙比相同的內(nèi)部穩(wěn)定土體，滲蝕曲線的起始階段和發(fā)展階段幾乎重合，而從過渡階段起，不同細(xì)粒質(zhì)量分?jǐn)?shù)土體的滲蝕曲線有所偏離，表現(xiàn)為細(xì)粒質(zhì)量分?jǐn)?shù)越大，延續(xù)階段的滲蝕率越大。

表2 滲蝕數(shù)值模型計(jì)算參數(shù)

(a) G4-10組土樣滲蝕率隨時間的變化(數(shù)值模擬結(jié)果)；(b) 被侵蝕出的細(xì)顆粒質(zhì)量隨時間的變化(室內(nèi)試驗(yàn)結(jié)果[3])

1—G2-10；2—G3-10；3—G4-10；4—G5-10；5—G6-10。

1—G2-10；2—G3-10；3—G4-10；4—G5-10；5—G6-10。

1—G2-10；2—G2-20；3—G2-30；4—G3-10；5—G3-20；6—G3-30。

3.2 數(shù)值試樣顆粒滲蝕過程

圖8所示為G6-10組土樣在滲蝕過程中顆粒位移云圖。根據(jù)滲蝕發(fā)展曲線可知：4 s時滲蝕開始啟動，局部細(xì)顆粒有向右移動的趨勢，6 s時，滲蝕開始加速發(fā)展，在加速發(fā)展過程中，內(nèi)部細(xì)顆粒穿過粗顆粒的孔隙并緩慢移動，逐漸形成一條條滲蝕通道；從20 s到40 s，越來越多的細(xì)顆粒聚集在大顆粒間的孔隙中，此時，左側(cè)的顆粒繼續(xù)穿過粗顆?？紫队鲎?，因此顆粒的滲蝕速率慢慢降低，為滲蝕過渡階段；40 s后進(jìn)入滲蝕的延續(xù)階段，易滲蝕土體的細(xì)顆粒由于間隙比大，仍能從粗顆?？紫吨辛鞒觯€(wěn)定土體的細(xì)顆粒則逐漸堵在粗顆粒周圍，趨于穩(wěn)定。

時間/s: (a) 6；(b) 20；(c) 40；(d) 120

3.3 滲蝕過程中細(xì)觀參數(shù)的變化

5組土樣在滲蝕過程中細(xì)顆粒平均位移隨時間的變化如圖9所示。由圖9可見：施加滲流場初期細(xì)顆粒位移快速增加，土體的間隙比越大，細(xì)顆粒的平均位移越大。隨著時間的推移，內(nèi)部穩(wěn)定土體中細(xì)顆粒的平均位移逐漸趨于穩(wěn)定，而易滲蝕土體中細(xì)顆粒平均位移在延續(xù)階段仍有增大趨勢。細(xì)顆粒平均速度隨時間的變化如圖10所示。在施加滲流場后，細(xì)顆粒在水流的拖曳下獲得某一速度，由于粗顆粒的阻礙，越來越多細(xì)顆粒撞上粗顆粒并聚集在粗顆粒周圍，因此，細(xì)顆粒的平均速度在滲蝕發(fā)展階段快速下降并趨于穩(wěn)定；而間隙比較大的土體，細(xì)顆粒更容易穿過粗顆?？紫?，形成孔隙通道。由圖10可知：到延續(xù)階段后，G2-10和G3-10這2組內(nèi)部穩(wěn)定土體細(xì)顆粒平均速度趨于0 mm/s，而G6-10和G5-10這2組易滲蝕土體平均速度仍保持某一常數(shù)。

1—G2-10；2—G3-10；3—G4-10；4—G5-10；5—G6-10。

3.4 滲蝕過程中分區(qū)參數(shù)隨時間的變化

為了更好地說明滲蝕過程中細(xì)顆粒的移動，G6-10組模型從左向右被分為5等份[27]，模型分區(qū)示意圖如圖11所示。圖12所示為G6-10組土樣模型各分區(qū)細(xì)顆粒數(shù)隨時間的變化曲線。將最初各分區(qū)的顆粒數(shù)歸零，顆粒數(shù)為負(fù)說明流失的顆粒多于流入的顆粒，顆粒數(shù)為正說明流入的顆粒多于流出的顆粒。從圖12可以看出：最左側(cè)I區(qū)的細(xì)顆粒由于沒有顆粒補(bǔ)充，所以，顆粒數(shù)一直在減小；由于有I區(qū)細(xì)顆粒的流入，II區(qū)的顆粒數(shù)在滲蝕前期基本保持不變，但是，隨著時間的推移，細(xì)顆?！叭氩环蟪觥保瑪?shù)量逐漸減少；III區(qū)和IV區(qū)流入和流失的顆粒數(shù)相當(dāng)，總細(xì)顆粒數(shù)基本保持不變，但是，隨著滲蝕的繼續(xù)發(fā)展，可以推測III區(qū)和IV區(qū)的顆粒數(shù)也會減?。籚區(qū)由于靠近出口邊界，細(xì)顆粒最容易被流體帶出，因此，滲蝕起始和發(fā)展階段流失的細(xì)顆粒數(shù)遠(yuǎn)遠(yuǎn)多于流入的顆粒數(shù)，過渡階段的細(xì)顆粒數(shù)反而有所增加，延續(xù)階段細(xì)顆粒數(shù)略微減少。

1—G2-10；2—G3-10；3—G4-10；4—G5-10；5—G6-10。

圖11 G6-10組土樣模型分區(qū)示意圖

3.5 不同初始流速下滲蝕率隨時間的變化

為了分析初始流速對滲蝕過程的影響，分別選取G6-10組易滲蝕土體和G3-10組內(nèi)部穩(wěn)定土體，在左側(cè)邊界施加不同的初始流速。圖13所示為不同初始流速下土體滲蝕曲線。從圖13(a)可以看出：易滲蝕土體能否發(fā)生持續(xù)滲蝕還取決于初始流速，初始流速越大，滲蝕起始階段的時間越短，發(fā)展階段的滲蝕速率越大，延續(xù)階段的開始時間點(diǎn)越遲，滲蝕率越大。同時，當(dāng)流速大于0.3 mm/s時，滲蝕率在延續(xù)階段仍有增長的趨勢，因此，可以認(rèn)為在相同顆粒級配下，存在某一臨界水流速度，超過該臨界水流速度，易滲蝕土體將發(fā)生持續(xù)滲蝕。而對于圖13(b)中的內(nèi)部穩(wěn)定土體，由于受顆粒級配的限制，流速并不影響延續(xù)階段滲蝕率趨于恒定，即在目前的研究階段下，內(nèi)部穩(wěn)定土體不存在臨界水流速度。

1—I區(qū)；2—II區(qū)；3—III區(qū)；4—IV區(qū)；5—V區(qū)。

初始流速/(mm?s?1)：1—0.1；2—0.2；3—0.3；4—0.4；5—0.5；6—0.3；7—0.5；8—0.7；9—0.9；10—1.1。

4 結(jié)論

1) 土體的間隙比越大，內(nèi)部越不穩(wěn)定，越容易被侵蝕，細(xì)顆粒流失量越多；易滲蝕土體的滲蝕曲線中起始階段時間較短，在加速發(fā)展的過程中滲蝕速率大，延續(xù)階段的滲蝕率仍有較大的增長趨勢，而內(nèi)部穩(wěn)定土體的滲蝕率逐漸趨于穩(wěn)定。

2) 內(nèi)部穩(wěn)定土體在間隙比相同的情況下，細(xì)顆粒質(zhì)量分?jǐn)?shù)不影響滲蝕曲線的起始階段和發(fā)展階段，而在延續(xù)階段，細(xì)粒質(zhì)量分?jǐn)?shù)越大，滲蝕率越大。

3) 初始流速是影響滲蝕發(fā)生的重要外部因素，流速越大，滲蝕曲線的起始階段越短，發(fā)展階段滲蝕速率越大，滲蝕率越大。對于易滲蝕土體，當(dāng)超過某一臨界水流速度時，易滲蝕土體將發(fā)生持續(xù)滲蝕；而對于內(nèi)部穩(wěn)定土體，流速并不影響延續(xù)階段的滲蝕率趨于穩(wěn)定。

[1] LUO Yulong, QIAO Liang, LIU Xingxing, et al. Hydro- mechanical experiments on suffusion under long-term large hydraulic heads[J]. Natural Hazards, 2013, 65(3): 1361?1377.

[2] CHANG Dongsheng. Internal erosion and overtopping erosion of earth dams and landslide dams[D]. Hong Kong: The Hong Kong University of Science and Technology. Department of Civil and Environment Engineering, 2012: 1?3.

[3] KAZUKI H, AKIHIRO T. Suffusion-induced change in spatial distribution of fine fractions in embankment subjected to seepage flow[J]. Soils and Foundations, 2015, 55(5): 1293?1304.

[4] KHAKSAR N E, ESLAMI A. Assessment of the likelihood of suffusion in alluvial soils: case history[J]. Bulletin of Engineering Geology and the Environment, 2015, 74(2): 611? 620.

[5] KEZDI A. Soil physics: selected topics (developments in geotechnical engineering)[M]. Amsterdam, Netherlands: Elsevier, 1979: 1?160.

[6] KENNEY T C, LAU D. Internal stability of granular filters[J]. Canadian Geotechnical Journal, 1985, 22(2): 215–225.

[7] KENNEY T C, LAU D. Internal stability of granular filters: reply[J]. Canadian Geotechnical Journal, 1986, 23(3): 420?423.

[8] MORACI N, MANDAGLIO M C, IELO D. Analysis of the internal stability of granular soils using different methods[J]. Canadian Geotechnical Journal, 2014, 51(9): 1063?1072.

[9] YANG K S, WANG J Y. E xperiment and statistical assessment on piping failures in soils with different gradations[J]. Marine Georesources & Geotechnology, 2017, 35(4): 512?527.

[10] MAROT D, BENDAHMANE F, BGUYEN H H. Influence of angularity of coarse fraction grains on internal erosion process [C]//Paris, France: ICSE6, 2012: 887?894.

[11] TERZAGHI K, PECK R B. Soil mechanics in engineering practice[M]. 2nd ed. London, UK: John Wiley and Sons, 1967: 149?150.

[12] SKEMPTON A W, BROGAN J M. Experiments on piping in sandy gravels[J]. Géotechnique, 1994, 44(3): 449?460.

[13] MOFFAT R, FANNIN R J. A hydromechanical relation governing internal stability of cohesionless soil[J]. Canadian Geotechnical Journal, 2011, 48(3): 413?424.

[14] LIANG Yue, YEH T C J, WANG Junjie. Onset of suffusion in upward seepage under isotropic and anisotropic stress conditions [J]. European Journal of Environmental and Civil Engineering, 2017: 1?15.

[15] ROCHIM A, MAROT D, SIBILLE L, et al. Effects of hydraulic loading history on suffusion susceptibility of cohesionless soils [J]. Journal of Geotechnical and Geoenvironmental Engineering, 2017, 143(7): 04017025-1?25.

[16] ZOU Yuhua, CHEN Qun, CHEN Xiaoqing, et al. Discrete numerical modeling of particle transport in granular filters[J]. Computers & Geotechnics, 2013, 47: 48?56.

[17] CUI Xilin, LI Jun, Chan A, et al. Coupled DEM-LBM simulation of internal fluidisation induced by a leaking pipe[J]. Powder Technology, 2014, 254(2): 299?306.

[18] WANG M, Feng Y T, Pand e G N, et al. Numerical modelling of fluid-induced soil erosion in granular filters using a coupled bonded particle lattice Boltzmann method[J]. Computers & Geotechnics, 2017, 82: 134?143.

[19] 常利營, 陳群, 葉發(fā)明. 均勻顆粒間接觸沖刷的顆粒流數(shù)值模擬[J]. 巖土工程學(xué)報(bào)，2016, 38(增刊2): 312?317. CHANG Liying, CHEN Qun, YE Faming. Particle flow simulation for contact erosion between uniform particles[J]. Chinese Journal of Geotechnical Engineering, 2016, 38(Suppl 2): 312?317.

[20] TSUJI Y, KAWAGUCHI T, TANATA T. Discrete particle simulation of two-dimensional fluidized bed[J]. Powder Technology, 1993, 77(1): 79?87.

[21] PFC 5.0 documentation[M]. Minneapolis, USA: Itasca Consulting Group Inc., 2016: 1?2.

[22] Felice R D. The voidage function for fluid-particle interaction systems[J]. International Journal of Multiphase Flow, 1994, 20(1): 153?159.

[23] Xu B H, Yu A B. Numerical simulation of the gas-solid flow in a fluidized bed by combining discrete particle method with computational fluid dynamics[J]. Chemical Engineering Science, 1997, 52(16): 2785?2809.

[24] Guyer J E, Wheeler D,Warren J A. FiPy: partial differential equations with python[J]. Computing in Science & Engineering, 2009, 11(3): 6?15.

[25] Shire T, O’SULLIVAN C. Micromechanical assessment of an internal stability criterion[J]. Acta Geotechnica, 2013, 8: 81?90.

[26] 劉先珊, 張林, 秦鵬偉. 考慮顆粒尺寸效應(yīng)的油井出砂數(shù)值模擬[C]//顆粒材料計(jì)算力學(xué)會議. 蘭州, 2014: 257?284. LIU Xianshan, ZHANG Lin, QIN Pengwei. Numerical simulation of sand production in oil well considering particle size effect[C]//Conference of Computational Mechanics of Granular Materials. Lanzhou, 2014: 257?284.

[27] TAO Hui, TAO Junliang. Quantitative analysis of piping erosion micro-mechanisms with coupled CFD and DEM method[J]. Acta Geotechnica, 2017, 12(3): 573?592.

Mesoscopicnumerical simulation for suffusion process of gap-graded sandy soil

CAI Yuanqiang1, 2, ZHANG Zhixiang1, CAO Zhigang1, YAN Shuhao1

(1. School of Civil Engineering and Architecture, Zhejiang University, Hangzhou 310000, China; 2. School of Civil Engineering and Architecture, Zhejiang University of Technology, Hangzhou 310000, China)

To investigate the meso-mechanism of suffusion between gap-graded sandy soil, the discrete element method coupled with computational fluid dynamics was used to simulate the suffusion process of nine groups of graded soils based on Kézdi’s criteria. The changes of the meso-parameters of the internal stable soil and the unstable soil during the erosion process were analyzed, and the influence of the initial flow rate on the erosion process was discussed. The results show that the eroded process of soil can be divided into four stages, i.e., initiation, development, transition and progression. The particle size distribution of soils directly affects the development of erosion. The larger the gap ratio is, the less stable the soil is, and it is easily eroded. When the gap ratio is greater than 4, the initiation stage of soil suffusion is short, the erosion rate in the development stage is large and has a steady growth trend in the progression stage. At the same gap ratio, the mass fraction of fine particles rarely affects the initial stage and development stage of the erosive curve, but the more the fine particles are, the greater the erosion rate is at the continuation stage. The fluid boundary condition is an important external factor. For easily erosive soils, when water flow rates exceed the critical value, the erosion rate still maintains an increasing trend in the progression stage. For the internal stable soil, the flow rate hardly affects the stable tend of the erosion rate in the progression stage.

suffusion; particle size distribution; gap ratio; mass fraction of fine particles; flow velocity

TU462

A

1672?7207(2019)05?1144?10

10.11817/j.issn.1672?7207.2019.05.018

2018?06?22；

2018?08?22

國家重點(diǎn)研發(fā)計(jì)劃項(xiàng)目(2016YFC0800203); 國家自然科學(xué)基金資助項(xiàng)目(51578500，51778571) (Project (2016YFC0800203) supported by the National Key Research and Development Program of China; Projects(51578500，51778571) supported by the National Natural Science Foundation of China)

曹志剛，博士，副教授，從事土動力學(xué)、路基工程和隧道工程研究；E-mail：caozhigang2011@zju.edu.cn

(編輯 伍錦花)