# A Fully Discrete Symmetric Finite Volume Element Approximation of Nonlocal Reactive Flows in Porous Media

Authors: Zhe Yin; Qiang Xu
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

We study symmetric finite volume element approximations for two-dimensional 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 inL2-norm. Numerical example is presented to illustrate the accuracy of our method.


---

## Body

## 1. Introduction

In this paper, we consider symmetric finite volume element discretizations of the following initial value problem for the operator equation foru=u(t):
(1)ut+𝒜u+∫0t‍ℬ(t,s)u(s)ds=f(t),u(0)=u0,
where 𝒜 is a strongly elliptic differential operator and ℬ is a second-order 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 Fickian-type 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], non-Newtonian 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 step-size 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 two-dimensional 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 L2-norm.Throughout this paper we useC (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: findu=u(x,t) such that
(2)ut-∇·(A∇u)-∫0t∇·(B∇u(s))ds=f,(x,t)∈Ω×(0,T],u(x,t)=0,(x,t)∈∂Ω×(0,T],u(x,0)=u0(x),x∈Ω,
where Ω is a bounded polygonal domain in R2 with boundary ∂Ω. A={aij(x)} is a 2×2 real-valued symmetric and uniformly positive definite matrix, B={bij(x,t,s)} is a 2×2 matrix. f=f(x,t) and u0(x) 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 𝒜:H01(Ω)→H-1(Ω) and ℬ:H01(Ω)→H-1(Ω):
(3)(𝒜u,v)=∫Ω‍A(x)∇u·∇vdx,(ℬu,v)=∫Ω‍B(x,t,s)∇u·∇vdx,
for all v∈H01(Ω) and for any t,s∈(0,T).We assume thatΩ is a convex polygonal domain. The domain Ω is split into triangular finite elements K. The elements K are considered to be a closed set, and the triangulation is denoted by 𝒯h. Then Ω¯=⋃K∈𝒯h‍K and Nh denotes all nodes or vertices:
(4)Nh={p:pisavertexofelement,K∈𝒯handp∈Ω¯}.
In order to accommodate Dirichlet boundary conditions, we will also need the set of vertices that are internal to Ω, denoted by Nh0, that is, Nh0=Nh∩Ω. For a given vertex xi, we define by Π(i) the index set of all neighbors of xi in Nh.For a given triangulation𝒯h, we construct a dual mesh 𝒯h* based upon 𝒯h 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 K∈𝒯h a point q is selected; similarly, on each of the three edges xixj¯ of K a point xij is selected; then q is connected with the points xij by straight lines (see Figure 1).Figure 1
A control volume.Thus, around each vertexxj∈Nh0, we associate a control volume Vj∈𝒯h*, which consists of the union of the subelements K∈𝒯h, which have xj as a vertex. Also, let γij denote the interface of two control volumes Vi and Vj:γij=Vi∩Vj,j∈∏‍(i). We call the partition regular or quasiuniform, if there exists a positive constant C>0 such that C-1h2≤meas(Vi)≤Ch2, for all Vi∈𝒯h*. Here, h is the maximal diameter of all elements K∈𝒯h. In the first (and most popular) control volume partition, the point q is chosen to be the medicenter (the center of gravity or centroid) of the finite element K and the points xij are chosen to be the midpoints of the edges of K. This type of control volume can be introduced for any finite element partition 𝒯h and leads to relatively simple calculations. Besides, if the finite element partition 𝒯h is locally regular, that is, there is a constant C such that ChK2≤meas(K)≤hK2,  hK=diam(K), for all elements K∈𝒯h, then the finite volume partition 𝒯h* is also locally regular.In this paper, we will also use the construction of the control volumes in which the pointq is the circumcenter of the element K. Then obviously, γij are the perpendicular bisectors of the three edges of K. 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 spaceSh:
(5)Sh={v∈C(Ω):v∣KislinearforallK∈Thandv∣∂Ω=0},
and its dual volume element space Sh*:
(6)Sh*={v∈L2(Ω):v∣VisconstantforallV∈Th*andv∣∂Ω=0{L2}}.
Obviously, Sh=span{ϕi(x):xi∈Nh0} and Sh*=span{χi(x):xi∈Nh0}, where ϕi are the standard nodal linear basis functions associated with the node xi and χi are the characteristic functions of the volume Vi. Let Ih:C(Ω)→Sh be the piecewise linear interpolation operator and let Ih*:C(Ω)→Sh* be the piecewise constant interpolation operator. That is,
(7)Ihu=∑xi∈Nh‍u(xi)ϕi(x),Ih*u=∑xi∈Nh‍u(xi)χi(x).The semidiscrete finite volume element approximationuh of (2) is a solution to the problem: find uh(t)∈Sh, for t>0 such that
(8)(uh,t,vh)+A(uh,vh)+∫0t‍B(t,s;uh(s),vh)ds=(f,vh),vh∈Sh*,uh(0)=u0,h∈Sh.Here, the bilinear formsA(u,v) and B(t,s;u,v) are defined by
(9)A(u,v)=-∑xi∈Nh‍vi∫∂Vi‍A(x)∇u·ndx,(u,v)∈Sh×Sh*,B(t,s;u,v)=-∑xi∈Nh‍vi∫∂Vi‍B(x,t,s)∇u·ndx,(u,v)∈Sh×Sh*,
where n denotes the outer-normal direction to the domain under consideration.It is more convenient to rewrite (8) in the following form, which we use in the exposition:
(10)(uh,t,Ih*vh)+A(uh,Ih*vh)+∫0t‍B(t,s;uh(s),Ih*vh)ds=(f,Ih*vh),vh∈Sh.Next, we define the fully discrete time stepping schemes. LetΔt>0 be a time-step size and tn=nΔt,  un=u(tn),  and  ∂tun=(un-un-1)/Δt.uh(0)∈Sh is an approximation of u0, such that ∥uh(0)-u0∥≤Ch2∥u0∥2.The backward Euler scheme is defined to be the solution ofuhn∈Sh such that
(11)(uhn-uhn-1Δt,Ih*vh)+A(uhn,Ih*vh)+∑k=0n-1‍ωn,kB(tn,tk;uhk,Ih*vh)=(fn,Ih*vh),
where ωn,j are the weights, and the quadrature error is given by for any smooth functions g and M and its error
(12)qn(g)=∫0tn‍M(tn,s)g(s)ds-∑k=0n-1‍ωn,jMn,jg(tj)
satisfies
(13)|qn(g)|≤CΔt∫0tn‍(|g|+|g′|)dt.Then we present fully discrete symmetric finite volume element scheme as the following: finduhn∈Sh, such that
(14)(∂tuhn,Ih*vh)+A-(uhn,Ih*vh)+∑k=0n-1‍ωn,kB(tn,tk;uhk,Ih*vh)=(fn,Ih*vh)+A-(uhn-1,Ih*vh)-A(uhn-1,Ih*vh),
where A-(u,v) is defined by
(15)A-(u,v)=-∑xi∈Nhvi∫∂ViA-(x)∇u·ndx,(u,v)∈Sh×Sh*,(16)A-∣K=A-K,A-K=1meas(K)∫K‍A(x)dx.
We note that A-(uhn-1,Ih*vh)-A(uhn-1,Ih*vh) is a modified term.


## 3. Some Auxiliary Results

Lemma 1 (see [17]).
If the matrixA=(aij)2×2 is constant over each element K∈𝒯h, then
(17)-∑xi∈Nh‍vi∫∂Vi‍A∇u·nds=∫Ω‍A∇u·∇vdx,∀u,v∈Sh.Proposition 2.
Consider the following:(18)-∑xi∈Nh‍vi∫∂Vi‍A-∇uhn·nds=∫Ω‍A∇uhn·∇vdx,∀uh,v∈Sh.Proof.
By the definition of (16), A- is piecewise constant on 𝒯h. Apply Lemma 1, and note that ∇uh,∇v are constants on K, we can get (18).It can be deduced from Proposition2 that the coefficient matrix of (14) is symmetric.Remark 3.
For eachK∈𝒯h, if we select a point xK∈K and replace (16) by A-K=A(xK),  for all  K∈𝒯h,  A-∣K=A-K from the analysis given in the following we can get the same error estimates as in Theorem 12.
We define a spaceVh={v∈L∞:v∣K= constant, for all K∈𝒯h}, then introduce Ph0: C(Ω-)→Vh by Ph0v=v(qK) or Ph0v=(1/|K|)∫K‍vdx,  for all K∈𝒯h.Lemma 4 (see [21]).
DefineAh=(Ph0aij)2×2, ah(uh,v)=∫Ω‍Ah∇uh·∇vdx. Then there exist two positive constants C1 and β, such that
(19)|ah(w,v)|≤C1∥w∥1∥v∥1,∀w,v∈Sh,ah(v,v)≥β∥v∥12,∀v∈Sh.Next we define some discrete norms onSh:
(20)∥uh∥0,h2=(uh,uh)0,h=(Ih*uh,Ih*uh),|uh|1,h2=∑xi∈Nh‍∑xj∈Π(i)‍meas(Vi)(ui-ujdij)2,∥uh∥1,h2=∥uh∥0,h2+|uh|1,h2,|||uh|||2=(uh,Ih*uh),
where dij=d(xi,xj), the distance between xi and xj. Obviously, these norms are well defined for uh∈Sh* as well and ∥uh∥0,h=|||uh|||.Lemma 5 (see [20]).
There exist two positive constantsC2,C3, independent of h, such that
(21)C2∥v∥0,h≤∥vh∥≤C3∥v∥0,h,∀vh∈Sh,C2|||vh|||≤∥vh∥≤C3|||vh|||,∀vh∈Sh,C2∥vh∥1,h≤∥vh∥1≤C3∥v∥1,h,∀vh∈Sh.Lemma 6 (see [20]).
Assume that the jumps (if any) of coefficient matricesA(x) and B(x,t,s) are aligned with the finite element partition 𝒯h, and over each element K∈𝒯h their entries are W1,∞(K)-functions. Then,(a)
there are positive constantsh0 and c0,c1, independent of h and u, such that for all 0<h≤h0,
(22)|A(uh,Ih*vh)|≤c1|uh|1,h|vh|1,h,∀uh,vh∈Sh,A(vh,Ih*vh)≥c0∥vh∥1,h2,∀vh∈Sh;(b)
forT>0 is fixed there is a constant C=C(T)>0 independent of  h and u, such that for 0<t≤T,
(23)|B(t,s;uh,Ih*vh)|≤C|uh|1,h|vh|1,h,∀uh,vh∈Sh.


## 4.L2-Norm Error Estimate

Now we introduce linear functionslij(u), which are used in the error analysis of the finite volume element method:
(24)lij(u)=-∫γij‍A∇(Ihu-u)·ndx,γij=Vi∩Vj.
The following estimate is a simple consequence of the Bramble-Hilbert lemma.Lemma 7 (see [20]).
Ifu∈H2(Ω), then there is a positive constant C>0, independent of  h, such that for eij=∪{K∣K∩γij≠∅,K∈𝒯h},
(25)|lij(u)|≤Ch∥A∥0,∞|u|2,eij.Lemma 8 (see [20]).
(a) Ifu∈H2(Ω), then there exists a positive constant C>0, independent of h and u, such that
(26)|A(u-Ihu,Ih*vh)|≤Ch∥u∥2|vh|1,h,vh∈Sh.
(b) Ifu(·)∈L∞(H2), then, for T>0 fixed there is a constant C=C(T)>0, independent of h and u, such that for 0<t≤T(27)|B(t,s;u-Ihu,Ih*vh)|≤Ch∥u∥2|vh|1,h,vh∈Sh.
For any fixed t>0, one can define the Ritz projection function u(x,t) and the operator Rh:H01∩H2→Sh, such that
(28)A(u-Rhu,Ih*vh)=0,∀vh∈Sh.Remark 9.
If the partition𝒯h* is regular (quasiuniform) and u is H2-regular, then
(29)∥u-Rhu∥1≤Ch∥u∥2.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 Ritz-Volterra projection and has been used in the analysis of the finite element method for integrodifferential equations.We define the Ritz-Volterra projectionVhu of a function u(x,t) defined on the cylinder Ω×[0,T]. The Ritz-Volterra projection Vh:L∞(H01∩H2)→L∞(Sh) is defined for t≥0 by
(30)A(u-Vhu,Ih*vh)+∫0t‍B(t,s;u(s)-Vhu(s),Ih*vh)ds=0,∀vh∈Sh.Lemma 10 (see [20]).
Assume thatu∈L∞(H01∩H2) for T>0 fixed, there is a constant C=C(T)>0 independent of  h and u, such that for all 0<t≤T,
(31)∥(u-Vhu)(t)∥1≤Ch(∥u(t)∥2+∫0t‍∥u(s)∥2ds).Now we consider anL2-estimate for the Ritz-Volterra projection. This estimate is optimal with respect to the order of convergence, but requires W3,p-regularity of the solution. Therefore, it is suboptimal with respect to the regularity of the solution and can be useful for p close to 1. Namely, we have the following result.Lemma 11 (see [20]).
Assume that, for somep>1, u∈L∞(W3,p(Ω)). Then for T>0 fixed there exists a positive constant C=C(T)>0 independent of h and u, such that, for all 0<t≤T,
(32)∥(u-Vhu)(t)∥≤Ch2(∥u(t)∥3,p+∫0t‍∥u(s)∥3,pds).Then we prove the followingL2-norm error estimate.Theorem 12.
Letu and uhn be the solution of (2) and (14), respectively. Then we have a constant C>0 independent of h, Δt, and u, such that
(33)∥un-uhn∥≤Ch2(∥u0∥3,p+∫0tn‍∥ut∥3,pds)+CΔt∫0tn‍(∥utt∥+∥ut∥1)ds.Proof.
Letθn=uhn-Vhun,ρn=un-Vhun, from (2), (14), and (30), we have
(34)(∂tθn,Ih*vh)+A-(θn,Ih*vh)+∑k=0n-1‍wn,kB(tn,tk;θk,Ih*vh)=(utn-∂tun,Ih*vh)+(∂tρn,Ih*vh)+[A-(θn-1,Ih*vh)-A(θn-1,Ih*vh)]-Δt[A-(∂tVhun,Ih*vh)-A(∂tVhun,Ih*vh)]+[{∑k=0n-1‍}∫0tn‍B(tn,s;Vhun(s),Ih*vh)ds-∑k=0n-1‍wn,kB(tn,tk;Vhuk,Ih*vh)].
Let τn=utn-∂tun+∂tρn,
(35)LA--A(θn-1,Ih*vh)=[A-(θn-1,Ih*vh)-A(θn-1,Ih*vh)],LA--A(∂tVhun,Ih*vh)=[A-(∂tVhun,Ih*vh)-A(∂tVhun,Ih*vh){A-}],qn(Vhu,Ih*vh)=[{∑k=0n-1‍}∫0tn‍B(tn,s;Vhun(s),Ih*vh)ds-∑k=0n-1‍wn,kB(tn,tk;Vhuk,Ih*vh)].
Seting vh=θn in (34), we get
(36)(∂tθn,Ih*θn)+A-(θn,Ih*θn)=(τn,Ih*θn)+LA--A(θn-1,Ih*θn)-ΔtLA--A(∂tVhun,Ih*θn)+qn(Vhu,Ih*θn)-∑k=0n-1‍wn,kB(tn,tk;θk,Ih*θn).
We estimate the terms on the right-hand side:
(37)|LA--A(θn-1,Ih*θn)|≤Ch∥A∥1,∞∫Ω‍|∇θn-1||∇θn|dx|LA--A(θn-1,Ih*θn)|≤Ch∥A∥1,∞∥∇θn-1∥∥∇θn∥|LA--A(θn-1,Ih*θn)|≤C∥θn-1∥2+ε|θn|12,|ΔtLA--A(∂tVhun,Ih*θn)|≤ChΔt∥∇∂tVhun∥∥∇θn∥|ΔtLA--A(∂tVhun,Ih*θn)|≤CΔt∥∇∂tVhun∥∥θn∥|ΔtLA--A(∂tVhun,Ih*θn)|≤CΔt∫tn-1tn‍∥ut∥12ds+C∥θn∥2,|qn(Vhu,Ih*θn)|≤ε|θn|12+C(Δt∫0tn‍|DtVhu(s)|1ds)2|qn(Vhu,Ih*θn)|≤ε|θn|12+C(Δt∫0tn‍∥ut(s)∥1ds)2,|∑k=0n-1‍wn,kB(tn,tk;θk,Ih*θn)|≤ε|θn|12+C∑k=0n-1‍Δt|θk|12,|(τn,Ih*θn)|≤C∥θn∥2+CΔt∫tn-1tn‍∥utt(s)∥2ds|(τn,Ih*θn)|+Ch41Δt∫tn-1tn‍∥ut(s)∥3,p2ds.
Let ε=C0/6, and use numerical quadrature error estimates to get
(38)∥θn∥2-∥θn-1∥22Δt+C0|θn|12≤C02|θn|12+C∑k=0n-1‍Δt|θk|12+C(Δt∫0tn‍∥ut(s)∥1ds)2+C∥θn-1∥2+C∥θn∥2+CΔt∫tn-1tn‍∥ut∥12dt+CΔt∫tn-1tn‍∥utt(s)∥2ds+Ch41Δt∫tn-1tn‍∥ut(s)∥3,p2ds.
Thus, multiplying each term on both sides of (38) by Δt, summing over n, and employing Gronwall’s lemma, we obtain
(39)∥θn∥≤C∥θ0∥+Ch2∫0tn‍∥ut∥3,pds+CΔt∫0tn‍(∥utt∥+∥ut∥1)ds.
And ∥θ0∥≤Ch2∥u0∥3,p; 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 Ω=(0,1)×(0,1), T=1, Th be Delaunay triangulation generated by EasyMesh [22] over Ω with mesh size h as shown in Figure 2 and time step be Δt.Figure 2
Triangulation and circumcenter dual mesh.We consider the case ofa11=x1+x2+3, a12=x1+x2+4, a12=a21=-(x1+x2), b11=1, b12=1/2, b12=b21=0, the exact solution u(x,t)=(t2+sin(πt))sin(πx1)sin(πx2), and the initial function u(x,0)=0. Denote the numerical solution of u(x,T) by uh and set γh=∥u(x,T)-u2h∥/∥u(x,T)-uh∥ with fixed ratio Δt/h2. The numerical results are presented in Tables 1 and 2. It is observed that the results support our theory.Table 1
Error behavior forΔt/h2=3.90625e-4.
Δ
t
h
∥
u
(
x
,
T
)
-
u
h
∥
γ
h
1
/
10
1
/
16
6.853070480441
e
-
3
1
/
40
1
/
32
1.690155085494
e
-
3
4.0
1
/
160
1
/
64
4.225622073870
e
-
4
4.0
1
/
640
1
/
128
1.052607231633
e
-
4
4.0Table 2
Error behavior forΔt/h2=1.953125e-4.
Δ
t
h
∥
u
(
x
,
T
)
-
u
h
∥
γ
h
1
/
20
1
/
16
4.671190948063
e
-
3
1
/
80
1
/
32
1.150735715485
e
-
3
4.0
1
/
320
1
/
64
2.880323203978
e
-
4
4.0
1
/
1280
1
/
128
7.164754698408
e
-
5
4.0



---

*Source: 175904-2013-02-13.xml*