[发明专利]一种地下水环境的建模及数值模拟方法

权利要求书

1.一种地下水环境的建模及数值模拟方法，其包括以下步骤：

1)前处理：设置两个属性文件，一是模拟区土壤属性文件，包括研究区边界、X与Y方向网格间距、区域土壤渗透系数、土壤有效孔隙度、污染物在此种土壤中的纵向和横向的弥散度、土壤非饱和带扩散系数和溶质吸附系数；二是地表高程与地下水水位属性文件；

应用IDL软件的空间插值函数KRIG2D()或TRI_SURF()对地表高程、地下水位进行空间网格插值，同时划分空间有限差分数值模拟网格，得到两个相同的差分网格矩阵：地表高程矩阵matrix1和地下水位矩阵matrix2；应用SIZE()函数确定研究区网格属性，包括单元网格数、单元网格行列数；应用随机数值生成函数FINDGEN()生成随机数，然后将网格计算单元的节点坐标生成序列数组；

2)二维水流模型计算

二维水流模型计算包括：溶质非饱和带迁移模拟的地下水埋深矩阵计算以及饱和带二维地下水流场计算部分；

利用前处理得到的地表高程矩阵与地下水位矩阵相减得非饱和带厚度h₀，即地下水位埋深矩阵matrixh0：

h₀：matrixh0＝matrix1-matrix2

该矩阵直接用于非饱和带溶质迁移模拟计算；

饱和带二维溶质迁移模型的水动力条件是地下水流场的分布，该流场由二维水流模型计算得到，应用二维水流模型计算地下水流场：

①计算单元节点流速

利用达西定律：

v=KIne]]>

其中，v为流速；K为渗透系数；I为水力坡度；n_e为有效孔隙度；

②计算单元边界面流速

求得节点流速，根据网格中相邻两个节点的流速，取两个节点流速的平均值，计算得到网格单元边界面上流速分量的IDL矩阵；

3)建立非饱和带一维溶质运移模型；

t1=ΔHvz---(1)]]>

其中，t₁为溶质在非饱和带中的运移总时间；ΔH为非饱和带厚度，ΔH＝h₀；v_z＝D/Δz表示假设非饱和带厚度为Δz＝1m时，污染物垂向扩散速度，D为垂向扩散系数，量纲为L²/T；

4)建立饱和带溶质迁移二维数值模型：

R∂C∂t1=∂∂x(Dxx∂C∂x)+∂∂x(Dxy∂C∂y)+∂∂y(Dyy∂C∂y)+∂∂y(Dyx∂C∂x)----(2)]]>

∂∂x(vxC)-∂∂y(vyC)+qsθCs-λ1C-λ2ρbθC‾]]>

其中，R为等温吸附系数，C为溶解浓度，单位为ML^-3，t₁为饱和带中溶质迁移数值模拟的时间，D_ij为弥散系数张量，单位为L²T^-1，x、y为计算距离L，v_x和v_y为地下水流速，单位为LT^-1，，q_s为源/汇处单位体积含水层流量(，单位为M³T^-1，θ为含水层孔隙度，C_s为源/汇项浓度，单位为ML^-1，λ₁为溶解相的反应速率常数，单位为T^-1，λ₂为吸附相的反应速率常数，单位为T^-1，ρ_b为空隙介质体积密度，单位为ML^-3，为吸附浓度，单位为MM^-1；

5)将公式(2)中的各项式进行近似表示：

在规则间距的节点中心网格中，公式(2)右边第一项在单元(i，j)近似为：

∂∂x(Dxx∂C∂x)≈1ΔxDxx(i,j+1/2)(Ci,j+1-Ci,j)-Dxx(i,j-1/2)(Ci,j-Ci,j-1)Δx---(3)]]>

公式(3)表示x方向上，由浓度梯度引起从x方向上进入单元(i，j)的净弥散通量；

公式(2)右边第二项近似为：

∂∂x(Dxy∂C∂y)]]>

≈1ΔxDxy(i,j+1/2)(Ci+1,j+1/2-Ci-1,j+1/2)-Dxy(i,j-1/2)(Ci+1,j-1/2-Ci-1,j-1/2)2Δy---(4)]]>

=14ΔxΔy[Dxy(i,j+1/2)(Ci+1,j+1+Ci+1,j-Ci-1,j+1-Ci-1,j)-]]>

Dxy(i,j-1/2)(Ci+1,j+Ci+1,j-1-Ci-1,j-Ci-1,j-1)]]]>

公式(4)表示由y方向的浓度梯度引起的从x方向进入单元(i，j)的净弥散通量；

同理，得到由y方向的浓度梯度引起的从y方向上进入单元(i，j)的净弥散通量，即公式(2)右边的第三项为：

∂∂y(Dyy∂C∂y)≈1ΔyDyy(i+1/2,j)(Ci+1,j-Ci,j)-Dyy(i-1/2,j)(Ci,j-Ci-1,j)Δy---(5)]]>

由x方向的浓度梯度引起的从y方向进入单元(i，j)的净弥散通量，即公式(2)右边的第四项，写为：

∂∂y(Dyx∂C∂x)]]>

≈1ΔyDyx(i+1/2,j)(Ci+1/2,j+1-Ci+1/2,j-1)-Dyx(i-1/2,j)(Ci-1/2,j+1-Ci-1/2,j-1)2Δx---(6)]]>

=14ΔyΔx[Dyx(i+1/2,j)(Ci+1,j+1+Ci,j+1-Ci+1,j-1-Ci,j-1)-]]>

Dyx(i-1/2,j)(Ci,j+1+Ci-1,j+1-Ci,j-1-Ci-1,j-1)]]]>

公式(2)右边第五项为从x方向进入单元(i，j)的净对流通量，近似为：

∂∂x(vxC)≈1Δx{vx(i,j+1/2)[(1-α)Ci,j+αCi,j+1]----(7)]]>

vx(i,j-1/2)[(1-α)Ci,j-1+αCi,j]}]]>

第六项为从y方向进入单元(i，j)的净对流通量，近似为：

∂∂y(vyC)≈1Δy{vy(i+1/2,j)[(1-α)Ci,j+αCi+1,j]----(8)]]>

vy(i-1/2,j)[(1-α)Ci-1,j+αCi,j]}]]>

其中，α为空间加权因子；中心加权，α等于0.5；采用上风加权方案：

公式(2)右边第七项是由源/汇引起的流入或流出单元(i，j)净质量通量；

公式(2)右边第八项和第九项是单元(i，j)内一级化学反应引起的质量损失或增加量；吸附浓度表示为溶解浓度的函数，具体取决于所用的吸附等温线；

公式(2)的左边是储存于单元(i，j)的质量变化率，近似为：

R∂C∂t≈RCi,jn+1-Ci,jnΔt---(9)]]>

其中延迟因子R由吸附等温线确定；

6)对以上各分量差分计算过程详细实现过程如下：

①弥散系数D：由边界流速分量及公式(10)，利用IDL矩阵运算计算弥散系数分量；

Dxx(i,j±1/2)=αLvx(i,j±1/2)2vi,j±1/2+αTvy(i,j±1/2)2vi,j±1/2+D*]]>

Dyy(i±1/2,j)=αLvy(i±1/2,j)2vi±1/2,j+αTvx(i±1/2,j)2vi±1/2,j+D*---(10)]]>

Dxy(i,j±1/2)=(αL-αT)vx(i,j±1/2)vy(i,j±1/2)vi,j±1/2]]>

Dyx(i±1/2,j)=(αL-αT)vx(i±1/2,j)vy(i±1/2,j)vi±1/2,j]]>

其中：D^*为有效分子扩散系数；

v*,*=vx(*,*)2+vy(*,*)2---(11)]]>

②弥散项：设置初始溶质浓度场，即确定污染事件发生位置及该物质的饱和溶解度作为模型的初始浓度分布矩阵C；由弥散项计算公式(3)、公式(4)、公式(5)及公式(6)，结合弥散系数计算结果矩阵，利用IDL矩阵运算计算各弥散项；

③对流项：由对流项计算公式(7)、公式(8)，结合单元格边界流速分量、土壤弥散度属性，利用IDL矩阵运算计算各对流项；

④公式(2)的第七项，源汇项：由于快速处理模型、源汇项设置为0；

⑤公式(2)的第八、第九项，化学反应项：根据各化学物质的等温吸附曲线获取；

⑥单元节点浓度：由公式(9)及各个分量矩阵计算结果，利用IDL矩阵运算计算单元节点后一时刻浓度场；

7)返回步骤5)中的计算过程，即实现一个时间段后，区域流场及溶质浓度场的计算；在步骤5)前加一个总时段的循环控制语句，即实现各时段区域流场及溶质浓度场的计算；

8)结果输出：经过上述步骤实时输出地表高程矩阵matrix1、地下水位矩阵matrix2、地下水埋深矩阵matrixh0、地下水流场、溶质浓度场和时间项。

2.如权利要求1所述的一种地下水环境的建模及数值模拟方法，其特征在于：所述步骤4)中，根据中心加权方案，从节点浓度获得界面浓度：

Ci+1,j+1/2=Ci+1,j+1+Ci+1,j2---(12)]]>