黃海均,嚴(yán)新軍*,毛海濤,王正成

(1.新疆農(nóng)業(yè)大學(xué)水利與土木工程學(xué)院，烏魯木齊 830052；2.河北工程大學(xué)水利水電學(xué)院，邯鄲 056002；3.重慶三峽學(xué)院土木工程學(xué)院，萬州 404100)

堤(壩)基中常出現(xiàn)互層的粗粒土透水層，是滲透薄弱區(qū)域，極易發(fā)生滲透破壞，具有代表性的粗粒土互層主要由砂礫石層、細(xì)砂層和砂礫石層組成，滲透時各土層微觀的顆粒運移規(guī)律對于揭示堤(壩)基滲透變形和破壞機(jī)理至關(guān)重要[1]。因此，需要深入研究層狀粗粒土滲透變形顆粒運移規(guī)律。

近年來，中外學(xué)者針對管涌發(fā)生機(jī)理開展了大量的研究工作。劉昌軍等[2]、梁越等[3]利用室內(nèi)砂槽模型試驗對單、雙層堤(壩)基管涌的發(fā)生及發(fā)展過程進(jìn)行了模擬，將管涌破壞的發(fā)生及發(fā)展過程劃分為4個階段；賈愷等[4]利用自行設(shè)計試驗裝置對雙層堤(壩)基管涌通道擴(kuò)展機(jī)制進(jìn)行試驗，得出了雙層堤(壩)基管涌通道是否擴(kuò)展取決于通道內(nèi)沖刷水流、邊壁滲流力和砂層顆粒條件等因素。Shamy等[5]運用顆粒流方法，對上游水位上升時層狀堤(壩)基發(fā)生管涌導(dǎo)致水工建筑物出現(xiàn)較大位移和沉降的過程進(jìn)行了模擬；文獻(xiàn)[6-7]給出了單一砂層堤(壩)基上臨界水頭的理論解答。不難發(fā)現(xiàn)，目前中外對此類問題的研究主要集中于管涌臨界滲透破壞條件的判別和單雙層堤(壩)基管涌通道發(fā)展機(jī)理方向，而對此類堤(壩)基管涌發(fā)展過程中顆粒運移特點研究極少，因此很有必要對層狀覆蓋層開展相關(guān)研究。

針對典型堤(壩)基的層狀粗粒土進(jìn)行研究，利用顆粒流軟件PFC3D建立對應(yīng)數(shù)值模型，模擬堤(壩)基顆粒運移特點及滲透變形發(fā)展過程并分析顆粒流失量、下沉量等參數(shù)的變化，為進(jìn)一步認(rèn)識層狀覆蓋層堤(壩)基滲透破壞提供參考。

1 顆粒流程序計算原理

1.1 流體相方程

在固液兩相體系中，相對于固體顆粒間孔隙體積變形量，液體的體積變形量幾乎可以忽略不計，因此可將顆?？紫吨械牧黧w視為密度不變的不可壓縮流體，在求解滲流場的每一時步中，滿足流體運動的連續(xù)性方程和平均Navier-Stokes[8]方程：

(1)

(2)

式中，n為顆?？紫堵?；t為時間；為拉普拉斯算子，為流體的速度矢量；ρf為流體密度；p為壓力梯度；τ為流體黏性應(yīng)力張量；g為重力加速度矢量；fint為單位體積內(nèi)顆粒與流體的相互作用力矢量。

1.2 固體相方程

對于固體顆粒，考慮流體作用時，其主要受到流體對顆粒、顆粒間相互作用及外應(yīng)力等作用力，并通過牛頓第二定律來實現(xiàn)模擬顆粒的運動，此時顆粒運動方程為

(3)

(4)

式中，mp為顆粒的質(zhì)量；vp為顆粒速度矢量；fg為重力加速度矢量；fc為接觸c(c=1，2，…，n)處顆粒間的接觸力；fd為顆粒在流體作用下受到拖曳力，包含浮力及固液之間的相互作用力；Ip為顆粒的轉(zhuǎn)動慣量；wp為顆粒的轉(zhuǎn)動速度矢量；rc為顆粒間接觸處指向顆粒中心的方向矢量。

1.3 流固耦合方程

固體顆粒在單個流體網(wǎng)格內(nèi)的流動情況如圖1所示。

圖1 流體單元滲流示意圖Fig.1 Seepage diagram of fluid unit

流體單元網(wǎng)格體積為ΔV=ΔxΔyΔz，假定流體單元網(wǎng)格內(nèi)包含np個顆粒，則流體單元網(wǎng)格內(nèi)的孔隙率

(5)

式(5)中，dpi(i=1，2，…，n)為顆粒粒徑。

考慮流體單元網(wǎng)格內(nèi)土體顆粒的受力平衡條件，可得單個固體顆粒受到的作用力

(6)

式(6)中，fdij為單個顆粒所受作用力分量；pj為流體壓力梯度分量。

對于穩(wěn)定無分叉流，假設(shè)顆粒流體間的相互作用力只來源于壓力梯度，即：

fintj=npj

(7)

將式(7)代入式(6)，可得單個固體顆粒受到作用力的一般表達(dá)式：

(8)

當(dāng)雷諾數(shù)為1～10的層流時，German將水力半徑理論公式代入Darcy公式得壓力梯度方程為

(9)

對于雷諾數(shù)大于10的紊流，Ergun[9]基于水力半徑理論提出了Ergun公式：

(10)

1.4 流固耦合過程

在顆粒流程序中，整個流體區(qū)域被劃分為若干個固定大小的流體計算單元，采用基于半隱式算法(SIMPLE法)結(jié)構(gòu)的計算流體動力學(xué)(CFD)計算技術(shù)對液相流動方程進(jìn)行求解。在計算過程中，首先根據(jù)力-位移定理計算得到顆粒間的接觸力，而后將該力與流體對顆粒的拖曳力一起代入下一時步顆粒運動方程中，以此確定顆粒新的位置、速度及接觸力等。同時顆粒的運動又將引起流體單元內(nèi)壓力、流速、孔隙率的變化，進(jìn)而影響到顆粒間接觸力及流體對顆粒拖曳力的變化，可見流固耦合模擬過程為一個動態(tài)循環(huán)平衡過程，其計算過程如圖2所示。

圖2 PFC3D流固耦合循環(huán)計算過程Fig.2 PFC3D fluid-solid coupling cycle computation process

2 數(shù)值模型建立

2.1 建模依據(jù)

中國西南和西北的河谷覆蓋層在縱向結(jié)構(gòu)上多具有強(qiáng)、弱透水層交替層疊的特點，各層巖性、物理力學(xué)特性等差異較大，較有代表性是砂卵礫石和細(xì)砂互層結(jié)構(gòu)如圖3(a)，各層取樣如圖3(b)、圖3(c)所示。

圖3 典型粗粒土及粗細(xì)相間的層狀覆蓋層Fig.3 Typical coarse grained soil and coarse and fine layered overburden

典型砂卵礫石層和細(xì)砂層的基本特征及透水性等級如表1所示。

表1 粗細(xì)粒土基本特征及透水性等級Table 1 Basic characteristics and permeability grade of coarse and fine-grained soil

根據(jù)上述典型的層狀粗粒土層，運用離散元軟件PFC3D來建立滲流水-土耦合模型，厘清顆粒運移規(guī)律。

2.2 數(shù)值模型生成及細(xì)觀參數(shù)標(biāo)定

2.2.1 數(shù)值模型生成

為了更準(zhǔn)確地對實際堤(壩)基中層狀粗粒土結(jié)構(gòu)的滲透變形過程進(jìn)行模擬，模型生成的顆粒必須能反映實際土體顆粒級配關(guān)系，實際土體中顆粒粒徑范圍(0.075～50 mm)較寬，如按照實際土體顆粒級配建立模型，將生成數(shù)量十分巨大的顆粒，導(dǎo)致計算效率降低甚至無法完成計算。因此，需對實際土體級配進(jìn)行一定的近似處理，本文數(shù)值模型中顆粒最大最小粒徑分別取為10 mm和0.25 mm，由工程地質(zhì)圖3及表1可知，上部土層為間斷級配的砂礫石土層，以2 mm作為粗細(xì)顆粒的區(qū)分粒徑，同時保持生成顆粒質(zhì)量百分比與實際土體粗細(xì)顆粒百分比基本一致，近似處理后的土體級配如表2所示。應(yīng)注意的是，該顆粒粒徑處理方法可能對試驗精度造成一定的誤差，但對滲透變形中顆粒運移過程的研究是可行的。此外，為降低顆粒模型中總的自由度，這里運用土工離心機(jī)試驗中的相似性原理，建立長×寬×高為120 mm×60 mm×60 mm的顆粒流模型如圖4所示。同時，為使建立的模型與實際土體的工程地質(zhì)分層情況基本一致，利用PFC3D內(nèi)置的fish語言編寫程序以生成不同土層，模型從上至下依次為：上部砂礫石層厚26 mm，中間細(xì)砂層厚6 mm，下部砂礫石層厚26 mm。此外，為防止上覆土層與模型上部墻體發(fā)生接觸沖刷，除管涌口外，在墻體上部邊界生成顆粒粒徑為2 mm規(guī)則排列的固定顆粒邊界。

表2 PFC數(shù)值模型土樣參數(shù)Table 2 Soil sample parameters of PFC numerical model

圖4 堤(壩)基中多層粗粒土結(jié)構(gòu)顆粒流模型Fig.4 Particle flow model of multi-layer coarse-grained soil structure in embankment (dam) foundation

2.2.2 細(xì)觀參數(shù)標(biāo)定

為使所建模型能更有效的反映實際滲透變形過程，需要準(zhǔn)確找出宏觀材料與細(xì)觀顆粒之間物理力學(xué)屬性的聯(lián)系。在利用PFC3D來模擬滲流破壞時，顆粒運動規(guī)律主要受到細(xì)觀參數(shù)摩擦系數(shù)的影響[10]，因此在數(shù)值模擬之前，必須對顆粒的摩擦系數(shù)進(jìn)行標(biāo)定，但顆粒宏-細(xì)觀參數(shù)之間存在較大的差異。研究表明，在散粒體材料中，材料的內(nèi)摩擦角是其摩擦特性的綜合體現(xiàn)[11]，因此采用“反演模擬法”[12]來對顆粒的細(xì)觀摩擦系數(shù)進(jìn)行標(biāo)定，室內(nèi)試驗測得砂礫石和細(xì)砂的宏觀內(nèi)摩擦角分別為35°10′、26°30′，然后利用PFC3D進(jìn)行三軸數(shù)值模擬試驗(圖5)，通過不斷調(diào)整顆粒細(xì)觀摩擦系數(shù)來獲取顆粒的細(xì)觀內(nèi)摩擦角，將得到的細(xì)觀內(nèi)摩擦角與宏觀內(nèi)摩擦角進(jìn)行對比，直到兩者接近或相等為止，最終得到粗礫和細(xì)砂的摩擦系數(shù)為1.25與0.79。模型具體細(xì)觀參數(shù)如表3所示。

圖5 三軸數(shù)值模擬試驗Fig.5 Triaxial numerical simulation test

表3 PFC模型具體細(xì)觀參數(shù)Table 3 Specific meso-parameters of the model

參照表2及表3中模型細(xì)觀參數(shù)，利用PFC3D中fish語言編寫函數(shù)生成相應(yīng)土層顆粒，各土層之間賦予線性接觸模型，采用半徑擴(kuò)大法生成土顆粒，先將顆粒粒徑縮小，再將顆粒粒徑放大到指定的孔隙率，而后對土顆粒施加重力并進(jìn)行循環(huán)計算，直至土顆粒間不平衡力較小或消除為止，以確保土層達(dá)到穩(wěn)定狀態(tài)。

2.3 流體網(wǎng)格劃分及監(jiān)測單元設(shè)置

PFC3D中運用固定粗糙網(wǎng)格法對流體進(jìn)行處理，需將流體區(qū)域劃分為若干流體單元網(wǎng)格，且應(yīng)確保所劃分的流體網(wǎng)格中均包含一定數(shù)量的顆粒。故考慮將整個模型沿長、寬、高方向分別劃分為8、4、4份，單個流體單元網(wǎng)格尺寸為15 mm×15 mm×15 mm。再進(jìn)行流體計算之前，需對邊界條件進(jìn)行設(shè)定，模型左側(cè)設(shè)定為施加水頭邊界，水頭壓力大小為100 kPa，上部、底部及側(cè)壁均設(shè)定為剛性不透水非滑移邊界，將模型上部X=12 mm、Y=39 mm、Z=24 mm處設(shè)置為零壓力流體出流口，出流口長寬為3 mm×3 mm[見圖6(a)]，該流體單元網(wǎng)格模型圖詳見圖6(b)。應(yīng)注意的是，在左側(cè)邊界施加水頭之前，為避免流域內(nèi)形成紊亂的流場，應(yīng)首先將顆粒固定，然后施加水頭，直至流場趨于穩(wěn)定，此后釋放顆粒，讓固體顆粒與流體進(jìn)行作用。

圖6 計算邊界及流體單元Fig.6 Computational boundary and fluid unit

在模型開始計算之前，為清晰記錄各土層不同部位顆粒流失的情況，分別在各土層設(shè)置以下監(jiān)測區(qū)域，上部砂礫石層及細(xì)砂層從壓力上游端至下游端分別編號為S1、S2、S3、S4及Z1、Z2、Z3、Z4，如圖7所示，由于細(xì)砂層對下部砂礫石層起阻擋作用，導(dǎo)致下部砂礫石層顆粒很難發(fā)生流失，因此不對下部砂礫石層中顆粒流失進(jìn)行監(jiān)測。

圖7 模型不同監(jiān)測區(qū)域劃分Fig.7 Division of different monitoring areas in the model

3 數(shù)值試樣結(jié)果分析

3.1 上部砂礫石層滲透破壞過程

3.1.1 上部砂礫石層細(xì)顆粒遷移過程

該粗粒土層狀結(jié)構(gòu)的上部砂礫石層中細(xì)顆粒遷移過程如圖8所示。

從圖8可以看出，因上部砂礫石為管涌型土，粗顆粒之間存在一定的孔隙，在初始水頭壓力作用下，上部砂礫石層中細(xì)顆粒從上游端逐漸向管涌口運移，上部砂礫石層中細(xì)顆粒在管涌口陸續(xù)流失，但并未出現(xiàn)顆粒整體流失及細(xì)砂層顆粒侵入上部砂礫石層中的現(xiàn)象 [見圖8(a)]。在計算時間步達(dá)到30萬步時，上部砂礫石層粗顆粒中細(xì)顆粒持續(xù)從上游端向管涌口處移動，導(dǎo)致在上游水頭S1區(qū)出現(xiàn)較大孔隙，在細(xì)顆粒移動過程中，有部分細(xì)顆粒因遇到較小孔隙而被阻擋形成阻塞區(qū)，此時，細(xì)砂層Z3和Z4區(qū)已有部分顆粒進(jìn)入上部砂礫石層中，且由于S1區(qū)有較大孔隙出現(xiàn)，細(xì)砂層Z1區(qū)部分顆粒也開始進(jìn)入上部砂礫石層中[見圖8(b)]。隨著計算時間步達(dá)到50萬步時，上部砂礫石層中細(xì)顆粒沿著粗顆粒間孔隙不斷從上游端向管涌口處移動，細(xì)砂層中顆粒進(jìn)入上部砂礫石層中含量進(jìn)一步增加，有部分顆粒已進(jìn)入上部砂礫石層中部，以至于細(xì)砂層出現(xiàn)較大變形[見圖8(c)]。

圖8 上部砂礫石層中細(xì)顆粒遷移過程Fig.8 Migration process of fine particles in upper gravel layer

3.1.2 上部砂礫石層不同區(qū)域細(xì)顆粒流失量

上部砂礫石層中不同區(qū)域細(xì)顆粒流失量隨計算時間步的變化曲線如圖9所示。顆粒流失量定義為該區(qū)域細(xì)顆粒流失的總體積與水頭壓力施加前該區(qū)域細(xì)顆粒總體積之比。

圖9 上部砂礫石層不同區(qū)域細(xì)顆粒流失量情況Fig.9 Fine particle loss in different areas of upper gravel layer

由圖9可知，在數(shù)值試驗過程中，上部砂礫石層上游端S1區(qū)細(xì)顆粒流失量最大，約20%，導(dǎo)致S1區(qū)粗顆粒間有較大孔隙形成。S2、S3區(qū)細(xì)顆粒流失量相對較少，當(dāng)時間步增加至20萬步左右時，S2、S3區(qū)顆粒流失量分別增至3.5%、3%，而后直至計算結(jié)束，流失量分別降低至2.5%、1%，這是由于細(xì)顆粒從上游端流向管涌口處時逐漸在S2和S3區(qū)形成阻塞區(qū)所致。當(dāng)時間步計算至25萬步時，S4區(qū)細(xì)顆粒流失量增至5%左右，之后隨時間步的增加，S4區(qū)顆粒流失量幾乎沒有變化，這是因為在試驗初期管涌口顆粒不斷流失，隨著時間步的增加，上游端細(xì)顆粒持續(xù)流向管涌口，同時，細(xì)砂層顆粒陸續(xù)進(jìn)入上部砂礫石層中，導(dǎo)致S4區(qū)細(xì)顆粒流失量逐漸減少直至幾乎不變。

通過對上部砂礫石層中細(xì)顆粒遷移現(xiàn)象及流失曲線分析可知，在滲透力作用下，上部砂礫石層中細(xì)顆粒在粗顆?？紫堕g移動而后被帶出土體外，逐漸形成穩(wěn)定的管涌通道，屬于典型的管涌破壞。

3.2 細(xì)砂層滲透破壞過程

3.2.1 細(xì)砂層顆粒遷移過程

該層狀結(jié)構(gòu)中的細(xì)砂層中顆粒遷移過程如圖10所示。

圖10 細(xì)砂層中顆粒遷移過程Fig.10 Particle migration in fine sand

從圖10(a)可清晰觀察到，在施加水頭壓力初期，細(xì)砂層顆粒幾乎沒有發(fā)生。如圖10(b)所示，隨著計算時間步增加至30萬步，細(xì)砂層中上游端Z1區(qū)及管涌口附近Z3及Z4區(qū)顆粒已部分侵入上部砂礫石層中，甚至在管涌口正下方細(xì)砂層Z4區(qū)已有少許顆粒進(jìn)入上部砂礫石層中部，這是由于上部砂礫石層中細(xì)顆粒的逐漸流失，給予了細(xì)砂層進(jìn)入上部砂礫石層的空間。當(dāng)時間步進(jìn)一步增加至50萬步過程中，細(xì)砂層Z1、Z3及Z4區(qū)處顆粒進(jìn)入上部砂礫石層中越來越多，且在細(xì)砂層Z4區(qū)已有一些顆粒流失至管涌口附近，隨著顆粒持續(xù)流失，導(dǎo)致細(xì)砂層逐漸出現(xiàn)小范圍的變形[見圖10(c)]。

3.2.2 細(xì)砂層不同區(qū)域顆粒流失體積

細(xì)砂層不同區(qū)域進(jìn)入上部砂礫石層顆粒流失體積隨計算時間步的變化曲線如圖11所示。

圖11 細(xì)砂層不同區(qū)域顆粒流失體積變化Fig.11 Variation of particle loss volume in different areas of fine sand bed

由圖11可知，在數(shù)值試驗計算過程中，細(xì)砂層在管涌口下方Z4區(qū)處顆粒持續(xù)流失并且流失顆粒含量最多，約3.5×10-8m3，其余區(qū)域顆粒流失量相對較少。在數(shù)值試驗初期，僅有曲線S4隨時間步的增加而增加，而其余曲線無明顯變化，表明細(xì)砂層最先在管涌口正下方Z4區(qū)處出現(xiàn)顆粒流失。隨著計算時間步的增加至30萬步，曲線Z1、Z2及Z3均隨著計算時間步的增加而增加，說明Z1、Z2及Z3區(qū)顆粒都已發(fā)生流失，且在此過程中，可見Z2區(qū)顆粒是最后發(fā)生流失的，這是由于上部砂礫石層S2區(qū)有阻塞區(qū)形成，阻礙了細(xì)砂層顆粒的進(jìn)入。當(dāng)計算時間步增加至50萬步時，曲線Z1隨時間步緩緩增加，對應(yīng)細(xì)砂層Z1區(qū)顆粒流失量逐漸增多，這是因為上部砂礫石層中上游端S1區(qū)細(xì)顆粒逐漸流失，以至形成較大孔隙，給予了Z1區(qū)顆粒進(jìn)入上部砂礫石層空間。

通過對中間細(xì)砂層中顆粒遷移現(xiàn)象及流失曲線分析可知，細(xì)砂層中顆粒最先在管涌口正下方Z4區(qū)發(fā)生流失，其余區(qū)域顆粒流失相對較晚，且顆粒流失量均隨著計算時間步的增加而增加，導(dǎo)致細(xì)砂層出現(xiàn)小范圍的變形。

3.3 上部砂礫石層下沉量

如圖12所示為上部砂礫石層下沉量隨計算時間步的變化曲線，下沉量定義為上部砂礫石層顆粒垂直向下方向運動距離總和。從圖12中可以看出，隨著計算時間的增加，上部砂礫石層的下沉量逐漸增加至0.012 m，這是由于細(xì)砂層中顆粒不斷流失所致，這種情況下當(dāng)上部砂礫石層細(xì)顆粒含量流失達(dá)到一定程度，隨著管涌通道的貫通而發(fā)生破壞，將對上部建筑物產(chǎn)生重大危害。通過顆粒流模擬，上部砂礫石層及細(xì)砂層所表現(xiàn)出的滲透破壞現(xiàn)象，這與劉杰[13]對間斷級配土滲透破壞機(jī)理試驗研究所得結(jié)論相吻合。

圖12 上部砂礫石層下沉量Fig.12 Subsidence of upper gravel layer

4 結(jié)論

在PFC3D細(xì)觀參數(shù)標(biāo)定時，采用“反演模擬法”，較準(zhǔn)確地建立了材料宏觀參數(shù)與顆粒細(xì)觀參數(shù)間聯(lián)系，有效的模擬了多層粗粒土層滲透變形的發(fā)展過程，獲得了堤(壩)基滲透變形中顆粒流失及下沉量等情況，為此類地基的滲透破壞的研究提供了一定參考，并得到以下結(jié)論。

(1)滲透破壞主要發(fā)生在上部砂礫石層中，隨著滲透破壞的持續(xù)發(fā)生，逐漸影響下部土層；作為滲流出口的上部砂礫石層，該層中細(xì)顆粒在粗顆?？紫堕g移動而后逐漸被帶出土體外，逐漸形成穩(wěn)定的管涌通道，具有明顯的溯源性。

(2)隨著計算時間的增加，上部砂礫石層的下沉量逐漸增加，當(dāng)上部砂礫石層細(xì)顆粒含量流失達(dá)到一定程度，堤(壩)基發(fā)生破壞，將對上部建筑物產(chǎn)生重大危害。

(3)中間細(xì)砂層在上部砂礫石層管涌破壞后，其顆粒最先在管涌口正下方Z4區(qū)發(fā)生流失，其余區(qū)域顆粒流失相對較晚，且顆粒流失量均隨著計算時間步的增加而增加，導(dǎo)致細(xì)砂層出現(xiàn)小范圍的變形出現(xiàn)。

(4)由于幾何條件不滿足，在上部土層不發(fā)生明顯破壞的情況下，不會對深層的砂礫石層顆粒運移造成影響。

綜上，上部強(qiáng)透水層是滲透破壞的關(guān)鍵區(qū)域，在實際工程中應(yīng)重點處理；深層次的強(qiáng)透水層對堤壩基滲透破壞的貢獻(xiàn)較小，在上部土層不發(fā)生破壞的情況下，深部土層一般不會發(fā)生顆粒流失現(xiàn)象。