非饱和土坡中平均有效饱和度分布模拟

许长辉,孙树林,李 方,2,张德恒,3

(1.河海大学 地球科学与工程学院,江苏 南京 211100；2.南京工程高等职业学校,江苏 南京 210035；3.南京工程学院 建筑学院,江苏 南京 211167)

非饱和土坡中平均有效饱和度分布模拟

许长辉1,孙树林1,李 方1,2,张德恒1,3

(1.河海大学 地球科学与工程学院,江苏 南京 211100；2.南京工程高等职业学校,江苏 南京 210035；3.南京工程学院 建筑学院,江苏 南京 211167)

由于直接测定土水特征曲线存在成本高、繁杂、费时等缺点，采用经验公式法预测土水特征曲线越来越受到重视。在非饱和多孔介质中，流体的运动特征主要表现为流体进入和填充孔隙过程中,水的传输和存储量的变化。可以运用COMSOLMultiphysics软件中的Richard方程接口，解决二维非饱和流问题。用该软件模拟非饱和土中的平均有效饱和度分布，预测土坡内和传感器周围的平均有效饱和度，并与实际情况进行对比分析，其结果具有一定的工程实践研究意义。

土-水特征曲线;Richrad方程;平均有效饱和度

影响土质边坡稳定的因素很多, 其中水是诱发滑坡最常见的外因。地球表面分布的土由于气候原因, 大部分是非饱和土, 非饱和土与饱和土最主要的区别在于土中存在负孔隙水压力, 负的孔隙水压力在土体中会产生基质吸力。基质吸力的存在加大了土颗粒间的压力, 使土的抗剪强度增大, 这对保持土坡的稳定起着重要的作用[1]。然而在非饱和土中的基质吸力很容易受外界环境变化的影响, 其中影响最大的是含水量的变化。

非饱和土中的基质吸力或压力水头和含水率之间有一定的对应关系，这种关系称为土水特征曲线[2]。由于直接测定土水特征曲线存在成本高、繁杂、费时等缺点，采用经验公式法预测土水特征曲线越来越受到重视。在土水分特征曲线经验模型中，Brooks-Corey(BC)(BrooksandCorey，1964)模型[3]和vanGenuchten(VG)(vanGenuchten，1980)模型[4]是广泛应用的两种。在描述非饱和多孔介质中流体的运动特征时，一个主要的问题是确定在流体进入和填充孔隙过程中，水的传输和存储量的变化，而这些数据的获得往往是很困难的[5]。此外，在土的吸水和脱水的过程中，往往出现明显的非线性。Richard方程接口提供了一个界面，该界面可自动给出在渗流过程中基于VG或BC方程的材料参数和土水特征曲线[6]。本文运用COMSOLMultiphysics软件模拟非饱和土中的平均饱和度分布，来预测土体内和传感器周围的平均饱和度。

1 土水特征曲线模型简介

土-水特征曲线描述了土体吸力与含水量之间的关系，它主要受土的质地、结构、干容重、温度等影响。此外, 土的膨胀收缩、吸附性离子的种类和数量等因素也影响土水分特征曲线。由于影响因素较多, 且关系复杂，目前尚不能从理论上推求土的基质势与含水量之间的关系, 一般常用经验公式或简单模型表示[7-10]。

BrooksandCorey模型可表示为：

(1)

VanGenuchten模型可表示为：

(2)

式中：θ为体积含水量(cm3/cm3)；θr为滞留含水量(cm3/cm3)；θs为饱和含水量(cm3/cm3)，h为土的吸力(cm)，hd为土的进气吸力(cm)；λ、a、m、n为拟合参数，m=1-1/n。由式(1)和式(2)可知,当土处于饱和状态时,式(1)中土的吸力h等于进气吸力,式(2)中土的吸力h等于零。从这种意义上讲,Brooks-Corey模型更符合脱湿曲线,而vanGenuchten模型更符合吸湿曲线。但由于两种模型在实际应用中各具特点,如Brooks-Corey模型形式简单,便于推求描述土水运动模型和确定土的水分运动参数的简单方法,而vanGenuchten模型适用土的质地范围比较宽,同时可以使饱和土的吸力为0,符合吸湿过程中土的吸力变化特点[11]。因此两模型在实际应用中没有严格区别脱湿和吸湿过程。同时Brooks-Corey模型与vanGenuchten模型也可以进行简单转化,如果省略vanGenuchten模型等式右边分母项中1,式(2)就可变为式(1)。

2 分析模型的定义

COMSOLMultiphysics是一个专业的有限元数值分析软件，是对基于偏微分方程的多物理场模型进行建模和仿真计算的交互式环境开发系统。亚利桑那大学水文与水资源学院的AndrewHinnell、AlexFurman和TyFerre三位学者[12-13]是最早利用COMSOLMultiphysics软件中的PDE接口，成功模拟灌溉过程中土中饱和度的分布和变化过程的学者。在COMSOLMultiphysics软件中创建一个二维非饱和流物理场，选择Richard方程接口，将两个土样参数分别输入模型中，开始模拟便可得到在相同条件下，两个不同土坡的有效饱和度、压力水头和流速的分布。通过软件模拟可以分别得到每个时间点的压力、流速和有效饱和度的分布情况图，这里我们选择15分钟时来分析。

问题设置如下，两个不同土样组成的均质非饱和土坡，截面积均为直角梯形，底边边长分别为1m和2m，高为2m。土样1的水力学特性符合vanGenuchten公式，土样2的材料特性满足Brooks-Corey公式。每个土坡内有三个不透水棒，半径均为0.1m。杆之间的间距为0.5m，杆的轴线平行于土样的中心线，如图1所示。不透水棒设置完毕后，压力水头仍是均匀的，水在压力作用下开始向下渗流。由于所有的垂直切片是相同的，所以可选取一个横截面，模拟15min内流场的变化情况。

通过调研，很多工科高职高专学院没有开设机械检验检测技术专业，没有专任测量实训实验指导教师，机械产品测量实训教学一般是由《互换性与测量技术》的任课教师担任。其教学过程是实验教师根据某一实验项目，首先介绍测量方法及原理、测量任务及要求，其次面对学生做演示，最后根据测量设备的数量分组，学生操作及完成实验报告。其教学组织模式如图1所示。

3 控制方程与模型数据

从图2和图3中可以看出，棒周围的流场发生了明显的变化，而其他部分变化基本上是均匀的。土坡1和土坡2的棒周围压力在15分钟后的分布有明显的不同。土坡1中棒上部压力较大，但土坡2中棒周围的压力更为集中，表明土坡2的渗透能力更差一些。

(3)

方程的第一项说明在饱和-非饱和状态下，储水量会随时间变化而变化。当孔隙完全被水充满时，就会对流体和空隙空间产生轻微压力。容积比湿度可以表示为

式中，压力水头Hp(m)为因变量，C表示容积比湿度(m-1)，Se是有效饱和度，S是储水系数(m-1)，t为时间，K为渗透系数(m/s)，D表示垂直方向的坐标(m)。方程式不显含水的体积分量θ，但它是一个依赖于Hp的数值。因为C、Se和K伴随Hp和θ的变化而变化，因此该过程是一个非线性的过程。

图4显示了棒-土边界处，在0、60、300、600、900s时，有效饱和度的分布。角坐标(-180°,0°)和(0°,180°)分别对应于边界的上、下两部分。在该图中，实线和虚线分别表示土坡1和土坡2的分析结果。有效饱和度的图形初始值为一条水平线。可见，棒上侧的土体比棒下侧的土体要湿润一些。随着时间的推移上侧的有效饱和度逐渐增大，同时下侧部分区域的有效饱和度也在增大。因此，该棒周围边界处最终会完全湿润。

(4)

它描述了含水量体积随水压力的变化关系。储水系数表示饱和状态下，随孔隙的压缩和膨胀，储水量的变化。为了模拟储水系数，本比例使用了比存储选项

S=ρfg(xp+θxf)

(5)

式中，ρf表示流体密度(kg/m3)，g是重力加速度，χp和χf表示固体颗粒和流体的压缩率(m·s2/kg)。VG和BC公式可以描述C、Se、K、θ随Hp变化而变化的关系，前提条件是θs和θr，以及其他定义岩土介质性质的常数α、n、m和l已知。在VG方程中假定当流体的压力等与大气压力时，土为饱和状态，为了得到该问题的解答，必须指定初始条件和边界条件。最初，试样有相同的压力水头Hp0。边界条件如表1所示。其中n是垂直于边界的单位向量。

该公司最大的4.0智库之一是其去年10月推出的SOPHIA服务平台。此平台的基本功能是发出实时信息和数据，优化机器和系统的性能和生产率，利用收集到的信息对生产过程和机器性能进行详细分析，发现失常情况，在维修、订购配件、避免故障发生等方面给客户以协助。Biesse 2018—2020计划的目标之一是进一步改进其服务，从“快速排除故障”到“有预见性的主动服务”。Biesse已具有创造数字化、自动化、具有紧密联系的工厂的能力。

我把碗搁在床头柜上，在二丫的后背塞了个枕头。去拉二丫的手，我像烫着了一样又缩回来——这手太像莲米的手了，也是一把干柴！我说：“二丫，细婶儿给你炖了鸡汤，吃口呵。”

模型其他参数如表2所示。

表1 边界条件设定

表2VG模型和BC模型参数

4 土坡内平均饱和度分布的数值模拟

设想在两个非饱和土坡中，土体均匀但彼此类型不同，土中各插入三个不透水传感器。要获得传感器周围的水量重分布情况和整个区域的饱和度分布，并判断与实际情况的差异。

因为土中的气压等于大气压力，所以只使用Richard方程就可以描述土中的饱和-非饱和流问题。模型的控制方程为

第三，鸦片贸易合法化，并明确规定了鸦片进口关税税率。《中英通商章程善后条约》第五款规定：“向来洋药、铜钱、米谷、豆石、硝磺、白铅等物，例皆不准通商，现定稍宽其禁，听商遵行纳税贸易。洋药准其进口，议定每百觔纳税银叁拾两”“洋药”即鸦片，并规定鸦片进口实行从量税，每百斤纳税三十两白银。

图5为对棒边界上土的平均饱和度与整个土试样的平均饱和度进行对比的结果图。图中散点表示棒周围不同位置处的有效饱和度。其中，圆圈散点区域对应于土坡1，方框散点区域对应于土坡2。图中上半部分的深色实线，表示土坡1中棒周围的平均有效饱和度；而下半部分的深色实线，则表示土坡2中棒周围的平均有效饱和度。对应的两条浅色水平线则分别表示土坡的平均饱和度。由此可见，棒周围图的平均有效饱和度是逐渐增大的，而土坡的平均饱和度却保持不变。对于土坡1，不透水棒对饱和度的影响更为明显。这是因为土坡1的渗透系数更高，在短时间内含水量变化更快，因而其饱和度变化比土坡2更为明显。

11月22日消息，财政部发文，自2018年11月7日起至2021年11月6日止，对境外机构投资境内债券市场取得的债券利息收入暂免征收企业所得税和增值税。上述暂免征收企业所得税的范围不包括境外机构在境内设立的机构、场所取得的与该机构、场所有实际联系的债券利息。

受我国大气污染形势依然严峻及碳排放总量的限制，燃煤电厂发展空间日益受限，清洁能源逐步成为未来能源发展新的增长点。作为清洁低碳化石能源，天然气发电的环境效益及其在能源发展格局中的重要作用越来越受到关注。根据规划，2020年天然气在我国一次能源消费中的比重要达到10%以上，利用量将达到3600亿立方米，成为未来我国能源革命的重要引擎。随着天然气占中国一次能源消费比重不断提高，以及发展天然气发电的重要意义日益凸显，未来我国天然气发电前景广阔。

5 结论

1)通过本文的建模分析可以发现，在相同条件下，不同的非饱和土中饱和度的变化与其渗透系数是有关的。

2)渗透系数大的土中，有效饱和度随时间的变化比渗透系数小的土要大。

3)本文模型还发现在模拟过程中，土坡的平均有效饱和度一直保持不变，而棒周围的平均有效饱和度是随时间逐渐增大的。这是因为在棒的周围发生了流场的变化，棒周围的压力场也在随时间变化，导致棒周围的平均有效饱和度也随时间增大。

[1]李荣建,于玉贞.非饱和土边坡稳定分析方法探讨[J].岩土力学，2007，28(10):2061-2064.

[2]陈 辉,韦昌富.非饱和土土-水特征曲线预估方法研究[J].岩土力学,2013,34(1):128-132.

[3]BROOKSRH,COREYAT.Hydraulicpropertiesofporousmedia[J].ColoradoStatesUniversityHydrol.Paper，1964(3)：1-27.

[4]VANGENUCHTENR.Predictingthehydraulicconductivityofunsaturatedsoils[J].SoilSci.Soc.Am.J.，1980(44): 892- 898.

[5]武海霞,张 铮,王洪义.东北黑土区非饱和土壤水分运动数值模拟研究[J].河北工程大学学报:自然科学版,2010,27(2):60-62.

[6]MUKHLISINM,BAIDILLAHMR.EffectofsoilhydraulicPropertiesmodelonslopestabilityanalysisbasedonstrengthreductionmethod.[J].JournaloftheGeologicalSocietyofIndia，2014，83(5):586-594.

[7]范严伟，邓 燕，王波雷.土壤水分特征曲线VG模型参数求解对比研究[J].人民黄河，2008，30(5)：49-50.

[8]彭建平，邵爱军.基于Matlab方法确定VG模型参数[J].水文地质工程地质，2006，6(7):25-28.

[9]伊盼盼，牛圣宽，柳燕子.非饱和土的土水特征曲线测试方法研究[J].中国农村水利水电，2016(1)：125-128.

[10]李 峰，缴锡云.田间土壤水分特征曲线参数反演[J].河海大学学报:自然科学版，2009，37(4):374-377.

[11]刘小文，叶云雪.不同和影响因素下非饱和红土土-水特征曲线的实验研究[J].水文地质工程地质，2015,42(2)：97-104.

[12]TH,M,VANGENUCHTEN.Aclosed-formequationforpredictingthehydraulicofconductivityofunsaturatedsoils.[J].SoilSci.Soc.Am.J.，1980(44)： 892-898.

[13]BROOKSR.H,COREYA.T.PROPERRropertiesofporousmediaaffectingfluidflow[J].J.Irrig.DrainageDiv.,ASCEProc,1966，72(2):61-88.

(责任编辑 李军)

Simulation of average effective saturation distribution in unsaturated soil slope

XU Changhui1, SUN Shulin1, LI Fang1，2, ZHANG Deheng1，3

(1.College of Earth Sciences and Engineering, Hohai University, Jiangsu Nanjing 211100, China;2. Jiangsu Nanjng Engineering Vocational College, Jiangsu Nanjing 210035, China;3.Department of Architecture Civil Engineering, Nanjing Institute of Technology, Jiangsu Nanjing 211167,China)

Thecommanmethodstodirectlydeterminatethesoilwatercharacteristiccurvehavetheshortcomingsofhighcost,complicatedandtime-consuming.Sotheempiricalformulamethodtopredictsoilwatercharacteristiccurveismoreandmoreaccounted.Themainmovementcharacteristicsof?fluidarethechangesofwatertransportandstorageintheprocessofenteringandfillingtheporeinunsaturatedporousmedia.ThisproblemcanbesoledbyusingRichardequationinCOMSOLMultiphysics.Thedistributionofunsaturatedsoilaverageeffectivesaturationwassimulatedandtheaverageeffectivesaturationinthesoilandaroundsensorwasforecasted.Thesimulationiscomparedwiththeactualsituationandtheresultshavepracticalsignificance.

soil-watercharacteristiccurve;Richradequation;averageeffectivesaturation

1673-9469(2016)04-0017-05doi:10.3969/j.issn.1673-9469.2016.04.005

2016-05-24

河海大学国家重点实验室开放研究基金项目(2005408911)；留学回国人员科研基金资助项目(20071108)

许长辉(1990-)，男，福建莆田人，硕士，研究方向为环境工程地质。

TU

A