一种对流-扩散传质过程的数值模拟方法

摘要

本发明涉及一种对流-扩散传质过程的数值模拟方法，属于地下水污染问题数值模拟技术领域。该方法以流网为计算单元进行地下水溶质运移的数值分析，根据流网的特殊性质，利用流网的流线与等势线对求解域进行离散。在流网单元中，由两条流线构成一个封闭的管状曲面称为流管；流管内、外的物质交换通过浓度梯度作用下的分子扩散发生，因此在对流－扩散传质过程中，流管之间的物质交流将不再受对流的影响，而仅依靠浓度梯度下的分子扩散，即不再需要弥散系数这一难以准确测得的参数。本发明提供的一种对流-扩散传质过程的数值模拟方法，以流网为计算单元进行对流-扩散传质过程的数值模拟，解决了已有对流－扩散问题数值求解方法中所需要的水动力弥散参数难以准确确定这一难题。

说明书

技术领域

本发明属于地下水污染问题数值模拟技术领域，涉及一种对流-扩散传质过程的数值模拟 方法。

背景技术

经济的快速发展伴随着日益严重的水污染问题。垃圾填埋场滤出液渗漏、海水入侵、核 废物及生产生活废水严重威胁人们的用水安全。污染物在对流与扩散共同作用下随地下水运 动是污染物迁移的主要途径。地下水污染预测与防治需要建立描述污染物运移过程的数学模 型，其中弥散参数及渗流场的确定是污染物运移模型建立与求解的关键。然而，由于渗流场 及污染物运移过程的复杂性，水动力弥散参数难以精确测定。水动力弥散过程由分子扩散与 机械弥散组成，与分子扩散系数这一可测定的物质固有参数不同，室内或野外测量得到的机 械弥散系数受多种因素影响，具有明显的尺寸效应，参数的准确性及工程应用的广泛性非常 有限。研究人员曾深入探索了影响水动力弥散参数不同因素，尤其是尺寸效应对水动力弥散 参数的影响，但成果尚未具有普遍性的指导意义。数值模拟是进行对流－扩散问题求解的主 要方法，而水动力弥散参数是对流－扩散过程数值求解的关键参数，在水动力弥散参数难以 准确测定的前提下，对流－扩散问题数值求解结果的精度及工程应用受到很大局限。

发明内容

有鉴于此，本发明的目的在于提供一种对流-扩散传质过程的数值模拟方法，该方法根据 流网的特殊性质，以流网为计算单元进行地下水溶质运移的数值分析，在分析时不再需要弥 散系数，从而解决了已有对流－扩散问题数值求解方法中所需要的水动力弥散参数难以准确 确定这一难题。

为达到上述目的，本发明提供如下技术方案：

一种对流-扩散传质过程的数值模拟方法，该方法包括以下步骤：

S1：通过以下公式，求解稳态平面渗流场势函数，

$\frac{\partial \left(k_x\frac{\partial H}{\partial x}\right)}{\partial x}+\frac{\partial \left(k_y\frac{\partial H}{\partial y}\right)}{\partial y}=0,$

其中，k_x，k_y为求解域中x，y方向的渗透系数；H为水头势函数；

S2：通过以下公式，求解稳态平面渗流场流函数，

ψ＝∫u_xdy-u_ydx

其中，u_x、u_y为x、y方向的渗透流速；

S3：利用插值法，得到水头势函数与流函数的等值线，并求解等值线的交点坐标；

S4：利用直线段将相邻的等值线交点相连，将求解域划分为四边形网格，所述四边形网 络即为流网单元格；

S5：建立流网单元格内的质量守恒方程；

对单元P，其周围相邻的四个单元分别为W、S、E、N，单元P内的溶质质量守恒方程 为

$a_P{C}_P=a_W{C}_W+a_E{C}_E+a_S{C}_S+a_N{C}_N+a_P^0{C}_P^0+S_u,$

其中，C_P、C_W、C_E、C_S、C_N为t+Δt时刻控制体积P、W、E、S、N内的溶液浓度；为 t时刻控制体积P内的污染物浓度；S_u为与源汇项有关的系数；a_P、a_W、a_E、a_S、a_N、为离散方程系数；

S6：计算单元内溶质质量守恒方程的系数a_W、a_E、a_S、a_N；采用中心差分时各系数为

$\left(\begin{array}{c}{a}_W=D_w+\frac{E_w}{2};a_E=D_e-\frac{F_e}{2}\\ a_S=D_s+\frac{F_s}{2};a_N=D_n-\frac{F_n}{2}\end{array}\right),$

其中，F为通过控制容积边界界面的对流量；D为界面上扩散阻力的系数；根据流网的 性质，n、s边界的对流量为0；且w、e边界上流体的流量都等于流管内流体的流量q；各界 面上的F与D的表达式为：

F_w＝qC_w；F_e＝qC_e；F_n＝F_s＝0

$D_w=\frac{D_m{A}_w}{\delta x_{\mathrm{WP}}};D_e=\frac{D_m{A}_e}{\delta x_{\mathrm{EP}}};D_s=\frac{D_m{A}_s}{\delta y_{\mathrm{SP}}};D_n=\frac{D_m{A}_n}{\delta y_{\mathrm{NP}}}$

△F＝F_e-F_w+F_n-F_s＝q(C_e-C_w)

其中，q为流管内的流量通量；D_m为扩散系数；A为界面在垂直于扩散方向上的投影面 积；C_w，C_e为控制体积P的w、e两界面上的溶液浓度；δx、δy为控制体积节点与相邻控制 体积节点的距离；

S7：通过以下公式，计算单元内溶质质量守恒方程的系数

$a_P^0=\frac{1}{2\mathrm{\Delta t}}\underset{i=0}{\overset{3}{\Sigma }}(x_i{y}_{i+1}-x_{i+1}{y}_i),$

其中，x_i、y_i为控制容积P四个角点的坐标；Δt为数值模拟的时间步长；

S8：通过以下公式，计算单元内溶质质量守恒方程的系数S_u、S_P，

$\bar{S}\mathrm{\Delta V}=S_u+S_p{C}_P,$

其中，为源汇项在控制容积中的平均值；

S9：通过以下公式，计算单元内溶质质量守恒方程的系数a_P，

$a_P=a_W+a_E+a_S+a_N+a_P^0+\mathrm{\Delta F}-S_P;$

S10：对每个计算单元都列出溶质质量守恒方程，并求解由各单元质量守恒方程组成的线 性方程线。

进一步，所述S2中的u_x、u_y的值通过求解水头分布及渗透系数得到。

进一步，所述S4中流网单元格的疏密由S3中等值线的数目控制。

进一步，所述S6中的C_w，C_e的值利用线性插值法获得。

本发明的有益效果在于：本发明提供的一种对流-扩散传质过程的数值模拟方法，该方法 根据流网的特点，以流网为计算单元进行地下水溶质运移的数值分析具有以下优势：①由两 条流线构成一个封闭的管状曲面称为流管，由于流线不会相交，因此流管内、外的流体不会 穿越管壁，即流管内的流体不与外界有流体物制交流，因此，流管内的流动为沿着流线方向 的一维流动。在溶质运移过程中，流管之间的物质交流将不再受对流的影响，而仅依靠浓度 梯度下的分子扩散。因此在分析时不再需要弥散系数。②流管内流量恒定，因此高流速的位 置流网单元较密，低流速处流网单元疏。高流速处一般都是溶质运移分析时离散单元需要加 密的位置，而低流速处疏松的单元也可以节省计算资源，即流网单元格可以较好的自适应地 层及渗流场非均质对单元疏密的要求。

附图说明

为了使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明作进一步的 详细描述，其中：

图1为本发明所述方法的流程图；

图2为流网单元示意图；

图3为利用本方法进行的某渗流场基于流网单元的网格划分；

图4为利用本方法进行的数值模拟结果与室内试验结果的对比。

具体实施方式

下面将结合附图，对本发明的优选实施例进行详细的描述。

针对对流－扩散问题数值求解过程中所需要的水动力弥散参数难以准确确定这一难题， 本发明提出一种复杂地质条件下以流网为计算单元对对流-扩散传质过程进行计算机数值模 拟的方法。该方法根据流网的特殊性质，利用流网的流线与等势线对求解域进行离散。在流 网单元中，由两条流线构成一个封闭的管状曲面称为流管。流管内、外的物质交换通过浓度 梯度作用下的分子扩散发生，因此在对流－扩散传质过程中，流管之间的物质交流将不再受 对流的影响，而仅依靠浓度梯度下的分子扩散，即不再需要弥散系数这一难以准确测得的参 数。

如图1所示，该方法包括以下步骤：

S1：求解稳态平面渗流场势函数。势函数由求解一定边值条件下的稳态渗流场方程得到， 如下式：

$\frac{\partial \left(k_x\frac{\partial H}{\partial x}\right)}{\partial x}+\frac{\partial \left(k_y\frac{\partial H}{\partial y}\right)}{\partial y}=0,$

式中，k_x，k_y为求解域中x，y方向的渗透系数；H为水头势函数。

S2：求解稳态平面渗流场流函数。在得到水头分布后，流函数可以根据下式求得：

ψ＝∫u_xdy-u_ydx，

式中，u_x、u_y为x、y方向的渗透流速，可由求解得到的水头分布及渗透系数求得。

S3：利用插值法，得到水头势函数与流函数的等值线，并求得等值线的交点坐标。

S4：利用直线段将相邻的等值线交点相连，将求解域划分为四边形网格，所述四边形网 络即为流网单元格；流网单元格的疏密由S3中等值线的数目控制。

S5：建立流网单元格内的质量守恒方程，对单元P，其周围相邻的四个单元分别为W、S、 E、N，那么单元P内的溶质质量守恒方程为：

$a_P{C}_P=a_W{C}_W+a_E{C}_E+a_S{C}_S+a_N{C}_N+a_P^0{C}_P^0+S_u,$

式中，C_P、C_W、C_E、C_S、C_N为t+Δt时刻控制体积P、W、E、S、N内的溶液浓度；为t时刻控制体积P内的污染物浓度；S_u为与源汇项有关的系数；a_P、a_W、a_E、a_S、a_N、 为离散方程系数。

S6：计算单元内溶质质量守恒方程的系数a_W、a_E、a_S、a_N，采用中心差分时各系数为：

$\left(\begin{array}{c}{a}_W=D_w+\frac{E_w}{2};a_E=D_e-\frac{F_e}{2}\\ a_S=D_s+\frac{F_s}{2};a_N=D_n-\frac{F_n}{2}\end{array}\right),$

式中，F为通过控制容积边界界面的对流量；D为界面上扩散阻力的系数。根据流网的 性质，n、s边界的对流量为0；且w、e边界上流体的流量都等于流管内流体的流量q。各界 面上的F与D的表达式为：

F_w＝qC_w；F_e＝qC_e；F_n＝F_s＝0

$D_w=\frac{D_m{A}_w}{\delta x_{\mathrm{WP}}};D_e=\frac{D_m{A}_e}{\delta x_{\mathrm{EP}}};D_s=\frac{D_m{A}_s}{\delta y_{\mathrm{SP}}};D_n=\frac{D_m{A}_n}{\delta y_{\mathrm{NP}}}$

△F＝F_e-F_w+F_n-F_s＝q(C_e-C_w)

式中，q为流管内的流量通量；D_m为扩散系数；A为界面在垂直于扩散方向上的投影面 积；C_w，C_e为控制体积P的w、e两界面上的溶液浓度，利用线性插值法获得；δx、δy为 控制体积节点与相邻控制体积节点的距离。

S7：计算单元内溶质质量守恒方程的系数

$a_P^0=\frac{1}{2\mathrm{\Delta t}}\underset{i=0}{\overset{3}{\Sigma }}(x_i{y}_{i+1}-x_{i+1}{y}_i),$

x_i、y_i为控制容积P四个角点的坐标；Δt为数值模拟的时间步长。

S8：计算单元内溶质质量守恒方程的系数S_u、S_P：

$\bar{S}\mathrm{\Delta V}=S_u+S_p{C}_P$

式中，为源汇项在控制容积中的平均值。

S9：计算单元内溶质质量守恒方程的系数a_P：

$a_P=a_W+a_E+a_S+a_N+a_P^0+\mathrm{\Delta F}-S_P.$

S10：对每个计算单元都列出溶质质量守恒方程，并求解由各单元质量守恒方程组成的线 性方程线。

图2本方法利用流网划分的单元示意图。在渗流场中，由于地层的不均匀性与各向异性， 流线与等势线可能正交也可能斜角。图2为流网中等势线与流线划分网格局部的一般形式。 其中控制单元为P，与其相邻的4个控制单元为W、S、E、N。控制单元为P的4个角点， 即等势线与流线的交点，分别为ω_P、η_P、λ_P、ε_P，4个角点在直角坐标系中的坐标分别为(x₀,y₀)、 (x₁,y₁)、(x₂,y₂)、(x₃,y₃)。控制体积P的节点坐标为：

$\left(\begin{array}{c}{c}_{\mathrm{Px}}=\frac{\underset{i=0}{\overset{3}{\Sigma }}(x_i^P+x_{i+1}^P)(x_i^P{y}_{i+1}^P-x_{i+1}^P{y}_i^P)}{3\underset{i=0}{\overset{n-1}{\Sigma }}(x_i^P{y}_{i+1}^P-x_{i+1}^P{y}_i^P)}\\ c_{\mathrm{Py}}=\frac{\underset{i=0}{\overset{3}{\Sigma }}(y_i^P+y_{i+1}^P)(x_i^P{y}_{i+1}^P-x_{i+1}^P{y}_i^P)}{3\underset{i=0}{\overset{n-1}{\Sigma }}(x_i^P{y}_{i+1}^P-x_{i+1}^P{y}_i^P)}\end{array}\right)$

其它控制体积的节点坐标可以用相同方法得到，因此可以得到相邻控制体积节点间的距 离，如：

$\delta x_{\mathrm{WP}}=\sqrt{(c_{\mathrm{Px}}^2-c_{\mathrm{Wx}}^2)+(c_{\mathrm{Py}}^2-c_{\mathrm{Wy}}^2)}$

式中，(C_Wx,C_Wy)为控制体积W的节点坐标。控制体积节点与其它相邻控制体积节点间的 距离可用上式类似形式得到。

为求得垂直与扩散方向的界面面积，如A_w，首先令直线的斜率分别为k₁、k₂， 即：

$k_1=\frac{c_{\mathrm{Py}}-c_{\mathrm{Wy}}}{c_{\mathrm{Px}}-c_{\mathrm{Wx}}};k_2=\frac{y_3-y_0}{x_3-x_0}$

则直线间的夹角为：

${\alpha }_{\mathrm{WP}}=\arctan \frac{|k_1-k_2|}{1+k_1{k}_2}$

因此可得：

$A_w=\sin {\alpha }_{\mathrm{WP}}·L_{\overline{{ϵ}_P{\omega }_P}}=\frac{|k_1-k_2|·\sqrt{(x_0^2-x_3^2)+(y_0^2-y_3^2)}}{\sqrt{1+k_1^2{k}_2^2+k_1^2+k_2^2}}$

其余各界面投影面积可用上式类似形式获得。

图3为利用本方法进行的某渗流场基于流网单元的网格划分。矩形渗流场长度为30cm， 高度为15cm。上下边界为隔水边界，左边界为定水头入水边界，右边界为定水头出水边界。 在定水头差作用下渗流场发生稳定渗流。利用渗流场流网单元对渗流场进行离散，离散后的 有限体积边长相等，取Δx＝Δy＝5.0mm，即每根流管的宽度为5.0mm。

图4为图3所示模型的计算实例。在左右边界施加定水头使模型内发生一维稳定渗流， 渗透流速为4.2×10^-4cm/s。在入水边界注入相对浓度为1.0的氯化钠溶液，取时间步长Δt＝1.0s， 利用图3所示流网单元及本项发明所述方法对氯化钠在渗流场中的运移过程进行模拟。沿溶 质运移方向均匀设置3个测点监测氯化钠溶液浓度的变化，并与室内模型试验结果相对比， 可以看出本项发明所述方法与模型试验结果非常吻合。

最后说明的是，以上优选实施例仅用以说明本发明的技术方案而非限制，尽管通过上述 优选实施例已经对本发明进行了详细的描述，但本领域技术人员应当理解，可以在形式上和 细节上对其作出各种各样的改变，而不偏离本发明权利要求书所限定的范围。