基于van Genuchten-Mualem模型的饱和-非饱和介质流动随机数值分析
不同区域土壤水下渗在四水转化中作用研究
不同区域土壤水下渗在四水转化中作用研究[摘要] 根据水循环机理及水平衡机理，现有水存在形式主要为大气水、地表水、土壤水、地下水，四种水存在形式的转化具有复杂的联系。
其中土壤水在转化过程中占据中心位置。
土壤水是指地面以下至潜水面以上土壤层中的水分，是联系地表水、地下水的纽带，同时也是物质传输和运移的载体。
本实验基于灌溉下渗实验模型，主要利用时域反射仪（TDR）连续测量土壤剖面水分含量的方法，分别观测北京通州区土壤和昆明土壤在降水入渗与再分配过程。
根据采集的数据经过图像处理与HYDRUS-1D软件模拟图像进行对比得出不同种土壤水入渗过程中不同深度初始含水率以及下渗率、剖面含水率等对实验的影响，据此得出相关结论。
[关键字]土壤下渗灌溉下渗模型时域反射仪初始含水率剖面含水率1 概述“四水”转化是一个庞大的系统工程，对于其中的每一方面而言，都是一个相对独立的子系统，每个子系统各自独立、自成体系，但它们又相互关联、相互依存、相互转化和相互制约，形成了一个独特的循环转化系统。
2 试验条件与方法研究区北京市通州区通州区位于北京市东南部，通州区地处永定河、潮白河冲积洪积平原，地势平坦，自西北向东南倾斜。
其土质多为潮黄土、两合土、沙壤土，土壤肥沃，质地适中。
云南省昆明市位于北纬亚热带，素以”春城”之称而享誉中外。
四是冬无严寒，日照充足，天睛少雨。
两地从地形、气候、降水、土壤等多方面具有不同特点，差异较大，以两地为研究对象进行对比研究实验结果将相对明显。
本实验的土柱灌溉下渗实验模型是由北京土和昆明土的两个土柱以及时域反射仪组成。
从试验地选取优质土进行分装，实验装置如图1所示。
土柱装填在有机玻璃柱中，有机玻璃柱高度65.0cm，内径28.1cm，在不同的深度剖面上水平埋设时域反射仪TDR）探针，每10分钟自动测定一次土壤水分含水率。
时域反射仪的基本原理是，高频电磁脉冲沿传输线在土壤中传播的速度依赖于土壤的介电特性，在一定的电磁波频率范围内（50M～10GHz），矿物质、空气和水的介电特性为常数，因此土体的介电常数主要依赖于土壤容积含水量，这样可以建立土壤容积含水量与土壤介电常数的经验方程，TDR通过测量高频电磁脉冲在土壤中的传播速度求得土壤的介电常数从而计算出土壤的含水量[2]。
中国农业大学科技成果——饱和-非饱和介质中水分和溶质不规则迁移机理与模拟
中国农业大学科技成果——饱和-非饱和介质中水分和溶质不规则迁移机理与模拟成果简介饱和-非饱和介质中水分与溶质迁移转化理论，是农业水资源可持续利用与水土环境保护的重要理论基础。
本项目主要以认识含水介质中水分与溶质不规则迁移规律，发展相应的模拟理论与方法为目标，以室内外试验、理论分析和数值模拟为手段，系统开展了非饱和土壤水分动态随机与分形模拟理论、含水介质中溶质非费克迁移理论及模拟模型、含水层中抽水井附近非达西流理论与模型、区域地下水动态模拟的理论与方法等方面的研究。
技术内容（1）提出了土壤水力特性参数求解的新方法，建立了非饱和土壤水分动态随机模拟的新理论；（2）探讨了介质非均质性对溶质非费克迁移的影响机理，建立了基于分数微分理论的吸附性溶质迁移的模拟模型，发展了模型求解的新方法；（3）研究了非达西流动的机理，建立了抽水条件下非达西流动的模拟模型，提出了抽水井水力学计算的新方法；（4）结合我国北方地区的实际情况，建立了地下水动态模拟的人工神经网络（ANN）与FEFLOW耦合模型。
该项目解决了含水介质中水分与溶质不规则迁移的一些关键性理论问题，提高了土壤与含水层中水分与溶质迁移模拟与预测的精度和可靠性，揭示了饱和与非饱和介质中水分与溶质不规则迁移的内在机理，深化了水分与溶质迁移不规则现象的认识，丰富了农业资源与环境、农业水利学科的内容，对相关学科的发展有重要促进和推动作用。
该项目在非饱和水分随机运动、非达西流动、非费克迁移的微观机理和宏观规律的研究方面具有明显的创新性，国际著名水文学家Fred Molz教授等认为我们的有关成果达到了国际先进水平。
研究成果在农业水资源高效利用和水土环境保护等领域中具有潜在的应用价值。
投资预算及效益分析研究成果在农业水资源高效利用和水土环境保护等领域中具有潜在的应用价值。
封丘地区土壤传递函数的研究
*国家重点基础研究发展规划项目(G1999011803)和ACIAR 项目(LWR1/96/164)资助收稿日期:2000-12-28;收到修改稿日期:2002-03-24封丘地区土壤传递函数的研究*朱安宁 张佳宝(中国科学院南京土壤研究所,南京 210008) 陈效民 陈德立(南京农业大学资源与环境学院,南京 210095) (澳大利亚墨尔本大学)摘 要 以黄淮海平原封丘地区的潮土和风沙土为研究对象,试图寻求解决土壤水动力学特性参数问题的更加实用可行的方法。
根据大量的土壤基本物理性质和土壤持水数据,利用多元逐步回归分析方法,分别建立了van Genuchten 模型的参数H r 、H s 、a 、n 的土壤传递函数模型,并通过统计分析和数学模拟进行模型验证。
研究结果表明,van Genuchten 模型参数与土壤的质地、有机质含量及容重等基本物理性质有一定的线性相关性。
van Genuchten 模型的参数估计模型对粘性土壤的拟合效果较好,而拟合砂性土壤时误差较大。
关键词 土壤物理性质,水分特征曲线,van Genuchten 模型,土壤传递函数,数学模拟中图分类号 S152172当前,区域土壤水分和溶质运移的研究已成为土壤物理学研究的一大热点,然而,用计算机数值模拟土壤水和溶质运移所遇到的主要问题之一是难于获取土壤水分运动参数,如土壤水分特征曲线h (H )、饱和导水率K s 、非饱和导水率K (H )等。
Childs 早在1940年就注意到土壤物理性质影响土壤水分特性,特别是70年代末以来,许多学者在研究土壤水分运动参数与土壤基本理化性质关系方面做了大量的工作[1,2],从而建立了研究土壤水分运动参数的一种新方法)))PTFs 法(Pedo -Transfer Functions,土壤传递函数)。
PTFs 法是用统计模型由土壤基本性质预测土壤水分运动参数的方法[3,4],即由已知的土壤水分运动参数和与之相对应的土壤基本理化性质建立多元回归方程,用此回归方程结合土壤基本理化性质来预测土壤的水分运动参数。
基于HYDRU S模型的盐碱地土壤水盐运移模拟
基于HYDRU S模型的盐碱地土壤水盐运移模拟潘延鑫;罗纨;贾忠华;井思媛;李山;武迪【摘要】To determine the movement of salt and water in the saline-alkali flatland of Lubotan,Shaanxi province, based on the saturated-unsaturated soil water and solute transport theory,field monitoring data of local water and salt for many years was applied to simulate the rules of local soil water and salt movement,using the HYDRUS 1D numerical model .The soil water and salt changes were analyzed,and the reasonable field irrigation quota was proposed .The results showed that during the whole reproductive period,soil water content had a similar variation trend under different irrigation quotas .Considering water saving and salt control,farmland irrigation quota of 500 m3·hm-2 is reasonable to control salt accumulation in soil .The simulated results of soil water and salt migration using HYDRUS model are basically consistent with the observed values in field experiment,and the results can be referred for farmland management of water and salt in semblable salinity regions .%为了解陕西卤泊滩盐碱地的水盐运移情况，基于当地2009—2013年田间水盐监测资料，应用饱和－非饱和土壤水分及溶质运移理论，利用HYDRUS－1 D数值模型对当地土壤水分、盐分运移规律进行数值模拟，分析了盐碱地的水盐变化状况，确定合理的田间灌水定额。
非饱和土壤导水率试验计算与模拟分析
非饱和土壤导水率试验计算与模拟分析胡钜鑫;虎胆·吐马尔白;穆丽德尔·托伙加;杨未静【摘要】以非饱和土壤导水率作为研究对象,用瞬时剖面法计算两种土壤非饱和土壤导水率,并与RETC中不同模型的模拟结果进行对比,研究瞬时剖面法计算结果的可靠性.结果表明:两种土壤的K-h与lgK-h模拟曲线和实测值均吻合较好,实测值和不同模型的模拟值均属于高度性相关,且K-θ实测曲线与各模型的模拟曲线变化规律相似,处于各模拟曲线之间.综上所述,瞬时剖面法计算结果与模拟结果相似,具有一定的准确性,可以直接使用在实际生产运用过程中.【期刊名称】《石河子大学学报（自然科学版）》【年(卷),期】2019(037)001【总页数】7页(P105-111)【关键词】非饱和土壤导水率;瞬时剖面法;van Genuchten模型;Mualem模型【作者】胡钜鑫;虎胆·吐马尔白;穆丽德尔·托伙加;杨未静【作者单位】新疆农业大学水利与土木工程学院,新疆乌鲁木齐市,830052;水文水资源与水利工程科学国家重点实验室,江苏南京,210098;新疆农业大学水利与土木工程学院,新疆乌鲁木齐市,830052;新疆农业大学水利与土木工程学院,新疆乌鲁木齐市,830052;新疆农业大学水利与土木工程学院,新疆乌鲁木齐市,830052【正文语种】中文【中图分类】S152.7非饱和土壤土导水率K 是土壤水分参数中的重要参数之一，它反⒊了土壤中的水分在非饱和状态下的运动规律。
非饱和土壤导水率的测定方法包括直接法和间接法，直接法又分为田间测定和室内测定。
田间测定方法包括结壳法[1]、圆盘入渗法[2-4]、双环法[5]等，室内测定方法包括瞬时剖面法、垂直下渗通量法、零通量法[6]等。
其中直接测量法通常耗时耗力，不易测量，因此大部分学者常选⒚间接方法求取非饱和导水率，包括土壤水分再分布法[7-8]，或者通过水分特征曲线C 和水平扩散度D 公式推求非饱和土壤导水率K[9]，另外通过模拟软件[10]，例如Hydrus 和RETC 通过土壤质地资料推求非饱和导水率[11-13]。
自由渗流面及渗流量推求的数值计算方法探讨
自由渗流面及渗流量推求的数值计算方法探讨秦甜甜;丁国辉;程勤波【摘要】自由渗流面具有复杂的非线性,较难确定.本文采用开源地下水数值模拟程序MODFLOW及SUTRA,分别运用MODFLOW模型中干湿单元转化技术、SUTRA模型中单元渗透矩阵调整法以及本文建立的缓变渗透系数矩阵法推求自由渗流面.对比其求解结果表明,采用MODFLOW运用干湿转化技术求解自由渗透面的方法稳定性最好、精度最高,而采用缓变渗透系数矩阵法的SUTRA程序,改善了传统单元渗透矩阵调整法的不稳定性,提高了数值计算精度,避免了MODFLOW必须矩形网格的局限性,是一种实用的计算自由渗流面及估算地下水与河流水量交换量的方法.【期刊名称】《地下水》【年(卷),期】2018(040)004【总页数】3页(P12-14)【关键词】MODFLOW;SUTRA;单元渗透矩阵调整法;缓变渗透系数矩阵法;自由渗流面【作者】秦甜甜;丁国辉;程勤波【作者单位】江苏省地质调查研究院,江苏南京 210018;南京市测绘勘察研究院股份有限公司,江苏南京 210005;河海大学水文水资源学院/河海大学水文水资源与水利工程科学国家重点实验室,江苏南京 210098【正文语种】中文【中图分类】P641.2自由渗流面与地下水潜水水面特征相似，具有复杂的非线性，是河流等地表水入渗补给地下水、大坝渗流等定量化计算的难点，目前，尚未发现确定自由渗流面的解析法，在实际工程应用中常用数值方法求解。
求解该类问题的数值方法可分为两类：变网格法和固定网格法。
由于变网格法操作复杂，容易使计算网格畸变，难以处理存在水平介质分层及各种复杂夹层的情况等缺点,因而研究较少。
为此，国内外学者提出了许多与固定网格法相关的处理技术，主要有剩余流量法、初流量法、节点虚流量法、单元渗透矩阵调整法、复合单元法及截止负压法等[1-6]。
但这些处理技术大多方法复杂，操作困难，应用不便。
本文利用开源地下水数值模拟软件MODFLOW[7]及SUTRA[8]，提出了采用干湿单元转化、缓变渗透系数矩阵法推求自由渗流面的方法，成功推求了自由渗流面，有效地提高了目前河流与地下水交互作用计算的求解精度。
不同出流方式下垂直流人工湿地数值模拟
不同出流方式下垂直流人工湿地数值模拟范思雨;龙天渝【摘要】建立了二级垂直流人工湿地的二维非恒定微生物动力学模型,针对淹没出流方式、不同的等间隔淹没出流方式和定水头出流方式,采用数值模拟的方法研究二级垂直流人工湿地处理重庆市径流雨水的效果;结果表明,在相同工况下,间隔较大的淹没出流方式能达到最大的有效体积比和水力效率,水流在三层填料中与各种类微生物接触更加充分、持久,在此出流方式下促进了溶解氧的传质,复氧后溶解氧浓度回升更快、浓度峰值更大、高浓度维持时间更长,因此促进了由微生物主导的生物化学降解反应,二级垂直流湿地对COD、氮和磷的去除效果也优于其他出流方式.【期刊名称】《重庆工商大学学报（自然科学版）》【年(卷),期】2017(034)006【总页数】9页(P110-118)【关键词】人工湿地;出流方式;CW2D;数值模拟【作者】范思雨;龙天渝【作者单位】重庆大学三峡库区生态环境教育部重点实验室,重庆,400045;重庆大学三峡库区生态环境教育部重点实验室,重庆,400045;重庆大学低碳绿色建筑国际联合研究中心,重庆,400045【正文语种】中文【中图分类】X703人工湿地污水处理技术因其显著的低能耗、低维护成本、处理范围广和良好的景观效应而被广泛应用[1]，又因其较强的抗冲击负荷能力而在蓄纳、处理城市径流雨水上有很好的应用前景[2-6]。
目前大多数人工湿地的设计与运行依然建立在经验的基础上，对其全运行过程中复杂的物理-化学-生物协同作用认知的“黑箱”现象依然存在，而这一作用与人工湿地内的水流规律密切相关[5]。
在模拟研究中将湿地的出流方式抽象成边界条件，边界条件的设定在流体力学的计算中极大地影响了水流规律，因此对人工湿地出流方式的探讨和研究对提高其处理效率、指导工程设计具有重要意义。
目前对于人工湿地进/出水方式的研究主要是两个方面：一是讨论穿孔管在进水处和出水处的配/集水的均匀性；二是讨论进/出水口在垂直方向上的位置(如上、中、下)及其组合搭配对配水均匀性和处理效果的影响。
求van Genuchten模型参数的AM-MCMC方法
求van Genuchten模型参数的AM-MCMC方法石晓蕾;徐绍辉;廖凯华【期刊名称】《土壤》【年(卷),期】2012(44)2【摘要】采用基于自适应采样算法的马尔科夫链蒙特卡罗方法(简称AM-MCMC)来估算描述土壤水分特征曲线的van Genuchten模型的参数,并推求出模型参数的后验分布,从而为模型参数的不确定性分析提供依据.结果表明,对于van Genuchten 模型而言,采用AM-MCMC算法能得到模型参数后验均值和方差的分布,并且能推求出模型参数的置信区间,所以用这种算法来求解van Genuchten方程的参数是行之有效的,为求解van Genuchten模型参数提供了一种新的思路.%This paper adopted the adaptive Metropolis algorithm of Markov chain Monte Carlo approach (AM-MCMC) to estimate the parameters of van Genuchten model, which describes the soil water characteristic curve. The posterior statistics of the van Genuchten model parameters were obtained to estimated the model uncertainty. It was found that the distribution of the posterior means and variances of the model parameters could be obtained. Confidence intervals of the model parameters could be ascertained. It was proved effective and novel in solving van Genuchten equation with this algorithm.【总页数】6页(P345-350)【作者】石晓蕾;徐绍辉;廖凯华【作者单位】青岛大学环境科学系,青岛266071;青岛大学环境科学系,青岛266071;南京大学水科学系,南京210093【正文语种】中文【中图分类】S152【相关文献】1.van Genuchten模型参数的物理意义 [J], 陈卫金;程东会;陶伟2.土壤持水曲线van Genuchten模型求参的一种新方法 [J], 刘贤赵;李嘉竹;张振华3.土壤持水曲线van Genuchten模型求参的Matlab实现 [J], 魏义长;刘作新;康玲玲;王云璋;时明立4.利用优化方法求算Van Genuchten方程参数 [J], 李春友;任理;李保国5.土壤水分特征曲线Van-Genuchten模型参数的土壤传输函数比选 [J], 李彬楠;樊贵盛因版权原因，仅展示原文概要，查看原文内容请购买。
温度梯度作用下非饱和黄土水分运移规律
1.1控制方程 土体蒸发过程主要受环境温度、风速、太阳辐射
及土体特性参数等因素的影响，同时土体热量传输也
是影响土体水分蒸发的重要因素。 在数值计算中,考
虑大气环境与土体水分、热量的交换，建立蒸发条件
下非饱和土体水分热量传输模型，并研究分析蒸发状
态时，非饱和黄土内部土体水分热量传输变化量。
(1) 土体水分传输方程在蒸发状态时，探究 非饱和土体水分运动变化规律, 仅考虑土体水分在
K(h)为土体的渗透系数(cm/s);S为植物根系的
水分吸收量，表层取值0o
(2) 土体热量传输方程建立土体热量传输方
程时，不考虑气态水的扩散过程，只考虑液态水对土
体热量传输的影响，其方程式为
dT d A (3) T Cq J,
dz
dZ
(2)
其中:C(3)为介质的比热容(J・g 1・C 1)Cw为
液体比热容(J・g 1・C 1)；(3)为土体导热率(J・
(a)距温顶度端梯距度离5〜/1c5m°C
(b)距温顶度端梯距度离5〜/2c0m°C
Fig. 2
(c)距温顶度端梯距度离5〜/2c5m°C
* lh
亠5h
(d)距温顶度端梯距度离5〜/3c0m°C
20 h - 24 h
图2不同温度梯度作用下土样温度场动态变化
Dynamic change of soil temperature field under different temperature gradients
等方法，对土体中水热运移耦合规律进行了研究；李 彦龙等［57］采用模型模拟试验,对温度梯度条件下的 水分运移规律、非饱和土的基质势梯度和结合水特征 进行了研究分析；石兰君等［］基于室外实际模拟，用 HYDRUS-1D软件对浅层包气带水分运动状态建立 了数值计算模拟,并研究分析了不同时间土体不同深 度处含水量的变化规律；闫亚景等［9］在对边坡水分迁 移研究中，采用了高密度电法;周宏［0］分析了干旱气 候条件下土壤水分迁移的驱动力与能量之间关系；林 宗泽等［11］采用理论与现场验证试验的手段，对土体 水分蒸发过程进行了分析。
HYDRUS 1D使用手册-以亭子口水库模拟为例
一、软件介绍
Hydrus-1D 是美国盐土实验室开发的，计算包气带水分、溶质运移规律的软 件， 用它可以计算在不同边界条件和初始条件下的数学模型。若将坐标原点选在 地面，取 z 轴向下为正，则一维饱和—非饱和带水分运移基本方程为：
其中，如前所述，有些边界条件是时变边界，需要输入时变的记录。 9、上、下边界条件数据输入
则上边界需要输入不同时刻的水 此例代表上下边界均选择变水头边界条件， 头；下边界需要输入不同时刻的地下水位（实际也是水头） 。 加入选择的是定水头或定流量等边界，则仅需要输入水头或流量的值。特别 的，隔水边界是一种流量为 0 的边界。 有些边界带有模型参数，如 Deep Drainage 边界。 10、土壤剖面-图像编辑器 剖面-图像编辑器可将模拟土层的离散划分、土壤质地划分、土壤亚区划分、 土壤初始的含水量（负压）等条件进行编辑，并用图像的形式显现出来。 在此步骤中选择菜单栏中的 condition，其下有 6 个选项。 Material Distribution 剖面离散 Root Distribution 根系分布 Scaling Factors 尺度因素 Initial Conditions 初始条件（初始的负压含水量分布等） Subregion Distribution 子区划分 Observation Nodes 设置观察点 其中在 2.模拟土壤的几何属性中，已经输入了土质的种类、土壤亚区的划分 数目等，这里还要指定土壤整个剖面中，亚区的起止范围及对应的质地。点击 Edit condition，然后在所示土柱中拖动鼠标进行编辑。 剖面离散这一项根据自己的需要，将土壤划分有限元计算单元。里面对 Edit Node 进行编辑，可以编辑节点数目，节点划分密度等。 设置观察点时，直接用鼠标在土柱需要观察的点进行左击。 设置初始条件和根系分布时，需要设定范围边界和边界点的值，再使值线性 分布。此处设置往往比较麻烦，一般下一步还可以设置。
间接推求非饱和土壤导水参数的方法
间接推求非饱和土壤导水参数的方法【摘要】本文介绍了四种间接推求非饱和土壤导水参数的方法，包括土壤水分再分布过程法、土壤水分特征曲线拟合模型法、简单入渗法和土壤转换函数法，并对四种方法出归纳，提出了进一步研究的方向。
【关键词】非饱和土壤导水参数；水分再分布；水分特征曲线；简单入渗法；PTFs确定非饱和土壤导水参数的方法有两大类：直接法和间接法。
鉴于直接法，耗时多，价格昂贵，所得精度难以保证，空间差异性大，间接法成为许多专家学者研究的焦点。
本文通过介绍推求非饱和土壤导水参数的主要几种间接方法，做出归纳，以期为非饱和土壤导水参数间接法的进一步研究起到推动作用。
1.土壤水分再分布过程法邵明安[1]提出的一种土壤导水参数模型，将一维垂直和水平土壤水分再分布过程结合起来，根据土壤湿润锋湿度与湿润剖面平均湿度之间的函数关系，推出非饱和导水率K(θ)的解析表达式。
这样，只要利用土壤水分再分布的湿润过程，确定表达式中的常数，就可直接算出土壤的非饱和导水率。
根据湿润锋湿度与平均湿度拟合关系的不同，推出K(θ)有不同的表达式，本文采用湿润锋湿度与平均湿度为幂函数关系：K(θ)=（θ-θinm■■-Z■■-n■m■■■-X■■）上式：θ为体积含水率，θi为初始含水率，H为水深，Z、X分别为垂直和水平方向水分运移的距离，Z■和X■为所加其水面刚刚消失时初始湿润深度，m、n、m1、n1、a、b均为拟合参数。
2.土壤水分特征曲线拟合模型法2.1 Mualem模型Mualem模型关于相对导水率的表达式如下：Kr(Θ)=Θ1/2■■/■■■（1）式中：Kr为相对导水率，即K/Ks，h为水头值，Θ为无量纲含水量。
无量纲含水量表达式Θ=■（2）求解方程（1）时需要建立无量纲含水量与水头值的关系式，采用幂函数形式：Θ=（ａｈ）－ｂ（3）将（3）式代入方程（1）积分得：Kr(Θ)＝Θ５／２＋２／ｂ（4）2.2 Burdine模型Burdine建立的相对导水率表达式如下：Kr(Θ)＝Θ２■■／■■（5）将（3）式代入方程（5）积分得：Kr(Θ)＝Θ３＋２／ｂ（6）2.3 Van Genuchten模型[2]Van Genuchten水分特征曲线模型为：Θ＝［１＋（ａｈ）ｎ］－ｍ（7）式中：α是标定参数，与土壤平均孔隙半径成正比，n和m是土壤水分特征曲线的形状参数或孔隙分布指数,且有m=1-1/n。
一种非饱和土相对渗透系数的试验数值联合估计法
一种非饱和土相对渗透系数的试验数值联合估计法作者：赵晓龙邱秀梅卞汉兵邱庆泰来源：《南水北调与水利科技》2015年第06期摘要：由于非饱和土相对渗透系数对饱和度的变化十分敏感，因此通过试验测定非饱和土相对渗透系数并不容易。
以广义达西定律和van Genuchten模型为理论依据，可通过简单的蒸发试验和数值分析，获得非饱和状态下土壤的相对渗透系数。
模拟结果和试验数据的对比分析表明，该方法得到的非饱和土相对渗透系数具有一定的可靠性，在有限元反分析模型参数方面具有一定的创新性，但其实际精度还有待通过试验对比来进一步验证。
关键词：非饱和土；饱和度；蒸发试验；数值分析；相对渗透系数中图分类号：TU43 文献标志码：A 文章编号：1672-1683（2015）06-1114-04Abstract：The relative hydraulic conductivity of unsaturated soil is sensitive to the change of saturation，so it is hard to determine the relative hydraulic conductivity through experiments.Based on the generalized Darcy′s Law and van Genuchten′s model，a simplified numerical model for predicting the relative hydraulic conductivity of unsaturated soil was derived.Through the simple evaporation experiment and numerical analysis，the relative hydraulic conductivity with function of water saturation was obtained.The numerical results were compared with the observed experimental data，which suggested that this method has certain reliability and it has certain innovation in the aspects of the finite element back analysis of model parameters.However，the practical accuracy of the new method needs further verification through experimental comparison.Key words：unsaturated soil；saturation；evaporation experiment；numerical analysis；relative hydraulic conductivity土的渗透系数及相对渗透系数是土的主要特性指标之一，是描述孔隙水在土壤孔隙中流动性质一种度量[1]，是反映水在非饱和土壤孔隙中迁移的重要参数，在雨水引起的边坡稳定性、路基变形、土石坝黏土心墙渗流控制、岩土工程施工、农业土壤环境研究中起着重要的作用，因此是国内外非饱和土研究中的一个重要课题。
推求土壤非饱和运动参数的方法
推求土壤非饱和运动参数的方法硕士生：景为学科专业名称：土壤学研究方向：土壤水分动力学指导教师：邵明安研究员准确获取能代表田间土壤条件的土壤水分运动参数（土壤水分特征曲线（或比水容重C）、土壤导水率K和土壤水分扩散率D）是模拟土壤中水和溶质运动的基础。
三个参数中，以预测非饱和导水率最为困难，原因之一在于直接测定困难。
对土壤水分运动参数空间变异性认识的加深将有助于预报田间水分和溶质迁移过程，也有助于完善参数确定的方法，使之更具普遍性。
在以往的研究中，已有许多直接测定或间接推求这些参数的方法。
本文选取了其中的三种代表方法，以实测水分特征曲线作为标准进行比较，评价各自的优缺点及适应范围。
三种方法是：(1)实测土壤水分特征曲线；(2)用简单入渗法推求van Genuchten水分特征曲线模型中的参数α和n,通过实测饱和导水率Ks，结合导水率模型而获得非饱和导水率K；(3) 根据土壤水分水平和垂直再分布过程直接推求非饱和导水率K和扩散率D。
研究结果表明：1．四种非扰动土壤饱和导水率具很大的差别，其半方差随间距加大而增加，但很快达到一个稳定值,此值即为其变异性的空间尺寸，沙土、黑垆土的空间尺寸为2m，黄绵土的为2.24m，娄土的则更小。
2．土壤水分再分布实验表明，用三种函数拟合湿润锋湿度与平均湿度的关系时，以指数函数拟合计算的比水容重值与实测值最为吻合，尤其是沙土、黄绵土、娄土。
3．利用简单入渗法估计van Genuchten水分特征曲线模型模型中的参数时，α和n值推求的准确度就主要取决于S值测定的准确度，而S的准确测定较易实现，由此可断定简单入渗法的准确性较高。
4．在三种推求导水参数的方法中，水分再分布方法准确性较差，但它无需测定水分特征曲线即可直接得到土壤导水参数K和D，是一种非常简便的方法，尤其适宜于黄土高原沙土导水参数的测定；由简单入渗法获得的水分特征曲线与实测值吻合最好，随着质地变细，拟合效果更好，适合于黄土高原黄绵土、黑垆土和娄土导水参数的测定，而且还解决了Van Genuchten模型中参数不唯一的问题，实验简便，省时（约需2天），计算简单，结果准确，具有很大的优越性。
土体初始饱和度对降雨入渗的影响
土体初始饱和度对降雨入渗的影响薛松;童富果;郝霜;王军【摘要】Based on the theory of gas-liquid two-phase flow, the finite element method was used to simulate rainfall infiltration in loam soil with different degrees of saturation to investigate the effects of the initial degree of saturation on rainfall infiltration. Numerical results show that the pore air pressure increases gradually in the initial stage of rainfall infiltration;then, the infiltration rate decreases to a stable value because of the jacking force caused by pore air. The stable infiltration rate is mainly determined by the matrix suction and permeability characteristics of soil, varying with the initial degree of saturation. With the variation of the initial degree of saturation, the infiltration rate reaches its maximum or minimum values. The maximum infiltration rate usually appears in the stage of residual or full saturation, while the minimum infiltration rate occurs in the unsaturated stage, depending on soil unsaturated characteristics. Affected by the matrix suction of soil and the relative permeability of water and air, the infiltration rate gradually decreases with the increase of the degree of saturation at first. When the infiltration rate drops to the minimum, it starts to rise with the increase of the degree of saturation. According to the results of numerical simulation, a formula is proposed for calculation of the relatively stable infiltration rate for different initial degrees of saturation.%基于水、气两相流理论，采用有限单元法计算了壤土在不同饱和度下的降雨入渗过程，以研究初始饱和度变化对降雨入渗的影响。
两相流_中文
流体流动模型: 两相流两相流简介以下例子分析多孔介质中的两相流. 描述流体如何不发生混合同时通过一个孔隙可以回答很多环境和工业问题. 不幸的是, 多相分析很复杂需要求解一些未知量以及多个相关变量. 多相彼此之间的液体属性与每个液相压力和饱和级有关.这个问题验证两相流按照美国环境保护署的实验结果(参考文献1). 这个试验表明实验和数值计算的两相流结果相符. 通过这些实验, 科研人员估计流体对(空气-水, 空气-石油, 和石油-水)变化的流动, 然后使实验结果和计算机使用保持和渗透解析表达式数值计算的结果相匹配. 这个讨论根据Lincoln土壤和使用Mualem (参考文献2)和van Genuchten(参考文献3)给出的流体属性解析公式展开.图2-42: 在Hopmans和其他(参考文献1)两相流实验几何结构.这个例子分多步. 第一部分建立水和空气的两相流模型; 方程求解压力. 饱和度随解而变化. 需要强调假设在整个分析时间里土体中存在一些残余空气和水. 模型追踪空气前锋移动过湿液体通过观察饱和度而不是假设离散界面. 第二部分修改空气－水属性对空气－石油和空气－水系统进行仿真.模型定义在空气和水的实验中, 插入空气到充满水和沙的实验柱体中. 进入的空气 (对该流体对为非湿相)强迫水(湿相)向柱体基础出口流动. 在入口, 空气压力随时间逐步增加, 并且在柱体顶端没有水存在. 在向出口流动, 水通过一个空气流动不可渗透的圆盘. 空气和水都不能通过柱体垂向壁面. 出口处的水压随时间变化, 对应着滴定管中液体高度增加. 柱体总体高度8.34 cm, 半径6cm, 并且圆盘厚度0.74 cm. 实验时间170 hours.以下讨论给出两相流仿真计算的步骤. 首先观察控制方程和定义流体保持和渗透的解析关系(参考文献2和参考文献3). 然后给出一些详细的执行过程和模型数据表. 后面给出结果. 最后在图形用户界面逐步建立模型.控制方程和边界条件多孔介质两相流满足湿(w)和非湿(nw)流体的独立方程:(2-4)(2-5)其中θs 为总气隙或者饱和体积率; Se 为有效饱和函数; t 为时间(s); κint 为多孔介质固有渗透率(m 2); k r 为给定液体的相对渗透函数; η为流体动力粘度(kg/(m·s)); p 为压力(kg/(m·s 2)); ρ为流体密度(kg/m 3); g 为重力加速度; 和D 为垂向坐标(例如, x , y , 或者z )(m).如果流体分布为连续的, 流体在整个分析时间内未完全填充土壤, 给定湿相容积率, θw , 和非湿项, θnw . 对湿相, θ从0或者很小的残余值θr 到总体孔隙率, θs . 有效饱和度, Se , 源自比例系数θ关于θs 和θr 从0到1变化.θ和Se 为系统中所有液体压力的函数. 定义毛细压力:(2-6).气孔可以在给定时间内被一种气体完全填充: (2-7).饱和度随毛细压力如何变化, 如下所示(2-8)其中C 为指定容积, 下表“p ”表示单位压力.使用方程2-6, 方程2-7, 和方程2-8在方程2-4和方程2-5简化数值模型. 控制方程变为:(2-9)(2-10).可以同时求解系统方程得到p w 和p nw . 这个例子中, 两种流体不可压缩, 但是无需如此假设.初始, 柱体中的水和空气服从静力分布. 边界条件允许水仅存在于土壤柱体的基础中. 对于湿相, 边界条件为(2-11)其中n为边界法向的单位向量.因为空气从柱体顶端进入但是没有出来, 所以非湿项的边界条件为(2-12)保持和渗透关系可以从试验数据, 任意数学公式, 和模型中的其他方程内插建立两相流分析并定义系数θ, C, Se, k r, 和p c如何同时变化. 已有模型从参考文献2和参考文献3中的保持和渗透关系表示系数θ, C, Se, 和k r作为p c函数的变化. 因为p c较大并且由于θ, C, Se,和k r的变化较小, 这些表达式转化为与相同水柱等高的毛细压力或者毛细压力头如H c = p c/(ρwater g). 水压属性相对湿流体为(2-13)其中α, n, m, 和L为van Genuchten参数表示土壤属性. 注意两相流, van Genuchten-Mualem公式与H c 有关.对于非湿流体, 属性为(2-14)对湿相的定义同理可得.不同流体对在空气-水, 空气-油, 和油-水之间转换试验, 可以根据Leverett (参考文献4)灵活使用界面张力比例. 不考虑流体对Leverett比例调整柱体顶端的非湿相压力生成相同体积的湿流体从柱体底端流出. 使用With Leverett比例, 在流体对之间转换要求使用正确的流体对的流体属性ρ和η并根据以下公式调整边界和初始压力在这些方程中, σ表示不同流体之间的界面张力, 并且下标表示流体对. 这些值出现在这部分内容末尾的表格中. 例如, σao /σaw等于0.373, 和σwo/σaw等于0.534; 此外, 对空气－水系统第一项非湿相压力头(水以米为单位)为0.4 m, 空气－油系统为0.1 m, 和水－油系统为 0.2 m.因为多孔介质的相对渗透率和保持属性与流体流过该介质有关, 所以要求改变流体对的时候也改变模型的保持和渗透属性. 这个要求可以通过插入新试验数据或者调整数学公式实现. 这个模型中, 作者通过对解析公式拟合得到渗透率和保持率参数. 通过调整参数α, n, κs, 和θr为每种流体得到最佳拟合. 以下数据表显示不同流体对对应的比值α上面给出的σ比值近似相等.执行过程: 数值微分估计C这个例子使用解析表达式估计指定湿度, C. 因为C为q和H c的曲线斜率, 也可以通过使用COMSOL Multiphysics微分运算符计算得到C例如Cp,w(p w) = pdiff(theta_w,Hc)ρwater-1g-1因此可以书写任意表达式或者使用数据θ. 这个方程对流动体积率的微分考虑到毛细压力头. 此处除以水的重量将表达式转化为与压力变化相关. 相同类型的符号diff可进行微分运算. 两个运算符的差别在于diff可识别空间和时间的导数并且已经定义在COMSOL Multiphysics中(例如, p x, p y, p t), 然而pdiff为纯符号不能应用到与变量相关的链式法则中.执行过程: 逐步改变边界以下部分逐步减少通过使用内插函数定义入口处时变步进非湿相压力. 在COMSOL Multiphysics中内插直接可选. 打开选项菜单对话框, 建立时间和相应压力头表, 为插值函数定义名称, 并在其他需要使用函数的地方通过函数名调用, 在求解域设定对话框. 为了激活创建的函数, 输入带变量的函数名(例如, p_nw), 变量即括号中的时间t. 命令如下p_nw(t) = H_nw_t(t)*rho_water*g_w方程中出现水的密度因为参考文献1利用水柱高度定义边界压力.数据模型中使用的数据对应参考文献1给出的Lincoln 沙中空气－水的试验结果:不同流体对的van Genuchten 参数如下所示:根据下表入口处的空气压力头随时间增加.变量 描述 单位表达式g rAcceleration due to gravity m/s 2 9.82 ρfw Fluid density, water kg/m 31000ηw Dynamic viscosity, water Pa·s1·10-3ρfg Fluid density, gas kg/m 31.28ηg Dynamic viscosity, gasPa·s1.81·10-5κint Intrinsic permeability, column m 2 2.48·10-12 κsPermeability, discm 21.33·10-14 q s,w Saturated volume fraction, column0.32 q s,wSaturated volume fraction, disc0.5 p nw,top Initial nonwetting phasepressure head, inletm water0.2变量 描述 单位 空气-水 空气-油 油-水q r,w Residualvolume fraction0.021 0.00001 0.0072α alpha parameter m-1 1.89 5.29 3.58n n parameter, column 2.811 3.002 3.1365 LL parameter, column 0.5 0.5 0.5 κsPermeability, discm -22.48·10-121.09·10-120.94·10-12压力头(m 水) 起始时间(小时)0.4 0 0.6 21.25 0.8 45.251.0 691.5 932 122.54 155在水的出口处, 滴定管的流体级随时间从0 m到0.1 m线性增加.结果和讨论图2-43给出试验柱体从COMSOL Multiphysics两相流解的早期图像. 阴影描述非湿相的有效饱和(空气), 同时箭头给出湿相(水)速度.图2-43: 两相流模型在0.1小时的解: 非湿相饱和(表面图), 湿相速度(箭头). 结果对应从US EPA (参考文献1)中Lincoln土壤的空气－水试验结果.图像给出非湿相流体进入土壤柱体并取代湿相流体. 因为多步压力变化压迫非湿相进入.图2-44给出沿着柱体不同高度的毛细压力在入口边界处的步进压力头. 指定后处理中输入时指定观察点. 两相流问题的解与参考文献1的结果完美一致.毛细压力头标定空气入口压力如图2-44轨迹一致说明试验成功. 为了获得渗透率和保持属性的高分辨率解, 作者在参考文献1设置压力步长足够大使土体中的冲击变为瞬时. 如图2-45所示, 渗透率在整个土体中瞬时变化.图2-44: 在Lincoln土壤(源自US EPA, 参考文献1)中空气－水流动的入口压力头(实线)和毛细压力头(虚线).图2-45: Lincoln土壤在 x = 0.03 m (源自US EPA, 参考文献1)处水的渗透函数(实线)和空气的渗透函数(虚线). 空气－油和油－水系统的两相流动解如图2-46和图2-47所示.图2-46: Lincoln土壤(源自US EPA, 参考文献1)空气－油的入口空气压力头(实线)和毛细压力头(虚线).图2-47: Lincoln土壤 (源自US EPA, 参考文献1)油－水的入口空气压力头(实线)和毛细压力头(虚线).使用COMSOL Multiphysics计算的空气－油和油－水两相流问题的结果证明与参考文献1的结果完美一致. 通过Leverett比例系数可以设置入口压力以便空气－油和油－水系统从空气－水试验生成体积流出率. 与空气－水系统相同, 空气－油试验的毛细压力头和空气入口压力为瞬态的. 然而, 对水－油系统, 非湿相和湿相之间有滞后.参考文献1. J.W. Hopmans, M.E. Grismer, J. Chen, and Y.P. Liu, Parameter estimation of two-fluid capillary pressure saturation and permeability functions, U.S. Environmental Protection Agency EPA/600/R-98/046, Cincinnati, Ohio, 1998.2. Y. Mualem, “A new model for predicting the hydraulic permeability of unsaturated porous media,” Water Res. Research, vol. 12, 1976, pp. 513–522.3. M.Th. van Genuchten, “A closed-form equation for predicting the hydraulic of conductivity of unsaturated soils,” Soil Sci. Soc. Am. J., vol. 44, 1980, pp. 892–898.4. M.C. Leverett, “Capillary behavior in porous solids,” Trans. AIME, vol. 142, 1941, pp. 152–169.模型库路径: Earth_Science_Module/Fluid_Flow/twophase_flow_air_water 使用图形用户界面建模: 空气-水系统模型导航视窗1 打开模型导航视窗下的新增页面, 在空间维度列表选择2D.2 从应用模式列表选择地球科学模块>流体流动>达西定律>压力分析>瞬态分析.3 在因变量编辑区输入p_w, 并在应用模式名编辑区输入w.4 单击多物理场按钮, 然后单击新增.5 再次在应用模式列表选择地球科学模块>流体流动>达西定律>压力分析>瞬态分析.6 重复以下步骤, 改变因变量为p_nw和应用模式名为nw. 单击多物理场按钮, 然后单击新增.7 单击确定.几何建模绘制两个矩形创建几何结构.1 从绘图菜单选择指定对象>矩形.2 输入以下设置;设置完成, 单击确定.参数值Width 0.06Height 0.0074x 0z 03 除了以下设置重复上面的过程:参数值Width 0.060.0834选项Height x 0 z4 在主工具条单击缩放至视窗大小按钮.1 选择选项>常数.2 输入以下名称, 表达式, 和描述(任意的); 设置完成, 单击确定.名称表达式描述SPH 1[h/s] Seconds per hour rho_water 1000[kg/m^3]Density, waterrho_w 1000[kg/m^3]Density, wetting fluidetaw 0.001[Pa*s]*SPH Dynamic viscosity, wetting fluid rho_nw1.28[kg/m^3]Density, nonwetting fluid etanw 0.0000181[Pa*s]*SPHDynamic viscosity, nonwetting fluid3 从选项菜单选择表达式>标量表达式. 输入以下名称和表达式(都在一条线上); 设定完成, 单击确定.名称表达式p_nw_top 0.2[m]*rho_water*g_nwp_nw_init p_nw_top+(8.34[cm]-y)*rho_nw*g_nw p_w_init -rho_w*g_w*yp_w0 0.1[m]*rho_water*g_w*(t/170[s]) Hc (p_nw-p_w)/(rho_water*g_w)Se_w (1+abs(alpha*Hc)^N)^(-M)*(Hc>0)+1*(Hc<=0) theta_w (theta_r+Se_w*(theta_s-theta_r))*(Hc>0)+theta_s*(Hc<=0) kr_w ((Se_w^L*(1-(1-Se_w^(1/M))^M)^2)+eps)*(Hc>0)+1*(Hc<=0) Cp ((alpha*M/(1-M)*(theta_s-theta_r)*Se_w^(1/M)*(1-Se_w^(1/M))^M))/(rho_water*g_w)*(Hc>0) Se_nw 1-Se_w theta_nw theta_s-theta_wkr_nw ((1-Se_w)^L*(1-Se_w^(1/M))^(2*M))*(Hc>0)+eps4 从选项菜单选择表达式>求解域表达式. 输入以下名称和表达式为两个求解域; 设置完成, 单击确定.名称求解域1求解域2theta_s 0.5 0.32 theta_r 0.0211 kap_s 1.34e-14 2.48e-12alpha1.89N 2.811M 1-1/NL 0.5现在使用内插方法为非湿相边界建立步进压力.5 从选项菜单, 选择函数.6 单击新增按钮.7 在函数名编辑区输入Hp_nw_t并选择内插法选项.8 单击确定.9 输入下表中的值; 完成输入, 单击确定.X Y0 0.421.20 0.421.25 0.645.20 0.645.25 0.868.95 0.869 1.092.95 1.093 1.5122.45 1.5122.5 2154.95 2155 4200 4物理量下面部分首先为每一个相建立材料属性, 初始条件, 和边界条件, 然后将两个方程建立联系.应用标量变量建立模型垂向和重力常数, 转到物理量菜单并选择标量变量. 按照如下参数修改默认值; 完成设定, 单击确定.名称表达式D_w yg_w 9.82[m/s^2]*SPH^2D_nw yg_nw 9.82[m/s^2]*SPH^2这步将模型基本时间单位从秒转化为小时. 因此, 任何时候用户交互界面出现的时间单位为秒—既有显式, 如m/s, 或隐式的, 如Pa (= kg/(m·s2))—代替小时.求解域设定—湿相1 从多物理场菜单选择达西定律(w).边界条件—湿相从物理量选择边界设定, 然后设置以下边界条件; 完成设置, 单击确定.求解域设定—非湿相边界条件—非湿相方程系统2 从物理量菜单选择求解域设定. 输入以下设定, 然后单击应用. 设定 求解域1 求解域2Storage term User defined User definedS0 Cp k Skap_s kap_s*kr_w ρfrho_w rho_w η etaw etaw3 单击初始选项卡. 同时选择求解域1和2, 并在p_w(t 0)编辑区为压力输入p_w_init . 单击确定.设定 边界2 边界1, 3, 5–7Boundary conditionPressure Zero flux/Symmetry p 0 p_w01 从多物理场菜单选择达西定律(nw).2 从物理量菜单选择求解域设定. 在求解域选择列表选择1, 然后清除启动本域的复选框.3 选择求解域2. 输入以下设定.设定 求解域2Storage term User definedSCp κSkap_s*kr_nw ρfrho_nw η etanw4 单击初始选项卡. 同时选择求解域1和2, 并在p_nw(t 0)编辑区为压力输入p_nw_init . 单击确定.1 从物理量菜单, 选择边界设定.2 设置以下边界条件; 完成设置, 单击确定.设置 边界3, 4, 7 边界5BoundaryconditionZero flux/Symmetry Pressure p 0 Hp_nw_t(t)*rho_water*g_nw1 从物理量菜单, 选择方程式系统>求解域设定.在对话框信息中, COMSOL Multiphysics转化输入值为系数求解有限元模型. 框中包含每种模型类型: 2D问题中的求解域, 边界, 和点.2 单击d a选项卡. 此处其中系数为压力变化率乘以时间. 4-元素矩阵原本对角线为: 湿相压力, p w, 对湿相方程; 和非湿相压力, p nw, 对非湿相方程. 因为两相问题有储存项对非湿相和湿相的方程, 可以给矩阵填充合适的项. 修改矩阵如下所示:PW PNWCp+eps -Cp+eps-Cp+eps Cp+eps注意eps为一个非常小数.3 单击确定.网格生成1 从网格菜单, 选择自由网格参数.2 单击边界选项卡. 使用Ctrl键选择边界4和5. 在最大单元尺寸编辑区输入0.001.3 选择边界2, 并在最大单元尺寸编辑区输入0.005.4 单击重化网格按钮, 然后单击确定.求解1 从求解菜单, 打开求解器参数对话框. 在求解器列表选择稳态(如果还未选中).2 在时间编辑区, 输入0,0:0.01:0.1,0.1:0.1:1,1:170. 单击确定.3 在主工具条单击求解按钮.后处理和可视化按照如下步骤生成图2-43:1 从后处理菜单, 打开绘图参数对话框.2 在绘图类型区, 选择表面和箭头复选框. 在不同时间的解列表选择0.1.3 单击表面选贤卡, 并在表达式编辑区输入Se_nw, 其表示名称为nw的应用模式计算的有效饱和(Se).4 单击箭头选项卡. 在内建物理参数列表选择Darcy定律(w)>速度场. 转到箭头位置区. 在点数列, x点和y点编辑区输入25. 转到箭头参数区. 单击颜色按钮, 改变颜色为黑色, 并单击确定. 清除自动复选框, 并在比例刻度编辑区输入0.5.5 单击确定.继续如下步骤生成图2-44:1 从后处理菜单, 打开剖面绘图参数对话框.2 单击点选项卡. 在表达式编辑区, 输入Hc. 在坐标区, x编辑区输入0.03 0.03 0.03, y编辑区输入0.02 0.04 0.06.3 单击线设定按钮. 设定线样式为虚线, 然后单击确定.4 单击应用.5 在表达式编辑区, 输入p_nw/(rho_water*g_nw). 在坐标区, x编辑区输入0.03, 在 y编辑区输入0.0834.6 单击线设定按钮. 设定线样式为实线, 然后单击确定.7 单击通用选项卡. 选择保留目前的图复选框, 然后单击应用.继续下面步骤生成图2-45:1 仍然在通用页面, 从绘图于列表选择新图.2 单击点选项卡. 在表达式编辑区输入kr_nw.3 在坐标区, x编辑区输入0.03 0.03 0.03, y编辑区输入0.02 0.04 0.06.4 单击线设定按钮. 设定线样式为虚线, 然后单击确定.5 单击应用.6 在表达式编辑区, 输入kr_w.7 单击线设定按钮. 设定线样式为实线, 然后单击确定.8 单击通用选项卡. 选择保留目前的图复选框, 然后单击确定.两相流: 改变流体对这是两相流动例子的第二部分. 这里可以修改例子在第一部分创建的空气－水模型文件以仿真计算空气－油和油－水两相系统的流动特征, 参见参考文献1. 在空气－油系统中, 油为湿相而空气为非湿相. 在油－水系统中, 水为湿相.将空气－水系统转变为空气－油或者油－水系统需要改变流体属性和几个多孔介质参数. 为了获得不同流体对生成湿相流出率其仿真过程与空气－水系统类似, 如前面讨论的在入口处界面张力随压力比例增加模型库路径: Earth_Science_Module/Fluid_Flow/twophase_flow_air_oil模型库路径: Earth_Science_Module/Fluid_Flow/twophase_flow_oil_water 使用图形用户界面建模: 空气－油和油－水以下部分打开文件twophase_flow_air_water.mph并修改为空气－油和油－水系统. 下面给出的表和说明描述了如何改变参数以满足两相系统仿真的要求.模型导航视窗1 打开模型导航视窗并单击模型数据库选项卡.2 在模型数据库列表中选择地球科学模块>流体流动>twophase flow air water.3 单击确定.选项物理量—非湿从多物理场菜单, 选择达西定律(nw).边界条件—非湿相1 从选项菜单, 打开常数对话框. 按下表所示修改密度, 粘度和界面张力; 完成设定, 单击确定. 约束 空气-油 油-水rho_w 800[kg/m^3]1000[kg/m^3] etaw 0.00392[Pa*s]*SPH 0.001[Pa*s]*SPH rho_nw 1.28[kg/m^3]800[kg/m^3] etanw 0.0000181[Pa*s]*SPH 0.00392[Pa*s]*SPH sigma_ao 0.0259sigma_aw 0.06810.0681 sigma_ow 0.0364 2 选择选项>表达式>标量表达式.3 修改以下表达式; 设定完成, 单击确定.项 空气-油 油-水p_nw_top 0.2[m]*rho_water*g_nw*sigma_ao/sigma_aw0.2[m]*rho_water*g_nw* sigma_ow/sigma_aw p_w0 0.1[m]*rho_w*g_w*(sigma_ao/ sigma_aw)*(t/170[s])0.1[m]*rho_w*g_w*(sigma_ow/sigma_aw)*(t/170[s])4 选择选项>表达式>求解域表达式.5 在求解域选择列表, 选择2.6 改变如下表达式的值; 完成设定, 单击确定.名称 表达式 kap_s 0.94e-12theta_r 0.0072alpha 3.58 N 3.13651 从物理量菜单, 选择边界设定.2 在边界选择列表选择5.3 设置如下条件; 完成设定, 单击确定.边界条件 变量 表达式Pressure(for Air/Oil p 0 Hp_nw_t(t)*rho_water*g_nw*sigma_ao/sigma_aw求解后处理和可视化重复生成图2-44的步骤生成图2-46和图2-47.本文档由lhpvibration 整理 model)Pressure(for Oil/Watermodel)p 0 Hp_nw_t(t)*rho_water* g_nw*sigma_ow/sigma_aw1 从求解菜单, 打开求解器参数对话框.2 在求解器列表, 选择稳态(如果还未选中).3 转到时间步长区. 在时间编辑区输入0,0:0.01:0.1,0.1:0.1:1,1:170. 单击确定.4 在主工具条单击求解按钮.。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
了更为全面的掌握土壤参数空间变异性对地下水运动的影响 还必须进一步探索多随机输入情况下的
随机模型分析方法
参考文献
杨金忠 蔡树英 黄冠华 叶自桐 多孔介质中水分及溶质运移的随机理论 北京 科学出版社
责任编辑 吕斌秀
上面第二个例子中 动量方法求解摄动方程大约 次 计算耗时 而谱方法仅求解了约 次 计
算耗时 很明显 基于 展开的谱方法的计算效率高于动量方法
饱和 非饱和土壤中的水流 模拟区域为宽
长
的剖面 剖分为
个单元 上边
界给定入渗率 ? 下边界为隔水边界 左右边界的下部分为定水头边界 水头分别为
和
上部分为无流量边界 其它参数为 饱和水力传导度的几何平均
法分解 把压力水头表示为多项式 通过摄动方法得到一系列关于水头展开式的偏微分方程 用有限差分法
进行求解 获得了压力水头的随机描述 并计算其均值和方差 应用本文的随机模型研究了二维非饱和以及饱和
非饱和介质流动的实例 结果与动量方法的计算结果一致 而且计算效率高于传统的动量方法
关键词 随机模型 饱和 非饱和流
均水头相同 图
两者计算的水头标准差存在一点偏差 图
计算水头标准差 谱方法只需
要展开前 项就可获得与动量方法近似的结果 项以后的展开项对结果没有太大影响
将上边界的入渗率提高到 ? 同时增大 和
的变异性 取
计
算结果如图 所示 可以看到 计算水头标准差需要 展开前 项 在该例中 谱方法和动量方法的
计算结果也有一点偏差 通过模拟 可以发现水头标准差在离开下边界后快速增加 到达重力占优区后
项都需要求解 只需
求解关于
项的摄动方程 这些性质使得数值计算具有较高的效率
实例分析
非饱和土壤中的入渗 本文考虑剖面二维流问题 模拟矩形区域为
下边界为定水头边界
左右边界为隔水边界 输入参数如下
何平均
?
?
上边界为入渗边界 即饱和水力传导度的几
图 两种方法计算的水头均值和标准差
图 为分别用谱方法和动量方法计算的水头统计矩 谱方法得到的零阶平均水头与动量方法的平
当
时 可以得到零阶平均水流方程 当
可以得到相应阶数的水流方程 这些方程中包
含有随机变量序列
和
结合随机微分方程的扰动展开 可以得到如下确定性偏微分方
程
式中 分别为 和 求解这些方程后就可得到水头的随机描述 水头的统计矩也可以通过计算得到
式
式 都是确定性的偏微分方程 可以用现有的任何一种数值方法求解 本文采用有限
分别为初始水头 一类边界水头和二类边界流量
非饱和水力传导度和容水度是水头和空间位置的函数 显然需要通过数学模型来表征它们的函数
关系 常用的经验模型有
模型和
模型 已有的随机分析大多采用相
对简单的
模型
模型在数学形式上较
模型复杂 但是
它能较好地拟合田间实验数据 目前在国际国内应用广泛 本文采用
模型
?
式中
为饱和水力传导度
为孔隙分布参数 为模型参数 有
?
?
为有效饱和度
分别是有效体积含水率 饱和体积含水率以及残余体积含水率 其中
从式 可以看出
因此式 是水力传导度的一个通式
容水度
?
由式 可得
?
?
式中 是贮水率
本文将对数饱和水力传导度
和对数孔隙分布参数
假定为服从正态
分布的随机场 并且协方差函数已知 这一假定与文献 一致 为使问题简化 本文仅研究 和
展开 数值分析
中图分类号
文献标识码
由于土壤和含水层的地质沉积过程的随机性 含水层的水文地质特征以及土壤的结构 构造和组成
等特性具有明显的空间变异性 在孔隙介质空间变异性的影响下 地下水运动和溶质运移具有不确定
性 近 多年来 已有大量的随机模型应用于水流和溶质运移的不确定性研究
研究方法主要
包括
模拟和动量方法
饱和 非饱和流的随机偏微分方程
饱和 非饱和介质中的水流运动满足以下控制方程
初始条件和边界条件如下
收稿日期
基金项目 国家自然科学基金项目资助
作者简介 李少龙
男 湖北人 博士生 主要从事地下水环境方面的研究
式中
为
流速
为压力水头
为源汇项
为容水度
导度 为空间向量 方向以向上为正
为边界上的单位外法线向量
为水力传 和
和流问题
本文在采用
非饱和水分特征模型的基础上 联合运用 展开 混沌多项式
展开和摄动方法 对饱和 非饱和流问题进行随机数值分析 首先 将土壤特性参数用 展开 同时
将水头表示为多项式展开 然后通过摄动方法得到一系列关于水头展开的偏微分方程 并用有限差分法
进行求解 水头的统计矩可通过它的随机描述直接计算得到
的次数与区域剖分的节点数目 成正比 而谱方法与 展开项数 成正比 较 明显要小很多
理论分析和计算结果都表明 谱方法比动量方法有更高的计算效率 因而也具有较好的实用性
与已有的一些随机模型不同 本文采用了
非饱和水分特征模型 本文的随机
模型考虑了饱和水力传导度 和孔隙分布参数 的空间变异性 但是忽略了模型参数 的变异性 为
头计算结果的影响不大 从图 和图 还看到 谱方法和动量方法的计算结果有一点差异 可能是由于
求解方程时的数值误差引起的 然而偏差很小 表明本文的随机模型是正确的
为了计算水头方差 动量方法在处理两随机场的情况下 必须在所有的节点上求解关于
的摄动方程 求解次数约为节点数目的 倍 而基于 展开的
谱方法 只需求解少量的关于水头展开系数的摄动方程 求解次数约为 展开截断项 的 倍 在
不相关的情况 当将 和
作为随机场时
就成为随机偏微分方程 它的解不再是
确定性的值 而是具有一定概率分布的随机函数
随机场的数学描述
随机场的
展开 空间随机函数
谱展开的基础上进行的 有界 对称 正定的
数构成了一个完备集 随机场的 展开为
的 展开是在
的协方差函数
具有相互正交的特征函数 这些特征函
式中
为随机空间函数的均值
基本保持不变 接近上边界时又增加 这一计算结果与动量方法的结果 是一致的
图 两种方法计算水头统计矩的比较
从这两个算例可以看到 随着输入方差和入渗率的加大 要求 展开 和
的项数也在同
时增加 当展开达到一定项数后 后面的展开项对结果的影响是很小的 这样 我们就可以用较少的展
开项获得较好的近似结果 本文还用谱方法将水头计算到二阶 结果表明水头的二阶展开项对平均水
水力坡度明显增大 最大的水头标准差也出现在弱透水层处
结论
本文联合运用
展开 混沌多项式展开和摄动方法 对饱和 非饱和流问题进行随机
图 饱和 非饱和介质中的平均流场和水头标准差
图 含弱透水层的饱和 非饱和介质中的平均流场和水头标准差
数值分析 将土壤特性参数作为随机场 采用 展开 水头用混沌多项式展开 通过摄动方法 推导
为零均值正交随机变量序列
和 分别为协方差
的特征函数和对应的特征值
与已有的随机模型的处理一样 本文假定 和
是二阶平稳随机场 其协方差函数为离散指
数型 在二维情况下有
式中 空间向量
为随机空间函数的方差 和 分别是两个方向上的相
关长度
对于矩形区域
可以获得特征函数和特征值的解析解
将和
表示为均值和扰动之和
对于 和
得到了一系列关于水头展开的确定性偏微分方程 最后用有限差分法求解 得到水头的随机描述后 利
用高斯正交随机变量的性质 统计矩可以从水头的随机描述中直接得到
通过非饱和流以及饱和 非饱和流的实例分析 并将模拟结果与动量方法进行比较 本文阐述了基
于 展开的随机模型的正确性和应用性 将统计矩计算到随机输入方差的一阶 动量方法求解方程
由于它
考虑到 和
的 展开 有
当式 和式 中的确定性系数计算出来后 就可以获得
淆的情况下 将
省略
的所有概率特征 后面在不引起混
摄动方程
为了数学上处理方便 利用变换
由式 可得
?
水头
的变异性依赖于 和
将这些函数展开为级数形式
同时将这些级数展开和式 代入式 式
的变异性 而
和
合并同阶项后可得
的变异性又依赖于 利用式
?
流动达到稳定状态时的模拟结果如图 所示 图 中 实线为总水头的等势线 带箭头的实线为
流线 在非饱和区 水流是垂直下渗的 在饱和区 由于左边界的总水头高于右边界的总水头 产生了从
左向右的平均流动 水头标准差如图 所示 在饱和区 水头标准差在边界处为 离开边界后逐渐
增加 在区域中心附近达到最大值 在非饱和区 由水面向上水头标准差迅速增大 到达重力占优区后
年月 文章编号
水利学报
第 卷第期
基于
模型的饱和 非饱和 介质流动随机数值分析
李少龙 杨金忠 蔡树英
武汉大学 水资源与水电工程科学国家重点实验室 湖北 武汉
摘要 基于
非饱和水分特征模型 联合运用
展开法 混沌多项式展开以及摄动
方法 对饱和 非饱和流问题进行随机数值分析 将土壤特性参数假定为协方差已知的随机函数 并按
基本保持不变 最后在上边界处取得最大值 在饱和 非饱和介质中 流动矩的空间变异性较单独的饱
和区或非饱和区都要复杂得多 本文的模拟结果和文献 是一致的
此外 本文还考察了介质非平稳特性对水头统计矩的影响 在上例的基础上 嵌入饱和水力传导度
的几何平均
? 的弱透水层 该土层厚宽其中心位于地面以下其余参数
不变 模拟结果如图 所示 可以看到介质的非平稳特性加剧了流动矩的空间变异性 弱透水层处的