在非饱和土渗流的应用

1、基于comsol的非饱和土渗流研究/comsol在岩土工程渗流的应用摘要：岩土工程的核心难点即解决地下水问题，一般岩土工程事故都是由于对地下水的影响重视不够而造成的，然而解决这一难点关键在于解决地下水渗流问题。目前对于非饱和土渗流研究的理论仍相对落后，本文结合非饱和土渗流场基本方程以及由水土特征曲线得到的相关渗流参数（渗透系数，体积含水量），阐明了如何解决渗透模型要求渗流场方程的连续性与现场实测数据的非连续性之间的矛盾，并利用comsol Multiphysics 软件对某工程中非饱和土渗流问题进行了模拟，并验证了Fredlund和xing（1994）土水特征曲线方程的正确性。这种解决非饱和土

2、渗流问题的思想可供学者参考。关键词：非饱和土；渗流场；渗流参数；连续性矛盾；Comsol MultiphysicsStudy on seepage of unsaturated soil seepage based on comsolAbstract: The core difficulty of geotechnical engineering is to solve groundwater problems, the general geotechnical engineering accidents are due to the impact of groundwater caused

3、by insufficient attention, however, the key to solve this difficult problem is to deal with the groundwater flow. At present, for the study of unsaturated soil seepage theory is still relatively backward, this paper combines basic equation of unsaturated soil seepage with soil-water characteristics

4、curve and obtains the relevant flow parameters (hydraulic conductivity, volumetric water content) from them, and illustrates how to solve the conflict between the seepage field penetration model requiring Equation of continuity and the measured data of non-continuity, and using the software comsol M

5、ultiphysics to simulate unsaturated soil seepage problems in one project and verified the right of Fredlund and xing (1994) soil-water characteristic curve equation. The idea of solving unsaturated soil seepage problems may be referred by similar projects.Key words: unsaturated soil; seepage field;

6、seepage parameters; continuous conflict; Comsol Multiphysics1引言岩土工程设计与施工的难点在于解决地下水问题，一般岩土工程事故都是由于对地下水的影响重视不够而造成的，像2003年7月14日上海轨道交通号线工程事故；2007年8月17日山东新汶煤矿透水事故；2008年11月15日杭州地铁工地塌陷事故以及2011年1月1日杭州余杭区-工地土方坍塌事故等等都是由于忽视地下水的影响而造成的。然而解决这一难题的关键在于解决地下水渗流问题。虽然众多国内外学者对土的渗流问题做了大量的研究，但是目前对于非饱和土渗流研究的理论以及实践应用仍相对落后。一

7、般来说，解决非饱和土渗流设计的问题以及与其相关的工程实践问题，可以归结于就具体的非饱和土渗流工程概况而建立渗流场基本方程，然后解这一渗流场基本方程，从而得出相关的渗流流线（水位）分布、水流渗流力矢量分布、流速矢量分布和相关的趋势，最后以此来指导实践施工。然而在求解非饱和土渗流场基本方程时，首要要解决两个重要未知参数，即体积含水量和渗透系数k，这两个参数在实际工程中是通过实验得到的，试验得到的是一系列孤立的点，然而这与渗流场基本方程建立于连续性模型相悖，这就给求解渗流场基本方程带来了很大的困难，于是国内外很多学者对此进行了大量的研究。为了解决非饱和土的体积含水量与渗透系数k测量的耗时以及模型的连

8、续性问题，Gardner， Brooks&Corey，Van Genuchten，McKee&Bumb和Fredlund&Xing等人1先后通过试验建立了一系列的土水特征方程即体积含水量与吸力的拟合连续性表达式；但是渗透系数k的连续性仍未解决，后来和H.Rahardjo等人通过大量的试验发现非饱和土的体积含水量和渗透系数k都与吸力存在密切的联系，并且依赖于土水特征曲线，于是建立了体积含水量与渗透系数k拟合关系。于此解决了求解非饱和土渗流场基本方程的首要问题。本文以以上内容为核心思想，利用Comsol Multiphysics软件，对具体工程中非饱和土渗流问题进行了数值模

9、拟，通过反馈验证了Fredlund和Xing4等人土水特征曲线方程的正确性。以上解决非饱和土渗流问题的思想可供学者参考。2模型的建立2.1非饱和土渗流场基本方程一般来说，解决非饱和土渗流设计的问题以及与其相关的工程实践问题，特别是数值分析，都归结于就具体的非饱和土渗流工程概况而建立渗流场基本方程，然后解这一渗流场基本方程，从而得出相关的渗流流线（水位）分布、水流渗流力矢量分布、流速矢量分布和相关的趋势，最后以此来指导实践施工。二维非饱和土渗流场基本方程为： （1）式中：H总水头；kxx方向的渗透系数（变量）；kyy方向的渗透系数（变量）；Q应用边界流；体积含水量；t时间。由于孔隙气压力通常为常

10、量ua，对体积含水量的改变没有影响，改变的体积含水量由下式决定： （2）式中：mwuw-关系曲线的斜率。H的表达式为： （3）式中：uw孔隙水压力；w水的重度；y海拔。（3）式可得： （4）由（1）（2）（4）得： （5）在求解上述非饱和土渗流场基本方程（1）需要考虑边界条件和初始条件，边界条件可以分为三类3：第一类边界条件为水头边界条件，即h|1=H1（x，y，t）；第二类边界条件为流量边界，即，其中qn为单位面积边界上穿过的已知流量；第三类边界条件为混合边界条件，即，其中，为参数。初始条件表示为：h（x，y，t）|t=0=h0（x，y，t）。2.2土水特征曲线方程在求解非饱和土渗流场基本方

11、程（1）时，要知道体积含水量和渗透系数k，而这两个参数在实际工程中是通过实验得到的，试验得到的是一系列孤立的点，然而这与渗流场基本方程建立于连续性模型相悖。为了解决非饱和土的体积含水量与渗透系数k测量的耗时以及模型的连续性问题，Gardner， Brooks&Corey，Van Genuchten，McKee&Bumb和Fredlund&Xing等人1先后通过试验建立了一系列的土水特征方程即体积含水量与吸力的拟合连续性表达式。但是没有形成统一的体系，后来E.C.Leong 和 H.Rahardjo通过大量实验对比在1中建议我们采用Fredlund 和 Xing（1994

12、）的土水特征方程，其为：或 （6）式中：w为体积含水量；为吸力，与孔隙气压力和孔隙水压力有关；s饱和体积含水量；r，a，b，c为未知量。通过试验测得非饱和土的体积含水量与吸力关系的数据点组，可以拟合成如（6）式的连续性表达式，从而解决的连续性。需要注意的是E.C.Leong 和 H.Rahardjo在1中提到Fredlund 和 Xing（1994）拟合效果最差，11个点的数据集，而不是6个点的数据集。2.3渗透参数试验方程上面解决了体积含水量非连续性的矛盾，渗透系数k的非连续仍未解决，后来很多学者人通过大量的试验发现非饱和土的体积含水量和渗透系数k都与吸力存在密切的联系，并且依赖于土水特征曲

13、线，于是建立了体积含水量与渗透系数k拟合关系。E.C.Leong 和 H.Rahardjo在2中提到的拟合关系如下： （7）式中：kr相对渗透系数，kr=kw/ks，kw为渗透系数，ks为饱和渗透系数；标准体积含水量，=（w-r）/（s-r），r为残余体积含水量；p为常量（一般在2.524.5之间）。由土水特征曲线得到饱和渗透系数ks和残余体积含水量r，以及由（6）式从而解决了非饱和土渗透系数k的连续性。3工程实例模拟3.1工程概况某水库如图1所示，具体物理参数在图中已给出，其中有8m的透水层；试验测得的体积含水量以及渗透系数k与吸力关系的数据如表1，水的重度w=9.8kN/m3，初始水位高度

14、为，H0=10m。44m8m12m4m10mABCDFE图1 某水库的物理参数Fig. 1 Physical parameters of a reservoir表1 土水特征参数Table 1 Soil-water characteristic parameters 实验组号渗透参数123456吸力（kPa）020406080100体积含水量0.38250.35000.20030.09960.07150.0627渗透系数k（m/s）0.00790.00320.000520×10-61.3×10-61.0×10-63.2模拟参数的处理（1）体积含水量根据表1中测得的

15、数据绘制出吸力与体积含水量（-关系图）如图2所示。图2 实测和关系图Fig. 2 Measured and diagramsFredlund和xing（1994）提到（6）式中的系数C（）系数接近于1，简化计算，取C（）=1，饱和体积含水量s=0.383，根据图2，拟合的含水量连续性方程（6）式，可得（6）式中的参数为a=30.72；b=3.851；c=1.257（一般可以用matlab或origin软件来拟合，本文通过origin拟合）。（2）渗透参数根据表1和图2，（7）式中饱和渗透系数ks=0.008m/s；残余体积含水量r=0.062；通过拟合可得到p=4.52再由（6）式就可以得到了

16、渗透系数k的连续性方程。本工程虽然已测出了渗透系数k和吸力的关系，仅供验证的目的。（3）mw处理将已求出未知参数的含水量连续性方程（6）式中的吸力换成-uw，并对uw求一阶导数即得到mw。3.3 comsol 的模拟在comsol中利用非饱和土渗流场基本方程（5）式，输入已求出的连续性表达式，以及边界条件（图1中，BC为第一类边界条件；CD为第二类边界条件，DE为第三类边界条件；在comsol中需要转化相应的Dirichlet boundary condition和Neumann boundary condition或这两类边界条件的组合），本例的模拟过程和模拟结果如图3图11图3 水库划分网

17、格Fig. 3 Reservoir mesh图4 水库初始应力分布图Fig. 4 Initial stress distribution of the reservoir图5 一年后水库应力分布图Fig. 5 Reservoir stress distribution after one year图6 初始浸润线位置Fig. 6 The initial location of saturation line图7 一年后浸润线位置Fig. 7 Location of saturation line after one year图8 初始流速矢量Fig. 8 Initial velocity ve

18、ctor图9 一年后流速矢量Fig. 9 Velocity vector after one year图10 初始流线图Fig. 10 Initial flow chart图11 一年后流线图Fig. 11 Flow chart after one year4结论本文提供了一种解决非饱和土渗流问题的方法，即结合非饱和土渗流场基本方程以及由水土特征曲线得到的相关渗流参数（渗透系数，体积含水量），以此解决渗透模型要求渗流场方程的连续性与现场实测数据的非连续性之间的矛盾，并利用数值软件对具体工程中非饱和土渗流问题进行模拟，从而得出相关的渗流流线（水位）分布、水流渗流力矢量分布、流速矢量分布和相关的趋势，最后以此来指导施工。参考文献1 E.C.Leong and