comsol在非饱和土渗流的应用

基于comsol的非饱和土渗流研究/comsol在岩土工程渗流的应用摘要：岩土工程的核心难点即解决地下水问题，一般岩土工程事故都是由于对地下水的影响重视不够而造成的，然而解决这一难点关键在于解决地下水渗流问题。

目前对于非饱和土渗流研究的理论仍相对落后，本文结合非饱和土渗流场基本方程以及由水土特征曲线得到的相关渗流参数（渗透系数，体积含水量），阐明了如何解决渗透模型要求渗流场方程的连续性与现场实测数据的非连续性之间的矛盾，并利用comsol Multiphysics 软件对某工程中非饱和土渗流问题进行了模拟，并验证了Fredlund和xing（1994）土水特征曲线方程的正确性。

这种解决非饱和土渗流问题的思想可供学者参考。

关键词：非饱和土；渗流场；渗流参数；连续性矛盾；Comsol Multiphysics Study on seepage of unsaturated soil seepage based on comsol Abstract:The core difficulty of geotechnical engineering is to solve groundwater problems, the general geotechnical engineering accidents are due to the impact of groundwater caused 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 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 Multiphysics 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; seepage parameters; continuous conflict; Comsol Multiphysics1引言岩土工程设计与施工的难点在于解决地下水问题，一般岩土工程事故都是由于对地下水的影响重视不够而造成的，像2003年7月14日上海轨道交通４号线工程事故；2007年8月17日山东新汶煤矿透水事故；2008年11月15日杭州地铁工地塌陷事故以及2011年1月1日杭州余杭区-工地土方坍塌事故等等都是由于忽视地下水的影响而造成的。

然而解决这一难题的关键在于解决地下水渗流问题。

虽然众多国内外学者对土的渗流问题做了大量的研究，但是目前对于非饱和土渗流研究的理论以及实践应用仍相对落后。

一般来说，解决非饱和土渗流设计的问题以及与其相关的工程实践问题，可以归结于就具体的非饱和土渗流工程概况而建立渗流场基本方程，然后解这一渗流场基本方程，从而得出相关的渗流流线（水位）分布、水流渗流力矢量分布、流速矢量分布和相关的趋势，最后以此来指导实践施工。

然而在求解非饱和土渗流场基本方程时，首要要解决两个重要未知参数，即体积含水量θ和渗透系数k，这两个参数在实际工程中是通过实验得到的，试验得到的是一系列孤立的点，然而这与渗流场基本方程建立于连续性模型相悖，这就给求解渗流场基本方程带来了很大的困难，于是国内外很多学者对此进行了大量的研究。

为了解决非饱和土的体积含水量θ与渗透系数k 测量的耗时以及模型的连续性问题，Gardner ， Brooks&Corey ，V an Genuchten ，McKee&Bumb 和Fredlund&Xing 等人[1]先后通过试验建立了一系列的土水特征方程即体积含水量θ与吸力Ψ的拟合连续性表达式；但是渗透系数k 的连续性仍未解决，后来E.C.Leong 和H.Rahardjo 等人通过大量的试验发现非饱和土的体积含水量θ和渗透系数k 都与吸力Ψ存在密切的联系，并且依赖于土水特征曲线，于是建立了体积含水量θ与渗透系数k 拟合关系。

于此解决了求解非饱和土渗流场基本方程的首要问题。

本文以以上内容为核心思想，利用Comsol Multiphysics 软件，对具体工程中非饱和土渗流问题进行了数值模拟，通过反馈验证了Fredlund 和Xing[4]等人土水特征曲线方程的正确性。

以上解决非饱和土渗流问题的思想可供学者参考。

2模型的建立2.1非饱和土渗流场基本方程一般来说，解决非饱和土渗流设计的问题以及与其相关的工程实践问题，特别是数值分析，都归结于就具体的非饱和土渗流工程概况而建立渗流场基本方程，然后解这一渗流场基本方程，从而得出相关的渗流流线（水位）分布、水流渗流力矢量分布、流速矢量分布和相关的趋势，最后以此来指导实践施工。

二维非饱和土渗流场基本方程为：x y H H k k Q x x y y t θ⎛⎫∂∂∂∂∂⎛⎫++= ⎪ ⎪∂∂∂∂∂⎝⎭⎝⎭ （1）式中：H —总水头；k x —x 方向的渗透系数（变量）；k y —y 方向的渗透系数（变量）；Q —应用边界流；θ—体积含水量；t —时间。

由于孔隙气压力通常为常量u a ，对体积含水量θ的改变没有影响，改变的体积含水量由下式决定：w w m u θ∂=∂ （2）式中：m w —u w -Θ关系曲线的斜率。

H 的表达式为：ww u H y γ=+ （3）式中：u w —孔隙水压力；γw 水的重度；y —海拔。

（3）式可得：()w w u H y γ=- （4）由（1）（2）（4）得：x y w w H H H k k Q m x x y y t γ⎛⎫∂∂∂∂∂⎛⎫++= ⎪ ⎪∂∂∂∂∂⎝⎭⎝⎭ （5）在求解上述非饱和土渗流场基本方程（1）需要考虑边界条件和初始条件，边界条件可以分为三类[3]：第一类边界条件为水头边界条件，即h |Γ1=H 1（x ，y ，t ）；第二类边界条件为流量边界，即22n q hH x y t n K Γ-∂==∂（，，），其中q n 为单位面积边界上穿过的已知流量；第三类边界条件为混合边界条件，即3hh n αβΓ∂+=∂，其中α，β为参数。

初始条件表示为：h （x ，y ，t ）|t =0=h 0（x ，y ，t ）。

2.2土水特征曲线方程在求解非饱和土渗流场基本方程（1）时，要知道体积含水量θ和渗透系数k ，而这两个参数在实际工程中是通过实验得到的，试验得到的是一系列孤立的点，然而这与渗流场基本方程建立于连续性模型相悖。

为了解决非饱和土的体积含水量θ与渗透系数k 测量的耗时以及模型的连续性问题，Gardner ， Brooks&Corey ，V an Genuchten ，McKee&Bumb 和Fredlund&Xing 等人[1]先后通过试验建立了一系列的土水特征方程即体积含水量θ与吸力Ψ的拟合连续性表达式。

但是没有形成统一的体系，后来E.C.Leong 和 H.Rahardjo 通过大量实验对比在[1]中建议我们采用Fredlund 和 Xing （1994）的土水特征方程，其为：ln 11100000ln 1lnr s w c b r e a ψψθθψψ⎡⎤⎛⎫+⎢⎥ ⎪⎝⎭⎢⎥=-⎢⎥⎛⎫⎧⎫⎡⎤⎪⎪⎛⎫+⎢⎥ ⎪+⎢⎥⎨⎬ ⎪⎢⎥⎝⎭⎣⎦⎝⎭⎢⎥⎪⎪⎣⎦⎩⎭或()lns w c b C e a θθψψ=⎧⎫⎡⎤⎪⎪⎛⎫+⎢⎥⎨⎬ ⎪⎝⎭⎢⎥⎪⎪⎣⎦⎩⎭ （6）式中：θw 为体积含水量；Ψ为吸力，与孔隙气压力和孔隙水压力有关；θs 饱和体积含水量；Ψr ，a ，b ，c 为未知量。

通过试验测得非饱和土的体积含水量θ与吸力Ψ关系的数据点组，可以拟合成如（6）式的连续性表达式，从而解决θ的连续性。

需要注意的是E.C.Leong 和 H.Rahardjo 在[1]中提到Fredlund 和 Xing （1994）拟合效果最差，11个点的数据集，而不是6个点的数据集。

2.3渗透参数试验方程上面解决了体积含水量θ非连续性的矛盾，渗透系数k 的非连续仍未解决，后来很多学者人通过大量的试验发现非饱和土的体积含水量θ和渗透系数k 都与吸力Ψ存在密切的联系，并且依赖于土水特征曲线，于是建立了体积含水量θ与渗透系数k 拟合关系。

E.C.Leong 和 H.Rahardjo 在[2]中提到的拟合关系如下：p r k =Θ （7）式中：k r —相对渗透系数，k r =k w /k s ，k w 为渗透系数，k s 为饱和渗透系数；Θ—标准体积含水量，Θ=（θw -θr ）/（θs -θr ），θr 为残余体积含水量；p 为常量（一般在2.5~24.5之间）。

由土水特征曲线得到饱和渗透系数k s 和残余体积含水量θr ，以及由（6）式从而解决了非饱和土渗透系数k 的连续性。

3工程实例模拟3.1工程概况某水库如图1所示，具体物理参数在图中已给出，其中有8m 的透水层；试验测得的体积含水量θ以及渗透系数k 与吸力Ψ关系的数据如表1，水的重度γw =9.8kN/m 3，初始水位高度为，H 0=10m 。