馮 梅，王 剛，，張甲波，白玉川，黃 哲

(1.河北省海洋地質(zhì)資源調(diào)查中心，河北 秦皇島 066000；2.河北省海洋岸線(xiàn)生態(tài)修復(fù)與智慧海洋工程研究中心，河北 秦皇島 066000；3.天津大學(xué)水利工程仿真與安全國(guó)家重點(diǎn)實(shí)驗(yàn)室，天津 300350)

1 研究背景

滲流導(dǎo)致的土體內(nèi)部泥沙顆粒運(yùn)動(dòng)是壩體、河堤失穩(wěn)破壞的重要原因[1-2]。根據(jù)泥沙運(yùn)動(dòng)形式的不同，滲流輸沙分為管涌、流土、接觸流土、接觸沖刷等多種形式。據(jù)統(tǒng)計(jì)，世界上約有50%的壩體破壞是由于滲流所致的泥沙流失引起的[3-4]。因此研究滲流作用下的泥沙運(yùn)動(dòng)規(guī)律具有重要的實(shí)際意義。

管涌是滲流輸沙的一種常見(jiàn)機(jī)制，常發(fā)生于洪水期的河堤及土石壩內(nèi)部[5]。由于其嚴(yán)重的破壞作用，研究者們針對(duì)其內(nèi)在機(jī)制和發(fā)展過(guò)程開(kāi)展了大量的研究[6-7]。本文所研究的滲流輸沙特指管涌現(xiàn)象。由于泥沙運(yùn)動(dòng)的過(guò)程是一個(gè)復(fù)雜的水土相互耦合作用的過(guò)程[8]，從而很難用數(shù)學(xué)語(yǔ)言對(duì)其進(jìn)行描述[9]。因此，在以往的研究中主要通過(guò)室內(nèi)試驗(yàn)或是現(xiàn)場(chǎng)測(cè)試進(jìn)行研究，并且利用試驗(yàn)結(jié)果建立經(jīng)驗(yàn)?zāi)Ｐ?，用?lái)表達(dá)輸沙率與各種參數(shù)之間的關(guān)系[7]。例如，Howard和McLane[10]基于土體受力平衡的理論，認(rèn)為滲流的輸沙量是由斜坡坡角與臨界坡角的相對(duì)關(guān)系決定的。以此為基礎(chǔ)，F(xiàn)ox等[11]提出了滲流剪應(yīng)力的計(jì)算方法，并且利用試驗(yàn)數(shù)據(jù)提出滲流輸沙量和滲流剪應(yīng)力的經(jīng)驗(yàn)關(guān)系。Chu-Agor等[12]提出了垂直河岸的滲流輸沙模型，認(rèn)為輸沙量是和滲流流量相關(guān)的函數(shù)。梁越等[9]在其流固耦合模型中假設(shè)了輸沙量與滲流流速之間的關(guān)系。另外，水力梯度是滲流中的一個(gè)重要參數(shù)，滲流輸沙的起動(dòng)條件也常用水力梯度描述[13-14]。比如Cividini等[15]通過(guò)抽水井周?chē)鷿B流場(chǎng)中的泥沙輸運(yùn)試驗(yàn)，建立了土體細(xì)顆粒的密度隨水力梯度的變化關(guān)系。再如Sterpi[16]利用取自意大利米蘭的黏土進(jìn)行研究，得到了不同滲流梯度情況下泥沙可動(dòng)顆粒的流失規(guī)律。Fox等[17]利用室內(nèi)的水平滲流試驗(yàn)結(jié)果，得出了滲流輸沙量與水力梯度以及臨界水力梯度的關(guān)系?？偨Y(jié)已有的經(jīng)驗(yàn)?zāi)Ｐ桶l(fā)現(xiàn)，滲流輸沙率常被描述為以下關(guān)系

q=Kx(X-Xcr)a

(1)

式中，q為輸沙率；Kx和a為與土體相關(guān)的參數(shù)；X為土體的滲流場(chǎng)參數(shù)，根據(jù)上文，X可表示滲流切應(yīng)力、滲流流量、滲流梯度、滲流速度等；Xcr為發(fā)生泥沙運(yùn)動(dòng)時(shí)X的臨界值。

近年來(lái)，隨著研究手段的發(fā)展和研究思想的創(chuàng)新，新的研究方法也被利用來(lái)研究滲流輸沙問(wèn)題[9，18-19]，但是多數(shù)的數(shù)學(xué)模型依然通過(guò)經(jīng)驗(yàn)方法計(jì)算輸沙量[12，16]。該類(lèi)數(shù)學(xué)模型能夠較完整的模擬出滲流輸沙過(guò)程以及土體結(jié)構(gòu)響應(yīng)規(guī)律[9，20]。但是，由于輸沙過(guò)程仍是由經(jīng)驗(yàn)公式計(jì)算而來(lái)，缺少對(duì)顆粒運(yùn)動(dòng)機(jī)制和運(yùn)動(dòng)狀態(tài)的定量描述。因此，其研究成果的普適性還存在問(wèn)題[21]，經(jīng)驗(yàn)公式中的參數(shù)確定還存在諸多困難。此外，研究者也試圖采用顆粒流方法模擬滲流輸沙過(guò)程[22-23]，模擬過(guò)程中能較好的捕捉泥沙顆粒運(yùn)動(dòng)軌跡和運(yùn)動(dòng)速度[24]。但是該類(lèi)方法計(jì)算效率不高[25]，所采用的土體條件和水力條件過(guò)于理想化，模型的驗(yàn)證以及對(duì)實(shí)際管涌現(xiàn)象的模擬還存在缺陷[21]。在這樣的背景下，亟需一種滲流輸沙的理論計(jì)算模型來(lái)彌補(bǔ)管涌模擬模型的不足。根據(jù)河流動(dòng)力學(xué)的基本理論，土顆粒運(yùn)動(dòng)狀態(tài)是精準(zhǔn)把握輸沙率的重要指標(biāo)參數(shù)[26]。而滲流作用下土體內(nèi)部顆粒的運(yùn)動(dòng)狀態(tài)是極為復(fù)雜的，而且存在很大的隨機(jī)性[27]。隨著管涌的不斷發(fā)展，顆粒運(yùn)動(dòng)的數(shù)量和強(qiáng)度均會(huì)發(fā)生明顯變化[28]。因此，如何精準(zhǔn)描述土顆粒的運(yùn)動(dòng)狀態(tài)成為滲流輸沙研究中重要的科學(xué)問(wèn)題之一。

本文從泥沙顆粒的受力狀態(tài)出發(fā)，推導(dǎo)泥沙顆粒在滲流場(chǎng)中運(yùn)動(dòng)的控制方程。通過(guò)分析滲流流速的概率分布模式，將單顆粒泥沙運(yùn)動(dòng)規(guī)律擴(kuò)展到整個(gè)滲流場(chǎng)，并由此計(jì)算得到泥沙運(yùn)動(dòng)概率和輸沙規(guī)律。通過(guò)和試驗(yàn)結(jié)果的對(duì)比，表明該模型具有合理性，可以對(duì)一維滲流作用下的泥沙運(yùn)動(dòng)規(guī)律進(jìn)行描述。

2 材料和方法

2.1 動(dòng)力學(xué)方程

如圖1所示，滲流場(chǎng)中單顆粒受力包括滲流力和水下重力[13]。其中滲流力包括2部分，即由于壓力差引起的滲流梯度力以及由于流動(dòng)引起的水流拖曳力[29-30]。

圖1 土體內(nèi)部可動(dòng)顆粒的受力分析

顆粒所受滲流力F和水下重力W′表達(dá)式為

F=FI+FD

(2)

(3)

式中，F(xiàn)D為水流拖曳力；FI為滲流梯度力；D為顆粒直徑；ρ′為顆粒的水下密度。

FD[26]和FI[31]的表達(dá)式分別為

(4)

(5)

式中，α為顆粒形狀修正系數(shù)，對(duì)于球形取1，對(duì)于泥沙顆粒取1.7；CD為拖曳力系數(shù)；ρw為水的密度；ui為顆粒運(yùn)動(dòng)速度；ubi為孔隙內(nèi)水流速度；I為水力梯度。

根據(jù)牛頓第二定律

(6)

將式(2)～式(5)帶入式(6)可得

(7)

式中拖曳力系數(shù)與雷諾數(shù)有關(guān)，滲流作用下的顆粒雷諾數(shù)Re的計(jì)算公式為

(8)

式中，v為水的動(dòng)力粘滯系數(shù)；u為流速，考慮到顆粒運(yùn)動(dòng)，u應(yīng)是水流與顆粒的相對(duì)速度，即u=ubi-ui。由于滲流流速較小，一般小于10-3m/s，顆粒大小一般為毫米級(jí)，數(shù)量級(jí)也為10-3m。因此，滲流雷諾數(shù)一般小于1。此時(shí)拖曳力系數(shù)與雷諾數(shù)的關(guān)系由斯托克斯定律確定[26]，即

(9)

將式(8)、式(9)帶入式(7)可得

(10)

式(10)即為一維滲流作用下泥沙顆粒運(yùn)動(dòng)的控制方程。

2.2 臨界起動(dòng)流速

滲流梯度和滲流流速的關(guān)系服從達(dá)西定律

ud=KI

(11)

式中，K為滲透系數(shù)；ud為滲流流速。

引入平均孔隙流速，滲流流速與平均孔隙流速的關(guān)系為

(12)

式中，n為土體孔隙率。

(13)

式(13)為顆粒運(yùn)動(dòng)的臨界孔隙流速，其與臨界水力梯度的關(guān)系為

(14)

可由式(13)和式(14)計(jì)算不同粒徑的臨界水力梯度。根據(jù)劉杰的研究成果[32]，管涌是部分顆粒的運(yùn)動(dòng)，選取d5作為特征粒徑進(jìn)行臨界水力梯度的計(jì)算是比較合理的。因此式(13)中，D=d5，d5為小于該粒徑的土顆粒質(zhì)量占總質(zhì)量的5%。

2.3 顆粒運(yùn)動(dòng)狀態(tài)分析

對(duì)于式(10)，左右兩邊對(duì)t積分，并考慮到t=0時(shí)，ui=0，可得

(15)

式(15)即為孔隙內(nèi)單個(gè)土顆粒的運(yùn)動(dòng)速度表達(dá)式。

2.4 輸沙率的確定

根據(jù)泥沙推移質(zhì)運(yùn)動(dòng)理論，垂直滲流作用下，單位面積單位時(shí)間內(nèi)的輸沙率可由下式計(jì)算：

(16)

式中，q為單位時(shí)間單位面積的滲流斷面上的輸沙率；N為單位體積內(nèi)可動(dòng)顆粒的個(gè)數(shù)；V為可動(dòng)顆粒的體積；P為可動(dòng)顆粒運(yùn)動(dòng)概率。

N和V的表達(dá)式為

(17)

(18)

式中，PV為可動(dòng)顆粒的體積比。

2.4.1P的確定方法

因?yàn)橥馏w孔隙有大有小，孔隙內(nèi)流速也呈大小不一的分布規(guī)律。只有孔隙內(nèi)的流速大于土顆粒的臨界起動(dòng)流速，土顆粒才會(huì)被水流帶走。但是如果顆粒運(yùn)動(dòng)軌跡中的下游孔隙小于顆粒直徑時(shí)，土顆粒依然不能被水流帶走。因此顆粒運(yùn)動(dòng)還應(yīng)該滿(mǎn)足一個(gè)條件，本文稱(chēng)為顆粒的輸移條件，即顆粒輸移必須滿(mǎn)足下游孔隙大于可動(dòng)顆粒的直徑。所以這里的P包括兩部分P1和P2，且P滿(mǎn)足

P=P1P2

(19)

式中，P1為土顆粒的起動(dòng)概率；P2為土顆粒的輸移概率。假設(shè)土體內(nèi)水流流動(dòng)符合管流泊肅葉理論，管流速度和圓管半徑的關(guān)系為[13]

(20)

這里假定公式(20)中的R服從均勻分布(0，2R)，R為平均圓管半徑，則可根據(jù)式(20)求得孔隙內(nèi)流速u(mài)bi的概率分布為

(21)

因此，P1的計(jì)算公式為

(22)

根據(jù)Hazen的研究[33]，土體孔隙直徑和顆粒級(jí)配有關(guān)：

D0=Bd10

(23)

式中，D0為孔隙直徑；B為參數(shù)，取值為0.21～0.25；d10為土體特征粒徑，表示小于該粒徑的土顆粒質(zhì)量占總質(zhì)量的10%。由此可得P2的計(jì)算方法為

P2=P(D

(24)

式中可根據(jù)土體級(jí)配的分布特征確定可動(dòng)顆粒D0的大小，然后根據(jù)顆粒的質(zhì)量比計(jì)算D小于D0的概率。

(25)

將式(15)、式(17)、式(18)、式(19)、式(25)帶入式(16)可求得單位時(shí)間單位面積的輸沙率。

2.5 滲流-泥沙場(chǎng)的耦合計(jì)算

前文已有論述，由于滲流輸沙的過(guò)程伴隨著滲流場(chǎng)的變化，因此對(duì)于輸沙率的計(jì)算需要耦合滲流場(chǎng)的變化。t時(shí)刻相對(duì)于t-1時(shí)刻，土體的孔隙率n、滲透系數(shù)K、可動(dòng)顆粒的體積分?jǐn)?shù)Pv會(huì)發(fā)生變化。而其他變量，比如臨界流速u(mài)cr、孔隙流速u(mài)b、起動(dòng)概率P、孔隙半徑R都可以由以上幾個(gè)變量求解得到。因此，在已知t-1時(shí)刻各參數(shù)的值的情況下，只需確定t時(shí)刻n、K和PV的計(jì)算方法，就可實(shí)現(xiàn)時(shí)間上的迭代計(jì)算，進(jìn)而實(shí)現(xiàn)滲流-泥沙場(chǎng)的耦合計(jì)算。

由輸沙率q的定義可得在t時(shí)刻，總輸沙量Qt的計(jì)算公式為

(26)

假設(shè)在輸沙過(guò)程中，沒(méi)有土體的變形，即輸沙的體積即為土體中孔隙增加的體積。則

(27)

根據(jù)劉杰的研究[32]，滲透系數(shù)的比值與孔隙率和d20有關(guān)。

(28)

假設(shè)水流帶走的土顆粒均小于(d20)0，而且在滲流輸沙的過(guò)程中，小于d20的顆粒粒徑是均勻的，此時(shí)可以計(jì)算輸沙Qt之后的(d20)t為

(29)

式中，(dx)0為可動(dòng)顆粒的最大粒徑；(Pm)t為粒徑小于最大可動(dòng)顆粒的土體所占的質(zhì)量分?jǐn)?shù)，其與輸沙量Qt有關(guān)，即

(30)

式中，m0為初始時(shí)刻的可動(dòng)顆粒的質(zhì)量，將式(29)、式(30)帶入式(28)可求得t時(shí)刻的滲透系數(shù)。

(PV)t為可動(dòng)顆粒的體積占總顆粒體積的比值，其與輸沙量的關(guān)系為

(31)

根據(jù)以上論述可通過(guò)式(27)、式(28)和式(31)實(shí)現(xiàn)滲流-泥沙場(chǎng)耦合計(jì)算，耦合計(jì)算流程如圖2所示。

圖2 耦合計(jì)算流程

3 模型驗(yàn)證

3.1 臨界水力梯度的驗(yàn)證

許多研究人員對(duì)于滲流的臨界水力梯度進(jìn)行了試驗(yàn)研究，目前已發(fā)表的有4個(gè)經(jīng)典的試驗(yàn)結(jié)果[14，28，34-35]，這些試驗(yàn)都以無(wú)黏性土為研究對(duì)象，符合本文計(jì)算方法的前提條件。Skempton試驗(yàn)研究了4種級(jí)配土體的臨界水力梯度，4組土體的級(jí)配如圖3所示。圖3中可以看出土樣A和土樣B級(jí)配較差，土體內(nèi)部不穩(wěn)，屬于管涌型土，而土樣C和土樣D屬于內(nèi)部穩(wěn)定性土體[36]。Ahlinhan試驗(yàn)中測(cè)定了5種土體在不同孔隙比情況下的臨界水力梯度，其土體級(jí)配如圖4所示，可以看出5組土樣中，A1和A2土樣級(jí)配較好，屬于流土型土；E2和E3土樣級(jí)配較差，屬于管涌型土，而E1屬于過(guò)渡性土[36]。Fleshman試驗(yàn)測(cè)定了天然沙的臨界水力梯度，其土體級(jí)配如圖5所示，其土樣均為穩(wěn)定性土體。陳亮試驗(yàn)測(cè)定了3種級(jí)配土體分別在4種孔隙率條件下的臨界水力梯度，其土體級(jí)配如圖6所示，可以發(fā)現(xiàn)土樣A和土樣B為管涌型土，土樣C為過(guò)渡性土。4組試驗(yàn)結(jié)果與計(jì)算值的對(duì)比，如圖7所示。本文計(jì)算方法與太沙基公式[37]的對(duì)比如圖8所示。

圖3 土樣級(jí)配曲線(xiàn)(Skempton試驗(yàn))

圖4 土樣級(jí)配曲線(xiàn)(Ahlinhan試驗(yàn))

圖5 土樣級(jí)配曲線(xiàn)(Fleshman試驗(yàn))

圖6 土樣級(jí)配曲線(xiàn)(陳亮試驗(yàn))

圖7 臨界水力梯度的驗(yàn)證

圖8 本文公式與太沙基方法的對(duì)比

根據(jù)圖7，總體看臨界水力梯度本文計(jì)算值與試驗(yàn)值接近(R2=0.79)。由圖7可知，只有極個(gè)別的數(shù)據(jù)點(diǎn)在置信區(qū)間為95%的曲線(xiàn)外(圖7中虛線(xiàn))，誤差較大的數(shù)據(jù)點(diǎn)來(lái)自陳亮試驗(yàn)，其試驗(yàn)中測(cè)定的土樣C在較小的孔隙率條件下臨界水力梯度達(dá)到2.4～2.6，這在其他試驗(yàn)或文獻(xiàn)中還未見(jiàn)到，一般當(dāng)水力梯度如此大的情況下，任何粒徑的土顆粒的重量都小于水力梯度力[23]。因此本文推測(cè)陳亮試驗(yàn)中，對(duì)應(yīng)的這兩組試驗(yàn)結(jié)果有測(cè)量誤差。其他試驗(yàn)中，Skempton試驗(yàn)結(jié)果略微大于計(jì)算值，而Fleshman的試驗(yàn)正好相反，Ahlinhan的試驗(yàn)對(duì)應(yīng)數(shù)據(jù)點(diǎn)分布在直線(xiàn)y=x兩側(cè)，與計(jì)算值較為接近。而從圖8可以看出，本文的計(jì)算值要優(yōu)于太沙基方法(R2=0.19)，太沙基的計(jì)算結(jié)果對(duì)于不同土體差別不大，在其公式中僅包含土體密度和孔隙率的參數(shù)，而不同的土顆粒，其密度和孔隙率均相差不大，其公式對(duì)于不同級(jí)配土體的考慮并沒(méi)有涉及。綜上，本文的計(jì)算公式適用性更強(qiáng)。

3.2 出沙量的驗(yàn)證

陳亮在試驗(yàn)中測(cè)定了滲流出沙和出沙后的滲透系數(shù)[28]，此外還列出了對(duì)于土樣A在4種不同孔隙率條件下的出沙量和滲透系數(shù)的變化以及土樣B在2種不同孔隙條件下的出沙量(土樣C的臨界水力梯度的試驗(yàn)存疑，此處并沒(méi)驗(yàn)證)，本文計(jì)算值與滲流輸沙量的對(duì)比結(jié)果如圖9所示。

圖9 累積出沙量的驗(yàn)證

由圖9可以看出，本文的計(jì)算值與試驗(yàn)值非常接近。出沙量隨時(shí)間逐漸增多，且呈現(xiàn)出兩個(gè)拐點(diǎn)。單位時(shí)間的出沙量，即輸沙率是隨著時(shí)間先增大，然后逐漸減小，最后基本為0。

3.3 滲透系數(shù)的驗(yàn)證

在土顆粒不斷流失的過(guò)程中，土體的孔隙率會(huì)逐漸變大，土顆粒也會(huì)逐漸變粗。該過(guò)程也會(huì)導(dǎo)致滲透系數(shù)的變化。在陳亮的試驗(yàn)中，同樣測(cè)定了不同時(shí)刻的土體滲透系數(shù)，在文獻(xiàn)[28]中，僅列出了土樣A在4組不同孔隙率條件下的滲透系數(shù)的變化過(guò)程(如圖10所示)。圖10中實(shí)線(xiàn)為本文的計(jì)算值，上下2條虛線(xiàn)為2倍的計(jì)算值和1/2計(jì)算值。

圖10 滲透系數(shù)的驗(yàn)證

由圖10可知，土體的滲透系數(shù)隨著時(shí)間是不斷增大的，這是土體孔隙率變大而且土體的顆粒也變粗的結(jié)果[38]。本文的計(jì)算結(jié)果大致與試驗(yàn)值是一致的，試驗(yàn)測(cè)量的滲透系數(shù)值在1/2～2倍的計(jì)算值范圍內(nèi)，其中土樣A1和土樣A3與試驗(yàn)結(jié)果驗(yàn)證較好。事實(shí)上，滲透系數(shù)受多種因素的影響，目前的計(jì)算方法使得滲透系數(shù)在1/3～3倍的真實(shí)值范圍內(nèi)就已經(jīng)非常困難[38-39]。另外，試驗(yàn)中精準(zhǔn)測(cè)量土體的滲透系數(shù)也是非常困難的[40]，因此本文得出的滲透系數(shù)的計(jì)算值與試驗(yàn)值吻合度較好。

4 結(jié)果分析

4.1 土顆粒的運(yùn)動(dòng)速度

式(15)列出了土體顆粒的運(yùn)動(dòng)狀態(tài)方程，根據(jù)方程可知當(dāng)滲流力大于重力時(shí)土顆粒做加速運(yùn)動(dòng)，土顆粒的運(yùn)動(dòng)速度是關(guān)于時(shí)間的函數(shù)，且與土體的顆粒大小有關(guān)。以文獻(xiàn)[28]中的土體A1為例，計(jì)算不同粒徑的顆粒由靜起動(dòng)的運(yùn)動(dòng)狀態(tài)，如圖11所示。由圖11可知，土顆粒由靜起動(dòng)時(shí)的加速運(yùn)動(dòng)過(guò)程，其加速度是逐漸減小，直至減小為0。而不同粒徑的顆粒的運(yùn)動(dòng)狀態(tài)有所區(qū)別。當(dāng)土顆粒粒徑小于0.1 mm時(shí)，土顆粒的加速運(yùn)動(dòng)進(jìn)程明顯加快，顆粒運(yùn)動(dòng)達(dá)到平衡狀態(tài)的時(shí)間較短，且達(dá)到平衡狀態(tài)時(shí)的末速度稍大，這是由于小顆粒較為容易被水流帶走而運(yùn)動(dòng)。當(dāng)土顆粒粒徑大于0.1 mm時(shí)，土顆粒的運(yùn)動(dòng)狀態(tài)隨粒徑的變化區(qū)別不大，粒徑從0.1 mm變化到2 mm的過(guò)程中顆粒的加速過(guò)程以及加速后的平衡狀態(tài)差別不大。

圖11 不同土顆粒的運(yùn)動(dòng)狀態(tài)

對(duì)于同種土體，當(dāng)土體孔隙率不同時(shí)，也可導(dǎo)致同粒徑的顆粒運(yùn)動(dòng)狀態(tài)不同，對(duì)于土體A1～A4，其差別在孔隙率和滲透系數(shù)的不同，這樣對(duì)于同樣的顆粒大小d10，不同土體中的運(yùn)動(dòng)狀態(tài)如圖12所示。由圖12可知，在孔隙性增大的情況下，土顆粒的運(yùn)動(dòng)速度也會(huì)增大。這主要是因?yàn)榭紫缎栽龃蠛髮?dǎo)致滲透系數(shù)增大，在水力梯度不變的情況下，滲流流速會(huì)增大，這就導(dǎo)致土顆粒的所受滲流力增大，因此可見(jiàn)隨著孔隙率的增加，幾條曲線(xiàn)的斜率也在增加，土顆粒加速度增大。而土顆粒加速運(yùn)動(dòng)后的平衡速度也是隨著孔隙率的增加而增大，這也是孔隙內(nèi)水流流速較大造成的。

圖12 不同孔隙性的土體顆粒的運(yùn)動(dòng)狀態(tài)

4.2 土顆粒的運(yùn)動(dòng)概率和輸沙率

2.4.1節(jié)中給出了土顆粒運(yùn)動(dòng)概率的計(jì)算方法，根據(jù)該方法可知在滲流輸沙的過(guò)程中土顆粒的運(yùn)動(dòng)概率也是不斷變化的。圖13為土樣A1～A4在滲流出沙的過(guò)程中，土顆粒運(yùn)動(dòng)概率的變化特征。由圖13可知，出土顆粒的運(yùn)動(dòng)概率隨著輸沙過(guò)程是逐漸增大的。隨著土體細(xì)顆粒的不斷流失，土體的孔隙率和滲透性逐漸變大。根據(jù)式(22)所示，孔隙流速將有更大的概率大于臨界流速，這時(shí)導(dǎo)致剩余土顆粒較為容易的發(fā)生起動(dòng)，即P1將會(huì)增大。另外，由于細(xì)顆粒的不斷流失，土體逐漸粗化，這時(shí)根據(jù)式(23)和式(24)土體的平均孔隙直徑將會(huì)逐漸變大，土顆粒的輸移概率P2也會(huì)增大。

圖13 土顆粒運(yùn)動(dòng)概率

此外，通過(guò)與圖9的對(duì)比分析可以發(fā)現(xiàn)，土顆粒的運(yùn)動(dòng)概率與細(xì)顆粒的累積流失量具有時(shí)間上的相關(guān)性。在顆粒流失的前20 min內(nèi)，土顆粒運(yùn)動(dòng)概率增長(zhǎng)速度很快，同時(shí)顆粒的流失速率也不斷增強(qiáng)。在20 min之后，土顆粒流失曲線(xiàn)出現(xiàn)拐點(diǎn)，顆粒的流失量逐漸減緩，主要是因?yàn)橥令w粒運(yùn)動(dòng)概率的增速減緩。通過(guò)輸沙率(單位時(shí)間內(nèi)的輸沙量，如圖14所示)的變化規(guī)律也能夠看出顆粒運(yùn)動(dòng)概率與輸沙過(guò)程的相關(guān)關(guān)系。滲流輸沙率呈現(xiàn)出先增大后減小的趨勢(shì)，這與梁越的理論模型的計(jì)算結(jié)果是一致的[9]。在輸沙進(jìn)行20 min內(nèi)，輸沙率逐漸增大，而在30 min后，除A3土體外，其余土樣的顆粒運(yùn)動(dòng)概率幾乎停滯增長(zhǎng)，輸沙率也開(kāi)始大幅減小而出現(xiàn)了“拖長(zhǎng)尾”的現(xiàn)象。根據(jù)式(16)中輸沙率的計(jì)算方法，此時(shí)雖然顆粒運(yùn)動(dòng)速度和概率都較大，但是可動(dòng)顆粒逐步流失，可動(dòng)顆粒的數(shù)量成為制約滲流輸沙強(qiáng)度的主要因素。

圖14 不同土體的輸沙率

4.3 初始孔隙率對(duì)輸沙過(guò)程的影響

本文以陳亮試驗(yàn)中的土樣為例討論孔隙率對(duì)出沙量的影響，在陳亮的試驗(yàn)中，孔隙比的大小在0.43～0.5之間，對(duì)應(yīng)的孔隙率在0.3～0.33之間，孔隙率相差并不大，因此在其試驗(yàn)中，出沙量也是比較接近的。再者4種土樣所受的水力梯度并不同，無(wú)法準(zhǔn)確獲得孔隙率這一單一要素對(duì)結(jié)果的影響。因此，本文以試驗(yàn)中的土樣為例，假設(shè)土體的初始孔隙率分別為0.3、0.4和0.5，水力梯度都為0.6，計(jì)算不同孔隙率條件下的出沙量的值如圖15所示。由圖15可知，孔隙率越小，最終的出沙量越大；但是單位時(shí)間的出沙量(輸沙率)即圖中曲線(xiàn)的斜率，是隨著孔隙率的減小而減小。這是因?yàn)椋瑢?duì)于等體積的土體，孔隙率較大導(dǎo)致單位土體中的土顆粒的質(zhì)量較小，可以被水流帶走的細(xì)顆粒也較少，因此最終的出沙量是小于孔隙率較小的土體。而對(duì)于輸沙率的變化規(guī)律則不同，由于土體孔隙率較大，細(xì)顆粒較為容易的從孔隙中流出，因此單位時(shí)間的出沙量就大于孔隙率較小的土體。另外，如前文所述，圖15中曲線(xiàn)也可以看到有兩個(gè)拐點(diǎn)。在水力梯度一定的情況下，影響出沙量的因素包括孔隙率和可動(dòng)顆粒的數(shù)量。而在沙量逐漸流失的過(guò)程中，孔隙率不斷增大，這能促進(jìn)輸沙；可動(dòng)顆粒的數(shù)量逐漸減少，這能抑制輸沙。圖15中的輸沙率先是增大然后減小，這說(shuō)明在滲流輸沙的初期孔隙率是影響輸沙率的主要因素，而在輸沙后期，可動(dòng)顆粒的數(shù)量是輸沙率變化的主要控制因素。

圖15 初始孔隙率對(duì)出沙量的影響

以上3種孔隙率的土體，在0.6的水力梯度作用下，其滲透系數(shù)的變化如圖16所示。圖中可以看出，隨著細(xì)顆粒的流失，土體的滲透系數(shù)是逐漸增大的。而隨著孔隙率的增大，不同土體的滲透系數(shù)的增加幅度逐漸變大，即不同孔隙率土體滲透系數(shù)的差值隨時(shí)間逐漸變大。這就是說(shuō)，滲流輸沙對(duì)孔隙率越大的土體的滲透性的影響也越大。另外，圖16與圖15表現(xiàn)出相同的規(guī)律，滲透系數(shù)隨時(shí)間的變化也呈現(xiàn)出兩個(gè)拐點(diǎn)，這也能證明孔隙率和可動(dòng)顆粒的數(shù)量對(duì)輸沙率的影響。

圖16 初始孔隙率對(duì)滲透系數(shù)的影響

4.4 水力梯度對(duì)輸沙過(guò)程的影響

同樣地，本文以陳亮試驗(yàn)中的土樣為例討論水力梯度的影響。假設(shè)土體孔隙率為0.3，而水力梯度分別為0.4、0.6、0.8和1.0，此時(shí)計(jì)算的出沙量和滲透系數(shù)變化分別如圖17和圖18所示。由圖17、18可知，對(duì)于水力梯度0.4的土體，在200 min內(nèi)，其出沙量要遠(yuǎn)遠(yuǎn)小于其他3個(gè)水力梯度，滲透系數(shù)的增長(zhǎng)幅度也非常有限。這是由于0.4的水力梯度小于該土體的臨界水力梯度，因此在較長(zhǎng)時(shí)間內(nèi)，滲流輸沙量很少。其余3個(gè)水力梯度下，200 min內(nèi)土體的最終輸沙量都相等，差別在土體達(dá)到最大輸沙量的時(shí)間上。當(dāng)水力梯度較大，對(duì)滲流輸沙的加速效應(yīng)顯著，輸沙所需時(shí)間明顯縮短。而且隨著時(shí)間的推移，這種加速效應(yīng)將逐漸增強(qiáng)。由圖18可知，當(dāng)水力梯度增加到1.0時(shí)，細(xì)顆粒將在30 min作用完全流失。在此過(guò)程中，較大水力梯度引起的輸沙量增長(zhǎng)迅速，最大能達(dá)到水力梯度為0.6情況下的3倍以上。同時(shí)伴隨著土體滲透系數(shù)的快速增加。該過(guò)程也可通過(guò)本文的理論模型進(jìn)行解釋。當(dāng)水力梯度較大，孔隙流速也更大，可動(dòng)顆粒所受滲流力(如式(4)和式(5)所示)更大，從而導(dǎo)致土顆粒具有更大的運(yùn)動(dòng)速度和運(yùn)動(dòng)概率，滲流輸沙的發(fā)展速度也將更快。

圖17 水力梯度對(duì)出沙量的影響

圖18 水力梯度對(duì)滲透系數(shù)的影響

5 結(jié) 論

本文針對(duì)滲流作用下泥沙輸運(yùn)過(guò)程，提出了泥沙顆粒運(yùn)動(dòng)的控制方程，并以此為基礎(chǔ)建立理論模型，提出了滲流作用下輸沙率的計(jì)算方法。利用試驗(yàn)數(shù)據(jù)驗(yàn)證了本理論模型的合理性。利用該模型對(duì)滲流作用下的泥沙運(yùn)動(dòng)狀態(tài)和影響因素進(jìn)行了分析，得出以下結(jié)論：

(1)本文提出由土顆粒單顆粒受力規(guī)律出發(fā)，推導(dǎo)了土顆粒的運(yùn)動(dòng)方程，該方程可用來(lái)計(jì)算土體的臨界水力梯度，通過(guò)與Terzaghi方法的對(duì)比，發(fā)現(xiàn)該計(jì)算方法可信度較高。

(2)土顆粒在由靜起動(dòng)的過(guò)程中做加速運(yùn)動(dòng)，加速歷程很短，其運(yùn)動(dòng)特征與土顆粒粒徑及孔隙率有關(guān)；滲流輸沙過(guò)程中顆粒運(yùn)動(dòng)概率不斷增加，主要是因?yàn)轭w?？紫吨饾u變大，可動(dòng)顆粒受力增加，輸運(yùn)的幾何容許度增加引起的。通過(guò)顆粒運(yùn)動(dòng)速度和運(yùn)動(dòng)概率可計(jì)算滲流輸沙率，輸沙過(guò)程與土顆粒運(yùn)動(dòng)概率具有很好的一致性。

(3)土體孔隙率和所處的水力梯度是影響滲流輸沙率的重要因素。孔隙率較大，輸沙率越大，總輸沙量越小。水力梯度對(duì)輸沙過(guò)程有明顯的加速效應(yīng)，且當(dāng)水力梯度低于臨界水力梯度時(shí)，土體顆粒流失明顯變慢。

(4)通過(guò)對(duì)出沙量和滲透系數(shù)的變化規(guī)律的分析，得出在輸沙的初始階段，孔隙率是控制輸沙率的主要因素，而在輸沙后期，可動(dòng)顆粒的數(shù)量是主要的限制因素。