滲流輸沙過程的顆粒運動狀態和影響因素分析

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

(1.河北省海洋地質資源調查中心，河北 秦皇島 066000；2.河北省海洋岸線生態修復與智慧海洋工程研究中心，河北 秦皇島 066000；3.天津大學水利工程仿真與安全國家重點實驗室，天津 300350)

1 研究背景

滲流導致的土體內部泥沙顆粒運動是壩體、河堤失穩破壞的重要原因[1-2]。根據泥沙運動形式的不同，滲流輸沙分為管涌、流土、接觸流土、接觸沖刷等多種形式。據統計，世界上約有50%的壩體破壞是由于滲流所致的泥沙流失引起的[3-4]。因此研究滲流作用下的泥沙運動規律具有重要的實際意義。

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

q=Kx(X-Xcr)a

(1)

式中，q為輸沙率；Kx和a為與土體相關的參數；X為土體的滲流場參數，根據上文，X可表示滲流切應力、滲流流量、滲流梯度、滲流速度等；Xcr為發生泥沙運動時X的臨界值。

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

本文從泥沙顆粒的受力狀態出發，推導泥沙顆粒在滲流場中運動的控制方程。通過分析滲流流速的概率分布模式，將單顆粒泥沙運動規律擴展到整個滲流場，并由此計算得到泥沙運動概率和輸沙規律。通過和試驗結果的對比，表明該模型具有合理性，可以對一維滲流作用下的泥沙運動規律進行描述。

2 材料和方法

2.1 動力學方程

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

圖1 土體內部可動顆粒的受力分析

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

F=FI+FD

(2)

(3)

式中，FD為水流拖曳力；FI為滲流梯度力；D為顆粒直徑；ρ′為顆粒的水下密度。

FD[26]和FI[31]的表達式分別為

(4)

(5)

式中，α為顆粒形狀修正系數，對于球形取1，對于泥沙顆粒取1.7；CD為拖曳力系數；ρw為水的密度；ui為顆粒運動速度；ubi為孔隙內水流速度；I為水力梯度。

根據牛頓第二定律

(6)

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

(7)

式中拖曳力系數與雷諾數有關，滲流作用下的顆粒雷諾數Re的計算公式為

(8)

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

(9)

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

(10)

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

2.2 臨界起動流速

滲流梯度和滲流流速的關系服從達西定律

ud=KI

(11)

式中，K為滲透系數；ud為滲流流速。

引入平均孔隙流速，滲流流速與平均孔隙流速的關系為

(12)

式中，n為土體孔隙率。

(13)

式(13)為顆粒運動的臨界孔隙流速，其與臨界水力梯度的關系為

(14)

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

2.3 顆粒運動狀態分析

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

(15)

式(15)即為孔隙內單個土顆粒的運動速度表達式。

2.4 輸沙率的確定

根據泥沙推移質運動理論，垂直滲流作用下，單位面積單位時間內的輸沙率可由下式計算：

(16)

式中，q為單位時間單位面積的滲流斷面上的輸沙率；N為單位體積內可動顆粒的個數；V為可動顆粒的體積；P為可動顆粒運動概率。

N和V的表達式為

(17)

(18)

式中，PV為可動顆粒的體積比。

2.4.1P的確定方法

因為土體孔隙有大有小，孔隙內流速也呈大小不一的分布規律。只有孔隙內的流速大于土顆粒的臨界起動流速，土顆粒才會被水流帶走。但是如果顆粒運動軌跡中的下游孔隙小于顆粒直徑時，土顆粒依然不能被水流帶走。因此顆粒運動還應該滿足一個條件，本文稱為顆粒的輸移條件，即顆粒輸移必須滿足下游孔隙大于可動顆粒的直徑。所以這里的P包括兩部分P1和P2，且P滿足

P=P1P2

(19)

式中，P1為土顆粒的起動概率；P2為土顆粒的輸移概率。假設土體內水流流動符合管流泊肅葉理論，管流速度和圓管半徑的關系為[13]

(20)

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

(21)

因此，P1的計算公式為

(22)

根據Hazen的研究[33]，土體孔隙直徑和顆粒級配有關：

D0=Bd10

(23)

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

P2=P(D

(24)

式中可根據土體級配的分布特征確定可動顆粒D0的大小，然后根據顆粒的質量比計算D小于D0的概率。

(25)

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

2.5 滲流-泥沙場的耦合計算

前文已有論述，由于滲流輸沙的過程伴隨著滲流場的變化，因此對于輸沙率的計算需要耦合滲流場的變化。t時刻相對于t-1時刻，土體的孔隙率n、滲透系數K、可動顆粒的體積分數Pv會發生變化。而其他變量，比如臨界流速ucr、孔隙流速ub、起動概率P、孔隙半徑R都可以由以上幾個變量求解得到。因此，在已知t-1時刻各參數的值的情況下，只需確定t時刻n、K和PV的計算方法，就可實現時間上的迭代計算，進而實現滲流-泥沙場的耦合計算。

由輸沙率q的定義可得在t時刻，總輸沙量Qt的計算公式為

(26)

假設在輸沙過程中，沒有土體的變形，即輸沙的體積即為土體中孔隙增加的體積。則

(27)

根據劉杰的研究[32]，滲透系數的比值與孔隙率和d20有關。

(28)

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

(29)

式中，(dx)0為可動顆粒的最大粒徑；(Pm)t為粒徑小于最大可動顆粒的土體所占的質量分數，其與輸沙量Qt有關，即

(30)

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

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

(31)

根據以上論述可通過式(27)、式(28)和式(31)實現滲流-泥沙場耦合計算，耦合計算流程如圖2所示。

圖2 耦合計算流程

3 模型驗證

3.1 臨界水力梯度的驗證

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

圖3 土樣級配曲線(Skempton試驗)

圖4 土樣級配曲線(Ahlinhan試驗)

圖5 土樣級配曲線(Fleshman試驗)

圖6 土樣級配曲線(陳亮試驗)

圖7 臨界水力梯度的驗證

圖8 本文公式與太沙基方法的對比

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

3.2 出沙量的驗證

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

圖9 累積出沙量的驗證

由圖9可以看出，本文的計算值與試驗值非常接近。出沙量隨時間逐漸增多，且呈現出兩個拐點。單位時間的出沙量，即輸沙率是隨著時間先增大，然后逐漸減小，最后基本為0。

3.3 滲透系數的驗證

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

圖10 滲透系數的驗證

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

4 結果分析

4.1 土顆粒的運動速度

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

圖11 不同土顆粒的運動狀態

對于同種土體，當土體孔隙率不同時，也可導致同粒徑的顆粒運動狀態不同，對于土體A1～A4，其差別在孔隙率和滲透系數的不同，這樣對于同樣的顆粒大小d10，不同土體中的運動狀態如圖12所示。由圖12可知，在孔隙性增大的情況下，土顆粒的運動速度也會增大。這主要是因為孔隙性增大后導致滲透系數增大，在水力梯度不變的情況下，滲流流速會增大，這就導致土顆粒的所受滲流力增大，因此可見隨著孔隙率的增加，幾條曲線的斜率也在增加，土顆粒加速度增大。而土顆粒加速運動后的平衡速度也是隨著孔隙率的增加而增大，這也是孔隙內水流流速較大造成的。

圖12 不同孔隙性的土體顆粒的運動狀態

4.2 土顆粒的運動概率和輸沙率

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

圖13 土顆粒運動概率

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

圖14 不同土體的輸沙率

4.3 初始孔隙率對輸沙過程的影響

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

圖15 初始孔隙率對出沙量的影響

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

圖16 初始孔隙率對滲透系數的影響

4.4 水力梯度對輸沙過程的影響

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

圖17 水力梯度對出沙量的影響

圖18 水力梯度對滲透系數的影響

5 結 論

本文針對滲流作用下泥沙輸運過程，提出了泥沙顆粒運動的控制方程，并以此為基礎建立理論模型，提出了滲流作用下輸沙率的計算方法。利用試驗數據驗證了本理論模型的合理性。利用該模型對滲流作用下的泥沙運動狀態和影響因素進行了分析，得出以下結論：

(1)本文提出由土顆粒單顆粒受力規律出發，推導了土顆粒的運動方程，該方程可用來計算土體的臨界水力梯度，通過與Terzaghi方法的對比，發現該計算方法可信度較高。

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

(3)土體孔隙率和所處的水力梯度是影響滲流輸沙率的重要因素?？紫堵瘦^大，輸沙率越大，總輸沙量越小。水力梯度對輸沙過程有明顯的加速效應，且當水力梯度低于臨界水力梯度時，土體顆粒流失明顯變慢。

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