土木建筑与环境工程  2013, Vol. 35 Issue (1): 40-45   PDF    
海塘渗透破坏数值模拟及渗透特性
胡云进1, 钟振1, 贺春雷2, 王国庆2    
1. 浙江大学 建筑工程学院, 杭州 310058;
2. 浙江省钱塘江管理局, 杭州 310016
收稿日期：2012-05-20
基金项目：国家自然科学基金资助项目（51279177，50809059）
作者简介：胡云进(1974-), 男, 副教授, 博士, 主要从事渗流水力学和岩体力学研究, (E-mail)huyunjin@zju.edu.cn。
摘要：为研究海塘的洪水渗透规律和渗透破坏机制，提出了可考虑土体渗透性随渗透破坏而变化的饱和非饱和/非稳定渗流分析模型，编制了相应的程序。选取实际海塘进行了洪水渗透过程的数值模拟，并将计算结果与现场监测试验结果进行了对比，表明所建立的分析模型和计算程序是正确的。工程实例分析结果表明，海塘洪水渗透过程是一个涉及土体渗透性变化的饱和非饱和/非稳定渗流过程，塘身浸润线变化经历4个阶段，海塘渗透破坏是大渗流比降下洪水长时间渗透所导致的。因此海塘的洪水渗透破坏涉及破坏区土体渗透性的变化以及饱和非饱和/非稳定渗流，按中国现行规范应用稳定渗流分析方法进行计算是存在问题的。
关键词：海塘    有限元模拟    渗透破坏    渗透特性    
Analysis of Numerical Simulation of Seepage Failure and Seepage Feature of Seawall
Hu Yunjin1, Zhong Zhen1, He Chunlei2, Wang Guoqing2    
1. College of Civil Engineering and Architecture, Zhejiang University, Hangzhou 310058, P. R. China;
2. The Qiantang River Administration of Zhejiang Province, Hangzhou 310016, P. R. China
Received: 2012-05-20
Abstract: In order to analyze seepage feature and seepage failure mechanisms of seawall under flood infiltration, a numerical model for saturated-unsaturated and unsteady seepage flow analysis considering soil permeability variation with seepage failure process was established. And the corresponding finite element computation program was made. Numerical simulation of flood seepage for typical section of real seawall was carried out by use of the program. The computed results were compared with the measured results of field test, which shows good agreement. Based on the numerical simulation results, it can be seen that flood seepage of seawall is a saturated-unsaturated and unsteady seepage flow process and seepage failure of seawall is caused by long time flood infiltration under large seepage gradient. Therefore, upon the flood infiltration, seepage failure of the seawall is a very complicated process, which involves soil permeability variation with seepage failure process and saturated-unsaturated/unsteady seepage flow. It is questionable to use the steady seepage analysis method to simulate seepage failure process of seawall.
Key Words: seawall    simulation by finite element method    seepage failure    seepage feature    

洪水渗透下海塘内的渗流是一种涉及土体渗透性变化的饱和非饱和/非稳定渗流。某些情况下，海塘内的渗流很难达到稳定渗流状态，渗透稳定性就较高。而在另外一些情况下，海塘内的渗流会较快达到稳定渗流状态，并在稳定渗流条件下长期工作，渗透稳定性就较低，易发生渗透破坏。因此，洪水渗透下海塘达到稳定渗流所需时间相差较大的话，采用稳定渗流分析方法得到相同的安全系数，海塘安全度实际上是不同的。而中国现行的《堤防工程设计规范》规定，采用稳定渗流分析方法计算出逸坡降，并以此评价海塘的渗透稳定性，这在很多情况下是不妥当的。为了准确评价海塘的渗透稳定性，需要研究海塘洪水渗透从非饱和/非稳定渗流达到饱和/稳定渗流的整个过程，以模拟海塘渗透破坏的发生和发展过程。而为了进行洪水渗透下海塘内渗流从非饱和/非稳定状态发展到饱和/稳定状态以及海塘渗透破坏动态发展过程的数值模拟，就必须采用饱和非饱和/非稳定渗流分析方法，同时考虑土体渗透性随海塘渗透破坏而变化的问题。从检索文献来看，在现有的海塘(或堤防)渗透破坏数值模拟方法中，大多采用增大渗透破坏区内土体渗透系数的方法来模拟渗透破坏发展过程，但只进行饱和渗流分析，未考虑从非饱和/非稳定渗流到饱和/稳定渗流的阶段，具体参见文献[1-6]。有进行饱和非饱和/非稳定渗流分析的文献，但渗流分析时没有考虑土体渗透性随海塘渗透破坏而变化的问题，具体参见文献[7-10]。因此为了获得更准确的洪水渗透下海塘渗透破坏模拟结果，需要建立可考虑土体渗透性随海塘渗透破坏而变化的饱和非饱和/非稳定渗流分析模型。

本文针对目前数值模拟方法的不足，首先建立可考虑土体渗透性随海塘渗透破坏而变化的饱和非饱和/非稳定渗流分析模型，并编制相应的有限元计算程序；然后选取实际海塘进行洪水渗透过程的数值模拟，并将计算结果与现场监测试验结果进行对比，验证分析模型和计算程序的正确性；最后根据数值模拟结果，总结归纳洪水渗透下海塘的渗透特征及规律。

1 饱和非饱和/非稳定渗流数学模型
1.1 控制方程

海塘为典型的线状建筑物，其渗流计算一般选典型断面按二维问题进行分析，控制方程如下[9]：

$ \begin{array}{l} \sum\limits_{i = 1}^2 {} \sum\limits_{j = 1}^2 {\frac{\partial }{{\partial {x_i}}}} [{k_r}\left( h \right){k_{ij}}\frac{\partial }{{\partial {x_i}}}\left( {h + {x_2}} \right)] - [C\left( h \right) + \\ \beta {S_S}]\frac{{\partial h}}{{\partial t}} -S = 0 \end{array} $ (1)

式中：kr(h)为相对渗透系数；kij为饱和渗透系数；C(h)为容水度；h为压力水头；kr和C均为压力水头的函数；β为选择参数，非饱和区为0，饱和区为1；SS为等效单位贮存量；S为源(汇)项；xi为坐标轴，其中x2为正向向上的铅直轴；t为时间。

有了符合海塘土质特性的h~θ、kr(h)~h及C(h)~h关系后(θ为含水量)，应用式(1)即可对洪水渗透下海塘的饱和非饱和/非稳定渗流场进行求解。其中h~θ关系可由试验实测确定，而kr(h)~h关系可根据h~θ关系及VG模型来推导。

1.2 初边值条件

初始条件可表示为：

$ h\left( {{x_i}, 0} \right) = {h_0}\left( {{x_i}, {t_0}} \right)\;\;\;i = 1, 2 $ (2)

式中：t0为初始时刻，h0为位置坐标的给定函数。

已知压力水头边界条件可表示为：

$ h({x_i}, t) = {h_c}({x_i}, t)\;\;\;\;\;i = 1, 2 $ (3)

式中：hc为位置坐标和时间的给定函数。

已知流量边界条件可表示为：

$ {k_r}\left( h \right)\sum\limits_{i = 1}^2 {} [\sum\limits_{j = 1}^2 {{k_{ij}}} \frac{{\partial h}}{{\partial {x_j}}} + {k_{i2}}]{n_i} = -q\left( {{x_i}, t} \right) $ (4)

式中：q为位置坐标和时间的给定函数；ni为边界面法向矢量的第i个分量。

此外，地表入渗边界条件和出逸边界条件及它们的其处理方法参见文献[11]。

2 渗透破坏模拟方法

运用达西渗流理论，采取扩大渗透破坏区渗透系数的方法，对渗透破坏区和未渗透破坏区统一进行渗流计算。引入1个土体渗透破坏的判别条件临界渗透坡降值，根据各单元的计算渗透坡降值与临界渗透坡降值比较结果，判断单元是否发生渗透破坏，对发生渗透破坏的单元，扩大其渗透系数值后重新进行渗流场的迭代计算；通过迭代计算，动态模拟海塘渗透破坏的发展过程。

2.1 渗透破坏单元的判别及其处理

由于渗透破坏区附近为渗流流态的急变区域，因此为了避免网格尺寸影响渗透破坏区附近单元的渗透坡降计算值，采用二次单元进行计算域剖分和渗流计算。在计算中选用各土层的允许坡降下限值作为临界渗透坡降值。渗透破坏单元的判别采用单元形心处计算渗透坡降是否大于临界渗透坡降来确定，当竖直向或水平向的计算渗透坡降值有1个大于临界渗透坡降时，即认为单元产生渗透破坏。实际分析计算时，首先判断单元是否处于渗透破坏区外边界，然后根据计算结果判别单元是否发生渗透破坏。如有渗透破坏单元，扩大渗透系数值后再重新进行渗流场计算，否则停止计算。

2.2 渗透破坏区渗透系数的取值

渗透破坏区可分为尖端过渡区、中间段和出口段3部分，其渗透系数是沿程变化的。其中中间段的渗透系数最大且不同位置处大体相等；出口段的含沙量较高，其渗透系数小于中间段的渗透系数；尖端过渡区的渗透系数最小值与土体相当，最大值与中间段相当。另外，渗透破坏区内渗流流态复杂、层紊流并存，渗透系数的确定非常麻烦。机时所限，网格剖分不可能太密，加之尖端过渡区渗透系数的剧烈变化也很难真实模拟，因此为了得到相对准确的结果，需要进行特殊处理。本文采取不区分出口段、中间段和过渡区，其渗透系数统一按渗透破坏前的单元渗透系数扩大100倍来选取。后面算例分析表明这样处理是比较合理的。

3 实例分析

基于上述饱和非饱和/非稳定渗流数学模型和渗透破坏模拟方法，编制了可考虑土体渗透性随渗透破坏而变化的饱和非饱和/非稳定渗流有限元计算程序，可进行洪水渗透下海塘的饱和非饱和/非稳定渗流分析。选择某海塘典型断面(该段海塘曾进行了现场监测试验[12])，采用该计算程序对其渗透破坏过程进行了数值模拟。计算断面及土层分布具体参见文献[13]；按海塘断面各土层渗透性的不同划分为11个渗透分区，各渗透分区的饱和渗透系数见表 1，非饱和水力参数按van-Genuchten模型拟定，其中的拟合参数按各渗透分区的土质类比确定，见表 1。初始条件参照现场渗透试验时的实际情况确定，即：初始塘身土体含水量按各土层实测含水量确定，堤基土体取饱和含水量，初始地下水位海塘内外侧齐平，均取为5.60 m；临水侧水头边界条件与现场原型渗透试验完全一致[12]。渗透薄弱土层(主要是含碎石、碎砖瓦的填土)的临界渗透坡降取0.15，其他土层的临界渗透坡降统一取0.40。

表 1 塘身、堤基土体饱和渗透系数及非饱和水力参数

将现场监测试验测压管水位实测值与本文有限元计算值的比较列于表 2。

表 2 测压管水位实测值与有限元计算值

由表 2可知，测压管水位实测值与本文计算值基本接近，最大差值均在0.29 m水头以下。说明本文所建立的数学模型和所编制的计算程序是正确的。

不同围堰水位下海塘渗透坡降及渗流场等势线分布图如图 1~3所示。图 1给出了洪水渗透4 d后的等势线及渗透坡降分布情况，由图可看出，不同土层分界面处等势线较密，坡脚处渗透坡降大约为0.12，远未达到土层的临界渗透坡降。图 2给出了洪水渗透14 d后的等势线及渗透坡降分布情况，坡脚和护塘地处等势线较密，坡脚处和护塘地的渗透坡降超过了土层的临界渗透坡降，发生了渗透破坏现象(渗透破坏区为图 2所示的两处阴影区域)。图 3给出了洪水回落2 d后的等势线及渗透坡降分布情况，不同土层分界面处等势线较密，各处渗透坡降值大大回落，渗透破坏区停止发展。

图 1 洪水渗透4 d后(临水侧水位8.0 m)的等势线及渗透坡降分布图

图 2 洪水渗透14 d后(临水侧水位9.38 m)的等势线及渗透坡降分布图

图 3 洪水回落2 d后(临水侧水位6.5 m)的等势线及渗透坡降分布图

根据不同围堰水位下海塘渗透坡降及渗流场等势线分布图(图 1~3)可知，洪水渗透下海塘塘身渗流是一个明显的饱和非饱和/非稳定渗流过程，塘身浸润线变化经历4个阶段。1)浸润线到达背坡脚以前阶段：此阶段海塘塘身的饱和区不断扩大，最终塘身背坡脚出现渗水(图 1)；2)出逸点沿坡面上升至最高点阶段：此阶段存在着坡脚渗漏，虽然最高水位已持续3 d，但海塘仍未出现渗透破坏；3)海塘开始渗透破坏阶段：此阶段塘身开始出现渗透破坏，从坡脚开始，逐渐往海塘内发展，至洪水开始回落时扩展至最大范围，主要是海塘背坡脚及附近的护塘地(图 2)；4)洪水回落阶段：出逸坡降逐渐减小，海塘渗透破坏区域停止发展，塘身内形成一个拱形浸润线，浸润线最高点向海塘内侧移动，塘身内侧长时间处于高浸润线、高饱和度的状态(图 3)。上述计算结果与现场监测试验结果基本吻合[12]。

4 结论

通过实际典型海塘的洪水渗透过程的数值模拟可得到以下结论：

1) 实例分析结果表明，本文的渗流分析模型和计算程序是正确的。洪水渗透下海塘塘身渗流是一个饱和非饱和/非稳定渗流过程，塘身浸润线变化经历4个阶段。(1)浸润线到达背坡脚以前阶段：该阶段饱和区不断扩大，直至背坡脚出现渗水；(2)出逸点沿坡面上升至最高点阶段：该阶段存在着坡脚渗漏，但海塘尚未出现渗透破坏；(3)海塘开始渗透破坏阶段：该阶段从坡脚开始渗透破坏，逐渐往海塘内发展，破坏区集中在海塘背坡脚及附近的护塘地；(4)洪水回落阶段：出逸坡降逐渐减小，破坏区域停止发展，塘身内形成一个拱形浸润线。

2) 由渗透破坏动态发展过程可知，海塘渗透破坏是大坡降下洪水长时间渗透所导致的。按现行规范进行海塘稳定渗流计算反映的是洪水渗透的最坏条件，难以模拟达到这种最坏条件所需的时间和整个过程。所以为了合理评价海塘的渗透稳定性，需要进行洪水渗透下从非饱和/非稳定渗流到饱和/稳定渗流的全过程模拟。

参考文献
[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] 周晓杰, 丁留谦, 姚秋玲, 等. 悬挂式防渗墙控制堤基渗透变形发展模型试验[J]. 水力发电学报, 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. (in Chinese)
[3] 张家发, 张伟, 袁耀宇. 防渗墙作用下堤防保护区地下水动态数值模拟分析Ⅱ:二元结构堤基条件下的数值模拟[J]. 长江科学院院报, 2007, 21(6): 47–50.
Zhang J F, Zhang W, Yuan Y Y. Numerical simulation on groundwater regime in land protected by dyke with cut-off wall Ⅱ:numerical simulation for dyke foundation of double layers[J]. Journal of Yangtze River Scientific Research Institute, 2007, 21(6): 47–50. (in Chinese)
[4] 丁留谦, 吴梦喜, 刘昌军, 等. 双层堤基管涌动态发展的有限元模拟[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. (in Chinese)
[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. DOI:10.3969/j.issn.1001-4179.2008.14.033 (in Chinese)
[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. DOI:10.3208/sandf1972.33.4_49
[8] Guarracino L, Quintana F. A constitutive model for water flow in unsaturated fractured rocks[J]. Hydrological Processes, 2009, 23: 697–701. DOI:10.1002/hyp.v23:5
[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. DOI:10.1016/j.tust.2009.06.002
[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] 李树巍, 吴雅峰, 王拥文. 钱塘江南岸西江塘海塘现场渗透试验观测成果分析[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. (in Chinese)
[13] 贺春雷. 海塘渗透变形模拟与影响因素分析[D]. 杭州: 浙江大学, 2010.
X