[发明专利]充气法治理承压含水层海水入侵的数值模拟方法

权利要求书

1.一种充气法治理承压含水层海水入侵的数值模拟方法，其特征是，包括如下步骤：

步骤一：建立地下水气-液二相流及溶质运移模型，包括基本控制方程及辅助方程，不考虑温度对系统的影响，具体如下：

模型的基本控制方程为：

∂Mκ∂t=-div(Fκ)+qκ---(1)]]>

式中，M^κ表示κ组分(纯水w、参考盐水b和空气a)的累积质量密度，F^κ为κ组分的流量密度，包括平流流量密度和扩散流量密度q^κ为组分κ的源汇项；

步骤二：模型求解：以TOUGH2/EOS7为工具，空间上采用积分形式的有限差分方法IFDM进行离散，时间上采用一阶向后差分的全隐式方法进行离散，模型的线性化采用Newton-Raphson迭代方法，最后得到大型稀疏系数矩阵的线性方程组；具体如下：

(1)空间上采用积分有限差分法IFDM进行离散

首先将计算域离散成子单元，其性质由形心点代表，分别对各个单元的质量平衡方程的积分形式进行空间离散，对于任意单元n，单元体积为V_n，边界面为Г_n，单元的质量平衡方程的积分形式如下：

ddt∫VnMκdVn=∫ΓnFκ·ndΓn+∫VnqκdVn---(2)]]>

式中，n为表面单元dГ_n的单位法向量，指向控制单元体内为正；

引入适当的体积平均值，有

∫VnMκdVn=VnMnκ---(3)]]>

∫VnqκdVn=qbκVn---(4)]]>

式中，为M^κ，q^κ在V_n上的平均值；

Г_n上的面积分可近似为其所包含的各个表面A_nm的面积分的平均值之和，有

∫ΓnFκ·ndΓn=ΣmAnmFnmκ---(5)]]>

式中，m为与单元n相邻的所有单元，A_nm是单元n和m相邻的交界面，是F^κ在面A_nm上的内法线方向的平均值；

将式(3)、(4)和(5)代入到式(2)中，得到一组关于时间的一阶微分方程组

dMnκdt=1VnΣmAnmFnmκ+qnκ---(6)]]>

(2)时间上采用一阶向后差分方法进行离散

对式(6)的时间微分采用一阶向后差分方法，得到任意单元的全隐式非线性方程组，见式(7)：

Rnκ,k+1=Mnκ,k+1-Mnκ,k-ΔtVn{ΣmAnmFnmκ,k+1+Vnqnκ,k+1}---(7)]]>

式中，引入了组分κ＝w，b，a的余量Δt为时间步长，上标k和k+1表示两相邻的时间步长指标；其中，分别表示k、k+1时刻M^κ在单元体积V_n上的平均值，表示k+1时刻F^κ在面A_nm上的内法线方向的平均值，表示k+1时刻q^κ在单元体积V_n上的平均值；右端的流量项和源汇项均采用新的时间步长值；

(3)Newton-Raphson迭代方法

运用Newton-Raphson迭代方法进行线性化，引入迭代指标p，对式(7)中的余量在迭代步p+1处进行泰勒级数展开，只保留一阶项，得到包含3×NEL即计算域内单元数个方程的线性方程组，并且以两迭代步的增量为未知量；最后得到大型稀疏系数矩阵的线性方程组，如式(8)：

-Σi∂Rnκ,k+1∂xi|p(xi,p+1-xi,p)=Rnκ,k+1(xi,p)---(8)]]>

步骤三：模型边界条件确定：模型计算中的边界条件包括Dirichlet边界条件和Neumann边界条件两种，其数学处理方法如下：

(1)Dirichlet边界条件

Dirichlet边界条件上，边界条件单元的主要变量在计算过程中保持不变，为此，设定边界条件单元的体积非常大，当边界单元的体积相对于土体单元很大时，与土体单元的流量交换将不会改变边界单元的主要变量值；

①对于空气边界，其边界条件单元的相态为仅有气相状态，主要变量为孔隙气压力p_g、参考盐水占相态的质量百分数X_b、空气占气相的质量分数和温度T，对于与大气接触的边界上p_g＝p_atm，对于有空气超压作用的边界上，如人工充气墙边界上，p_g＝p_atm+Δp，其中p_atm为大气压力，Δp为充气压力；X_b＝0.0；T为气温；

②对于已知水头边界，包括地下淡水水头边界和海水水头边界，其边界条件单元的相态均为仅有液相状态，主要变量为孔隙气压力p_g、参考盐水占液相的质量百分数X_b、空气占液相的质量分数和温度T，由于液相饱和状态下毛细压力为零，因此有p_g＝p_atm+ρ_lgh，其中ρ_l为淡水或海水的密度，h为边界上的水深；地下淡水水头边界上X_b等于零，海水边界上的X_b可根据参考盐水的特性以及海水的盐度计算出来；

(2)Neumann边界条件

Neumann边界条件描述的是系统与外界的流量交换情况，边界条件单元的单位流量对应式中的源汇项，流入为正，可以是常量，也可以随时间变化，对于不透水边界，是一类特殊的Neumann边界条件，边界上的流量为零，积分形式的有限差分法的处理很简单，就是不设置与之相邻的边界条件单元，则与不透水面单元的流量交换为零；

步骤四：模型验证：运用地下水气-液二相流及溶质运移模型，模拟分析承压含水层受咸、淡水共同作用的静力平衡情况，并与已有文献的计算结果进行对比验证；

步骤五：根据地下水气-液二相流及溶质运移模型，以咸-淡水静力平衡情况作为初始条件，施加充气作用，分析充气法治理承压含水层海水入侵的效果：

(1)盐度变化：充气结束时，相对盐度等值线与不透水底板的交点向海水方向撤退，海水入侵范围减小；

(2)水、气相压力及流场变化：水相压力和气相压力的变化规律基本一致，分为三个阶段：迅速增大至峰值、由峰值迅速减小、缓慢减小趋近于相对稳定值；另外，由于毛细压力的作用，观察点的水相压力略小于气相压力；在注入区及含水层顶部，水、气相压力增大的较多，形成指向海水侧的水力梯度，从而驱替入侵的咸水退出含水层；

(3)空气损失变化：充气初始时期的空气损失较小，然后迅速上升，之后逐渐趋于平稳，达到相对稳定值。

2.如权利要求1所述的充气法治理承压含水层海水入侵的数值模拟方法，其特征是，

(1)累积质量密度M^κ等于β相中组分κ的质量之和，其表达式为：

Mκ=ΣβXβκφSβρβ---(9)]]>

式中，β表示流体相(气相或液相)，液相为纯水和参考盐水的混合物，φ表示孔隙率，S_β为β相的饱和度，ρ_β为β相的密度，为κ组分占β相的质量百分数；

(2)平流流量的表达式为：

Fadvκ=ΣβXβκFβ---(10)]]>

式中，F_β为β相的平流流量，符合达西定律，其多相流形式的表达式为：

Fβ=-kρβkγβ(Sw)μβ(▿pβ-ρβg)---(11)]]>

其中，k为固有渗透系数，k_rβ为β相的相对渗透系数，μ_β为β相的粘滞性系数，p_β为β相的孔隙压力，g为重力加速度矢量；

液相的孔隙压力p_l为气相压力p_g和毛细压力p_c之和：

p_l＝p_g+p_c                     (12)

关于毛细压力～饱和度的关系，采用Leverett模型来表征：

p_c＝p₀γ[1.417(1-S_e)-2.120(1-S_e)²+1.263(1-S_e)³]        (13)

式中，p₀为进气值，γ为表面张力，有效水饱和度S_e的表达式为S_e＝(S_l-S_lr)/(S_ls-S_lr)，S_l为液相饱和度，S_lr为液相剩余饱和度，S_ls为液相最大饱和度；

关于相对渗透率～饱和度的关系吗，采用Fatt and Klikoff模型来表征：

液相：krl=Se3---(14)]]>

气相：k_rg＝(1-S_e)³       (15)

(3)扩散流量的表达式为：

Fdiffκ=-φτ0Σβτβ(Sβ)ρβdβκ▿Xβκ---(16)]]>

式中，为κ组分在β相中的分子扩散系数；τ₀τ_β(S_β)为迂曲度，τ₀是多孔介质特性参数，τ_β＝τ_β(S_β)是饱和度的函数；

分子扩散系数随压力p和温度T变化的计算式为：

式中，θ为温度系数，通常取值为1.8；

TOUGH2中提供了三种迂曲度模型：分子扩散系数τ₀一般取1.0

relative permeability model：τ₀τ_β＝τ₀k_rβ(S_β)

Millington-Quirk Model：τ0τβ=φ1/3Sβ10/3]]>

constant diffusivity：τ₀τ_β＝τ₀S_β              (18)

Relative permeability model为相对渗透率模型，Millington-Quirk Model是由Millington-Quirk提出的非饱和土壤导水率模型，constant diffusivity为常扩散系数模型。