一种膨胀性土壤水分运动过程模拟方法与流程

本发明涉及一种膨胀性土壤水分运动过程模拟方法，属于膨胀性土壤水分运动模拟技术领域。

背景技术：

当前，在非膨胀性土壤水分运动过程研究方面，国内外学者已经开展了大量室内外试验，并建立相关数学模型。然而，这些针对刚性土壤的研究成果不适用于膨胀性土壤水分运动过程描述。膨胀性土壤吸水膨胀，失水收缩，湿胀干缩过程极大地影响着土壤水分运动过程。世界上膨胀性土壤分布广泛，我国也是世界膨胀性土壤分布最为广泛的国家之一。膨胀性土壤变形不仅影响着土壤水文过程，也会给工程建设带来一系列问题，已成为工程地质学、水文学和土力学等学界关注的重点。土壤膨胀变形主要与初始含水量和上覆荷载等有关，膨胀力和膨胀变形随土壤增湿程度增加而增加。土壤吸水膨胀变形主要受膨胀力和自重应力的作用，其中膨胀力随土壤含水量变化而变化，自重应力随土壤深度变化而变化。随土壤深度的增加，土壤受力变形特征发生改变，进而引起土壤饱和含水量、饱和导水系数、饱和比容积等参数变化，变化的土壤水分运动参数改变了土壤剖面水分运动过程。然而，当前针对膨胀性土壤变形对土壤水分运动过程的影响研究仍处于探索阶段。

技术实现要素：

本发明所要解决的技术问题是：提供一种膨胀性土壤水分运动过程模拟方法，模拟膨胀性土壤水分运动过程，有助于完善土壤水分运动理论。

本发明为解决上述技术问题采用以下技术方案：

一种膨胀性土壤水分运动过程模拟方法，包括如下步骤：

步骤1，基于土壤变形受力分析，得到土壤变形过程，采用质量守恒定理，模拟膨胀性土壤水分运动过程，得到：

基质区水分运动过程方程为：

上边界：

下边界：

裂隙优先流区水分运动过程方程为：

上边界：

下边界：

其中，θ为土壤体积含水量，cm³/cm³；θ1为时段初的土壤含水量，cm³/cm³；e为孔隙度，cm³/cm³；e1为时段初的孔隙度，cm³/cm³；ke(ψ)为膨胀性土壤非饱和导水系数，cm/min；ψ为土壤水吸力，cm；wr为植物根系吸水量，cm³；we为两流区水量交换量，cm³；φ为土壤水势，cm；q为水分通量，cm/min；上述字母的下标f、j分别表示裂隙优先流区、基质区；t为时间，min；z表示坐标轴z轴方向；wf为裂隙优先流区面积比例，wj为基质区面积比例；

步骤2，采用改进的vangenuchten模型计算膨胀性土壤裂隙优先流区、基质区非饱和导水系数，非饱和导水系数计算公式为：

其中，ke(ψ)为膨胀性土壤非饱和导水系数，cm/min；ksh(e)为深度为h时膨胀性土壤饱和导水系数，cm/min；se为饱和度；m和n为参数，m＝1-1/n；

步骤3，分别计算膨胀性土壤裂隙优先流区、基质区的面积比例，其计算公式为：

wf＝dew

wj＝1-dew

其中，ew为由土壤吸水膨胀变形导致的孔隙度变化量，cm³/cm³，当土壤饱和时，wf等于0。

作为本发明的一种优选方案，步骤2所述深度为h时膨胀性土壤饱和导水系数ksh(e)，计算公式为：

其中，e0、k0分别为零压强下的孔隙度、饱和导水系数，cm³/cm³、cm/min；eh为深度为h时土壤饱和孔隙度，cm³/cm³；m′为与土壤质地有关的参数；e1为时段初的土壤孔隙度，cm³/cm³；ρd为土壤密度，g/cm³；a3为参数；α3为土壤膨胀特征曲线斜率；u为质量含水量，g/g；a和b均为参数；γ为土壤湿比重，n/cm³；h为土壤深度，cm。

作为本发明的一种优选方案，步骤2所述饱和度se，计算公式为：

其中，θ为土壤体积含水量，cm³/cm³；θr为残余含水量，cm³/cm³；θsh为深度为h时土壤饱和含水量，cm³/cm³；α、m、n均为参数，m＝1-1/n；ψ为土壤水吸力，cm。

作为本发明的一种优选方案，步骤3所述由土壤吸水膨胀导致的孔隙度变化量ew，表示为：

其中，e为孔隙度，cm³/cm³；e1为时段初的土壤孔隙度，cm³/cm³；ρsw为由土壤吸水膨胀导致的容重变化量，g/cm³；ρd为土壤密度，g/cm³。

作为本发明的一种优选方案，所述土壤深度为h时土壤饱和含水量θsh等于深度为h时土壤饱和孔隙度eh，cm³/cm³。

本发明采用以上技术方案与现有技术相比，具有以下技术效果：

1、本发明基于土壤变形受力分析，得到土壤变形过程，采用质量守恒定理，提出了膨胀性土壤水分运动过程模拟方法；并提出了膨胀性土壤水分运动参数计算方法，其中膨胀性土壤水分运动参数包括膨胀性土壤非饱和导水系数、比容积等。

2、本发明方法有助于完善土壤水分运动理论，对膨胀性土壤分布地区的水分管理与调控具有一定的指导价值。

附图说明

图1是本发明土壤膨胀变形示意图。

图2是本发明变形条件下计算单元土壤水量平衡示意图。

具体实施方式

下面详细描述本发明的实施方式，所述实施方式的示例在附图中示出。下面通过参考附图描述的实施方式是示例性的，仅用于解释本发明，而不能解释为对本发明的限制。

本发明提出了基于质量守恒原理，提出了膨胀性土壤水分运动过程计算方法，并考虑土壤膨胀性的土壤水分运动参数计算方法，相关研究有助于完善土壤水分运动理论，对膨胀性土壤水分管理与调控具有一定的指导价值。

1、膨胀性土壤水分运动过程计算步骤

如图1所示，受土壤膨胀力f和自重应力g影响，土壤发生了变形。p为土壤膨胀力与自重应力间的合力，z1和z2分别为微元变形前后垂向大小，vs和vs’分别为微元变形前后体积大小。为了便于分析，假设：(1)土壤变形只引起土壤孔隙度变化；(2)土壤变形是弹性变形，无滞后性。取一块长、宽和高为δx、δy和δz微小单元(如图2所示，u为微元)，其中δx×δy＝1×1。设研究土体单元内任意时刻t，在δt时间土体内含水量变化量由三部分构成，分别为：

单元体(δxδyδz)含水量变化：

δxδyδz(θ2-θ1)(1)

单元体变化量(δxδyδe)含水量变化：

δxδyδeθ1+δxδyδe(θ2-θ1)(2)

式中：e1为变形前土壤孔隙度(时段初)，cm³/cm³；e2为变形后土壤孔隙度(时段末)，cm³/cm³；δe＝e2-e1，为土壤孔隙度变化量，根据假设1，等于土壤体积变化量，cm³/cm³；θ1为时段初土壤含水量，cm³/cm³，θ2为时段末土壤含水量，cm³/cm³。

由于δx·δy·δe·(θ2-θ1)是高阶小量，可以忽略，则整个土体内在δt时间内含水量变化量可以表示为：

δx·δy·δz·(θ2-θ1)+δx·δy·δe·θ1(3)

一维条件下，在单位时间δt内，z轴方向由水势引起的水分通量可以表示为：

φ＝h0+z+ψ(6)

在单位时间δt内，土体内由植物根系吸水引起的水分通量可以表示为：

wr＝-sδxδyδzδt(7)

在单位时间δt内，土体内由两流区水量交换引起的水分通量可以表示为：

we＝-iδxδyδzδt(8)

则由质量守恒定律有：

考虑到两流区水分运动过程，其中基质区面积比例为wj，裂隙优先流区面积比例为wf，则方程进一步可以改为(基质区)：

上边界：

下边界：

优先流区为：

上边界：

下边界：

式中：x、y和z表示坐标轴；f为膨胀力，n；g为自重应力，n；θ为土壤体积含水量，cm³/cm³；ke(ψ)为土壤非饱和导水系数，cm/min；ψ为土壤水吸力，cm；φ为土壤水势，cm；h0为静水压力(当静水压力较小时，可以忽略)，cm；e为孔隙度，cm³/cm³；t为时间，min；e1为时段初的孔隙度，cm³/cm³；θ1为时段初的土壤含水量，cm³/cm³；vs为土壤体积，cm³；wr为植物根系吸水量，cm³；s为根系吸水强度，min^-1；we为两流区水量交换量，cm³；i为两流区水量交换强度，min^-1；q为水分通量，cm/min；下标f和j分别表示优先流区和基质区；wj为基质区面积比例；wf为裂隙优先流区面积比例，两者间面积比例可以采用膨胀特征曲线计算。

2、膨胀性土壤水分运动参数计算：优先流区和基质流区模型参数数值大小不同，因此应该需要单独计算。

1基质区模型水分运动参数计算：

1)膨胀力作用导致的孔隙度变化量

假设土壤膨胀变形是由土壤孔隙度的变化引起的，则当土壤吸水变形时，土壤膨胀力引起的孔隙度变化量可以表示为：

式中：ew为由土壤吸水膨胀导致的孔隙度变化量，cm³/cm³；ρd为土壤密度，g/cm³；e1为时段初孔隙度，cm³/cm³；ρsw为由土壤吸水膨胀导致的容重变化量，g/cm³。

土壤吸水后在膨胀力的作用下土壤自由变形，土壤变形各项同性，土壤膨胀变形量是土壤含水量的函数，可以采用三直线模型计算：

式中：ν为比容积，是土壤容重的倒数，cm³/g；u为质量含水量，g/g；α1、α2、α3为土壤膨胀特征曲线斜率；ua、ub、us分别为拐点处质量含水量，g/g；a1、a2、a3为参数。

于是，可得

i为三直线模型各分段(i＝1,2,3)。

2)自重应力作用导致的孔隙度变化量

同理，土壤自重应力引起的孔隙度变化量可以表示为：

式中：ep为由土壤自重应力导致的孔隙度变化量，cm³/cm³；ρsp为由土壤自重应力导致的容重变化量，cm³/g；可以采用土壤应力——应变关系曲线计算：

ρsp＝a+bln(γh)(20)

式中：γ为土壤湿比重，n/cm³；h为土壤深度，cm；a和b均为参数；其它符号意义同前。

3)合力作用导致的孔隙度变化量

土壤膨胀变形是膨胀力和自重应力共同作用的结果，则合力导致的土壤空隙变化量可以表示为：

de＝dew+dep＝δew+δep(21)

4)土壤饱和导水系数计算模型：受自重应力和膨胀力影响，土壤孔隙度随深度的变化而变化，导致土壤饱和导水系数随深度的变化而变化。土壤饱和导水系数随深度变化关系采用改进的lambe模型计算：

其中，e0、k0用恒体积法测定，为零压强下对应的孔隙度和饱和导水系数，cm³/cm³、cm/min；eh为深度为h时土壤饱和孔隙度，cm³/cm³；m′为与土壤质地有关的参数；ρd为土壤密度，g/cm³；a3为参数；α3为土壤膨胀特征曲线斜率；u为质量含水量，g/g；a和b均为参数；γ为土壤湿比重，n/cm³；h为土壤深度，cm。

5)土壤饱和含水量计算：当土壤达到饱和时，土壤饱和含水量等于孔隙度，则任一深度饱和含水量可以表示为：

θsh＝eh(24)

式中：θsh为土壤深度为h时土壤饱和含水量，cm³/cm³。

6)水分特征曲线计算模型：土壤水吸力ψ与土壤含水量θ间关系采用vangenuchten模型描述：

其中，θ为土壤体积含水量，cm³/cm³；θr为残余含水量，cm³/cm³(θr也随孔隙度的变化而变化，但土壤干燥时，土壤变形较小。同时，同一种土壤残余含水量相差较小，因此不同深度残余含水量可以互相代替)；θsh为深度为h时土壤饱和含水量，cm³/cm³；α、m、n均为参数(α、m、n随土壤孔隙度的变化而变化，但同一种土壤α、m、n相差较小，因此不同深度处可以互相代替)，m＝1-1/n；ψ为土壤水吸力，cm。

7)土壤非饱和导水系数计算模型：非饱和导水系数采用改进的vangenuchten模型计算：

式中：ke(ψ)为土壤非饱和导水系数，cm/min。

2裂隙优先流区模型水分运动参数计算：

由于优先流区和基质流区土壤变形受力特征一样，均由土壤膨胀力和自重应力作用引起孔隙度变化，进而导致水分运动参数变化。因此，裂隙优先流区模型水分运动参数计算方法同基质区。

3裂隙优先流区和基质区面积比例计算：

膨胀性土壤湿涨干缩，当土壤饱和时，土壤膨胀力作用导致的孔隙度ew达到最大，土壤裂隙闭合，则裂隙优先流面积比例wf为零，反之，当土壤未达到饱和时，由于裂隙的产生只与土壤膨胀力有关，土壤膨胀力作用导致的土壤变形各向同性(各个方向变形相等)，则水平方向裂隙优先流所占面积比例可以表示为：

wf＝dew(28)

wj＝1-dew(29)

式中：当土壤饱和时，wf等于0；ew为由土壤吸水膨胀变形导致的孔隙度变化量，cm³/cm³。

以上实施例仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明保护范围之内。