﻿ 利用熵格式计算地下水溶质运移方程
 文章快速检索
 同济大学学报(自然科学版)  2019, Vol. 47 Issue (8): 1175-1179.  DOI: 10.11908/j.issn.0253-374x.2019.08.014 0

### 引用本文

CHEN Rongsan, WANG Fanrui, ZOU Min. Entropy Scheme for Solute Transport Equation in Groundwater[J]. Journal of Tongji University (Natural Science), 2019, 47(8): 1175-1179. DOI: 10.11908/j.issn.0253-374x.2019.08.014

### 文章历史

Entropy Scheme for Solute Transport Equation in Groundwater
CHEN Rongsan , WANG Fanrui , ZOU Min
School of Mathematics and Physics, China University of Geosciences, Wuhan 430074, China
Abstract: The equation of solute transport in groundwater is divided into the convection equation and the diffusion equation by using the splitting method. The entropy scheme is applied to solve the conservation equation while the central difference scheme is applied to solve the diffusion equation. Numerical experiments show that the scheme does not produce excessive problems, and no physical oscillation occurs.
Key words: groundwater solute transport equation    entropy scheme    convection dominated

1 熵格式的描述

 $\frac{{\partial \rho }}{{\partial t}} + u\frac{{\partial \rho }}{{\partial x}} = {D_{\rm{L}}}\frac{{{\partial ^2}\rho }}{{\partial {x^2}}}$ (1)

 $\frac{{\partial U\left( \rho \right)}}{{\partial t}} + u\frac{{\partial U\left( \rho \right)}}{{\partial x}} = 2{D_{\rm{L}}}\rho \frac{{{\partial ^2}\rho }}{{\partial {x^2}}}$ (2)

 ${\rho _{i,n}} \approx \frac{1}{{\Delta x}}\int_{{x_{i - 1/2}}}^{{x_{i + 1/2}}} {\rho \left( {x,{t_n}} \right){\rm{d}}x}$ (3)

 ${U_{i,n}} \approx \frac{1}{{\Delta x}}\int_{{x_{i - 1/2}}}^{{x_{i + 1/2}}} {U\left( {\rho \left( {x,{t_n}} \right)} \right){\rm{d}}x}$ (4)

 $\left\{ \begin{array}{l} \frac{{\partial \rho }}{{\partial t}} + u\frac{{\partial \rho }}{{\partial x}} = 0\\ \frac{{\partial U\left( \rho \right)}}{{\partial t}} + u\frac{{\partial U\left( \rho \right)}}{{\partial x}} = 0 \end{array} \right.$ (5)

 $\left\{ \begin{array}{l} \frac{{\partial \rho }}{{\partial t}} = {D_{\rm{L}}}\frac{{{\partial ^2}\rho }}{{\partial {x^2}}}\\ \frac{{\partial U\left( \rho \right)}}{{\partial t}} = {D_{\rm{L}}}\rho \frac{{{\partial ^2}\rho }}{{\partial {x^2}}} \end{array} \right.$ (6)

(1) 重构. tn时刻第i个网格Ii内采用熵台阶重构，它是一个阶梯函数，如式(7)所示.

 $R\left( {x;{\rho _n},{U_n}} \right) = \left\{ \begin{array}{l} {\rho _{i,n}} - {d_{i,n}},\;\;\;\;{x_{i - 1/2}} < x \le {x_i}\\ {\rho _{i,n}} + {d_{i,n}},\;\;\;\;{x_i} < x \le {x_{i + 1/2}} \end{array} \right.$ (7)

 $\frac{1}{{\Delta x}}\int_{{x_{i - 1/2}}}^{{x_{i + 1/2}}} {R\left( {x;{\rho _n},{U_n}} \right){\rm{d}}x} = {\rho _{i,n}}$ (8)

 $\frac{1}{{\Delta x}}\int_{{x_{i - 1/2}}}^{{x_{i + 1/2}}} {U\left( {R\left( {x;{\rho _n},{U_n}} \right)} \right){\rm{d}}x} = {U_{i,n}}$ (9)

 ${d_{i,n,e}} = {\mathop{\rm sgn}} \left( {{C_{i + 1}} - {C_{i - 1}}} \right)\sqrt {{U_{i,n}} - {{\left( {{\rho _{i,n}}} \right)}^2}}$ (10)

di, n, e符号取为和(ρj+1, n-ρj-1, n)相同.另外，定义一个总变差下降(TVD)的台阶di, n, TVD

 ${d_{i,n,{\rm{TVD}}}} = \min \bmod \left( {{\rho _{i,n}} - {\rho _{i - 1,n}},{\rho _{i + 1,n}} - {\rho _{i,n}}} \right)$ (11)

 $\begin{array}{l} \min \bmod \left( {a,b} \right) = \\ \left\{ \begin{array}{l} s \cdot \min \left( {\left| a \right|,\left| b \right|} \right),如果\;{\mathop{\rm sgn}} \left( a \right) = {\mathop{\rm sgn}} \left( b \right) = s\\ 0,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;其他 \end{array} \right. \end{array}$ (12)

 ${d_{i,n}} = {\mathop{\rm sgn}} \left( {{\rho _{i + 1}} - {\rho _{i - 1}}} \right)\min \left( {\left| {{d_{i,n,e}}} \right|,\left| {{d_{i,n,{\rm{TVD}}}}} \right|} \right)$ (13)

(2) 发展.以重构函数R(x; ρn, Un)作为tn层的初值，求解以下初值问题:

 $\left\{ \begin{array}{l} {v_t} + u{v_x} = 0\\ v\left( {x,{t_n}} \right) = R\left( {x;{\rho _n},{U_n}} \right) \end{array} \right.$ (14)

(3) 网格平均.在t=tn+1时的数值解和数值熵分别计算为

 ${\rho _{i,n + 1}} = \frac{1}{{\Delta x}}\int_{{x_{i - 1/2}}}^{{x_{i + 1/2}}} {v\left( {x,{t_{n + 1}}} \right){\rm{d}}x}$ (15)
 ${U_{i,n + 1}} = \frac{1}{{\Delta x}}\int_{{x_{i - 1/2}}}^{{x_{i + 1/2}}} {U\left( {v\left( {x,{t_{n + 1}}} \right)} \right){\rm{d}}x}$ (16)

 $\begin{array}{l} {\rho _{i,n + 1}} = {\rho _{i,n}} + \frac{{2{D_{\rm{L}}}\Delta t}}{{\Delta {x^2}}} \cdot \\ \;\;\;\;\;\;\;\;\;\;\left( {{\rho _{i + 1,n}} - 2{\rho _{i,n}} + {\rho _{i - 1,n}}} \right) \end{array}$ (17)
 $\begin{array}{*{20}{c}} {{U_{i,n + 1}} = {U_{i,n}} + \frac{{2{D_{\rm{L}}}{\rho _{i,n}}\Delta t}}{{\Delta {x^2}}} \cdot }\\ {\left( {{\rho _{i + 1,n}} - 2{\rho _{i,n}} + {\rho _{i - 1,n}}} \right)} \end{array}$ (18)
2 数值试验和结果分析

 $\rho = \frac{{{\rho _0}}}{2}{\rm{erfc}}\left( {\frac{{x - ut}}{{2\sqrt {{D_{\rm{L}}}t} }}} \right) + \frac{{{\rho _0}}}{2}{{\rm{e}}^{\frac{{ux}}{{{D_{\rm{L}}}}}}}{\rm{erfc}}\left( {\frac{{x + ut}}{{2\sqrt {{D_{\rm{L}}}t} }}} \right)$ (19)

Peclet数定义为Pex/αPe是用来描述对流和弥散的相对大小，根据Pe的大小可以将对流-弥散方程分为两大类：当Pe>2时，称之为对流占优方程；当Pe≤2时，称之为弥散占优方程.计算区域为m，网格数取为30. 图 1~4是对流强度Pe分别为1、4、16、32，时间分别为1 d和4 d的数值解和解析解的质量浓度分布对比图，图中圆圈表示数值解，实线表示解析解.从图 1~4中不难看出，熵格式没有出现过量问题，也没有出现非物理振荡，所产生的数值弥散比较小.随着对流强度的增加，数值解和解析解吻合得更好，这说明熵格式在处理强对流问题方面具有很大优势.

 图 1 Pe=1时质量浓度分布 Fig.1 Mass concentration distribution at Pe=1
 图 2 Pe=4时质量浓度分布 Fig.2 Mass concentration distribution at Pe=4
 图 3 Pe=16时质量浓度分布 Fig.3 Mass concentration distribution at Pe=16
 图 4 Pe=32时质量浓度分布 Fig.4 Mass concentration distribution at Pe=32

 $\left\{ \begin{array}{l} \frac{{\partial \rho }}{{\partial t}} + \left( {x - 1} \right)\frac{{\partial \rho }}{{\partial x}} = \frac{1}{2}\left( {1 - x} \right)\frac{{{\partial ^2}\rho }}{{\partial {x^2}}} + 3\rho ,\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;x \in \left( {0,1} \right),t > 0\\ \rho \left( {x,0} \right) = x\left( {1 - x} \right),\;\;\;\;\;\;0 \le x \le 1\\ \rho \left( {0,t} \right) = \rho \left( {1,t} \right) = 0,\;\;\;\;\;t > 0 \end{array} \right.$ (20)

 图 5 t=0.6时质量浓度分布 Fig.5 Mass concentration distribution at t=0.6
3 结论

 [1] ROACHE P J. Computational fluid dynamics[M]. Albuquerque: Hermosa Publishing, 1976 [2] HEINRICH J C, HUYAKORN P S, ZIENKIEWICKY C, et al. An upwind finite element scheme for two dimensional convective transport equation[J]. International Journal for Numerical Methods in Engineering, 1977, 11: 131 DOI:10.1002/nme.1620110113 [3] SUN N Z, YEH W G. A proposed upstream weight numerical method for simulation pollutant transport in groundwater[J]. Water Resources Research, 1983, 19(6): 1489 DOI:10.1029/WR019i006p01489 [4] BROOKS A N, HUGHES T J R. Streamline upwind Petrov-Galerkin formulations for convection dominated flow with particular emphasis on the impressible Navier-Stokes equations[J]. Computer Methods Applied Mechanics and Engineering, 1982, 32: 199 DOI:10.1016/0045-7825(82)90071-8 [5] 梅一, 吴吉春. 地下水溶质运移数值模拟中减少误差的新方法[J]. 水科学进展, 2009, 20(5): 640 MEI Yi, WU Jichun. New method for reducing the numerical error in solving the problem of contaminant transport in groundwater[J]. Advances in Water Sciences, 2009, 20(5): 640 [6] LEISMAN H M, FRIND E O. A symmetric-matrix time integration scheme for the efficient solution of advection-dispersion problems[J]. Water Resources Research, 1989, 25(6): 1133 DOI:10.1029/WR025i006p01133 [7] 刘扬. 对流扩散方程的新型Crank-Nicholson差分格式[J]. 数学杂志, 2005, 25(4): 464 LIU Yang. A new Crank-Nicholson difference scheme for convection-diffusion equations[J]. Journal of Mathematics, 2005, 25(4): 464 [8] 成建梅, 胡进武. 饱和水流溶质运移问题数值解法综述[J]. 水文地质工程地质, 2003, 30(2): 99 CHENG Jianmei, HU Jinwu. Reviews on numerical method for solving solute transport problems in the saturated porous media[J]. Hydrogeology & Engineering Geology, 2003, 30(2): 99 DOI:10.3969/j.issn.1000-3665.2003.02.026 [9] 焦甜, 王军霞, 唐仲华, 等. 地下水溶质运移方程有限差分格式的实证研究[J]. 安全与环境工程, 2016, 23(3): 17 JIAO Tian, WANG Junxia, TANG Zhonghua, et al. Empirical study of the difference schemes of solute transport equation in groundwater[J]. Safety and Environmental Engineering, 2016, 23(3): 17 [10] CHENG Y, SHU C W. Superconvergence of local discontinuous galerkin methods for one-dimensional convection-diffusion equations[J]. Computers & Structures, 2009, 87: 630 [11] LI H X, WANG Z G, MAO D-K. Numerically neither dissipative nor compressive scheme for linear advection equation and its application to the Euler system[J]. Journal of Scientific Computing, 2008, 36: 285 DOI:10.1007/s10915-008-9192-x [12] CUI Y F, MAO D-K. Numerical method satisfying the first two conservation laws for the Kortewegde Vries equation[J]. Journal of Computational Physcics, 2007, 227: 376 DOI:10.1016/j.jcp.2007.07.031 [13] CUI Y F, MAO D-K. Error self-canceling of a difference scheme maintaining two conservation laws for linear advection equation[J]. Mathematics of Computation, 2012, 81: 715 [14] CHEN R S, MAO D-K. Entropy-TVD scheme for nonlinear scalar conservation laws[J]. Journal of Scientific Computing, 2011, 47: 150 DOI:10.1007/s10915-010-9431-9 [15] CHEN R S, MAO D-K. Improved Entropy-Ultra-bee scheme for the Euler system of gas dynamics[J]. Journal of Computational Mathematics, 2017, 35: 121 DOI:10.4208/jcm.1603-m2015-0338 [16] CHEN R S, ZOU M, XIAO L. Entropy-TVD scheme for the shallow water equations in one dimension[J]. Journal of Scientific Computing, 2017, 71: 822 DOI:10.1007/s10915-016-0322-6