comsol在非饱和土渗流的应用

基于comsol的非饱和土渗流研究/comsol在岩土工程渗流的应用摘要：岩土工程的核心难点即解决地下水问题，一般岩土工程事故都是由于对地下水的影响重视不够而造成的，然而解决这一难点关键在于解决地下水渗流问题。目前对于非饱和土渗流研究的理论仍相对落后，本文结合非饱和土渗流场基本方程以及由水土特征曲线得到的相关渗流参数（渗透系数，体积含水量），阐明了如何解决渗透模型要求渗流场方程的连续性与现场实测数据的非连续性之间的矛盾，并利用comsolMultiphysics软件对某工程中非饱和土渗流问题进行了模拟，并验证了Fredlund和xing（1994）土水特征曲线方程的正确性。这种解决非饱和土渗流问题的思想可供学者参考。关键词：非饱和土；渗流场；渗流参数；连续性矛盾；ComsolMultiphysicsStudyonseepageofunsaturatedsoilseepagebasedoncomsolAbstract:Thecoredifficultyofgeotechnicalengineeringistosolvegroundwaterproblems,thegeneralgeotechnicalengineeringaccidentsareduetotheimpactofgroundwatercausedbyinsufficientattention,however,thekeytosolvethisdifficultproblemistodealwiththegroundwaterflow.Atpresent,forthestudyofunsaturatedsoilseepagetheoryisstillrelativelybackward,thispapercombinesbasicequationofunsaturatedsoilseepagewithsoil-watercharacteristicscurveandobtainstherelevantflowparameters(hydraulicconductivity,volumetricwatercontent)fromthem,andillustrateshowtosolvetheconflictbetweentheseepagefieldpenetrationmodelrequiringEquationofcontinuityandthemeasureddataofnon-continuity,andusingthesoftwarecomsolMultiphysicstosimulateunsaturatedsoilseepageproblemsinoneprojectandverifiedtherightofFredlundandxing(1994)soil-watercharacteristiccurveequation.Theideaofsolvingunsaturatedsoilseepageproblemsmaybereferredbysimilarprojects.Keywords:unsaturatedsoil;seepagefield;seepageparameters;continuousconflict;ComsolMultiphysics1引言岩土工程设计与施工的难点在于解决地下水问题，一般岩土工程事故都是由于对地下水的影响重视不够而造成的，像2003年7月14日上海轨道交通４号线工程事故；2007年8月17日山东新汶煤矿透水事故；2008年11月15日杭州地铁工地塌陷事故以及2011年1月1日杭州余杭区-工地土方坍塌事故等等都是由于忽视地下水的影响而造成的。然而解决这一难题的关键在于解决地下水渗流问题。虽然众多国内外学者对土的渗流问题做了大量的研究，但是目前对于非饱和土渗流研究的理论以及实践应用仍相对落后。一般来说，解决非饱和土渗流设计的问题以及与其相关的工程实践问题，可以归结于就具体的非饱和土渗流工程概况而建立渗流场基本方程，然后解这一渗流场基本方程，从而得出相关的渗流流线（水位）分布、水流渗流力矢量分布、流速矢量分布和相关的趋势，最后以此来指导实践施工。然而在求解非饱和土渗流场基本方程时，首要要解决两个重要未知参数，即体积含水量θ和渗透系数k，这两个参数在实际工程中是通过实验得到的，试验得到的是一系列孤立的点，然而这与渗流场基本方程建立于连续性模型相悖，这就给求解渗流场基本方程带来了很大的困难，于是国内外很多学者对此进行了大量的研究。为了解决非饱和土的体积含水量θ与渗透系数k测量的耗时以及模型的连续性问题，Gardner，Brooks&Corey，VanGenuchten，McKee&Bumb和Fredlund&Xing等人[1]先后通过试验建立了一系列的土水特征方程即体积含水量θ与吸力Ψ的拟合连续性表达式；但是渗透系数k的连续性仍未解决，后来E.C.Leong和H.Rahardjo等人通过大量的试验发现非饱和土的体积含水量θ和渗透系数k都与吸力Ψ存在密切的联系，并且依赖于土水特征曲线，于是建立了体积含水量θ与渗透系数k拟合关系。于此解决了求解非饱和土渗流场基本方程的首要问题。本文以以上内容为核心思想，利用ComsolMultiphysics软件，对具体工程中非饱和土渗流问题进行了数值模拟，通过反馈验证了Fredlund和Xing[4]等人土水特征曲线方程的正确性。以上解决非饱和土渗流问题的思想可供学者参考。2模型的建立2.1非饱和土渗流场基本方程一般来说，解决非饱和土渗流设计的问题以及与其相关的工程实践问题，特别是数值分析，都归结于就具体的非饱和土渗流工程概况而建立渗流场基本方程，然后解这一渗流场基本方程，从而得出相关的渗流流线（水位）分布、水流渗流力矢量分布、流速矢量分布和相关的趋势，最后以此来指导实践施工。二维非饱和土渗流场基本方程为：xyHHkkQxxyyt（1）式中：H—总水头；kx—x方向的渗透系数（变量）；ky—y方向的渗透系数（变量）；Q—应用边界流；θ—体积含水量；t—时间。由于孔隙气压力通常为常量ua，对体积含水量θ的改变没有影响，改变的体积含水量由下式决定：wwmu（2）式中：mw—uw-Θ关系曲线的斜率。H的表达式为：wwuHy（3）式中：uw—孔隙水压力；γw水的重度；y—海拔。（3）式可得：wwuHy（4）由（1）（2）（4）得：xywwHHHkkQmxxyyt（5）在求解上述非饱和土渗流场基本方程（1）需要考虑边界条件和初始条件，边界条件可以分为三类[3]：第一类边界条件为水头边界条件，即h|Γ1=H1（x，y，t）；第二类边界条件为流量边界，即22nqhHxytnK（，，），其中qn为单位面积边界上穿过的已知流量；第三类边界条件为混合边界条件，即3hhn，其中α，β为参数。初始条件表示为：h（x，y，t）|t=0=h0（x，y，t）。2.2土水特征曲线方程在求解非饱和土渗流场基本方程（1）时，要知道体积含水量θ和渗透系数k，而这两个参数在实际工程中是通过实验得到的，试验得到的是一系列孤立的点，然而这与渗流场基本方程建立于连续性模型相悖。为了解决非饱和土的体积含水量θ与渗透系数k测量的耗时以及模型的连续性问题，Gardner，Brooks&Corey，VanGenuchten，McKee&Bumb和Fredlund&Xing等人[1]先后通过试验建立了一系列的土水特征方程即体积含水量θ与吸力Ψ的拟合连续性表达式。但是没有形成统一的体系，后来E.C.Leong和H.Rahardjo通过大量实验对比在[1]中建议我们采用Fredlund和Xing（1994）的土水特征方程，其为：ln11100000ln1lnrswcbrea或lnswcbCea（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]中提到的拟合关系如下：prk（7）式中：kr—相对渗透系数，kr=kw/ks，kw为渗透系数，ks为饱和渗透系数；Θ—标准体积含水量，Θ=（θw-θr）/（θs-θr），θr为残余体积含水量；p为常量（一般在2.5~24.5之间）。由土水特征曲线得到饱和渗透系数ks和残余体积含水量θr，以及由（6）式从而解决了非饱和土渗透系数k的连续性。3工程实例模拟3.1工程概况某水库如图1所示，具体物理参数在图中已给出，其中有8m的透水层；试验测得的体积含水量θ以及渗透系数k与吸力Ψ关系的数据如表1，水的重度γw=9.8kN/m3，初始水位高度为，H0=10m。图1某水库的物理参数Fig.1Physicalparametersofareservoir表1土水特征参数Table1Soil-watercharacteristicparameters实验组号渗透参数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中测得的数据绘制出吸力Ψ与体积含水量θ（Ψ-θ关系图）如图2所示。44m8m12m4m10mABCDFE图2实测Ψ和θ关系图Fig.2MeasuredΨ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）式就可以得到了渗透系数k的连续性方程。本工程虽然已测出了渗透系数k和吸力Ψ的关系，仅供验证的目的。（3）mw处理将已求出未知参数的含水量连续性方程（6）式中的吸力Ψ换成-uw，并对uw求一阶导数即得到mw。3.3comsol的模拟在comsol中利用非饱和土渗流场基本方程（5）式，输入已求出的连续性表达式，以及边界条件（图1中，BC为第一类边界条件；CD为第二类边界条件，DE为第三类边界条件；在comsol中需要转化相应的Dirichletboundarycondition和Neumannboundarycondition或这两类边界条件的组合），本例的模拟过程和模拟结果如图3~图11图3水库划分网格Fig.3Reservoirmesh图4水库初始应力分布图Fig.4Initialstressdistributionofthereservoir图5一年后水库应力分布图Fig.5Rese