Research Article  Open Access
A Fully Discrete Symmetric Finite Volume Element Approximation of Nonlocal Reactive Flows in Porous Media
Abstract
We study symmetric finite volume element approximations for twodimensional parabolic integrodifferential equations, arising in modeling of nonlocal reactive flows in porous media. It is proved that symmetric finite volume element approximations are convergent with optimal order in norm. Numerical example is presented to illustrate the accuracy of our method.
1. Introduction
In this paper, we consider symmetric finite volume element discretizations of the following initial value problem for the operator equation for : where is a strongly elliptic differential operator and is a secondorder elliptic differential operator in space. The operators and incorporate Dirichlet and Neumann boundary conditions. The problem (1) is an abstract form of an initial boundary value problem for a parabolic integrodifferential equation.
This model is very important in the transport of reactive and passive contaminats in aquifers, an area of active interdisciplinary research of mathematicians, engineers, and life scientists. From a mathematical point of view, the evolution of either a passive or reactive chemical within a velocity field exhibiting strong variation on many scales defies representation using classical Fickian theory. The evolution of a chemical in such a velocity field when modeled by Fickiantype theories leads to a dispersion tensor whose magnitude depends upon the time scales of observation. In order to avoid such difficulty, a new class of nonlocal models of transport have been derived. In this case, the constitutive relations involve either integrals or higher order derivatives, which take multiscales into consideration. We refer to [1, 2] for the derivation of the mathematical models and for the precise hypotheses and analysis.
Mathematical formulations of this kind also arise naturally in various engineering models, such as nonlocal reactive transport in underground water flows in porous media [3], heat conduction, radioactive nuclear decay in fluid flows [4], nonNewtonian fluid flows, or viscoelastic deformations of materials with memory (in particular polymers) [5], semiconductor modeling [6], and biotechnology. One very important characteristic of all these models is that they all express a conservation of a certain quantity (mass, momentum, heat, etc.) in any moment for any subdomain. This in many applications is the most desirable feature of the approximation method when it comes to the numerical solution of the corresponding initial boundary value problem.
This type of equations has been extensively treated by finite element, finite difference, and collocation methods in the last years [7–12].
The finite element method conserves the flux approximately; therefore, in the asymptotic limit (i.e., when the grid stepsize tends to zero) it produces adequate results. However, this could be a disadvantage when relatively coarse grids are used. Perhaps the most important property of the finite volume method is that it exactly conserves the approximate flux (heat, mass, etc.) over each computational cell. This important property combined with adequate accuracy and ease of implementation has contributed to the recent renewed interest in the method. The theoretical framework and the basic tools for the analysis of the finite volume element methods have been developed in the last decade (see, e.g., [13–19]). Reference [20] has given the finite volume element approximations of the problem (1). But in a general case the coefficient matrix of the linear system obtained from the finite volume element method is not symmetric. In this paper, we study symmetric modified finite volume element approximations for twodimensional parabolic integrodifferential equations, arising in modeling of nonlocal reactive flows in porous media. It is proved that symmetric modified finite volume element approximations are convergent with optimal order in norm.
Throughout this paper we use (without or with subscript or superscript) to denote a generic constant independent of the discretization parameter.
2. Fully Discrete Symmetric Finite Volume Element Scheme
Consider the following initial boundary value problem: find such that where is a bounded polygonal domain in with boundary . is a realvalued symmetric and uniformly positive definite matrix, is a matrix. and are known functions, which are assumed to be smooth so that problem (2) has a unique solution in a certain Sobolev space.
The problem (2) can be written in the form (1) by introducing the operators and : for all and for any .
We assume that is a convex polygonal domain. The domain is split into triangular finite elements . The elements are considered to be a closed set, and the triangulation is denoted by . Then and denotes all nodes or vertices: In order to accommodate Dirichlet boundary conditions, we will also need the set of vertices that are internal to , denoted by , that is, . For a given vertex , we define by the index set of all neighbors of in .
For a given triangulation , we construct a dual mesh based upon whose elements are called control volumes. In the finite volume methods, there are various ways to introduce the control volumes. Almost all approaches can be described in the following general scheme: in each triangle a point is selected; similarly, on each of the three edges of a point is selected; then is connected with the points by straight lines (see Figure 1).
Thus, around each vertex , we associate a control volume , which consists of the union of the subelements , which have as a vertex. Also, let denote the interface of two control volumes and . We call the partition regular or quasiuniform, if there exists a positive constant such that , for all . Here, is the maximal diameter of all elements . In the first (and most popular) control volume partition, the point is chosen to be the medicenter (the center of gravity or centroid) of the finite element and the points are chosen to be the midpoints of the edges of . This type of control volume can be introduced for any finite element partition and leads to relatively simple calculations. Besides, if the finite element partition is locally regular, that is, there is a constant such that , , for all elements , then the finite volume partition is also locally regular.
In this paper, we will also use the construction of the control volumes in which the point is the circumcenter of the element . Then obviously, are the perpendicular bisectors of the three edges of . This construction requires that all finite elements are triangles of acute type, which we will assume whenever such a triangulation is used.
We define the linear finite element space : and its dual volume element space : Obviously, and , where are the standard nodal linear basis functions associated with the node and are the characteristic functions of the volume . Let be the piecewise linear interpolation operator and let be the piecewise constant interpolation operator. That is,
The semidiscrete finite volume element approximation of (2) is a solution to the problem: find , for such that
Here, the bilinear forms and are defined by where denotes the outernormal direction to the domain under consideration.
It is more convenient to rewrite (8) in the following form, which we use in the exposition:
Next, we define the fully discrete time stepping schemes. Let be a timestep size and , , and is an approximation of , such that .
The backward Euler scheme is defined to be the solution of such that where are the weights, and the quadrature error is given by for any smooth functions and and its error satisfies
Then we present fully discrete symmetric finite volume element scheme as the following: find , such that where is defined by We note that is a modified term.
3. Some Auxiliary Results
Lemma 1 (see [17]). If the matrix is constant over each element , then
Proposition 2. Consider the following:
Proof. By the definition of (16), is piecewise constant on . Apply Lemma 1, and note that are constants on , we can get (18).
It can be deduced from Proposition 2 that the coefficient matrix of (14) is symmetric.
Remark 3. For each , if we select a point and replace (16) by , for all , from the analysis given in the following we can get the same error estimates as in Theorem 12.
We define a space constant, for all , then introduce : by or , for all .
Lemma 4 (see [21]). Define , . Then there exist two positive constants and , such that
Next we define some discrete norms on : where , the distance between and . Obviously, these norms are well defined for as well and .
Lemma 5 (see [20]). There exist two positive constants , independent of , such that
Lemma 6 (see [20]). Assume that the jumps (if any) of coefficient matrices and are aligned with the finite element partition , and over each element their entries are functions. Then,(a) there are positive constants and , independent of and , such that for all , (b) for is fixed there is a constant independent of and , such that for ,
4. Norm Error Estimate
Now we introduce linear functions , which are used in the error analysis of the finite volume element method: The following estimate is a simple consequence of the BrambleHilbert lemma.
Lemma 7 (see [20]). If , then there is a positive constant , independent of , such that for ,
Lemma 8 (see [20]). (a) If , then there exists a positive constant , independent of and , such that
(b) If , then, for fixed there is a constant , independent of and , such that for
For any fixed , one can define the Ritz projection function and the operator , such that
Remark 9. If the partition is regular (quasiuniform) and is regular, then
However, these estimates for the Ritz projection lead to suboptimal error estimates for the finite volume element solution of the integrodifferential equation. In order to obtain optimal order estimates, we need a projection that also takes into account the integral term. This type of projection was called by Cannon and Lin [7] the RitzVolterra projection and has been used in the analysis of the finite element method for integrodifferential equations.
We define the RitzVolterra projection of a function defined on the cylinder . The RitzVolterra projection is defined for by
Lemma 10 (see [20]). Assume that for fixed, there is a constant independent of and , such that for all ,
Now we consider an estimate for the RitzVolterra projection. This estimate is optimal with respect to the order of convergence, but requires regularity of the solution. Therefore, it is suboptimal with respect to the regularity of the solution and can be useful for close to 1. Namely, we have the following result.
Lemma 11 (see [20]). Assume that, for some , . Then for fixed there exists a positive constant independent of and , such that, for all ,
Then we prove the following norm error estimate.
Theorem 12. Let and be the solution of (2) and (14), respectively. Then we have a constant independent of , , and , such that
Proof. Let , from (2), (14), and (30), we have Let , Seting in (34), we get We estimate the terms on the righthand side: Let , and use numerical quadrature error estimates to get Thus, multiplying each term on both sides of (38) by , summing over , and employing Gronwall’s lemma, we obtain And ; hence (33) follows from the above analysis and Lemma 11.
5. Numerical Experiments
In this section, we present a numerical example for solving the problem (2) by using the symmetric modified finite volume element scheme presented in Section 2. Let , , be Delaunay triangulation generated by EasyMesh [22] over with mesh size as shown in Figure 2 and time step be .
We consider the case of , , , , , , the exact solution , and the initial function . Denote the numerical solution of by and set with fixed ratio . The numerical results are presented in Tables 1 and 2. It is observed that the results support our theory.


Acknowledgments
This work is supported by the Excellent Young and Middleaged Scientists Research Fund of Shandong Province (no. 2008BS09026); National Natural Science Foundation of China (no. 11171193); Natural Science Foundation of Shandong Province (no. ZR2011AM016).
References
 J. H. Cushman and T. R. Ginn, “Nonlocal dispersion in media with continuously evolving scales of heterogeneity,” Transport in Porous Media, vol. 13, no. 1, pp. 123–138, 1993. View at: Publisher Site  Google Scholar
 J. H. Cushman, B. X. Hu, and F. Deng, “Nonlocal reactive transport with physical and chemical heterogeneity: localization errors,” Water Resources Research, vol. 31, no. 9, pp. 2219–2237, 1995. View at: Publisher Site  Google Scholar
 G. Dagan, “The significance of heterogeneity of evolving scales to transport in porous formations,” Water Resources Research, vol. 30, no. 12, pp. 3327–3336, 1994. View at: Google Scholar
 V. V. Shelukhin, “A nonlocal in time model for radionuclides propagation in Stokes fluids, dynamics of fluids with free boundaries,” Siberian Branch of the Russian Academy of Sciences, Institute of Hydrodynamics, vol. 107, pp. 180–193, 1993. View at: Google Scholar
 M. Renardy, W. J. Hrusa, and J. A. Nohel, Mathematical Problems in Viscoelasticity, vol. 35, Longman Scientific and Technical, Essex, UK, 1987. View at: MathSciNet
 W. Allegretto, Y. Lin, and A. Zhou, “A box scheme for coupled systems resulting from microsensor thermistor problems,” Dynamics of Continuous, Discrete and Impulsive Systems, vol. 5, no. 1–4, pp. 209–223, 1999. View at: Google Scholar  Zentralblatt MATH  MathSciNet
 J. R. Cannon and Y. Lin, “Nonclassical ${H}^{1}$projection and Galerkin methods for nonlinear parabolic integrodifferential equations,” Calcolo, vol. 25, no. 3, pp. 187–201, 1988. View at: Publisher Site  Google Scholar  MathSciNet
 Y. P. Lin, V. Thomée, and L. B. Wahlbin, “RitzVolterra projections to finiteelement spaces and applications to integrodifferential and related equations,” SIAM Journal on Numerical Analysis, vol. 28, no. 4, pp. 1047–1070, 1995. View at: Google Scholar  Zentralblatt MATH  MathSciNet
 Y. P. Lin, “On maximum norm estimates for RitzVolterra projection with applications to some time dependent problems,” Journal of Computational Mathematics, vol. 15, no. 2, pp. 159–178, 1997. View at: Google Scholar  Zentralblatt MATH  MathSciNet
 I. H. Sloan and V. Thomée, “Time discretization of an integrodifferential equation of parabolic type,” SIAM Journal on Numerical Analysis, vol. 23, no. 5, pp. 1052–1061, 1986. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 T. Zhang, “The stability and approximation properties of RitzVolterra projection and application. I,” Numerical Mathematics. English Series, vol. 6, no. 1, pp. 57–76, 1997. View at: Google Scholar  Zentralblatt MATH  MathSciNet
 V. Thomée and L. B. Wahlbin, “Longtime numerical solution of a parabolic equation with memory,” Mathematics of Computation, vol. 62, no. 206, pp. 477–496, 1994. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 R. E. Bank and D. J. Rose, “Some error estimates for the box method,” SIAM Journal on Numerical Analysis, vol. 24, no. 4, pp. 777–787, 1987. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 Z. Q. Cai, “On the finite volume element method,” Numerische Mathematik, vol. 58, no. 7, pp. 713–735, 1991. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 Z. Q. Cai and S. McCormick, “On the accuracy of the finite volume element method for diffusion equations on composite grids,” SIAM Journal on Numerical Analysis, vol. 27, no. 3, pp. 636–655, 1990. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 W. Hackbusch, “On first and second order box schemes,” Computing, vol. 41, no. 4, pp. 277–296, 1989. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 J. Huang and S. Xi, “On the finite volume element method for general selfadjoint elliptic problems,” SIAM Journal on Numerical Analysis, vol. 35, no. 5, pp. 1762–1774, 1998. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 R. Li, Z. Chen, and W. Wu, Generalized Difference Methods for Differential Equations: Numerical Analysis of Finite Volume Methods, vol. 226, Marcel Dekker, New York, NY, USA, 2000. View at: MathSciNet
 I. D. Mishev, “Finite volume methods on Voronoi meshes,” Numerical Methods for Partial Differential Equations, vol. 14, no. 2, pp. 193–212, 1998. View at: Google Scholar  Zentralblatt MATH  MathSciNet
 R. Ewing, R. Lazarov, and Y. Lin, “Finite volume element approximations of nonlocal reactive flows in porous media,” Numerical Methods for Partial Differential Equations, vol. 16, no. 3, pp. 285–311, 2000. View at: Google Scholar  Zentralblatt MATH  MathSciNet
 S. Liang, X. Ma, and A. Zhou, “A symmetric finite volume scheme for selfadjoint elliptic problems,” Journal of Computational and Applied Mathematics, vol. 147, no. 1, pp. 121–136, 2002. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 EasyMesh web site, http://wwwdinma.univ.trieste.it/nirftc/research/easymesh/.
Copyright
Copyright © 2013 Zhe Yin and Qiang Xu. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.