充气法治理承压含水层海水入侵的数值模拟方法

摘要

本发明涉及防止海水入侵技术领域。为治理承压含水层海水入侵，本发明采用的技术方案是：充气法治理承压含水层海水入侵的数值模拟方法，包括如下步骤：步骤一：建立地下水气-液二相流及溶质运移模型，包括基本控制方程及辅助方程，不考虑温度对系统的影响；步骤二：模型求解；步骤三：模型边界条件确定：步骤四：模型验证：运用地下水气-液二相流及溶质运移模型，模拟分析承压含水层受咸、淡水共同作用的静力平衡情况，并与已有文献的计算结果进行对比验证；步骤五：分析充气法治理承压含水层海水入侵的效果。本发明主要应用于防止海水入侵。

主权项

1.一种充气法治理承压含水层海水入侵的数值模拟方法，其特征是，包括如下步骤：步骤一：建立地下水气-液二相流及溶质运移模型，包括基本控制方程及辅助方程，不考虑温度对系统的影响，具体如下：模型的基本控制方程为：$\frac{{\partial M}^{\kappa }}{\partial t}=-\mathrm{div}\left(F^{\kappa }\right)+q^{\kappa }---\left(1\right)$式中，M^κ表示κ组分(纯水w、参考盐水b和空气a)的累积质量密度，F^κ为κ组分的流量密度，包括平流流量密度和扩散流量密度q^κ为组分κ的源汇项；步骤二：模型求解：以TOUGH2/EOS7为工具，空间上采用积分形式的有限差分方法IFDM进行离散，时间上采用一阶向后差分的全隐式方法进行离散，模型的线性化采用Newton-Raphson迭代方法，最后得到大型稀疏系数矩阵的线性方程组；具体如下：(1)空间上采用积分有限差分法IFDM进行离散首先将计算域离散成子单元，其性质由形心点代表，分别对各个单元的质量平衡方程的积分形式进行空间离散，对于任意单元n，单元体积为V_n，边界面为Г_n，单元的质量平衡方程的积分形式如下：$\frac{d}{\mathrm{dt}}{\int }_{V_n}{M}^{\kappa }d{V}_n={\int }_{{\Gamma }_n}{F}^{\kappa }·\mathrm{nd}{\Gamma }_n+{\int }_{V_n}{q}^{\kappa }d{V}_n---\left(2\right)$式中，n为表面单元dГ_n的单位法向量，指向控制单元体内为正；引入适当的体积平均值，有${\int }_{V_n}{M}^{\kappa }d{V}_n=V_n{M}_n^{\kappa }---\left(3\right)$${\int }_{V_n}{q}^{\kappa }d{V}_n=q_b^{\kappa }{V}_n---\left(4\right)$式中，为M^κ，q^κ在V_n上的平均值；Г_n上的面积分可近似为其所包含的各个表面A_nm的面积分的平均值之和，有${\int }_{{\Gamma }_n}{F}^{\kappa }·\mathrm{nd}{\Gamma }_n=\underset{m}{\Sigma }{A}_{\mathrm{nm}}{F}_{\mathrm{nm}}^{\kappa }---\left(5\right)$式中，m为与单元n相邻的所有单元，A_nm是单元n和m相邻的交界面，是F^κ在面A_nm上的内法线方向的平均值；将式(3)、(4)和(5)代入到式(2)中，得到一组关于时间的一阶微分方程组$\frac{{\mathrm{dM}}_n^{\kappa }}{\mathrm{dt}}=\frac{1}{V_n}\underset{m}{\Sigma }{A}_{\mathrm{nm}}{F}_{\mathrm{nm}}^{\kappa }+q_n^{\kappa }---\left(6\right)$(2)时间上采用一阶向后差分方法进行离散对式(6)的时间微分采用一阶向后差分方法，得到任意单元的全隐式非线性方程组，见式(7)：$R_n^{\kappa ,k+1}=M_n^{\kappa ,k+1}-M_n^{\kappa ,k}-\frac{\mathrm{\Delta t}}{V_n}\{\underset{m}{\Sigma }{A}_{\mathrm{nm}}{F}_{\mathrm{nm}}^{\kappa ,k+1}+V_n{q}_n^{\kappa ,k+1}\}---\left(7\right)$式中，引入了组分κ＝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)：$-\underset{i}{\Sigma }\frac{\partial R_n^{\kappa ,k+1}}{\partial x_i}{|}_p(x_{i,p+1}-x_{i,p})=R_n^{\kappa ,k+1}\left(x_{i,p}\right)---\left(8\right)$步骤三：模型边界条件确定：模型计算中的边界条件包括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)空气损失变化：充气初始时期的空气损失较小，然后迅速上升，之后逐渐趋于平稳，达到相对稳定值。