Jump dominance on the contaminant transport residual error estimator

Abstracts

Considering the fact that the transport of contaminants occur in small advection regime, a residual estimator is used to evaluate the parabolic equation that describes the phenomena of advection-diffusion-reaction in the saturated porous medium. The correspondent numerical solution is obtained by the finite element method using the θ A-stable scheme and a Python code. The residual error estimator considers its component parts and enables analysis and comparisons of contributions to residual error.This analysis considers a problem sequence with a different number of elements in computational mesh. As a result of numerical simulations, there is a dominance of the jump residuals compared to other residual estimates and this dominance increases with both, the growth of elements number in the computational mesh and with time. Furthermore, the considered problem requires additional effort for the calculation of contributions associated with the L² projection of the contaminant source function on the finiteelement space.

residual error estimator; advection-diffusion-reaction equation; small advection regime; finite elements; contaminant transport in porous media

Considerando o fato de que o transporte de contaminantes ocorre em regime de pequena advecção, um estimador residual é usado para analisar a equação parabólica que descreve o fenômeno de advecção-difusão-reação em meio poroso saturado. A solução numérica correspondente é obtida pelo método de elementos finitos usando um esquema θ A-estável em código Python. O estimador de erro residual avalia separadamente as partes componentes do erro e permite a análise e comparação das contribuições para o erro residual. Essa análise considera uma sequência de problemas com diferentes números de elementos na malha computacional. As simulações numéricas indicam que o residual do salto é dominante em comparação com outros resíduos estimados e a dominância cresce com o número de elementos na malha e com o tempo. Adicionalmente, o problema considerado requer esforço adicional para o cálculo das contribuições associadas com a projeção L² da função fonte de contaminante no espaço de elementos finitos.

estimador de erro residual; equação de advecção-difusão-reação; regime de pequena advecção; elementos finitos; transporte de contaminantes em meios porosos

Jump dominance on the contaminant transport residual error estimator^*

J.P. MartinsI,^**; A. Firmiano^II; E. Wendland^I

^IDepartment of Hydraulics and Sanitary Engineering, EESC - USP, Av. Trabalhador Sancarlense, 400, 13566-590 São Carlos, SP, Brazil. E-mails: [e-mail]; [e-mail]

^IIAir Force Academy. E-mail: [e-mail]

ABSTRACT

Considering the fact that the transport of contaminants occur in small advection regime, a residual estimator is used to evaluate the parabolic equation that describes the phenomena of advection-diffusion-reaction in the saturated porous medium. The correspondent numerical solution is obtained by the finite element method using the θ A-stable scheme and a Python code. The residual error estimator considers its component parts and enables analysis and comparisons of contributions to residual error.This analysis considers a problem sequence with a different number of elements in computational mesh. As a result of numerical simulations, there is a dominance of the jump residuals compared to other residual estimates and this dominance increases with both, the growth of elements number in the computational mesh and with time. Furthermore, the considered problem requires additional effort for the calculation of contributions associated with the L² projection of the contaminant source function on the finiteelement space.

Keywords: residual error estimator, advection-diffusion-reaction equation, small advection regime, finite elements, contaminant transport in porous media.

RESUMO

Considerando o fato de que o transporte de contaminantes ocorre em regime de pequena advecção, um estimador residual é usado para analisar a equação parabólica que descreve o fenômeno de advecção-difusão-reação em meio poroso saturado. A solução numérica correspondente é obtida pelo método de elementos finitos usando um esquema θ A-estável em código Python. O estimador de erro residual avalia separadamente as partes componentes do erro e permite a análise e comparação das contribuições para o erro residual. Essa análise considera uma sequência de problemas com diferentes números de elementos na malha computacional. As simulações numéricas indicam que o residual do salto é dominante em comparação com outros resíduos estimados e a dominância cresce com o número de elementos na malha e com o tempo. Adicionalmente, o problema considerado requer esforço adicional para o cálculo das contribuições associadas com a projeção L² da função fonte de contaminante no espaço de elementos finitos.

Palavras-chave: estimador de erro residual, equação de advecção-difusão-reação, regime de pequena advecção, elementos finitos, transporte de contaminantes em meios porosos.

1 INTRODUCTION

Analytical solutions are far from encompassing the variety of phenomena present in contaminant transport and numerical methods are required to obtain an approximate solution of the corresponding mathematical model. Moreover, the reliability of computational methods depends on the discretization technique and the quality of the finite element mesh adopted. Although estimates of a priori are available, a posteriori estimates are fundamental for practical problems involving finite element method [13]. Once the numerical result is obtained, the a posteriori error estimator can be used to provide general or specific information about the quality of the numerical solution [3].

2 THE MODEL PROBLEM FOR THE CONTAMINANT TRANSPORT

The linear second order parabolic equation is used to describe the contaminant transport in a domain Ω ⊂

², with general unknown concentration solution C = C(x, y, t), data functions D, v, λ, f, g, C₀ and final time t_final is arbitrary, but is kept constant.

The data are real valued functions that may depend on space and time while the initial condition depends only on space [17]. Here, the equation describes the phenomena of advection-dispersion-reaction (ADR model) in a saturated porous media and, according to [1], can be presented by:

where Ω ⊂

² is a polygonal cross-section with a Lipschitz boundary Γ consisting of two disjoint parts Γ_D and Γ_N. The space dependent function C = C₀ in Ω for t = 0 is the initial condition while

C = C_D in Γ_D × (0, t_final] and n · D∇C = g in Γ_N × (0, t_final]

are the Dirichlet and Neumann boundary conditions. Assuming that the data satisfy the additional conditions [20]:

P1 - The dispersion D = D(x, y, t) is a continuously differentiable matrix-valued function and symmetric, uniformly positive definite and uniformly isotropic. Formally,

and

is moderate size constant.

P2 - The velocity v = v(x, y, t) = (υ_x(x, y, t), υ_y(x, y, t)) is a continuously differentiable vector-field and scaled such that

P3 - The reaction λ = λ (x, y, t) is a continuous non-negative scalar function.

P4 - There is a constant β such that λ - (1/2)div (v) > β for almost all (x, y) ∈ Ω and 0 < t < t_final. Moreover there is and a constant c_b > 0 of moderate size such that

P5 - The Dirichlet boundary Γ_D has positive measure (d - 1)-dimensional and includes the inflow boundary

With these additional assumptions, the contaminant transport regime can be classified into:

• dominant dispersion:

with constants of moderate size;

• regime of dominant reaction:

with a constant c_c of moderate size;

• regime of dominant advection: c_υ » ε.

A detailed discussion of these elements can be found in D. Praetorius [17] which is an extension of the results from Verfürth [20].

To derive the space-time discretization of (2.1), we consider a test function w ∈ where denotes the subspace of the Sobolev space H¹(Ω) = W^1,2(Ω), with functions that vanish on the Dirichlet boundary Γ_D. Multiply equation by a test function and use integration by parts to derive the weak form

Next, consider the partition I = {[t_{n - 1}, t_n]:1 < n < N_I} of the time interval [0, t_final] such that 0 = t₀ < t₁ ... < = t_final.

For every n with 1 < n < N_I denote by I_n = [t_{n - 1}, t_n] the n-th subinterval and τ_n = t_n - t_{n - 1} its length.

With every intermediate time t_n, 0 < n < N_I associate an admissible, affine equivalent, shape regular partition T_n of Ω and a corresponding finite element space X_n. In addition, the partitions I and T_n and the spaces X_n must satisfy the following assumptions:

• Ω ∪ Γ is the union of all elements in T_n;

• Affine equivalence, Admissibility and Shape-regularity;

• Non-degeneracy, transition condition and degree condition.

With a time discretization parameter θ ∈ [, 1] and the abbreviations Cⁿ = C(x, y, t_n), Dⁿ = D(x, y, t_n), vⁿ = v(x, y, t_n), λⁿ = λ(x, y, t_n), fⁿ = f(x, y, t_n), gⁿ = g(x, y, t_n) the finite element approximation with θ-A-stable-scheme is obtained by replacing the approximations in the weak form of the parabolic problem and is given by

Find Cⁿ ∈ Xⁿ, 0 < n < N_I, such that C⁰ = π₀C₀ and for n = 1, 2 ..., N_I

where Cⁿ = and w_n = .

In particular, θ = 1/2 gives the Crank-Nicolson scheme and θ = 1 gives implicit Euler scheme [20]. Thus, finite element formulation (2.7) can be rewritten as a(, w_n) = L(w_n). The term a(, w_n) is called the bilinear form and is defined by expression (2.8)

while the term L(w_n) is called linear form and defined by expression (2.9)

Although variational formulation (2.7) is the key in finite element method, forms (2.8) and (2.9) are essential for implementation using the FEniCS Project methodology [4]. These forms generate the linear system to be solved in each step of the simulation process. The solution in (n - 1)-th time step and the solution on the n-th time step are used to provide an error measure.

A detailed description of the residual estimator and assumptions presented here are in references [20, 19, 16, 17]. For finite elements method, a detailed discussion can be found in reference [2].

In the residual method, an element residuals R_K is defined by:

while an edge or face residual R_E is defined by:

where J is the jump operator, f_I(x, y, t) = π_n(θ f(x, y, t_n) + (1 - θ) f(x, y, t_{n - 1})) and g_I(x, y, t) = π_n(θ g(x, y, t_n) + (1 - θ) g(x, y, t_{n - 1})) the projection functions on the finite element space X_n [20]. According to Verfürth [20], if the transport regime is of small advection then the residual estimator is given by equation (2.12)

where

and

with weighting factors α_S = min{, -}, where S = {K, E} is an element or an edge/ face and = ∞ if β = 0.

3 IMPLEMENTATION

The Python numerical code considers the available methodology for the FEniCS Project. A complete description of this project can be found in references [8, 9, 10, 11, 12, 14] or at http://fenicsproject.org [4]. For graphical display of numerical solutions, Matplotlib/Scitools was used [15, 18]. Using Python language, the contaminant transport equation is implemented using the available tools from [4]. These tools provide conditions for setting the transport of contaminants directly through the use of bilinear and linear forms. Furthermore, mesh, initial conditions, dispersions, velocity field, projections and boundary conditions are defined using the available classes/tools described in documentation [4]. For example, mesh = UnitSquareMesh(n_x, n_y, 'crossed') gives a mesh with n_x, n_y triangular elements in each coordinate direction and crossed orientation. Reference [4] is a complete description of the available classes and tools.

After the solution in the n-th step is obtained, the residual error estimates (2.13) are implemented through the summation of the quantities defined in (2.10) and (2.11). Formally,

According to the designations adopted by Verfürth [20] and D. Praetorius [17], components ECⁿ, JCⁿ, BCⁿ, |||Cⁿ - C^{n - 1}|||² will be called, respectively: element, jump, boundary and time contributions for the n-th time step. These quantities allow us to rewrite equation (2.12) as a sum of spatial and temporal contributions, obtained in each step, weighted by time step. Formally,

From equation (3.2), the amount (ESⁿ)² = (ηⁿ)² + |||Cⁿ - C^{n - 1}|||² may be regarded as the residual error at each time step.

The code that implements the element, jump, boundary and time contributions follows an example available in [6]. This allows the calculation of η_I and enables comparisons of components or individual evaluation of contributions. This paper considers:

• Jump contribution versus element contribution , defined by

which provides information on the magnitude of the jumps on the elements;

• Time contributions weighted by time time step, defined by τ_n|||Cⁿ - C^{n - 1}|||² which provides information on time residual error.

4 RESULTS AND DISCUSSION

With an adaptation of a problem described in [7] and [17], this paper implements the contaminant transport equation in a two-dimensional rectangular domain Ω defined by points x = (x, y) ∈ ². A spatial adaptation of that problem consider Ω defined by points (x₀, y₀) = (0, 0), (x₁, y₀) = (80, 0), (x₁, y₁) = (80, 40) and (x₀, y₁) = (0, 40), inital condition C₀ = C(x, y, 0) = 0 and contaminant source f defined by

The physical parameters are the velocity field v = (υ_x, υ_y) = (0.864,0), the dispersion D = α_Dυ_x I_{2 × 2} with I_{2 × 2} the identity matrix of order two and α_D = 0.05m the diffusivity. Neumann boundary condition is defined, using bilinear and linear forms, by g = n · D∇C on Γ_N = {x₁} × (y₀, y₁) while Dirichlet boundaries are defined, using the available class, by C = 0 on Γ_D = Γ \ Γ_N.

For all t_n, 1< n < N_I, the numerical result of contaminant transport provides the distribution function Cⁿ ∈ ², which describes the contaminant concentration at each point of Ω. Figure (1) presents the numerical solution for t₇₀ = 70days, τ_n = 1.0day, lagrangian functions of order two, θ = 1.0 and finite element mesh with n_x = 200 = 2n_y triangular elements in each direction and left/rigth orientation.

For residual estimates, the conditions P1-P5 were verified and contaminant transport was classified in advection dominated [17]. However, the regime is the small advection since C_c = υ/ε is a constant of moderate size [5].

The numerical simulations consider finite element meshes with n_x = 2n_y triangular elements in each coordinate direction with left/right orientation, τ_n = 0.50, t_final = 100.0days, θ = 1.0 and lagrangian functions of order two. Figures (2)-a and (2)-b presents the results for temporal estimates and for ratio for a sequence of problems with n_x = [100, 200, 300, 400].

Numerical results provide that temporal residual decreases with time and as the number of elements in mesh increases. On the other hand, increases with time and as the number of elements in mesh increases. The increasing in with time is due to the contaminant front advances, which, in turn, reflects in the dominance of jump contributions under element contributions already in the coarse mesh. This dominance characterizes the jump contributions as the most relevant part of η_I because τ_n|||Cⁿ - C^{n - 1}|||²≈ 0 and g = n · D∇C] is, by definition, such that BCⁿ = 0.

With the results, η_I can be approximated by

where = 0 and (ESⁿ)² = (ηⁿ)² + |||Cⁿ - C^{n - 1}|||²≈ ECⁿ + JCⁿ ≈ JCⁿ.

This results and the approximation for η_I reveals that jump error is not a negligible quantity and complements the results obtained by Firmiano [5]. In fact the jump residual can be the most important part of the residual estimates.

5 CONCLUSION

The presence of residual estimates enables an analysis of the numerical solution, which provides information about the quality of the numerical results. Error estimator partition allows us to analyze and compare the component parts of residual error. This provides a better understanding of residual behavior for changing number of elements in mesh or time step. As a result of simulations, jump dominance manifests for all adopted meshes, which is due to the advances of contaminant front in the computational domain. As the number of finite elements increases, the magnitude of dominance becomes more significant, making elements and temporal contributions negligible when compared to jump residual. This dominance shows that jump residual is a important quantity in residual error estimator and can not be neglected.

Received on November 10, 2013

Accepted on April 9, 2014

• [1] J. Bear. Hydraulics of Groudwater, Dover Publications Inc., New York, (1979).
• [2] S.C. Brenner. The Mathematical Theory of Finite Element Methods, em "Texts in Applied Mathematics, v.15" (J.E. Marsden, L. Sirovich and S.S. Antman, eds.), Springer-Verlag, New York, (1994).
• [3] R.E. Ewing. A posteriori error estimation. Computer Methods in Applied Mechanics and Engineering, 82(1-3): 59-72, September (1990).
• [4] Fenics project documentation. http://fenicsproject.org/, 2012. Accessed in 05/2012.
• [5] A. Firmiano. "Um estimador de erro a posteriori para a equação do transporte de contaminantes em regime de pequena advecção". Tese de Doutorado, SHS/EESC/USP, (2010).
• [6] https://answers.launchpad.net/dolfin/+question/177108
  » link
• [7] G. Hofinger & F. Judex. Pollution in groundwater flow: definition of ARGESIM comparison C19. Simulation News Europe, 44/45: 51-52, (2005).
• [8] R.C. Kirby. Algorithm 839: Fiat, a new paradigm for computing finite element basis functions. ACM Transactions on Mathematical Software, 30(4) (2004), 502-516.
• [9] R.C. Kirby. Optimizing the evaluation of finite element matrices. SIAM Journal on Scientific Computing, 27(3) (2005), 741-758.
• [10] R.C. Kirby. A compiler for variational forms. ACM Transactions on Mathematical Software, 32(3), (2006).
• [11] R.C. Kirby. Efficient compilation of a class of variational forms. ACM Transactions on Mathematical Software, 33(3), (2007).
• [12] R.C. Kirby. Geometric optimization of the evaluation of finite element matrices. SIAM Journal on Scientific Computing, 29(2) (2007), 827-841.
• [13] D. Kuzmin. Goal-oriented a posteriori error estimates for transport problems. Mathematics andComputers in Simulation, 80(8): 1674-1683, April (2010).
• [14] A. Logg. Automating the finite element method. Archives of Computational Methods in Engineering, 14(2) (2007), 93-138.
• [15] Matplotlib. http://matplotlib.sourceforge.net/index.html, 2012. Accessed in 06/2012.
  » link
• [16] A. Papastavrou. A posteriori error estimators for stationary convection-diffusion problems: a computational comparison. Computer Methods in Applied Mechanics and Engineering, 189(2): 449-462, September (2000).
• [17] D. Praetorius. A space-time adaptive algorithm for linear parabolic problems, Asc report 07/2008, Institute for Analysis and Scientific Computing Vienna University of Technology-TU Wien, 2008, disponível em www.asc.tuwien.ac.at ISBN 978-3-902627-00-1.
• [18] Scitools: Python library for scientific computing. http://code.google.com/p/scitools/, 2012. Accessed in 06/2012.
• [19] R. Verfürth. A review of a posteriori error estimation techniques for elasticity problems. Computer Methods in Applied Mechanics and Engineering, 176(1-4): 419-440, July (1999).
• [20] R. Verfürth. Adaptive finite element methods lecture notes winter term 2007/08, 2007/2008, disponível em http://citeseerx.ist.psu.edu/viewdoc/summary-doi=10.1.1.155.4064

*

Paper presented at CMAC-SE 2013.

**

Corresponding author: João Paulo Martins dos Santos

Publication Dates

• Publication in this collection
  10 June 2014
• Date of issue
  Apr 2014

History

• Received
  10 Nov 2013
• Accepted
  09 Apr 2014