水位下降时岸坡渗流的多相流模型

水位下降时岸坡渗流的多相流模型

1饱和-非饱和渗流场如果岸壁外水位降低，则库水位变化区域内的土壤从饱和状态变为非饱和状态。水逐渐补充的初始孔空间逐渐被空气补充。在这种“气驱水”的渗透过程中，孔压的变化对岸壁的稳定起着决定性作用。因此，在研究水库水位下降时产生的滑动滑动机制时，研究岸壁的非稳定、饱和和渗透率非常重要。目前,岸坡在库水位下降情况下的边坡稳定性评价及设计大多是建立在饱和渗流的基础上,但饱和渗流无法考虑基质吸力对边坡稳定的贡献,即使在饱和-非饱和渗流场的研究中,大都假设孔隙气压力等于大气压力,非饱和土的基质吸力简化为负的孔隙水压力,在非饱和土的抗剪强度分析中只考虑孔隙水压力的变化,然而在水位下降时,饱和区的孔隙气压力的消散不可能瞬间完成,孔隙气压力的变化必将影响孔隙水压力及基质吸力的分布。本文根据水、空气的质量守恒定律和达西定律,结合多相流理论建立水-气二相流模型,对库水位下降时的非稳定渗流场进行分析,探索孔隙气压力的变化对非稳定渗流场的影响。2高水势电场数学模型近年来,多孔介质中多相流理论的研究和应用越来越广泛,例如石油工业、环境工程以及地下水的保护和水资源的管理等,建立了许多模拟多相流运动的数学模型。然而最常见的多相流系统是水-气二相流系统,岩土工程及水利工程的饱和非饱和渗流过程大多可看作为恒温的水-气二相流过程。2.1组分及含量测定相是指具有与相邻物质不同的性质及明确分界面的物质。非饱和土应认为是由4相组成,即固相(土粒)、气相、液相和水-气分界面(收缩膜)。假设固体骨架不变形,处于平衡状态的气相(g)和液相(l)在多孔介质孔隙中流动,满足达西定律,不计收缩膜的质量,考虑毛细作用力。组分是指各个相的化学组成成分,构成液相和气相的组分均为水(w)和空气(a)。液相包含液态水和溶解其中的空气,常见的岩土工程一般可视为恒温系统,空气的溶解性很小,符合Henry定律;气相包含干空气和水蒸气,均可看作理想气体,气相压力满足Dalton气体分压力定律,并假设水蒸气的分压力等于饱和水蒸汽压力。相态是指控制体积内水相、气相的存在与否情况,水-气二相流系统的相态结构转换图见图1所示。2.2控制方程式2.2.1单元内部法向量任意体积内空气和水所满足的质量守恒方程为水-气二相流模型的基本控制方程:式中:Vn为任意单元体;Γn为单元体的封闭边界;上标κ表示组分水(w)或空气(a);Fκ为κ的单位流量(kg/s),流入为正;n为dΓn上的法向量余弦,指向单元体内为正;qκ为κ的源汇项(kg/s);Mκ为κ在dVn上的质量(kg)。在多孔介质力学中,液相或气相的流动满足达西定律,水-气二相流中的达西定律的表达式为式中:下标β为液相l或气相g;uβ为β相的流速(m/s);k为固有渗透系数张量(m2);krβ为β相的相对渗透系数;µβ为β相的动力黏滞系数,(N⋅s/m2);pβ为β相的孔隙压力(Pa);g为重力加速度矢量(N/kg)。引入质量分数Xβκ来表示组分κ在β相中所占的质量分数,则Fκ=∑βXβκFβ,Mκ=φ∑βSβρβXβκ,其中,Fβ=ρβuβ为β相的单位流量(kg/s);φ为孔隙率,Sβ为β相的饱和度,ρβ为β相的密度(kg/m3)。2.2.2-饱和度关系曲线(1)毛细压力方程:pβ=pg+pcβ,其中,pg为气相压力;pcβ为β相毛细压力,以气相压力为参照,水相毛细压力为水相压力与气相压力的差值,气相毛细压力始终为0。(2)饱和度方程:∑βSβ=1。(3)质量分数方程:∑κXβκ=1。(4)毛细压力-饱和度关系曲线近年来,建立了许多表征毛细压力-饱和度关系曲线的模型,最常用的是vanGenuchten模型。式中:,其中,pcl为液相毛细压力;Se为有效水饱和度;Srw为剩余水水和度。p0、n和m为模型的拟合参数,需由试验数据拟合得到,其中p0为进气值。(5)相对渗透率-饱和度关系曲线表征相对渗透率-饱和度关系曲线的模型中最常用的是vanGenuchten-Mualem模型。式中:krl和krg分别为液相和气相的相对渗透系数;λ=n-1为拟合参数;Srg为剩余气饱和度;的定义为。2.3含nil个单元的计算域由上述可见,恒温下的水-气二相流系统中的模型变量较多,根据Gibbsian’s相律,自由度为:F=C+2-P=2,等于任意体积内的质量守恒方程数,其中C和P分别为组分数和相数。对于包含NEL个子单元的计算域,总自由度等于2×NEL。将所有变量分为:主要变量——可直接通过方程求解得到,次要变量——通过主要变量得到。因单元内的相态可能发生改变,每次迭代结果需根据相态准则判断,选择不同的主要变量来表征各种相态。表1列出所有相态所对应的主要变量及各相出现准则。在各种相态下,孔隙气压力均可作为主要变量,水-气二相流状态时,饱和度可作为主要变量,仅液相或气相时,质量分数Xla或Xga可作为主要变量。仅有液相时,气相出现的准则1:水蒸汽分压力与溶解的空气压力之和大于气相压力;仅有气相时,液相出现的准则2:水蒸汽分压力大于饱和水蒸汽压力。3模型的数值解和边界条件的处理3.1数值解算方法3.1.1单元表面anm的离散在水-气二相流模型的求解中,首先需对基本控制方程(1)这一高阶强非线性的复杂方程组进行空间离散,在目前多相流模型数值求解的研究中,积分有限差分法(IFDM)是最有效并广泛采用的方法。IFDM的核心是计算域中离散的子单元的性质由形心点代表,分别对各个单元的质量平衡方程的积分形式进行空间离散。引入适当的体积平均值,可以得到式中:Mnκ、qnκ分别为Mκ、qκ在nV上的平均值。Γn上的面积分可近似为其所包含的各个表面Anm的面积分的平均值之和。式中:m为与单元n相邻的所有单元;Anm为单元n和m相邻的交界面;Fnmκ为Fκ在面Anm上的内法线方向的平均值,离散后的表达式为式中:下标nm为单元n和m之间交界面的上游权重平均值;Dnm为节点n和m之间的距离;gnm为重力加速度在mn连线方向的分量。最后,可得到一组关于时间的一阶微分方程组:3.1.2两相邻时间步长的代数方程组对式(8)中的时间微分进行一阶向后差分,得到全隐式的非线性代数方程组为式中:Rnκ,k+1为κ的余量;∆t为时间步长;上标k和k+1为两相邻的时间步。右端的流量项和源汇项均采用新的时间步的值,计算表明,这种全隐式方法提高了数值稳定性。3.1.3计算域和各变量间的关系对式(9)应用Newton-Raphson迭代算法,得到式(10)所示的线性方程组。对于包含NEL个单元的计算域,所有的线性方程可写成如下矩阵形式:式中:因为κ=w,a,J为2×NEL阶的雅可比系数矩阵;∆X为独立变量在相邻两次迭代的增量列阵;B为已知向量列阵。J中各项∂Rn∂xi采用数值差分近似求得,由于与nV相邻的单元的数目有限,因此,雅克比矩阵J是稀疏矩阵,式(11)可通过求解大型稀疏系数矩阵的线性方程组直接求解。3.2空气边界单元水-气二相流模型的边界条件可根据液相或气相划分,但仍然分为两类:Dirichlet(第一类)边界条件和Neumann(第二类)边界条件。Dirichlet边界条件单元的主要变量始终保持不变,为了实现这一点,给边界单元的体积一个相对于土体单元很大的值,如1.0×1050m3,边界单元的主要变量初值不会因为与土体单元的流量交换而改变。对于已知水头边界,主要变量为pg=pl,Xla很小,为保证数值稳定,取1.0×10-10;对于空气边界,pg=patm,Xga很大,可取0.9999。Neumann边界条件描述的是系统与外界的流量交换情况,边界单元的单位流量对应式(1)中的源汇项,流入为正,可以是常量,也可以是随时间变化的,相态在计算过程中会随着流量交换情况而变化,而边界单元与土体之间的流量交换即为真正进入到土体的流量。水位变化范围内的边界条件是随时间变化的Dirichlet边界条件,目前通常是将水位变化过程等效为一系列的特征水位,即将非稳定渗流简化为一系列稳定渗流。然而这种简化在水-气二相流模型中不能很好地描述边界上相态的转化过程,在此将水位引起的卸载作用等效为边界单元有一个流出的流量,即随时间变化的Dirichlet边界条件通过Neumann边界条件来实现,边界单元上源汇项的计算式为式中:m为边界单元上的源汇项(kg/s),流出为负;b为拟合系数,与土体的性质有关;V为边界单元体积,1.0×1050m3;C为边界单元的压缩性,非常小,取1.0×10-6(1/Pa);∆P为水压力的变化,水位下降时为负;∆t为水压力变化所需的时间。4计算与分析4.1渗流逸低表面活性剂模型验证关于本文建立的恒温下的水-气二相流模型,目前还没有精确的理论解,但Muskat早在1937年给出了水流通过上、下游坝面垂直的简单土坝剖面的理论解,见图2(a)所示,得到下游面渗流逸出面的长度s=1.54m。后来许多学者对Muskat问题进行深入研究,得到土坝剖面形状、水位及渗流逸出面的长度之间的关系。在此,利用本文的水-气二相流模型来模拟Muskat问题,求得下游的渗流逸出面长度与理论解进行比较,从而验证本文所提出的水-气二相流模型的有效性。Muskat问题的土坝剖面网格见图2(b),计算得到的孔隙水压力及渗流逸出面长度和孔隙气压力分布图如图3、图4所示。由图3可知,孔隙水压力为0的等值线即为浸润线,下游渗流逸出面的长度约为1.67m。考虑到网格剖分、数值计算精度等误差,这个结果与Muskat给出的解析解1.54m具有一致性,证明了该模型的有效性。由图4所示的孔隙气压力分布可见,在自由面附近的孔隙气压力水头约为-0.01~0.01m,说明稳定渗流状态下,非饱和区基本与大气相通,与单相流模型中的假定是一致的。图4孔隙气压力水头分布(单位:m)Fig.4Poreairpressureheadcontourlines(unit:m)4.2初始条件假设某库岸边坡的横断面如图5(a)所示,坡比为1:3,坡外水位在26h内由HW=184.9m下降到LW=179.3m,降落的水位差为dh=5.6m,降落的速度为6.0×10-5m/s,土体孔隙率为0.4,网格剖分的横断面图见图5(b)。假设非饱和水力函数关系中的拟合参数为λ=0.457,n=1+λ,1/p0=0.5×10-4Pa-1,并且Srl=0.15,Ssl=0.9,Sgr=0.1。在稳定渗流情况下,上、下游水位以下的边界为已知水头边界,水位以上为空气边界。水位下降时,边界AB和MN为不透水边界(图5(a)),即边界上的流量为0,在本模型中只要去掉该边界上的边界单元即可;BC、CD和DE为随时间变化的Dirichlet边界条件,根据式(12),边界单元的源汇项m=-2.7E45kg/s。初始条件假设库岸边坡土体水相饱和,初始气压力为标准大气压。HW和LW水位时的稳定渗流的孔隙水压力分布图,以及库水位由HW降落到LW过程的浸润线变化过程见图6所示,由图6(c)可见,浸润线下降的速度明显滞后于库水位的下降速度。降落历时26h,水位由HW降为LW时的岸坡渗流状态分布图见图7所示。结果表明,浸润线以下为饱和区,孔隙中充满了重力水,这个区域内的孔隙水压力和孔隙气压力均为正值,且两者相等,毛细压力为0。然而在浸润线以上的区域,可为两部分,一部分为浸润线以上、毛细压力为0的分界线以下的区域,这个区域内的孔隙水压力和孔隙气压力均为负值,且两者相等,毛细压力为0,对应水相饱和,这个区域是毛细管水饱和带,水位下降后,在岸坡表面空气的趋替作用,孔隙中的水在重力作用下流向浸润线,由于气、液表面张力的作用,这部分水消散的非常慢;另一部分是毛细压力为负值的传统意义上的非饱和区,孔隙水压力和孔隙气压力均为负值,且孔隙水压力小于孔隙气压力,毛细压力为负值。5在非饱和排除系统中的应用本文通过建立水-气二相流模型,研究库水位下降的岸坡非稳定渗流场,得到如下几点结论。(1)将库水位下降下的岸坡非稳定渗流场作为水气二相流问题,采用水-气二相流模型进行研究是符合实际的物理意义的,计算结果合理,能够反映出气相对孔压的变化的影响;模型中所采用的各种边界条件的数值处理方法精确、有效,尤其是随时间变化的Dirichlet边界条件的处理方法,能够实现岸坡由水相饱和向非饱和状态的转变;将高效的积分有限差分法引入到渗流问题的求解中,有效减小了计算量,从而提高了计算速度。(2)水位下降时的岸坡可分为3部分,浸润线以下的饱和区、浸润线以上、毛细