海塘滲透破壞數(shù)值模擬及滲透特性

胡云進，鐘 振，賀春雷，王國慶

（1.浙江大學 建筑工程學院，杭州 310058；2.浙江省錢塘江管理局，杭州 310016）

洪水滲透下海塘內(nèi)的滲流是一種涉及土體滲透性變化的飽和非飽和／非穩(wěn)定滲流。某些情況下，海塘內(nèi)的滲流很難達到穩(wěn)定滲流狀態(tài)，滲透穩(wěn)定性就較高。而在另外一些情況下，海塘內(nèi)的滲流會較快達到穩(wěn)定滲流狀態(tài)，并在穩(wěn)定滲流條件下長期工作，滲透穩(wěn)定性就較低，易發(fā)生滲透破壞。因此，洪水滲透下海塘達到穩(wěn)定滲流所需時間相差較大的話，采用穩(wěn)定滲流分析方法得到相同的安全系數(shù)，海塘安全度實際上是不同的。而中國現(xiàn)行的《堤防工程設計規(guī)范》規(guī)定，采用穩(wěn)定滲流分析方法計算出逸坡降，并以此評價海塘的滲透穩(wěn)定性，這在很多情況下是不妥當?shù)?。為了準確評價海塘的滲透穩(wěn)定性，需要研究海塘洪水滲透從非飽和／非穩(wěn)定滲流達到飽和／穩(wěn)定滲流的整個過程，以模擬海塘滲透破壞的發(fā)生和發(fā)展過程。而為了進行洪水滲透下海塘內(nèi)滲流從非飽和／非穩(wěn)定狀態(tài)發(fā)展到飽和／穩(wěn)定狀態(tài)以及海塘滲透破壞動態(tài)發(fā)展過程的數(shù)值模擬，就必須采用飽和非飽和／非穩(wěn)定滲流分析方法，同時考慮土體滲透性隨海塘滲透破壞而變化的問題。從檢索文獻來看，在現(xiàn)有的海塘（或堤防）滲透破壞數(shù)值模擬方法中，大多采用增大滲透破壞區(qū)內(nèi)土體滲透系數(shù)的方法來模擬滲透破壞發(fā)展過程，但只進行飽和滲流分析，未考慮從非飽和／非穩(wěn)定滲流到飽和／穩(wěn)定滲流的階段，具體參見文獻［1－6］。有進行飽和非飽和／非穩(wěn)定滲流分析的文獻，但滲流分析時沒有考慮土體滲透性隨海塘滲透破壞而變化的問題，具體參見文獻［7－10］。因此為了獲得更準確的洪水滲透下海塘滲透破壞模擬結果，需要建立可考慮土體滲透性隨海塘滲透破壞而變化的飽和非飽和／非穩(wěn)定滲流分析模型。

本文針對目前數(shù)值模擬方法的不足，首先建立可考慮土體滲透性隨海塘滲透破壞而變化的飽和非飽和／非穩(wěn)定滲流分析模型，并編制相應的有限元計算程序；然后選取實際海塘進行洪水滲透過程的數(shù)值模擬，并將計算結果與現(xiàn)場監(jiān)測試驗結果進行對比，驗證分析模型和計算程序的正確性；最后根據(jù)數(shù)值模擬結果，總結歸納洪水滲透下海塘的滲透特征及規(guī)律。

1 飽和非飽和／非穩(wěn)定滲流數(shù)學模型

海塘為典型的線狀建筑物，其滲流計算一般選典型斷面按二維問題進行分析，控制方程如下［9］：

1.1 控制方程

式中：kr（h）為相對滲透系數(shù)；kij為飽和滲透系數(shù)；C（h）為容水度；h為壓力水頭；kr和C均為壓力水頭的函數(shù)；β為選擇參數(shù)，非飽和區(qū)為0，飽和區(qū)為1；SS為等效單位貯存量；S為源（匯）項；xi為坐標軸，其中x2為正向向上的鉛直軸；t為時間。

有了符合海塘土質(zhì)特性的h～θ、kr（h）～h及C（h）～h關系后（θ為含水量），應用式（1）即可對洪水滲透下海塘的飽和非飽和／非穩(wěn)定滲流場進行求解。其中h～θ關系可由試驗實測確定，而kr（h）～h關系可根據(jù)h～θ關系及VG模型來推導。

1.2 初邊值條件

初始條件可表示為：

式中：t0為初始時刻，h0為位置坐標的給定函數(shù)。

已知壓力水頭邊界條件可表示為：

式中：hc為位置坐標和時間的給定函數(shù)。

已知流量邊界條件可表示為：

式中：q為位置坐標和時間的給定函數(shù)；ni為邊界面法向矢量的第i個分量。

此外，地表入滲邊界條件和出逸邊界條件及它們的其處理方法參見文獻［11］。

2 滲透破壞模擬方法

運用達西滲流理論，采取擴大滲透破壞區(qū)滲透系數(shù)的方法，對滲透破壞區(qū)和未滲透破壞區(qū)統(tǒng)一進行滲流計算。引入1個土體滲透破壞的判別條件臨界滲透坡降值，根據(jù)各單元的計算滲透坡降值與臨界滲透坡降值比較結果，判斷單元是否發(fā)生滲透破壞，對發(fā)生滲透破壞的單元，擴大其滲透系數(shù)值后重新進行滲流場的迭代計算；通過迭代計算，動態(tài)模擬海塘滲透破壞的發(fā)展過程。

2.1 滲透破壞單元的判別及其處理

由于滲透破壞區(qū)附近為滲流流態(tài)的急變區(qū)域，因此為了避免網(wǎng)格尺寸影響滲透破壞區(qū)附近單元的滲透坡降計算值，采用二次單元進行計算域剖分和滲流計算。在計算中選用各土層的允許坡降下限值作為臨界滲透坡降值。滲透破壞單元的判別采用單元形心處計算滲透坡降是否大于臨界滲透坡降來確定，當豎直向或水平向的計算滲透坡降值有1個大于臨界滲透坡降時，即認為單元產(chǎn)生滲透破壞。實際分析計算時，首先判斷單元是否處于滲透破壞區(qū)外邊界，然后根據(jù)計算結果判別單元是否發(fā)生滲透破壞。如有滲透破壞單元，擴大滲透系數(shù)值后再重新進行滲流場計算，否則停止計算。

2.2 滲透破壞區(qū)滲透系數(shù)的取值

滲透破壞區(qū)可分為尖端過渡區(qū)、中間段和出口段3部分，其滲透系數(shù)是沿程變化的。其中中間段的滲透系數(shù)最大且不同位置處大體相等；出口段的含沙量較高，其滲透系數(shù)小于中間段的滲透系數(shù)；尖端過渡區(qū)的滲透系數(shù)最小值與土體相當，最大值與中間段相當。另外，滲透破壞區(qū)內(nèi)滲流流態(tài)復雜、層紊流并存，滲透系數(shù)的確定非常麻煩。機時所限，網(wǎng)格剖分不可能太密，加之尖端過渡區(qū)滲透系數(shù)的劇烈變化也很難真實模擬，因此為了得到相對準確的結果，需要進行特殊處理。本文采取不區(qū)分出口段、中間段和過渡區(qū)，其滲透系數(shù)統(tǒng)一按滲透破壞前的單元滲透系數(shù)擴大100倍來選取。后面算例分析表明這樣處理是比較合理的。

3 實例分析

基于上述飽和非飽和／非穩(wěn)定滲流數(shù)學模型和滲透破壞模擬方法，編制了可考慮土體滲透性隨滲透破壞而變化的飽和非飽和／非穩(wěn)定滲流有限元計算程序，可進行洪水滲透下海塘的飽和非飽和／非穩(wěn)定滲流分析。選擇某海塘典型斷面（該段海塘曾進行了現(xiàn)場監(jiān)測試驗［12］），采用該計算程序?qū)ζ錆B透破壞過程進行了數(shù)值模擬。計算斷面及土層分布具體參見文獻［13］；按海塘斷面各土層滲透性的不同劃分為11個滲透分區(qū)，各滲透分區(qū)的飽和滲透系數(shù)見表1，非飽和水力參數(shù)按van－Genuchten模型擬定，其中的擬合參數(shù)按各滲透分區(qū)的土質(zhì)類比確定，見表1。初始條件參照現(xiàn)場滲透試驗時的實際情況確定，即：初始塘身土體含水量按各土層實測含水量確定，堤基土體取飽和含水量，初始地下水位海塘內(nèi)外側齊平，均取為5.60m；臨水側水頭邊界條件與現(xiàn)場原型滲透試驗完全一致［12］。滲透薄弱土層（主要是含碎石、碎磚瓦的填土）的臨界滲透坡降取0.15，其他土層的臨界滲透坡降統(tǒng)一取0.40。

將現(xiàn)場監(jiān)測試驗測壓管水位實測值與本文有限元計算值的比較列于表2。

表1 塘身、堤基土體飽和滲透系數(shù)及非飽和水力參數(shù)

表2 測壓管水位實測值與有限元計算值 （m）

由表2可知，測壓管水位實測值與本文計算值基本接近，最大差值均在0.29m水頭以下。說明本文所建立的數(shù)學模型和所編制的計算程序是正確的。

不同圍堰水位下海塘滲透坡降及滲流場等勢線分布圖如圖1～3所示。圖1給出了洪水滲透4d后的等勢線及滲透坡降分布情況，由圖可看出，不同土層分界面處等勢線較密，坡腳處滲透坡降大約為0.12，遠未達到土層的臨界滲透坡降。圖2給出了洪水滲透14d后的等勢線及滲透坡降分布情況，坡腳和護塘地處等勢線較密，坡腳處和護塘地的滲透坡降超過了土層的臨界滲透坡降，發(fā)生了滲透破壞現(xiàn)象（滲透破壞區(qū)為圖2所示的兩處陰影區(qū)域）。圖3給出了洪水回落2d后的等勢線及滲透坡降分布情況，不同土層分界面處等勢線較密，各處滲透坡降值大大回落，滲透破壞區(qū)停止發(fā)展。

圖1 洪水滲透4d后（臨水側水位8.0m）的等勢線及滲透坡降分布圖

圖2 洪水滲透14d后（臨水側水位9.38m）的等勢線及滲透坡降分布圖

根據(jù)不同圍堰水位下海塘滲透坡降及滲流場等勢線分布圖（圖1～3）可知，洪水滲透下海塘塘身滲流是一個明顯的飽和非飽和／非穩(wěn)定滲流過程，塘身浸潤線變化經(jīng)歷4個階段。1）浸潤線到達背坡腳以前階段：此階段海塘塘身的飽和區(qū)不斷擴大，最終塘身背坡腳出現(xiàn)滲水（圖1）；2）出逸點沿坡面上升至最高點階段：此階段存在著坡腳滲漏，雖然最高水位已持續(xù)3d，但海塘仍未出現(xiàn)滲透破壞；3）海塘開始滲透破壞階段：此階段塘身開始出現(xiàn)滲透破壞，從坡腳開始，逐漸往海塘內(nèi)發(fā)展，至洪水開始回落時擴展至最大范圍，主要是海塘背坡腳及附近的護塘地（圖2）；4）洪水回落階段：出逸坡降逐漸減小，海塘滲透破壞區(qū)域停止發(fā)展，塘身內(nèi)形成一個拱形浸潤線，浸潤線最高點向海塘內(nèi)側移動，塘身內(nèi)側長時間處于高浸潤線、高飽和度的狀態(tài)（圖3）。上述計算結果與現(xiàn)場監(jiān)測試驗結果基本吻合［12］。

圖3 洪水回落2d后（臨水側水位6.5m）的等勢線及滲透坡降分布圖

4 結 論

通過實際典型海塘的洪水滲透過程的數(shù)值模擬可得到以下結論：

1）實例分析結果表明，本文的滲流分析模型和計算程序是正確的。洪水滲透下海塘塘身滲流是一個飽和

非飽和／非穩(wěn)定滲流過程，塘身浸潤線變化經(jīng)歷4個階段。（1）浸潤線到達背坡腳以前階段：該階段飽和區(qū)不斷擴大，直至背坡腳出現(xiàn)滲水；（2）出逸點沿坡面上升至最高點階段：該階段存在著坡腳滲漏，但海塘尚未出現(xiàn)滲透破壞；（3）海塘開始滲透破壞階段：該階段從坡腳開始滲透破壞，逐漸往海塘內(nèi)發(fā)展，破壞區(qū)集中在海塘背坡腳及附近的護塘地；（4）洪水回落階段：出逸坡降逐漸減小，破壞區(qū)域停止發(fā)展，塘身內(nèi)形成一個拱形浸潤線。

2）由滲透破壞動態(tài)發(fā)展過程可知，海塘滲透破壞是大坡降下洪水長時間滲透所導致的。按現(xiàn)行規(guī)范進行海塘穩(wěn)定滲流計算反映的是洪水滲透的最壞條件，難以模擬達到這種最壞條件所需的時間和整個過程。所以為了合理評價海塘的滲透穩(wěn)定性，需要進行洪水滲透下從非飽和／非穩(wěn)定滲流到飽和／穩(wěn)定滲流的全過程模擬。

［1］Yin J H.FE modelling of seepage in embankment soils with piping zone ［J］. Chinese Journal of Rock Mechanics and Engineering，1998，17（6）：679－686.

［2］周曉杰，丁留謙，姚秋玲，等.懸掛式防滲墻控制堤基滲透變形發(fā)展模型試驗［J］.水力發(fā)電學報，2007，26（2）：54－59.Zhou X J，Ding L Q，Yao Q L，et al.Laboratory model test for evolution of seepage deformation controlled by means of suspended cut－off wall in foundation of dike［J］.Journal of Hydroelectric Engineering，2007，26（2）：54－59.

［3］張家發(fā)，張偉，袁耀宇.防滲墻作用下堤防保護區(qū)地下水動態(tài)數(shù)值模擬分析Ⅱ：二元結構堤基條件下的數(shù)值模擬［J］.長江科學院院報，2007，（1）：27－31.Zhang J F，Zhang W，Yuan Y Y.Numerical simulation on groundwater regime in land protected by dyke with cut－off wall II： numerical simulation for dyke foundation of double layers［J］.Journal of Yangtze River Scientific Research Institute，2007，21（6）：47－50.

［4］丁留謙，吳夢喜，劉昌軍，等.雙層堤基管涌動態(tài)發(fā)展的有限元模擬［J］.水利水電技術，2007，38（2）：36－39.Ding L Q ，Wu M X，Liu C J，et al.Numerical simulation of dynamic development of piping in two－stra tum dike foundations ［J］. Water Resources and Hydropower Engineering，2007，38（2）：36－39.

［5］王國慶，王忠權，謝興華.錢塘江西江塘段海塘滲流分析及防滲方案研究［J］.人民長江，2008，39（14）：89－91.Wang G Q，Wang Z Q，Xie X H.Seepage analysis and impervious reinforcement study of the Xijiangtang bank of Qiantang river［J］.Yangtze River，2008，39（14）：89－91.

［6］Usama E S，Zeghal M.Coupled continuum discrete model for saturated granular soils ［J］.Journal of Engineering Mechanics，2005（4）：413－426.

［7］Kohgo Y，Nakano M，Miyazaki T.Theoretical aspects of constitutive modeling for unsaturated soils［J］.Soil and Foundations，1993，33（4）：49－63

［8］Guarracino L，Quintana F.A constitutive model for water flow in unsaturated fractured rocks ［J］.Hydrological Processes，2009，23：697－701.

［9］Wan C F，Robin F.Investigation of rate of erosion of soils in embankment dams［J］.Journal of Geotechnical and Geoenvironmental Engineering，2004，（4）：373－380.

［10］Fahimifar A，Zareifard M R.A theoretical solution for analysis of tunnels below groundwater considering the hydraulic－mechanical coupling ［J］. Tunnelling and Underground Space Technology，2009，24：634－646.

［11］Hu Y J，Su B Y，Zhan M L.Numerical analysis of saturated－unsaturated seepage flow in fractured rock mass due to surface infiltration ［J］.Journal of Hydrodynamics，2001，13（4）：28－33.

［12］李樹巍，吳雅峰，王擁文.錢塘江南岸西江塘海塘現(xiàn)場滲透試驗觀測成果分析［J］.城市道橋與防洪，2009（5）：93－95.Li S W，Wu Y F，Wang Y W.Analysis on the field permeability test results of the Xijiangtang at the south bank of Qiantang river［J］.Urban Roads Bridges ＆Flood Control，2009（5）：93－95.

［13］賀春雷.海塘滲透變形模擬與影響因素分析［D］.杭州：浙江大學，2010.