基于分形理论的非饱和土多相流相对渗透系数模型

曹 帅，宁金成，杨庆国

(1.河南交通职业技术学院 公路学院，郑州 450015； 2.重庆交通大学 土木建筑学院，重庆 400074)

1 研究背景

非饱和土的渗透破坏与暴雨诱发滑坡、高库大坝水力破坏、洪涝管涌破坏等工程问题[1-2]密切相关，而相对渗透系数则是表征非饱和土渗流规律最为关键的参数之一。因此，合理确定非饱和土的相对渗透系数对于准确把握非饱和土的渗流力学行为，进而为工程设计及灾害防治提供依据，具有重要的理论价值与工程意义。

众所周知，非饱和土的相对渗透系数受土体细观孔隙分布特征(如孔隙的孔径大小、不同孔径孔隙的体积分布、孔隙形状等)及饱和度影响极大[3-5]。然而，由于非饱和土在较高饱和度时，相对渗透系数对饱和度的变化非常敏感，饱和度的微小变化会使相对渗透系数发生相当大的改变;而在较低饱和度时，非饱和土的相对渗透系数很小，造成该系数测量难度大、费时长，使得直接通过物理测定获得该参数的准确值变得极为困难。因此，提出各种预测模型来预测非饱和土的相对渗透系数成为较为实际的方法。不少学者对该问题进行了深入研究[6-15]，如Burdine[6]和Mualem[7]通过提出理想的土壤孔隙模型建立了非饱和土液相渗流模型，该模型依赖于土水特征曲线试验数据来计算液相的相对渗透系数。Van Genuchten[8]通过提出经验的土水特征曲线模型,并分别与Mualem[6]和Burdine[7]的非饱和多孔介质广义模型相结合，得到了三参数非饱和土的相对渗透系数模型，该模型一定程度上反映了土壤持水与孔隙分布规律的关系，但模型参数为经验参数，无具体物理含义。Yang和Mohanty[9]结合Kosugi[10]的土水特征曲线函数，推导了多孔介质非湿润相相对渗透率的一般表达式，并使用该广义公式研究了Burdine[6]、Mualem[7]、Alexander等[11]、Skaggs[12]模型对非湿润相渗透率的预测精度。

土壤孔隙形状的不规则性和连通性使得准确估算非饱和土渗透系数仍是一个具有挑战性的课题，尤其在土壤多相流运移模型方面。为研究土壤多相流运移模型，本文将土壤孔隙简化为一簇服从分形定律的迂曲毛细管，通过建立液相与气相在多孔介质中的输运特征模型，结合Young-Laplace方程与Buckingham-Darcy定律得到了非饱和土的液相与气相的相对渗透系数模型，该模型仅包含2个参数，模型参数物理意义明确且均可通过试验测得。最后，通过8组试验数据验证了本文模型的合理性与适用性。

2 Young-Laplace方程

毛细管模型为理解空气-水-土交界面上的液体压力与蒸汽压之间的关系提供了物理解释。一般情况下，可用圆柱面、椭球面和球面等简单的几何形状来近似表示空气-水-土交界面。Young-Laplace方程[16]用双曲率模型(椭球面如图1所示)建立了基质吸力(ua-uw)与交界面几何形状关系的公式，其表达式为

图1 椭球体状气-液-固交界面Fig.1 Ellipsoid gas-liquid-solid interface

(1)

式中：ua与uw分别为气相和液相的压力；Ts为水表面张力；θ为接触角；r1、r2分别为椭圆的半长轴和半短轴。

为方便计算且不失一般性，本文对毛细管模型做出如下假设：

(1)土体孔隙可简化为横截面为圆形的迂曲毛细管。

(2)对于不同孔径的毛细管，接触角θ均相等且为常数。

(3)不考虑温度变化对水表面张力的影响，即水表面张力Ts为常数。

由r1=r2=r，ua-uw=rwgh，代入式(1)可得

(2)

式中：h为基质水头；ρw为液相密度；g为重力加速度。

3 分形尺度分析及模型建立

3.1 分形尺度分析

众所周知，非饱和土土体中分布着大量复杂的孔隙，其分布及自身尺寸则直接影响了土体的水力特性。为描述该特性，本文采用一簇具有不同半径的毛细管来代替这些孔隙，如图2所示。

图2 分形多孔介质孔隙简化模型Fig.2 Simplified pore modelof fractal porous media

在土中取一个具有代表性的圆柱体单元(半径为R0，长度为L)，假设土体孔隙可用一簇具有不同半径的毛细管代替，其半径r在一定范围内，即rminr的毛细管的累积数量近似服从以下分形定律[17-18]，即

(3)

式中：ε为分形标度，此处指土壤孔隙半径；Df为分形维数。

Yu[18]提出分形维数Df可由式(4)计算。

(4)

式中：dE为欧式空间维度数，在二维欧式空间里dE=2，在三维欧式空间里dE=3；φ为孔隙率。在二维欧式空间里，多孔介质的孔隙分形维数一般有1

将式(4)对r微分，可得到毛细管半径从r到r+dr的毛细管累积面积，即

(5)

由式(5)可知代表毛细管的累积数量随着孔径的增大而减少。

假设每根毛细管中的水都处于可流动状态，根据热力学局部平衡理论，在具有一定含水量的土体中，当土中水运动时(一般分为干燥过程和湿润过程)，存在一个临界半径R(R为毛细管恰好饱和状态时的临界孔隙半径)，使得R为控制着毛细管饱和状态和非饱和状态的临界半径，即毛细管半径r=R时恰好为饱和状态，其水头高度恰好为L，则当毛细管半径为rmin≤r≤R时，毛细管为饱和状态；当R

3.2 非饱和土多相流相对渗透系数模型建立

非饱和多孔介质中的多相流问题可以分为多个单相流问题。在具有代表性的圆柱体单元中，某一特定含水量下的非饱和多孔介质的孔隙气压力ua和孔隙水压力uw为常数。因此，根据Bear[19]提出的Navier-Stokes方程特解，单根孔隙的气相和液相体积流率可分别表示为：

(6)

(7)

式中：μa和μw分别为气相和液相的流体动力黏度；ua和uw分别为非饱和土多孔介质的孔隙气压力和孔隙水压力；ra和rw分别为液相和气相密度；z为垂直坐标；g为重力加速度。

结合式(5)和式(7)可得具有代表性的圆柱体单元横截面的液相总体积流率Qw，即

(8)

结合式(2)和式(8)，液相的总体积流率可用基质水头表示为

(9)

由Buckingham-Darcy定律[20]可知，非饱和多孔介质的液相总流率还可表示为

(10)

式中kw为非饱和多孔介质液相的渗透率。结合式(9)与式(10)可得

(11)

同理可得气相(R～rmax)的渗透率ka可表示为

(12)

由式(11)和式(12)可知，当h=hmax时，土壤为干燥的，此时土壤孔隙被气相填满；当h=hmin时，土壤孔隙被水充满，土壤是饱和的。因此，饱和土壤的渗透率ksat可表示为

(13)

结合式(11)—式(13)可得非饱和土的液相和气相的相对渗透系数，即：

(14)

(15)

式中：kr,w和kr,a分别为非饱和土液相和气相的相对渗透系数；hmin和hmax分别为最小和最大基质水头。土壤的饱和度Se可由饱和土定义求得，即

式中:Vw为土中孔隙水的体积；Vv为孔隙总体积。

因此，联立式(2)与式(16)可得饱和度与基质水头关系的表达式，即

结合式(14)、式(15)和式(17)可得非饱和多孔介质液相与气相的相对渗透系数与饱和度关系的表达式，即：

式中ξ=hmax/hmin。

式(18)和式(19)即为本文提出的非饱和土多相流(气相与液相)相对渗透系数模型，该模型仅包含Df和ξ2个参数，且Df和ξ可分别通过计盒法和轴平移技术测得。

4 模型验证与讨论

本部分采用不同土类的相对渗透系数试验数据验证所提出模型的可行性与合理性，这些土包括Sable de Riviere土、Gllat 砂质壤土和Rubicon 砂质壤土[21]；Grimsel Test Site土[22]；Mixed 砂和Oakley 砂[23]。采用MatLab全局搜索功能拟合试验数据并将得到的模型参数列于表1。

表1 不同土类的相对渗透系数模型参数Table 1 Parameters of relative permeability coefficientmodel for different soil types

图3为Sable de Riviere土、Gllat 砂质壤土和Rubicon 砂质壤土液相的相对渗透系数试验数据与本文模型模拟曲线的对比。由图3可知，本文模型模拟曲线与试验数据吻合良好。其中，对于Sable de Riviere土的液相相对渗透系数试验数据，当土体饱和时土壤的相对渗透系数为1，即土壤渗透率为饱和渗透率。当土壤的饱和度<1时，土壤的相对渗透率也随之不断减小，直至饱和度为0时，相对渗透系数也变为0，此时土壤气相的渗透率达到最大值。Gllat 砂质壤土和Rubicon 砂质壤土的相对渗透系数试验数据也表现出了类似的特征。

图3 不同非饱和土液相的相对渗透率试验数据与本文模型模拟曲线对比Fig.3 Comparison between test data of liquid-phaserelative permeability coefficient of different unsaturatedsoils and the simulated curve by the present model

图4为Grimsel Test Site土、Mixed砂和Oakley砂的气相相对渗透系数与本文模型模拟曲线对比。

图4 不同非饱和土气相的相对渗透率试验数据与本文模型模拟曲线对比Fig.4 Comparison between test data of gas-phaserelative permeability of different unsaturated soils andthe simulated curve by the present model

由图4可知本文模型可较好地反映非饱和多孔介质的气相渗流规律。对于Grimsel Test Site土气相的相对试验数据，当土壤完全干燥时，土壤的气相渗透率最大且其值随着饱和度的增大而不断减小。值得注意的是，当土壤的气相传导能力接近0时，土壤的饱和度并未到达1，这主要是因为当土壤含水量达到一定值以后会将土壤孔隙堵住,从而限制了气相的输运。Mixed砂和Oakley砂的气相相对渗透系数试验数据也表现出了与Grimsel Test Site土类似的渗流特征。

孔隙分布规律对非饱和土的渗流规律有着不可忽视的影响。图5表明了分形维数Df对非饱和土多相渗流的影响。由图5可知，在同一饱和度下分形维数越大，土壤液相的相对渗透系数越小，气相的相对渗透系数越大。这是因为分形维数越大表明土壤的较小孔隙占比越大，由于孔隙水表面张力的存在，土壤的导水能力越差，而气相的输运能力相对就越强。

图5 孔隙尺寸分布对非饱和土相对渗透系数的影响Fig.5 Influence of pore size distribution on relativepermeability coefficient of unsaturated soil

由上述分析可知，本文提出的多相相对渗透系数模型可较好地模拟非饱和土的多相渗流规律。由于本文在重构土壤孔隙系统时是将土壤孔隙简化为一簇光滑的毛细管，因此忽略了孔隙的粗糙度、连通性及不规则性对土壤渗流规律的影响，造成本文模型在预测土壤的相对渗透系数时仍存在一些误差，这些因素对土壤的相对渗透系数的影响也是笔者进一步研究的目标。

5 结 论

本文在可用毛细管系统近似代替土孔隙的假设上，基于Young-Laplace方程和分形理论，探讨了非饱和土的多相渗流特征，并得到了以下结论：

(1)将土壤孔隙简化为等截面迂曲毛细管，假设土壤孔隙尺寸分布服从分形定律并重构了土壤孔隙系统模型。

(2)结合Young-Laplace方程与Buckingham- Darcy定律得到了非饱和土的液相与气相的相对渗透系数模型，该模型仅包含2个参数且模型参数物理意义明确。采用6组不同土类的试验数据对本文提出的模型进行了验证，结果显示本文提出的多相相对渗透系数模型与试验数据吻合良好，表明了本文模型可较好地模拟非饱和土的多相相对流渗透系数。

(3)讨论了分形维数Df对非饱和土液相与气相的渗流特性的影响。由分析可知同一饱和度下，分形维数Df越大，土壤液相的相对渗透系数越小，气相的相对渗透系数越大。