CN104537232A - 考虑Lisse现象的浅层地下水位预测方法 - Google Patents

Abstract

本发明涉及预测浅层地下水位波动技术领域，为模拟分析降雨条件下，考虑Lisse现象的地下水位的波动情况。为此，本发明采取的技术方案是，考虑Lisse现象的浅层地下水位预测方法，包括如下步骤：步骤一：建立非饱和土的水-气二相流模型，包括基本控制方程及本构关系，假定系统一直处于恒温状态；步骤二：模型求解；步骤三：模型边界条件确定；步骤四：利用非饱和土的水-气二相流模型，以稳定渗流情况作为降雨入渗的初始条件，模拟分析降雨条件下，考虑Lisse现象的浅层地下水位的波动情况：本发明主要应用于预测浅层地下水位波动。

Description

考虑Lisse现象的浅层地下水位预测方法

技术领域

本发明涉及预测浅层地下水位波动技术领域；具体讲，涉及考虑Lisse现象的浅层地下水位预测方法。

技术背景

地下水作为水资源的重要组成部分，以水质较好、分布广泛以及便于开发等优点，成为理想的供水水源之一，所以采用观测井对地下水位进行实时监测也越来越重要。然而对于降雨条件下，地下水位变化较大较快的情形，观测井的监测数据可能存在偏差。因此，采用数学模型模拟地下水位，更能精细地反映出水位的波动情况，且避免了不必要的人力和物力。

在传统意义上，地下水位与土壤中饱和区的顶部相一致，然而近年来，随着对饱和-非饱和渗流研究的不断深入，越来越多的研究发现，地下水位以上存在饱和区，地下水位以下存在非饱和区，因此，根据饱和度来定义地下水位是不合理的。本文采用孔隙水压力等于零的线(面)，即浸润线，来代表地下水位，这也是土壤中水分流入观测井的必要条件。

Lisse现象即为，在强降雨条件下，观测井中水位快速上升，上升高度为湿润峰深度的若干倍，且水位衰退须持续几小时至几天不等，而此时地下水并未获得真正的补给。其内在原因为，雨水封闭了土体表面，湿润峰以下的气体受到压缩，导致非饱和区中孔隙气压力增加，由此井中水位受到压迫而上升。因此，不能准确地识别出Lisse现象，可能会过量估计地下水量，进而影响地下水的开采工作。

发明内容

为克服现有技术的不足，模拟分析降雨条件下，考虑Lisse现象的地下水位的波动情况。为此，本发明采取的技术方案是，考虑Lisse现象的浅层地下水位预测方法，包括如下步骤：

步骤一：建立非饱和土的水-气二相流模型，包括基本控制方程及本构关系，假定系统一直处于恒温状态，具体如下：

模型的基本控制方程为：

$\frac{\∂M^{\mathrm{\κ}}}{\∂t}=-\mathrm{div}{F}^{\mathrm{\κ}}+q^{\mathrm{\κ}}---\left(1\right)$

式中，M^κ表示包括空气a和水w的κ组分的累积质量密度；F^κ为κ组分的平流流量；q^κ为κ组分的源汇项，t表示时间；

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

式中，β表示液相l或气相g，其中液相包括液态孔隙水及溶解的空气，气相包括干燥空气和水蒸气；为κ组分占β相的质量百分数；φ为孔隙率；S_β为β相的饱和度；ρ_β为β相的密度；

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

$F^{\mathrm{\κ}}=\underset{\mathrm{\β}}{\mathrm{\Σ}}{X}_{\mathrm{\β}}^{\mathrm{\κ}}{F}_{\mathrm{\β}}---\left(3\right)$

式中，F_β为β相的平流流量，遵循达西定律，其表达式为：

$F_{\mathrm{\β}}=-K\frac{{\mathrm{\ρ}}_{\mathrm{\β}}{k}_{\mathrm{r\β}}\left(S_l\right)}{{\mathrm{\μ}}_{\mathrm{\β}}}(\▿p_{\mathrm{\β}}-{\mathrm{\ρ}}_{\mathrm{\β}}g)---\left(4\right)$

其中，K为固有渗透率；k_rβ为β相的相对渗透率，与液相饱和度S_l有关；μ_β为β相的粘滞性系数；p_β为β相的压力；为β相的压力梯度；g为重力加速度矢量；

毛细压力p_c为液相压力p_l与气相压力p_g之间的差值：

p_c＝p_l-p_g               (5)

模型的本构关系包括毛细压力-饱和度关系与相对渗透率-饱和度关系，其中关于毛细压力-饱和度关系曲线采用van Genuchten模型简称VG模型表示：

p_c＝-p₀[(S^*)^-1/λ-1]^1-λ(-p_max≤p_c≤0)   (6)

式中，p₀为进气值；λ为模型拟合参数，与土体均匀程度有关；S^*为有效饱和度，表示为S^*＝(S_l-S_lr)/(S_ls-S_lr)，S_l为液相饱和度，S_lr为残余液相饱和度，S_ls为饱和液相饱和度；

关于相对渗透率-饱和度的关系，其中液相相对渗透率k_rl采用van Genuchten-Mualem模型简称VG-M模型表征：

$k_{\mathrm{rl}}=\left\{\begin{array}{cc}\sqrt{S^{*}}{[1-{(1-{\left(S^{*}\right)}^{1/\mathrm{\&lambda;}})}^{\mathrm{\&lambda;}}]}^2& (S_l<S_{\mathrm{ls}})\\ 1& (S_l\&GreaterEqual;S_{\mathrm{ls}})\end{array}\right.---\left(7\right)$

气相相对渗透率k_rg采用科里(Corey)提出的表达式：

$k_{\mathrm{rg}}=\left\{\begin{array}{cc}1-k_{\mathrm{rl}}& (S_{\mathrm{gr}}=0)\\ {(1-\hat{S})}^2(1-{\hat{S}}^2)& (S_{\mathrm{gr}}>0)\end{array}\right.---\left(8\right)$

式中，S_gr为残余气相饱和度；

步骤二：模型求解：以TOUGH2/EOS3为工具，将变量分为主要变量和次要变量，其中，次要变量可通过主要变量求得。空间上采用积分形式的有限差分方法(IFDM)进行离散，时间上采用一阶向后差分的全隐式方法进行离散，模型的线性化采用牛顿-拉斐逊(Newton-Raphson)迭代方法，最后得到大型稀疏系数矩阵的线性方程组；

步骤三：模型边界条件确定：边界条件包括狄利克雷(Dirichlet)边界条件和黎曼(Neumann)边界条件，其数学处理方法如下：

(1)Dirichlet边界条件

将边界条件单元的体积设为1×10⁵⁰m³。包括空气边界和已知水头边界：

①对于空气边界，仅有气相，其主要变量为p_g，和T，其中p_g为气相压力，为空气占气相的质量百分数，取值为0.9999，T为系统温度；

②对于已知水头边界，仅有液相状态，其主要变量为p_l，和T，其中p_l为液相压力，为液相中空气所占的质量百分数，取值为1.0×10^-10，T为系统温度；

(2)Neumann边界条件

Neumann边界条件描述的是系统与外界的流量交换情况，流入为正，流出为负；

①对于降雨入渗条件，可通过m_r(t)＝ρ_wA_eQ_w(t)，设置一个适当大小的源汇项来实现，其中，源汇项m_r表示单位时间内通过土体表面法向的水分质量，流入为正，流出为负；A_e为土体的有效面积；Q_w为降雨强度，ρ_w：水的密度；

②对于不透水边界或不入流边界，看作一类特殊的Neumann边界条件，边界上的流量设为零；

步骤四：利用非饱和土的水-气二相流模型，以稳定渗流情况作为降雨入渗的初始条件，模拟分析降雨条件下，考虑Lisse现象的浅层地下水位的波动情况：

(1)渗流场变化：降雨过程中，非饱和区的气体受到湿润峰下移的压缩，孔隙气压力(p_g-p_atm)/ρ_wg增加；接近地表处的孔隙水压力(p_l-p_atm)/ρ_wg大于零，形成暂态饱和区；在湿润峰下方，由于孔隙气压力的增加，孔隙水压力偏离了初始值，浸润线即p_l＝0的线或面上升，降雨结束之后，非饱和区中的气体可以自由溢出，孔隙气压力减小为零，湿润峰下方的孔隙水压力也恢复到初始值，当湿润峰到达初始地下水位时，初始地下水位上方形成暂态饱和区，浸润线再次上升；

(2)地下水位变化：降雨过程中，地下水位首先快速增加，之后缓慢增加，并在降雨结束时达到最大值，而此时湿润峰并未到达初始地下水位处，说明水位的快速上升是由Lisse现象引起的；当降雨结束后，浸润线快速回到初始地下水位，之后逐渐上升，说明此时地下水获得真正的补给，因此，地下水位的波动有两部分原因：Lisse现象和真正的补给，在地下水开采工作中，没能正确地识别Lisse现象，可能会造成地下水的过量开采。

步骤二进一步具体为：

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

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

$\frac{d}{\mathrm{dt}}{\&Integral;}_{V_n}{M}^{\mathrm{\&kappa;}}d{V}_n={\&Integral;}_{{\mathrm{\&Gamma;}}_n}{F}^{\mathrm{\&kappa;}}\&CenterDot;\mathrm{nd}{\mathrm{\&Gamma;}}_n+{\&Integral;}_{V_n}{q}^{\mathrm{\&kappa;}}d{V}_n---\left(9\right)$

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

引入体积平均值：

${\&Integral;}_{V_n}{M}^{\mathrm{\&kappa;}}d{V}_n=V_n{M}_n^{\mathrm{\&kappa;}}---\left(10\right)$

${\&Integral;}_{V_n}{q}^{\mathrm{\&kappa;}}d{V}_n=q_n^{\mathrm{\&kappa;}}{V}_n---\left(11\right)$

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

表面单元Γ_n上的面积分可近似为其所包含的各个表面A_nm上面积分的平均值之和：

${\&Integral;}_{{\mathrm{\&Gamma;}}_n}{F}^{\mathrm{\&kappa;}}\&CenterDot;\mathrm{nd}{\mathrm{\&Gamma;}}_n=\underset{m}{\mathrm{\&Sigma;}}{A}_{\mathrm{nm}}{F}_{\mathrm{nm}}^{\mathrm{\&kappa;}}---\left(12\right)$

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

将式(10)、(11)和(12)代入到式(9)中，得到一组关于时间的一阶微分方程组

$\frac{d{M}_n^w}{\mathrm{dt}}=\frac{1}{V_n}\underset{m}{\mathrm{\&Sigma;}}{A}_{\mathrm{nm}}{F}_{\mathrm{nm}}^{\mathrm{\&kappa;}}+q_n^{\mathrm{\&kappa;}}---\left(13\right)$

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

采用一阶向后差分方法，对式(13)进行时间上的离散，得到任一单元的全隐式非线性方程组：

$R_n^{\mathrm{\&kappa;},k+1}=M_n^{\mathrm{\&kappa;},k+1}-M_n^{\mathrm{\&kappa;},k}-\frac{\mathrm{\&Delta;t}}{V_n}\{\underset{m}{\mathrm{\&Sigma;}}{A}_{\mathrm{nm}}{F}_{\mathrm{nm}}^{\mathrm{\&kappa;},k+1}+V_n{q}_n^{\mathrm{\&kappa;},k+1}\}---\left(14\right)$

式中，引入了组分κ的余量上标k和k+1分别表示两相邻的时间步长指标；△t为时间步长；

(3)Newton-Raphson迭代方法

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

$-\underset{i}{\mathrm{\&Sigma;}}\frac{\&PartialD;R_n^{\mathrm{\&kappa;},k+1}}{\&PartialD;x_i}{|}_p(x_{i,p+1}-x_{i,p})=R_n^{\mathrm{\&kappa;},k+1}\left(x_{i,p}\right)---\left(15\right).$

与已有技术相比，本发明的技术特点与效果：

(1)建立非饱和土的水-气二相流模型，模拟了降雨条件下，考虑Lisse现象的浅层地下水位的波动情况。

(2)本发明避免了野外监测的高费用以及结果偏差问题，模拟结果合理可信，避免了由于Lisse现象造成的地下水过量开采的后果。

附图说明

图1土体几何形状。

图2(a)毛细压力-饱和度关系曲线；(b)相对渗透率-饱和度关系曲线。

图3降雨过程中不同时刻的土体剖面(a)液相饱和度；(b)毛细压力；

(c)孔隙气压力；(d)孔隙水压力的分布情况。

图4降雨结束后不同时刻的土体剖面(a)液相饱和度；(b)毛细压力；

(c)孔隙气压力；(d)孔隙水压力的分布情况。

图5(a)B点压力水头和气体流速；(b)C点压力水头和液相饱和度；(c)D点压力水头随时间的变化情况。

图6Lisse现象原理图。

图7地下水位随时间的变化情况

图8本发明的模拟方法流程图

具体实施方式

本发明旨在基于非饱和土中多孔介质的多相流理论，在同时考虑水相和气相流动规律及其相互作用的基础上，建立了非饱和土的水-气二相流模型，来模拟分析降雨条件下，考虑Lisse现象的地下水位的波动情况。

本发明采用的技术方案是：

步骤一：建立非饱和土的水-气二相流模型，包括基本控制方程及本构关系，假定系统一直处于恒温状态。具体如下：

模型的基本控制方程为：

$\frac{\&PartialD;M^{\mathrm{\&kappa;}}}{\&PartialD;t}=-\mathrm{div}{F}^{\mathrm{\&kappa;}}+q^{\mathrm{\&kappa;}}---\left(1\right)$

式中，M^κ表示κ组分(空气a和水w)的累积质量密度；F^κ为κ组分的平流流量；q^κ为κ组分的源汇项。

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

式中，β表示液相(l)或气相(g)，其中液相包括液态孔隙水及溶解的空气，气相包括干燥空气和水蒸气；为κ组分占β相的质量百分数；φ为孔隙率；S_β为β相的饱和度；ρ_β为β相的密度。

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

$F^{\mathrm{\&kappa;}}=\underset{\mathrm{\&beta;}}{\mathrm{\&Sigma;}}{X}_{\mathrm{\&beta;}}^{\mathrm{\&kappa;}}{F}_{\mathrm{\&beta;}}---\left(3\right)$

式中，F_β为β相的平流流量，遵循达西定律，其表达式为：

$F_{\mathrm{\&beta;}}=-K\frac{{\mathrm{\&rho;}}_{\mathrm{\&beta;}}{k}_{\mathrm{r\&beta;}}\left(S_l\right)}{{\mathrm{\&mu;}}_{\mathrm{\&beta;}}}(\&dtri;p_{\mathrm{\&beta;}}-{\mathrm{\&rho;}}_{\mathrm{\&beta;}}g)---\left(4\right)$

其中，K为固有渗透率；k_rβ为β相的相对渗透率，与液相饱和度S_l有关；μ_β为β相的粘滞性系数；p_β为β相的压力；为β相的压力梯度；g为重力加速度矢量。

毛细压力p_c为液相压力p_l与气相压力p_g之间的差值：

p_c＝p_l-p_g             (5)

模型的本构关系包括毛细压力-饱和度关系与相对渗透率-饱和度关系。其中关于毛细压力-饱和度关系曲线采用van Genuchten模型(VG模型)表示：

p_c＝-p₀[(S^*)^-1/λ-1]^1-λ(-p_max≤p_c≤0)    (6)

式中，p₀为进气值；λ为模型拟合参数，与土体均匀程度有关；S^*为有效饱和度，可表示为S^*＝(S_l-S_lr)/(S_ls-S_lr)，S_l为液相饱和度，S_lr为残余液相饱和度，S_ls为饱和液相饱和度。

关于相对渗透率-饱和度的关系，其中液相相对渗透率k_rl采用van Genuchten-Mualem模型(VG-M模型)表征：

$k_{\mathrm{rl}}=\left\{\begin{array}{cc}\sqrt{S^{*}}{[1-{(1-{\left(S^{*}\right)}^{1/\mathrm{\&lambda;}})}^{\mathrm{\&lambda;}}]}^2& (S_l<S_{\mathrm{ls}})\\ 1& (S_l\&GreaterEqual;S_{\mathrm{ls}})\end{array}\right.---\left(7\right)$

气相相对渗透率k_rg采用Corey(科里)提出的表达式：

$k_{\mathrm{rg}}=\left\{\begin{array}{cc}1-k_{\mathrm{rl}}& (S_{\mathrm{gr}}=0)\\ {(1-\hat{S})}^2(1-{\hat{S}}^2)& (S_{\mathrm{gr}}>0)\end{array}\right.---\left(8\right)$

式中，S_gr为残余气相饱和度。

步骤二：模型求解：以TOUGH2/EOS3为工具，将变量分为主要变量和次要变量，其中，次要变量可通过主要变量求得。空间上采用积分形式的有限差分方法(IFDM)进行离散，时间上采用一阶向后差分的全隐式方法进行离散，模型的线性化采用Newton-Raphson(牛顿-拉斐逊)迭代方法，最后得到大型稀疏系数矩阵的线性方程组；具体如下：

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

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

$\frac{d}{\mathrm{dt}}{\&Integral;}_{V_n}{M}^{\mathrm{\&kappa;}}d{V}_n={\&Integral;}_{{\mathrm{\&Gamma;}}_n}{F}^{\mathrm{\&kappa;}}\&CenterDot;\mathrm{nd}{\mathrm{\&Gamma;}}_n+{\&Integral;}_{V_n}{q}^{\mathrm{\&kappa;}}d{V}_n---\left(9\right)$

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

引入体积平均值：

${\&Integral;}_{V_n}{M}^{\mathrm{\&kappa;}}d{V}_n=V_n{M}_n^{\mathrm{\&kappa;}}---\left(10\right)$

${\&Integral;}_{V_n}{q}^{\mathrm{\&kappa;}}d{V}_n=q_n^{\mathrm{\&kappa;}}{V}_n---\left(11\right)$

式中，为M^κ，q^κ在单元体积V_n上的平均值。

表面单元Γ_n上的面积分可近似为其所包含的各个表面A_nm上面积分的平均值之和：

${\&Integral;}_{{\mathrm{\&Gamma;}}_n}{F}^{\mathrm{\&kappa;}}\&CenterDot;\mathrm{nd}{\mathrm{\&Gamma;}}_n=\underset{m}{\mathrm{\&Sigma;}}{A}_{\mathrm{nm}}{F}_{\mathrm{nm}}^{\mathrm{\&kappa;}}---\left(12\right)$

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

将式(10)、(11)和(12)代入到式(9)中，得到一组关于时间的一阶微分方程组

$\frac{d{M}_n^w}{\mathrm{dt}}=\frac{1}{V_n}\underset{m}{\mathrm{\&Sigma;}}{A}_{\mathrm{nm}}{F}_{\mathrm{nm}}^{\mathrm{\&kappa;}}+q_n^{\mathrm{\&kappa;}}---\left(13\right)$

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

采用一阶向后差分方法，对式(13)进行时间上的离散，得到任一单元的全隐式非线性方程组：

$R_n^{\mathrm{\&kappa;},k+1}=M_n^{\mathrm{\&kappa;},k+1}-M_n^{\mathrm{\&kappa;},k}-\frac{\mathrm{\&Delta;t}}{V_n}\{\underset{m}{\mathrm{\&Sigma;}}{A}_{\mathrm{nm}}{F}_{\mathrm{nm}}^{\mathrm{\&kappa;},k+1}+V_n{q}_n^{\mathrm{\&kappa;},k+1}\}---\left(14\right)$

式中，引入了组分κ的余量上标k和k+1分别表示两相邻的时间步长指标；△t为时间步长；式子右端的平流流量项和源汇项均采用新的时间步长值，这种全隐式方法提高了模型求解的数值稳定性。

(3)Newton-Raphson迭代方法

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

$-\underset{i}{\mathrm{\&Sigma;}}\frac{\&PartialD;R_n^{\mathrm{\&kappa;},k+1}}{\&PartialD;x_i}{|}_p(x_{i,p+1}-x_{i,p})=R_n^{\mathrm{\&kappa;},k+1}\left(x_{i,p}\right)---\left(15\right)$

步骤三：模型边界条件确定：边界条件包括Dirichlet(狄利克雷)边界条件和Neumann(黎曼)边界条件，其数学处理方法如下：

(1)Dirichlet边界条件

在Dirichlet边界条件的数学处理方法中，一般将边界条件单元的体积看作非常大，即认为边界条件单元与相邻土体单元的流量交换，不影响边界条件单元的状态，因此可将边界条件单元的体积设为1×10⁵⁰m³。包括空气边界和已知水头边界：

①对于空气边界，仅有气相，其主要变量为p_g，和T，其中p_g为气相压力，为空气占气相的质量百分数，取值为0.9999，T为系统温度；

②对于已知水头边界，仅有液相状态，其主要变量为p_l，和T，其中p_l为液相压力，为液相中空气所占的质量百分数，取值为1.0×10^-10，T为系统温度。

(2)Neumann边界条件

Neumann边界条件描述的是系统与外界的流量交换情况，流入为正，流出为负，可以是常量，也可以随时间变化。

①对于降雨入渗条件，可通过m_r(t)＝ρ_wA_eQ_w(t)，设置一个适当大小的源汇项来实现，其中，源汇项m_r表示单位时间内通过土体表面法向的水分质量，流入为正，流出为负；A_e为土体的有效面积；Q_w为降雨强度。

②对于不透水边界或不入流边界，可看作一类特殊的Neumann边界条件，边界上的流量设为零。

步骤四：利用非饱和土的水-气二相流模型，以稳定渗流情况作为降雨入渗的初始条件，模拟分析降雨条件下，考虑Lisse现象的浅层地下水位的波动情况：

(1)渗流场变化：降雨过程中，非饱和区的气体受到湿润峰下移的压缩，孔隙气压力(p_g-p_atm)/ρ_wg增加；接近地表处的孔隙水压力(p_l-p_atm)/ρ_wg大于零，形成暂态饱和区；在湿润峰下方，由于孔隙气压力的增加，孔隙水压力偏离了初始值，浸润线(p_l＝0的线或面)上升。降雨结束之后，非饱和区中的气体可以自由溢出，孔隙气压力减小为零，湿润峰下方的孔隙水压力也恢复到初始值。当湿润峰到达初始地下水位时，初始地下水位上方形成暂态饱和区，浸润线再次上升。

(2)地下水位变化：降雨过程中，地下水位首先快速增加，之后缓慢增加，并在降雨结束时达到最大值。而此时湿润峰并未到达初始地下水位处，说明水位的快速上升是由Lisse现象引起的。当降雨结束后，浸润线快速回到初始地下水位，之后逐渐上升，说明此时地下水获得真正的补给。因此，地下水位的波动有两部分原因：Lisse现象和真正的补给，在地下水开采工作中，没能正确地识别Lisse现象，可能会造成地下水的过量开采。

下面结合附图，针对模拟降雨条件下，考虑Lisse现象的浅层地下水位的波动情况进行具体说明，其步骤如下：

(1)建立数学模型：模型选取均质各向同性的土体作为研究对象，模拟范围长100m，宽100m，厚度为2m，其剖面图见图1。地下水位位于地表以下0.6m，观测点B、C、D分别位于地表以下0.35m，0.58m和0.60m。土体材料特性参数取值见表1。毛细压力-饱和度、相对渗透率-饱和度关系曲线分别采用VG模型和VG-M模型，如图2所示。假定含水层系统处于恒温状态，T＝25℃。

表1土体材料特性参数取值

模型的基本控制方程为：

$\frac{\&PartialD;M^{\mathrm{\&kappa;}}}{\&PartialD;t}=-\mathrm{div}{F}^{\mathrm{\&kappa;}}+q^{\mathrm{\&kappa;}}---\left(1\right)$

式中，M^κ表示κ组分(空气a和水w)的累积质量密度；F^κ为κ组分的平流流量；q^κ为κ组分的源汇项。

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

式中，β表示液相(l)或气相(g)，其中液相包括液态孔隙水及溶解的空气，气相包括干燥空气和水蒸气；为κ组分占β相的质量百分数；φ为孔隙率；S_β为β相的饱和度；ρ_β为β相的密度。

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

$F^{\mathrm{\&kappa;}}=\underset{\mathrm{\&beta;}}{\mathrm{\&Sigma;}}{X}_{\mathrm{\&beta;}}^{\mathrm{\&kappa;}}{F}_{\mathrm{\&beta;}}---\left(3\right)$

式中，F_β为β相的平流流量，遵循达西定律，其表达式为：

$F_{\mathrm{\&beta;}}=-K\frac{{\mathrm{\&rho;}}_{\mathrm{\&beta;}}{k}_{\mathrm{r\&beta;}}\left(S_l\right)}{{\mathrm{\&mu;}}_{\mathrm{\&beta;}}}(\&dtri;p_{\mathrm{\&beta;}}-{\mathrm{\&rho;}}_{\mathrm{\&beta;}}g)---\left(4\right)$

其中，K为固有渗透率；k_rβ为β相的相对渗透率，与液相饱和度S_l有关；μ_β为β相的粘滞性系数；p_β为β相的压力；为β相的压力梯度；g为重力加速度矢量。

毛细压力p_c为液相压力p_l与气相压力p_g之间的差值：

p_c＝p_l-p_g                (5)

模型的本构关系包括毛细压力-饱和度关系与相对渗透率饱和度关系。其中关于毛细压力-饱和度关系曲线采用van Genuchten模型(VG模型)表示：

p_c＝-p₀[(S^*)^-1/λ-1]^1-λ(-p_max≤p_c≤0)    (6)

式中，p₀为进气值；λ为模型拟合参数，与土体均匀程度有关；S^*为有效饱和度，可表示为S^*＝(S_l-S_lr)/(S_ls-S_lr)，S_l为液相饱和度，S_lr为残余液相饱和度，S_ls为饱和液相饱和度。

关于相对渗透率-饱和度的关系，其中液相相对渗透率k_rl采用van Genuchten-Mualem模型(VG-M模型)表征：

$k_{\mathrm{rl}}=\left\{\begin{array}{cc}\sqrt{S^{*}}{[1-{(1-{\left(S^{*}\right)}^{1/\mathrm{\&lambda;}})}^{\mathrm{\&lambda;}}]}^2& (S_l<S_{\mathrm{ls}})\\ 1& (S_l\&GreaterEqual;S_{\mathrm{ls}})\end{array}\right.---\left(7\right)$

气相相对渗透率k_rg采用Corey(科里)提出的表达式：

$k_{\mathrm{rg}}=\left\{\begin{array}{cc}1-k_{\mathrm{rl}}& (S_{\mathrm{gr}}=0)\\ {(1-\hat{S})}^2(1-{\hat{S}}^2)& (S_{\mathrm{gr}}>0)\end{array}\right.---\left(8\right)$

式中，S_gr为残余气相饱和度。

(2)模型求解：为了保证模型计算结果的精确度，竖直方向上，模型网格尺寸在饱和区为0.05m，非饱和区为0.01m；水平方向上，网格尺寸为1m。以TOUGH2/EOS3为工具，空间上采用积分形式的有限差分方法(IFDM)进行离散，时间上采用一阶向后差分的全隐式方法进行离散，见式(9)；模型的线性化采用Newton-Raphson(牛顿-拉斐逊)迭代方法，最后得到大型稀疏系数矩阵的线性方程组，见式(10)：

$R_n^{\mathrm{\&kappa;},k+1}=M_n^{\mathrm{\&kappa;},k+1}-M_n^{\mathrm{\&kappa;},k}-\frac{\mathrm{\&Delta;t}}{V_n}\{\underset{m}{\mathrm{\&Sigma;}}{A}_{\mathrm{nm}}{F}_{\mathrm{nm}}^{\mathrm{\&kappa;},k+1}+V_n{q}_n^{\mathrm{\&kappa;},k+1}\}---\left(9\right)$

式中，引入了组分κ的余量 分别为M^κ，q^κ在单元体积V_n上的平均值；m为与单元n相邻的所有单元；A_nm为单元n和m相邻的交界面；为F^κ在面A_nm上内法线方向的平均值；上标k和k+1分别表示两相邻的时间步长指标；△t为时间步长。

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

$-\underset{i}{\mathrm{\&Sigma;}}\frac{\&PartialD;R_n^{\mathrm{\&kappa;},k+1}}{\&PartialD;x_i}{|}_p(x_{i,p+1}-x_{i,p})=R_n^{\mathrm{\&kappa;},k+1}\left(x_{i,p}\right)---\left(10\right)$

(3)模型边界条件确定：首先计算稳定渗流情况，作为降雨入渗的初始条件。在稳定渗流情况下，其边界条件为：地下水位以下的侧面边界为已知水头边界，即p_l＝p_atm+ρ_wg(1.8-Z)，T＝25℃；地表的边界条件为空气边界，即p_g＝p_atm，T＝25℃；区域底部及其它边界为不透水边界。模型的初始条件设为土体水相饱和，主要变量p_l＝1.013×10⁵pa，T＝25℃，在此定解条件下，运行模型直至达到毛细压力-重力平衡状态，即为降雨模型的初始条件。

以稳定渗流状态作为降雨入渗模拟的初始条件，研究区域的侧面与底部为不透水边界，地表为空气边界。区域顶部根据源汇项公式m_r(t)＝ρ_wA_eQ_w(t)，施加45mm/h的雨强，持续时间为1h。

(4)根据非饱和土的水-气二相流模型，分析降雨条件下，非饱和区中的渗流场及考虑Lisse现象的浅层地下水位的波动情况：

由于研究土体为砂土，因此其毛细管厚度很小，可以忽略不计。暂态饱和区是一种由非饱和状态过渡到完全饱和状态的区域，此时毛细压力仍然存在；孔隙水压力((p_l-p_atm)/ρ_wg)和孔隙气压力((p_g-p_atm)/ρ_wg)均为正值，且孔隙水压力小于孔隙气压力。图3为降雨过程中不同时刻的液相饱和度、毛细压力和孔隙压力的分布情况。降雨开始后，湿润峰下移，并在降雨结束时达到0.19m。毛细压力的变化趋势与饱和度变化相似，在湿润区负毛细压力减小，湿润区下方仍保持初始值。由于入渗水流下移，湿润峰下方的空气受到压缩，因此非饱和区的孔隙气压力增加，并在降雨结束时达到最大值，0.13m。接近地表处的孔隙水压力显著增加，且大于零，说明此处形成了暂态饱和区；在湿润峰下方，由于孔隙气压力的增加，孔隙水压力偏离了初始值，且偏离值的大小等于相应的孔隙气压力的变化。

图4为降雨结束后不同时刻的液相饱和度、毛细压力和孔隙压力的分布情况。降雨结束之后，湿润峰继续下移，在6h至10h之间到达初始地下水位。在湿润锋到达初始地下水位之前，由于非饱和区中的气体可以溢出地表，因此初始地下水位以上的孔隙气压力减小为零，且湿润区下方的孔隙水压力恢复至初始值；而由于蒸发作用，接近地表处的孔隙水压力小于零，暂态饱和区恢复到初始的非饱和状态。当湿润峰到达初始地下水位时，其上方的孔隙气压力开始逐渐增加；孔隙水压力开始偏离初始值；同时由于封闭气泡，毛细压力仍然存在，即在地下水位上方，形成暂态饱和区。

图5(a)为观测点B处孔隙压力、毛细压力和气体流速随时间的变化情况。在降雨过程中，孔隙气压力一开始迅速增大，之后由于地下水的阻碍缓慢增加，并在降雨结束时达到最大值0.13m。孔隙水压力的变化趋势与孔隙气压力基本一致。由于湿润峰未到达B点，因此其毛细压力和饱和度均保持不变。降雨一结束，孔隙气压力和孔隙水压力均迅速减小，分别接近于零和初始值。之后孔隙水压力与毛细压力的变化趋势一致。在降雨开始后的4小时，B点的负毛细压力开始减小，说明此时湿润峰到达该处，饱和度增加。由B点竖直方向的气体流速随时间的变化情况来看，其中气体向上为正，向下为负，降雨一开始，由于雨水下渗，B点气体受到压缩，向下进入土体内部，导致孔隙气压力迅速增加；之后由于受到地下水的限制，气体流速快速减小，因此孔隙气压力缓慢增加；当降雨结束时，气流反转，气体溢出地表，造成B点孔隙气压力减小为零；当湿润峰到达B点时，气流向上，这是因为气体受到水流驱替。在整个入渗过程中，B点的孔隙水压力均为负值，说明B点一直处于非饱和状态。

图5(b)观测点C处孔隙压力、毛细压力和液相饱和度随时间的变化情况。降雨过程中，C点孔隙气压力与B点变化情况相似。其孔隙水压力由初始的负值逐渐增加，并在降雨结束时达到最大值0.106m，说明C点水分获得气流提供的能量，克服了毛管力，能够自由流动并传递水压力。而C点的饱和度和毛细压力在降雨过程中基本保持不变。当降雨结束时，C点孔隙气压力和孔隙水压力均迅速减小，分别接近于零和初始值。在7h时，C点饱和度增加，说明湿润锋到达该点，其孔隙水压力和孔隙气压力也开始逐渐增加，但孔隙水压力始终小于孔隙气压力，毛细压力存在。在18h时，孔隙水压力变为正值，浸润线上升。当湿润锋到达C点之后，且其孔隙水压力大于零时，C点处于暂态饱和状态。

图5(c)观测点D处孔隙压力和毛细压力随时间的变化情况。由于D点位于初始地下水位处，因此其初始的孔隙气压力和孔隙水压力均为零，且在整个入渗过程中，毛细压力保持为零，孔隙水压力等于孔隙气压力。降雨过程中，D点孔隙气压力和孔隙水压力增加，且在降雨结束时达到最大值，0.13m。降雨结束之后，孔隙气压力和孔隙水压力均恢复到零。在7h时，由于湿润锋到达该点，两者开始逐渐增加。在整个入渗过程中，D点始终处于饱和状态，孔隙水压力的增加有两点原因：降雨过程中的气流推动和获得真正补给。

图6为Lisse现象的原理解释，地下水位以下的土体和观测井可看作连通器。观测点D位于土体中的地下水位处，点E位于井中，且与D点相平。根据连通器原理，D点的液相压力等于E点的液相压力即：

$p_l^D=p_l^E---\left(11\right)$

在降雨开始之前，井中的水位应与地下水位相平，D点和E点的液相压力为：

$p_l^D=p_l^E=p_{\mathrm{atm}}---\left(12\right)$

在降雨过程中，由于初始地下水位和湿润峰之间的空气受到压缩，因此D点的气相压力增加(如图5c)。对于任一时刻，D点增加的气相压力为水分提供了额外的能量，导致浸润线上升。此时，D点的水相压力等于其气相压力：

$p_l^D=p_g^E=p_{\mathrm{atm}}+\mathrm{\&Delta;}{p}_g^D---\left(13\right)$

根据连通器原理，土体中的孔隙水将会流入井中，直至井中水位与浸润线相平。此时，假设井中水位上升高度为△H，则E点的水相压力为：

$p_l^E=p_{\mathrm{atm}}+{\mathrm{\&rho;}}_w\mathrm{g\&Delta;H}---\left(14\right)$

由方程(11)，水位上升高度△H可表示为：

$\mathrm{\&Delta;H}=\mathrm{\&Delta;}{p}_g^D/{\mathrm{\&rho;}}_{w}g---\left(15\right)$

此时，湿润峰并未到达初始地下水位，水位上升则表示发生了Lisse现象。而且，Lisse现象是由地下水位处的孔隙气压力控制的。当湿润峰到达初始地下水位时，地下水获得真正的补给，浸润线再次上升。

图7为浸润线上升高度随时间的变化情况，其中初始地下水位为基准线。浸润线的变化趋势和初始地下水位处的孔隙气压力变化趋势一致：降雨过程中，首先快速增加，之后缓慢增加，并在降雨结束时达到最大值0.13m。而此时湿润峰并未到达初始地下水位处，说明水位的快速上升是由Lisse现象引起的。当降雨结束后，浸润线快速回到初始地下水位，之后在7h处逐渐上升，说明此时地下水获得真正的补给。因此，当用浸润线代表地下水位时，水位波动包括两部分原因：Lisse现象和真正补给。

综上所述，对本实施例总结如下：

(1)针对降雨条件下，考虑Lisse现象的浅层地下水位的波动，采用孔隙水压力等于零的线(面)来代表地下水位，并利用数值模拟的方法研究了水位的变化情况。首先利用非饱和土的水-气二相流模型，模拟了稳定渗流情况，作为降雨入渗的初始条件，之后分析了降雨条件下，非饱和土中渗流场及水位的变化情况。

(2)对非饱和区中渗流场的分析表明，降雨过程中，非饱和区的气体受到湿润峰下移的压缩，孔隙气压力增加；接近地表处的孔隙水压力大于零，形成暂态饱和区；在湿润峰下方，由于孔隙气压力的增加，孔隙水压力偏离了初始值。降雨结束之后，非饱和区中的气体可以自由溢出，孔隙气压力减小为零，湿润峰下方的孔隙水压力也恢复到初始值。当湿润峰到达初始地下水位时，初始地下水位上方形成暂态饱和区，浸润线上升。

(3)降雨过程中，地下水位首先快速增加，之后缓慢增加，并在降雨结束时达到最大值。而此时湿润峰并未到达初始地下水位处，说明水位的快速上升是由Lisse现象引起的。当降雨结束后，浸润线快速回到初始地下水位，之后逐渐上升，说明此时地下水获得真正的补给。因此，地下水位的波动有两部分原因：Lisse现象和真正的补给，在地下水开采工作中，没能正确地识别Lisse现象，可能造成地下水过量开采的后果。然而，由于土体材料特性、降雨强度和地形特点的不同，Lisse现象的形式多种多样，水位回退也不尽相同，因此需要相应开展很多的研究工作。

Claims (2)

1.一种考虑Lisse现象的浅层地下水位预测方法，其特征是，包括如下步骤：

步骤一：建立非饱和土的水-气二相流模型，包括基本控制方程及本构关系，假定系统一直处于恒温状态，具体如下：

模型的基本控制方程为：

$\frac{{\&PartialD;M}^{\mathrm{\&kappa;}}}{\&PartialD;t}=-\mathrm{div}{F}^{\mathrm{\&kappa;}}+q^{\mathrm{\&kappa;}}---\left(1\right)$

式中，M^κ表示包括空气a和水w的κ组分的累积质量密度；F^κ为κ组分的平流流量；q^κ为κ组分的源汇项，t表示时间；

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

式中，β表示液相l或气相g，其中液相包括液态孔隙水及溶解的空气，气相包括干燥空气和水蒸气；为κ组分占β相的质量百分数；φ为孔隙率；S_β为β相的饱和度；ρ_β为β相的密度；

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

$F^{\mathrm{\&kappa;}}=\underset{\mathrm{\&beta;}}{\mathrm{\&Sigma;}}{X}_{\mathrm{\&beta;}}^{\mathrm{\&kappa;}}{F}_{\mathrm{\&beta;}}---\left(3\right)$

式中，F_β为β相的平流流量，遵循达西定律，其表达式为：

$F_{\mathrm{\&beta;}}=-K\frac{{\mathrm{\&rho;}}_{\mathrm{\&beta;}}{k}_{\mathrm{r\&beta;}}\left(S_l\right)}{{\mathrm{\&mu;}}_{\mathrm{\&beta;}}}({\&dtri;p}_{\mathrm{\&beta;}}-{\mathrm{\&rho;}}_{\mathrm{\&beta;}}g)---\left(4\right)$

其中，K为固有渗透率；k_rβ为β相的相对渗透率，与液相饱和度S_l有关；μ_β为β相的粘滞性系数；p_β为β相的压力；为β相的压力梯度；g为重力加速度矢量；

毛细压力p_c为液相压力p_l与气相压力p_g之间的差值：

p_c＝p_l-p_g    (5)

模型的本构关系包括毛细压力-饱和度关系与相对渗透率-饱和度关系，其中关于毛细压力-饱和度关系曲线采用van Genuchten模型简称VG模型表示：

$p_c={-p}_0{[{\left(S^{*}\right)}^{-1/\mathrm{\&lambda;}}-1]}^{1-\mathrm{\&lambda;}},({-p}_{\max }\&le;p_c\&le;0)---\left(6\right)$

式中，p₀为进气值；λ为模型拟合参数，与土体均匀程度有关；S^*为有效饱和度，表示为S^*＝(S_l-S_lr)/(S_ls-S_lr)，S_l为液相饱和度，S_lr为残余液相饱和度，S_ls为饱和液相饱和度；

关于相对渗透率-饱和度的关系，其中液相相对渗透率k_rl采用van Genuchten-Mualem模型简称VG-M模型表征：

$k_{\mathrm{rl}}=\left\{\begin{array}{cc}\sqrt{S^{*}}{[1-{(1-{\left(S^{*}\right)}^{1/\mathrm{\&lambda;}})}^{\mathrm{\&lambda;}}]}^2& (S_l<S_{\mathrm{ls}})\\ 1& (S_l\&GreaterEqual;S_{\mathrm{ls}})\end{array}\right.---\left(7\right)$

气相相对渗透率k_rg采用科里(Corey)提出的表达式：

$k_{\mathrm{rg}}=\left\{\begin{array}{cc}1-k_{\mathrm{rl}}& (S_{\mathrm{gr}}=0)\\ {(1-\hat{S})}^2(1-{\hat{S}}^2)& (S_{\mathrm{gr}}>0)\end{array}\right.---\left(8\right)$

式中， $\hat{S}=(S_l-S_{\mathrm{lr}})/(1-S_{\mathrm{lr}}-S_{\mathrm{gr}}),$ S_gr为残余气相饱和度；

步骤二：模型求解：以TOUGH2/EOS3为工具，将变量分为主要变量和次要变量，其中，次要变量可通过主要变量求得，空间上采用积分形式的有限差分方法(IFDM)进行离散，时间上采用一阶向后差分的全隐式方法进行离散，模型的线性化采用牛顿-拉斐逊(Newton-Raphson)迭代方法，最后得到大型稀疏系数矩阵的线性方程组；

步骤三：模型边界条件确定：边界条件包括狄利克雷(Dirichlet)边界条件和黎曼(Neumann)边界条件，其数学处理方法如下：

(1)Dirichlet边界条件

将边界条件单元的体积设为1×10⁵⁰m³，包括空气边界和已知水头边界：

①对于空气边界，仅有气相，其主要变量为p_g，和T，其中p_g为气相压力，为空气占气相的质量百分数，T为系统温度；

②对于已知水头边界，仅有液相状态，其主要变量为p_l，和T，其中p_l为液相压力，为液相中空气所占的质量百分数，T为系统温度；

(2)Neumann边界条件

Neumann边界条件描述的是系统与外界的流量交换情况，流入为正，流出为负；

①对于降雨入渗条件，通过m_r(t)＝ρ_wA_eQ_w(t)，设置一个适当大小的源汇项来实现，其中，源汇项m_r表示单位时间内通过土体表面法向的水分质量，流入为正，流出为负；A_e为土体的有效面积；Q_w为降雨强度，ρ_w为水的密度；

②对于不透水边界或不入流边界，看作一类特殊的Neumann边界条件，边界上的流量设为零；

步骤四：利用非饱和土的水-气二相流模型，以稳定渗流情况作为降雨入渗的初始条件，模拟分析降雨条件下，考虑Lisse现象的浅层地下水位的波动情况：

(1)渗流场变化：降雨过程中，非饱和区的气体受到湿润峰下移的压缩，孔隙气压力(p_g-p_atm)/ρ_wg增加；接近地表处的孔隙水压力(p_l-p_atm)/ρ_wg大于零，形成暂态饱和区；在湿润峰下方，由于孔隙气压力的增加，孔隙水压力偏离了初始值，浸润线(p_l＝0的线或面)上升，降雨结束之后，非饱和区中的气体可以自由溢出，孔隙气压力减小为零，湿润峰下方的孔隙水压力也恢复到初始值，当湿润峰到达初始地下水位时，初始地下水位上方形成暂态饱和区，浸润线再次上升；

(2)地下水位变化：降雨过程中，地下水位首先快速增加，之后缓慢增加，并在降雨结束时达到最大值，而此时湿润峰并未到达初始地下水位处，说明水位的快速上升是由Lisse现象引起的；当降雨结束后，浸润线快速回到初始地下水位，之后逐渐上升，说明此时地下水获得真正的补给，因此，地下水位的波动有两部分原因：Lisse现象和真正的补给，在地下水开采工作中，没能正确地识别Lisse现象，可能会造成地下水的过量开采。

2.如权利要求1所述的考虑Lisse现象的浅层地下水位预测方法，其特征是，步骤二进一步具体为：

步骤二进一步具体为：

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

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

$\frac{d}{\mathrm{dt}}{\&Integral;}_{V_n}{M}^{\mathrm{\&kappa;}}{\mathrm{dV}}_n={\&Integral;}_{{\mathrm{\&Gamma;}}_n}{F}^{\mathrm{\&kappa;}}\&CenterDot;\mathrm{nd}{\mathrm{\&Gamma;}}_n+{\&Integral;}_{V_n}{q}^{\mathrm{\&kappa;}}{\mathrm{dV}}_n---\left(9\right)$

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

引入体积平均值：

${\&Integral;}_{V_n}{M}^{\mathrm{\&kappa;}}{\mathrm{dV}}_n=V_n{M}_n^{\mathrm{\&kappa;}}---\left(10\right)$

${\&Integral;}_{V_n}{q}^{\mathrm{\&kappa;}}{\mathrm{dV}}_n=q_n^{\mathrm{\&kappa;}}{V}_n---\left(11\right)$

式中，M^κ，q^κ在单元体积V_n上的平均值；

表面单元Γ_n上的面积分可近似为其所包含的各个表面A_nm上面积分的平均值之和：

${\&Integral;}_{{\mathrm{\&Gamma;}}_n}{F}^{\mathrm{\&kappa;}}\&CenterDot;\mathrm{nd}{\mathrm{\&Gamma;}}_n=\underset{m}{\mathrm{\&Sigma;}}{A}_{\mathrm{nm}}{F}_{\mathrm{nm}}^{\mathrm{\&kappa;}}---\left(12\right)$

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

将式(10)、(11)和(12)代入到式(9)中，得到一组关于时间的一阶微分方程组

$\frac{{\mathrm{dM}}_n^w}{\mathrm{dt}}=\frac{1}{V_n}\underset{m}{\mathrm{\&Sigma;}}{A}_{\mathrm{nm}}{F}_{\mathrm{nm}}^{\mathrm{\&kappa;}}+q_n^{\mathrm{\&kappa;}}---\left(13\right)$

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

采用一阶向后差分方法，对式(13)进行时间上的离散，得到任一单元的全隐式非线性方程组：

$R_n^{\mathrm{\&kappa;},k+1}=M_n^{\mathrm{\&kappa;},k+1}-M_n^{\mathrm{\&kappa;},k}-\frac{\mathrm{\&Delta;t}}{V_n}\{\underset{m}{\mathrm{\&Sigma;}}{A}_{\mathrm{nm}}{F}_{\mathrm{nm}}^{\mathrm{\&kappa;},k+1}+V_n{q}_n^{\mathrm{\&kappa;},k+1}\}---\left(14\right)$

式中，引入了组分κ的余量上标k和k+1分别表示两相邻的时间步长指标；△t为时间步长；

(3)Newton-Raphson迭代方法

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

$-\underset{i}{\mathrm{\&Sigma;}}\frac{{\&PartialD;R}_n^{\mathrm{\&kappa;}.k+1}}{{\&PartialD;x}_i}{|}_p(x_{i,p+1}-x_{i,p})=R_n^{\mathrm{\&kappa;},k+1}\left(x_{i,p}\right)---\left(15\right).$