王 霜，陳建生，鐘啟明

(1.河海大學水利水電學院，江蘇 南京 210098；2.南京鐵道職業(yè)技術學院鐵道工程系，江蘇 南京 210031；3.河海大學土木與交通學院，江蘇 南京 210098；4.南京水利科學研究院巖土工程研究所，江蘇 南京 210029)

2020年汛期，我國南方地區(qū)發(fā)生多輪強降雨，造成多地發(fā)生嚴重洪澇災害，江西贛江大堤、安徽河段堤防等多地出現管涌險情。管涌是江河大堤在汛期的常見險情之一，是造成大壩失事的一種很重要的滲透破壞形式[1-3]，工程意義上的堤基管涌通常將接觸沖刷、流土等滲透變形現象也歸結其中。

多年來，眾多學者采用試驗、數值模擬等方法對管涌問題進行了廣泛深入的研究。劉杰等[4-5]在室內進行了由滲流引起的堤基管涌試驗，論證了管涌的發(fā)生以及發(fā)展過程對堤壩安全穩(wěn)定性的影響。李廣信等[6]采用彩色砂作為示蹤材料，模擬了二元堤基結構及不同防滲墻深度下堤基內部管涌的發(fā)生與發(fā)展過程。羅玉龍等[7]對管涌發(fā)展機理提出了新認識：管涌是涉及孔隙水滲流、可動顆粒運移、多孔介質變形等眾多復雜行為的多相多場耦合過程。Rochim等[8]在水力梯度控制或流量控制條件下，對4種不同的無黏性土進行了滲流管涌試驗。研究表明，水力加載歷史對管涌臨界水力梯度有顯著影響。明攀等[9]利用聲發(fā)射技術實時監(jiān)測堤防的管涌過程，發(fā)現管涌過程的水力參數和聲發(fā)射參數具有相同的分布規(guī)律，可通過聲發(fā)射參數對管涌過程進行判別。

隨著計算機應用的迅速發(fā)展，數值模擬漸漸成為管涌研究的重要方法。殷建華[10]在國內首次采用有限元方法來研究管涌問題，計算了管涌發(fā)生區(qū)域的長度和滲透系數對堤身滲流的影響。倪曉東等[11]運用顆粒流模擬揭示了非穩(wěn)定流作用下管涌發(fā)展的細觀機制。Hu等[12]采用計算流體力學和離散元耦合的方法，對缺乏中間粒徑的管涌型土和級配良好的土體進行了管涌侵蝕的微觀和宏觀研究。賈愷等[13]考慮了雙層堤基管涌通道斷面擴展和上溯與通道尖端坡降之間的相互影響，通過編制有限元程序模擬研究了管涌通道發(fā)展趨勢。

目前，關于管涌的研究對象大多集中在雙層堤基以及堤基內部土體。事實上，堤防總是沿江河而建，江河兩岸由于歷史上的河流沖刷、淤積和改道等，堤基結構存在很大差異，大體上可以概括為單層、雙層、三層和多層4種類型[14-15]。對于不同的堤基結構，管涌的發(fā)生和發(fā)展過程存在很大差異。本文針對多層堤基結構，利用有限元分析軟件對堤基管涌的動態(tài)發(fā)展過程進行模擬。顆粒流失影響局部土體滲透性，進而改變滲流場的分布以及堤基內各區(qū)域的水力梯度，繼而又會影響土體顆粒流失，模擬過程中考慮了土顆粒流失和土體滲透性及其相互影響、相互耦合的過程。

1 管涌臨界水力梯度確定

通過有限元數值模擬得到不同上游水頭下堤基內部的水力梯度分布情況，將堤基內部水力梯度的分布與該區(qū)域土體發(fā)生破壞的臨界水力梯度進行比較，從而確定在該上游水頭下堤基內部是否有顆粒侵蝕，若無侵蝕，繼續(xù)逐級增加上游水頭，最終得到堤基內部土體發(fā)生侵蝕的上游水頭，也就是管涌發(fā)生的臨界上游水頭，同時可以得到該水頭下堤基內部顆粒侵蝕的位置與區(qū)域。

本文以細砂夾層埋入砂礫石深度為1 m的多層堤基為例，各層土體從上到下依次為黏土層、砂礫石層、細砂夾層、砂礫石層。砂礫石層以及細砂夾層的土料級配如圖1所示，土體的物理參數如表1所示，上覆黏土的滲透系數為5×10-6m/s。

表1 堤基土體物理參數

圖1 堤基內土體級配曲線

多層堤基中，由于滲透系數較小的細砂夾層的存在，將堤基的砂礫石層分成了細砂夾層上部與細砂夾層下部兩部分。堤基滲透破壞發(fā)生的標志主要有兩種形式：細砂夾層上部，管涌口下端砂礫石發(fā)生豎向侵蝕破壞，只有管涌口處顆粒首先沿豎向發(fā)生侵蝕形成空洞，才能為水平向的顆粒運移提供足夠的空間，侵蝕才能沿水平向逐漸向上游發(fā)展；細砂夾層下部，細砂夾層沖破，發(fā)生豎向流土破壞，同樣細砂夾層豎向流土破壞為水平向發(fā)展提供了空間。多層堤基管涌發(fā)生的標志主要取決于這兩種滲透破壞形式哪一種的臨界水力條件先得到滿足。

對于管涌口下端砂礫石豎向侵蝕破壞，其水力梯度Jc[16]為

(1)

式中：d為管涌時被沖蝕細粒粒徑，根據試驗觀測可取質量分數小于3%顆粒對應的粒徑d3，cm；k為砂礫石的滲透系數，cm/s；n1為砂礫石的孔隙率。將砂礫石土料的各參數代入，得到砂礫石顆粒發(fā)生豎向侵蝕的臨界水力梯度為0.36。

(2)

細砂夾層及其上部土顆粒發(fā)生豎向流土破壞的臨界水力梯度J為

(3)

將堤基土體的物理參數代入式(2)計算可得，當細砂夾層埋深為1 m時，細砂夾層及上部土體被沖破所需的細砂層底部豎向水力梯度為2。

本文采用逐級提升水頭的方式，模擬管涌集中滲漏通道的發(fā)生與發(fā)展過程。當堤基管涌口下端砂礫石的水力梯度小于0.36或者細砂夾層底部水力梯度小于2時，滲透破壞未發(fā)生，顆粒不發(fā)生流失，該種情況文中未做分析。

2 多層堤基管涌發(fā)展過程的模擬

2.1 上游水頭H=14 m

采用有限元軟件對細砂夾層埋深為1 m的多層堤基內部滲流場進行模擬，模型以向堤內方向為x軸正方向、向上為y軸正方向建立坐標系，堤基向兩邊無限延伸，模擬的是半無限邊界，上游側水頭的邊界條件作用于模型左側，下游側水頭的邊界條件作用于模型右側的黏土層表面，假定為零勢面。本文不考慮上部堤防滲流對堤基內部滲流場產生的影響，假定堤防為不透水體，同時假定距堤腳26 m處為黏土層薄弱處，薄弱區(qū)域直徑為4 m，即在模型中將x軸坐標為84～88 m的黏土層定為管涌口，滲透系數設定為0.5 m/s,計算結果如圖2和圖3所示。

圖2 多層堤基管涌口形成后水頭等值線(單位：m)

由圖2可以看出，細砂夾層底部水頭分布較均勻，斷面流量較小，管涌口下端的細砂夾層承擔了較大的水頭，細砂夾層底部y方向的水力梯度最大，最大水力梯度達到2.14(圖3(a))，超過了計算所得細砂夾層沖破的臨界水力梯度2，細砂夾層沖破，發(fā)生管涌破壞。而管涌口下端黏土層與砂礫石層界面砂礫石所承擔的水頭差非常小，豎向水力梯度幾乎為零(圖3(b)y=28 m處)，小于顆粒豎向運動的允許水力梯度，從而認定管涌口下端黏土層與砂礫石層界面豎向侵蝕破壞并未發(fā)生。因而，對于細砂夾層埋深為1 m的多層堤基，管涌破壞發(fā)生的標志為細砂夾層沖破，臨界上游水頭為14 m。

圖3 不同區(qū)域y方向水力梯度

由于細砂夾層的顆粒粒徑較小，管涌發(fā)生時，假設管涌口下端的細砂夾層顆粒被沖散，全部帶出管涌口。細砂夾層下部砂礫石的沖蝕深度同樣根據豎向水力梯度來確定，由圖3(b)可以看出，細砂夾層下部離細砂夾層越近，豎向水力梯度越大，細顆粒侵蝕粒徑也就越大，因而可以假設細砂夾層沖破時夾層下部1 m深度的砂礫石和上部的砂礫石被沖散，細顆粒被全部帶出，粗顆粒沉積下來。

管涌口在涌砂過程中，由于細砂夾層顆粒被沖散，帶出管涌口，上部砂礫石層顆粒下沉，管涌口下部形成與細砂夾層相等厚度的1 m深空洞，空洞的滲透系數與管涌口保持一致，設為0.5 m/s。細砂夾層上部被沖散的砂礫石下沉填充了原細砂夾層的位置，與細砂夾層下部1 m深度范圍內的沖蝕砂礫石形成了一個侵蝕區(qū)域，侵蝕區(qū)域如圖4所示，在管涌口下方y(tǒng)=25～27 m粗實線區(qū)域內，侵蝕區(qū)域內細顆粒流失，涌出的細顆粒引起相應區(qū)域孔隙率發(fā)生改變，則

(4)

式中：nt為管涌侵蝕區(qū)的孔隙率；Vvt為涌砂體積；Vv為管涌侵蝕區(qū)初始孔隙體積；V為管涌侵蝕區(qū)總體積；n0為初始孔隙率。

根據Kozeny-Carman方程[17]，土體滲透系數與孔隙率有如下關系：

(5)

土體未發(fā)生侵蝕前，土體的初始孔隙率為n0，初始滲透系數為k0，則

(6)

(7)

根據上述假設，管涌口下端砂礫石中的細顆粒全部流失，同時結合土體的初始顆粒級配，計算可得顆粒發(fā)生流失后的土體孔隙率nt以及土體等效粒徑Dht，代入式(7)得管涌口下端侵蝕區(qū)域砂礫石的滲透系數為2.8×10-2m/s。

圖4為細砂夾層沖破后x方向的水力梯度等值線圖，可以看出，水力梯度最大的區(qū)域發(fā)生在管涌口前端細砂夾層底部附近。為了進一步研究多層堤基管涌發(fā)生后管涌通道的位置以及管涌通道是如何發(fā)展的，需要對不同區(qū)域土體顆粒發(fā)生運移的臨界水力梯度進行分析和確定。

圖4 細砂夾層沖破后x方向水力梯度等值線

對于堤基內部砂礫石顆粒在水平滲流下所引起的水平管涌，其臨界水力梯度JH[16]為

(8)

式中：α為形狀系數，根據明茲[18]的試驗資料，對于各種砂粒α=1.16～1.17，對于銳角顆粒α=1.50～1.67；d為管涌時被沖動的細粒粒徑，根據試驗觀測可取d3；tanφ為砂礫石的內摩擦系數；γs1為砂礫石的容重;d3K為等效粒徑，可采用下式計算[18]：

(9)

式中：A為層流阻力系數，A=5.1[18]。將式(9)代入式(8)并采用一些常用數據γw=9.81 kN/m3、γs1=26 kN/m3、g= 981 cm/s2、ν= 0.011 6 cm2/s、α=1.165，可得管涌臨界水力梯度的近似公式：

(10)

將砂礫石土料的各參數代入式(10)，可得砂礫石顆粒發(fā)生水平向管涌的臨界水力梯度為0.24。

圖5(a)為黏土層與砂礫石層界面x方向水力梯度分布圖，可見水力梯度在管涌口邊界x=84 m處有所減小，該處水力梯度的減小主要受管涌口的影響，管涌口處的滲透系數很大，因而水頭差較小，水力梯度很?。辉诠苡靠谇岸?，由于管涌口處的水頭為0，水頭在管涌口前端發(fā)生突變，水力梯度最大，達到了0.29，大于顆粒發(fā)生水平向管涌的臨界水力梯度0.24，因而在黏土層與砂礫石層界面顆粒發(fā)生侵蝕，管涌通道在此界面最終發(fā)展到x=79 m處。將管涌口前端的水力梯度0.29代入式(10)反算，可得到顆粒粒徑小于0.15 mm的顆粒被侵蝕。為了簡化計算，假設管涌通道上不同位置處的顆粒侵蝕均相同，侵蝕顆粒均為粒徑小于0.15 mm的顆粒。將顆粒侵蝕后的土體孔隙率nt以及等效粒徑Dht代入式(7)，同時結合土體未被侵蝕時的初始參數，最終可計算出管涌通道的滲透系數為1.68×10-3m/s。

圖5 細砂夾層沖破后不同區(qū)域x方向水力梯度

圖5(b)為細砂夾層底部與砂礫石層界面x方向水力梯度分布圖，可見細砂夾層底部水力梯度大于臨界水力梯度的區(qū)域為x=79～84 m，因而管涌通道在此界面形成，且通道在細砂夾層底部發(fā)展到x=79 m處。為了簡化計算，假設細砂夾層底部和上部的砂礫石顆粒侵蝕相同，細砂夾層底部的通道滲透系數同樣采用1.68×10-3m/s。在該水頭下，細砂夾層上部和下部的砂礫石均發(fā)生侵蝕，也就是深層管涌破壞與淺層管涌破壞同時作用。

對于堤基內部細砂夾層在水平滲流下所引起的流土破壞，細砂顆?；旧峡烧J為是均一粒徑的，故式(8)中的d3K=d，于是得到水平向流土破壞的臨界水力梯度JB[16]計算公式：

(11)

式中：tanφ2為細砂的內摩擦系數，細砂的內摩擦角φ2可取28°。將細砂的各參數代入式(11)可得細砂夾層發(fā)生水平向流土破壞的臨界水力梯度為0.54；根據細砂夾層水平向流土破壞的臨界水力梯度與細砂內摩擦角之間的關系，同樣可以計算出細砂夾層發(fā)生豎向流土破壞的臨界水力梯度為1.0。

細砂夾層沖破后，細砂夾層底部的水力梯度最大(圖5)，細砂夾層如果發(fā)生流土破壞，發(fā)生在細砂夾層底部的可能性最大，然而，由圖5(b)細砂夾層底部x方向水力梯度分布可知，細砂夾層底部水力梯度小于發(fā)生水平向流土破壞的臨界水力梯度，因此，在該上游水頭下細砂夾層不會發(fā)生水平向的流土破壞。

2.2 上游水頭H=18 m

在上游水頭為14 m時，管涌通道形成，管涌通道在細砂夾層上部和下部均發(fā)展到x=79 m處，通道的滲透系數為1.68×10-3m/s，管涌口下端被侵蝕的砂礫石的滲透系數為2.8×10-2m/s，保持上游水頭14 m不變，管涌通道不繼續(xù)發(fā)展，此時提升上游水頭至18 m，得到管涌通道形成后在該水頭下x方向的水力梯度等值線如圖6所示。由于通道內顆粒侵蝕的粒徑較小，滲透系數增大的幅度不大，因而水力梯度最大的區(qū)域仍發(fā)生在管涌口前端細砂夾層底部附近，而不是管涌通道前端。

圖6 管涌通道形成后x方向水力梯度等值線

土體破壞區(qū)域的分析與上述分析過程相同，此處不再贅述。圖7(a)為黏土層與砂礫石層界面x方向水力梯度分布圖，可見水力梯度在管涌通道前端x=79 m處減小，然后在管涌通道區(qū)域內又會有小幅增大，在管涌口邊界水力梯度再次減小，由此可知，水力梯度的變化規(guī)律與該界面各區(qū)域土體滲透系數的不同有關，水力梯度在滲透系數增大處發(fā)生減??；在18 m上游水頭下，管涌通道在此界面繼續(xù)發(fā)展到x=70 m。

圖7(b)為細砂夾層底部與砂礫石層界面x方向水力梯度分布圖，同樣地，水力梯度在管涌通道前端和管涌口邊界減??；在18 m上游水頭下，細砂夾層底部砂礫石中的管涌通道繼續(xù)發(fā)展到x=73 m，細砂夾層未發(fā)生水平向的流土破壞。

圖7 管涌通道形成后不同區(qū)域x方向水力梯度

由于管涌通道內水力梯度分布并不相同，因而顆粒侵蝕的粒徑不同，最終使得通道內各處的滲透系數不同，為了簡化計算，假設土體侵蝕后，管涌通道內的滲透系數保持1.68×10-3m/s不變；同樣地，管涌口下端的砂礫石中有顆粒不斷地被帶出管涌口，也有水平方向來的顆粒不斷補充，其滲透系數是復雜多變的，為了簡化計算，假設管涌口下端砂礫石的滲透系數為2.8×10-2m/s不變，在下文的數值模擬中均保持管涌通道的滲透系數和管涌口下端砂礫石的滲透系數不變，不再重復說明。

2.3 上游水頭H=22 m

保持上游水頭18 m不變，管涌通道在細砂夾層上部和下部分別發(fā)展到x=70 m和x=73 m后不再繼續(xù)發(fā)展，此時提升上游水頭至22 m，圖8為該上游水頭下x方向的水力梯度等值線圖，水力梯度最大的區(qū)域仍然發(fā)生在管涌口前端細砂夾層底部附近。由圖9(a)可知，水力梯度同樣在管涌通道前端和管涌口邊界減小，在管涌通道內增加；通過水力梯度分布可以看出管涌通道繼續(xù)在黏土層與砂礫石層界面發(fā)展到x=42 m處。

圖8 管涌通道發(fā)展過程中x方向水力梯度等值線

圖9 管涌通道發(fā)展過程中不同區(qū)域x方向水力梯度

圖9(b)為細砂夾層底部與砂礫石層界面x方向水力梯度分布圖，同樣地，水力梯度變化規(guī)律與細砂夾層上部類似；然而，在22 m上游水頭下，細砂夾層底部的水力梯度最大值達到了0.54，正好等于由式(9)計算出的細砂發(fā)生水平向流土破壞的臨界水力梯度，因此細砂夾層在x=82 m處發(fā)生水平向的流土破壞，假設發(fā)生流土破壞的厚度為細砂夾層厚度的一半，同時，細砂夾層底部砂礫石中的管涌通道繼續(xù)發(fā)展到x=48 m。

當細砂夾層底部發(fā)生流土破壞后，上部的細砂和砂礫石發(fā)生沉降，使得黏土層與砂礫石層之間形成間隙，設間隙處的滲透系數與管涌口一致，取0.5 m/s，間隙的產生使得黏土層與砂礫石層界面變得相對薄弱，更易于侵蝕的發(fā)生。細砂夾層上部和下部管涌通道分別發(fā)展到x=42 m和x=48 m，且細砂夾層發(fā)生流土破壞后，保持上游水頭22 m不變，得到間隙產生后x方向水力梯度等值線如圖10所示，可見最大水力梯度發(fā)生在流土破壞后的細砂夾層前端。

圖10 間隙產生后x方向水力梯度等值線

由圖11可以看出，細砂夾層上部和下部砂礫石層中的管涌通道繼續(xù)發(fā)展最終貫通。在本文的數值模擬中，假定相同粒徑的顆粒被侵蝕，且管涌通道形成后，其滲透系數不變，由于侵蝕的顆粒粒徑較小，通道的滲透系數相對較小，管涌通道貫通后還存在較大的水力梯度，在水流的作用下，通道中的土體繼續(xù)侵蝕，管涌通道逐漸擴寬。

圖11 間隙產生后不同區(qū)域x方向水力梯度

圖12為間隙前端D截面y方向水力梯度，可以看出，間隙下部細砂夾層處y方向的水力梯度最大，達到了1.15，大于細砂夾層發(fā)生豎向流土破壞的臨界水力梯度1.0，因此，細砂夾層在豎向發(fā)生流土破壞，黏土層與砂礫石層之間的間隙逐漸向深部和上游發(fā)展。

圖12 間隙前端D截面y方向水力梯度

因而，在多層堤基管涌發(fā)展過程中，細砂夾層的存在將土體分為了多個區(qū)域，每個區(qū)域破壞發(fā)生的條件不同，在細砂夾層埋深為1 m的多層堤基結構中，通道發(fā)展主要經歷了細砂夾層上部和下部砂礫石層的水平向管涌破壞、細砂夾層水平向流土破壞和豎向流土破壞等破壞形式，是深層管涌破壞與淺層管涌破壞同時作用，產生的管涌通道深度很大，對大壩安全產生重大影響。

3 結 論

a.多層堤基結構存在滲透系數較小的細砂夾層，對于細砂夾層埋深為1 m的結構，管涌發(fā)生的標志為管涌口下端細砂夾層沖破，發(fā)生豎向流土破壞，顆粒大量流失，堤基內部滲流場重新分布。

b.管涌發(fā)展過程中，堤基內部土顆粒不斷流失，局部土體滲透性隨之發(fā)生變化，滲透性的變化影響滲流場的分布以及堤基內各區(qū)域的水力梯度，繼而又會影響土顆粒的流失，因而土顆粒流失和土體滲流是一個相互影響、相互耦合的過程。

c.管涌發(fā)展過程中存在多種類型的破壞形式，主要經歷了細砂夾層上部和下部砂礫石層的水平向管涌破壞、細砂夾層水平向流土破壞和豎向流土破壞等破壞形式，是深層管涌破壞與淺層管涌破壞同時作用，產生的管涌通道深度很大，對大壩安全產生重大影響。