- Email: [email protected]
On the stable solution of transient convection–diffusion equations
Journal of Computational and Applied Mathematics 280 (2015) 275–293
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
On the stable solution of transient convection–diffusion equations N.R. Bayramov ∗ , J.K. Kraus Johann Radon Institute for Computational and Applied Mathematics, Austrian Academy of Sciences, Altenberger Straße 69, 4040 Linz, Austria
highlights • • • • •
We consider transient convection–diffusion equation with dominating convection. The model has an application in optimal control of the asymmetric flow field-flow fractionation process. For a stable discretization we propose a monotone edge-averaged finite element (EAFE) scheme. EAFE is generalized to a time-dependent case and a new error estimate is proved. We present numerical results for comparison with a popular SUPG method.
article
info
Article history: Received 5 September 2013 Received in revised form 30 November 2014 Keywords: Transient convection–diffusion equation Monotone finite element scheme Dominating convection Error estimates for time-dependent problem
abstract A transient convection–diffusion equation is considered, which particularly arises in optimization problems for the asymmetric flow field-flow fractionation (AF4) process. A timedependent generalization of the monotone edge-averaged finite element (EAFE) scheme is used to obtain a stable discretization in the convection-dominated regime. New error estimates are proved for this scheme and a comparison with the popular SUPG method is presented. Numerical experiments demonstrate that the EAFE method is more suitable for problems where boundary layers are formed. © 2014 Elsevier B.V. All rights reserved.
1. Introduction The dynamics of the concentration of small particles ruled by an external velocity flow field plays an important role in various applications in natural sciences and industry. Typically such problems are described mathematically by convection–diffusion equations coupled with Stokes or Navier–Stokes equations, see e.g., [1,2]. It is well known that in the convection-dominated regime standard Galerkin approximations suffer from a loss of monotonicity, which may result in spurious oscillations and instabilities [3,2]. These difficulties have already been observed in context of finite difference approximations in the early 1960s. Pioneering works on the analysis of stable numerical schemes for one-dimensional convection-dominated convection–diffusion equations has been conducted by Samarskii [4] in context of finite difference and by Christie et al. [5] in context of finite element discretizations. The Streamline-Upwind Petrov–Galerkin (SUPG) or Streamline-Diffusion (SD) finite element method, introduced by Hughes and Brooks [6], see also [7–9], can be viewed as a multidimensional generalization of such optimal one-dimensional upwind schemes. By adding a consistent diffusion term in streamline direction, this method provides a stable discretization
∗
Corresponding author. E-mail addresses: [email protected] (N.R. Bayramov), [email protected] (J.K. Kraus).
http://dx.doi.org/10.1016/j.cam.2014.12.001 0377-0427/© 2014 Elsevier B.V. All rights reserved.
276
N.R. Bayramov, J.K. Kraus / Journal of Computational and Applied Mathematics 280 (2015) 275–293
tool for convection dominated transport equations. However, SUPG approximations may still develop boundary or internal layers in the crosswind direction (orthogonal to the flow direction) resulting in numerical crosswind layers [10]. Several so-called spurious oscillations at layers diminishing (SOLD) methods have been proposed since the mid 1980s with the aim to remove this drawback [11,12], for an overview see also [13]. Other approaches to a stable numerical solution of (evolutionary) convection–diffusion–reaction equations include discontinuous Galerkin (DG) finite element methods, see, e.g., [14,15] and the references therein, local projection stabilization (LPS) schemes [16,17], which are based on a scale separation and local projections into a larger-scale space, and finite element flux corrected transport schemes (FEM-FCT) [18,19]. Unlike in SUPG, SOLD, DG, or LPS schemes, FEM-FCT achieves the stabilization by modifying the system matrix and the right hand side vector on the algebraic level. For an overview and performance comparison of the above mentioned approaches see [20,21]. There are many applications, for example in semi-conductor device modeling [22,23] or in field-flow fractionation processes [24], where it is desirable that the numerical scheme is monotone, which means that it preserves the maximum principle. Here we consider a transient convection–diffusion problem arising in the optimal control of the asymmetric flow field-flow fractionation (AF4) process, see [25]. The simulation of the AF4 process plays an important role in a number of different applications in medicine, biology, and chemistry. We study the edge-averaged finite element (EAFE) scheme, which has been proposed in [26] and generalized to problems with tensor coefficients in [27], for transient convection–diffusion problems with a main focus on the aforementioned application. The main contribution of this work concerns the application of this monotone scheme to time-dependent problems, the derivation of a new error estimate (for the transient case), and an experimental comparison with the SUPG method. The remainder of the paper is organized as follows. In Section 2 we introduce a mathematical description of the problem. Two stable discretizations (SUPG and EAFE schemes) of the state equation are discussed in Section 3, where we also present a new error estimate for the EAFE scheme and comment on some more general convection–diffusion problems. Section 4 is devoted to numerical experiments. 2. Mathematical model 2.1. Problem statement We consider a flow that satisfies the following Stokes equation:
△V⃗ − ∇ p = 0, ∇ · V⃗ = 0, ⃗ V = g⃗ (u(t )),
in Ω , in Ω , on ∂ Ω .
(1)
Here V⃗ is the velocity field, p is the pressure, g⃗ (u) is a certain prescribed boundary velocity function (depending on a vector parameter u) and Ω is a bounded polygonal domain in Rd . One can see that, determined by the Dirichlet boundary data g⃗ (u(t )), the velocity flow V⃗ (·) = V⃗ (·; u(t )) is a function of a time-dependent control parameter u = u(t ). We study a mass-conserving convection–diffusion equation for the concentration c (x, t ) of particles or a substance driven by the velocity flow V⃗ (x; u(t )) satisfying the system (1), that is,
∂ ∂ t c + ∇ · V⃗ c − ε∇ c = 0, c (x, 0) = c0 (x), ε∇ c − V⃗ c · n⃗ = 0,
(x, t ) ∈ Ω × (0, T ) x∈Ω (x, t ) ∈ ∂ Ω × (0, T )
(2)
⃗ is the outward unit vector normal to the boundary ∂ Ω . where ε is a diffusion coefficient and n According to the boundary condition (last equation in (2)) an appropriate choice of the space for a variational formulation of problem (2) is V = H 1 (Ω ). In other words we are looking for a function c (·, t ) ∈ V which satisfies the system (2) at any time moment t ∈ (0, T ). Multiplying the first equation in (2) by an arbitrary test function ϕ ∈ V , integrating over Ω , using the divergence theorem, and taking into account the incompressibility of the flow V⃗ (i.e., ∇ · V⃗ = 0) yields 0 = ⟨ ∂∂t c , ϕ⟩L2 (Ω ) + ⟨∇ · (V⃗ c − ε∇ c ), ϕ⟩L2 (Ω )
=
∂ ⟨c , ϕ⟩L2 (Ω ) ∂t
+ ⟨ϕ, (V⃗ c − ε∇ c ) · n⃗⟩L2 (∂ Ω ) − ⟨V⃗ c − ε∇ c , ∇ϕ⟩L2 (Ω ) .
Further, by using the boundary condition we obtain
∂ ⟨c (t ), ϕ⟩L2 (Ω ) ∂t
c (0) = c0 ,
+ ε∇ c (t ) − V⃗ c (t ), ∇ϕ L2 (Ω ) = 0 ∀ϕ ∈ V ,
(3)
where the first equation in (3) should be understood in the sense of distributions in (0, T ), that is, for a.e. t ∈ (0, T ), and the initial condition is regarded as an equivalence of two elements of L2 (Ω ).
N.R. Bayramov, J.K. Kraus / Journal of Computational and Applied Mathematics 280 (2015) 275–293
277
Fig. 1. Domain Ω of the system and support of initial condition (as considered in applied optimization problem).
The need to consider the introduced state equation arises, for instance, in the problem of optimal control of a steady accumulation of particles within a subdomain ω ⊂ Ω . A simple form of this abstract problem is given by
c (x, T ) dx −→ min u(t )∈U Ω \ω c (x, T ) dx ≤ κi , i = 1, . . . , M
(4)
ωi
where {κi }M i=1 are given constants, U is a control set and ω = 1≤i≤M ωi . Note that problem (4) in general is rather difficult to solve because of the presence of both, control and state constraints. It can be modified in order to obtain a simpler optimal control problem without state constraints (by adding a penalty term in the functional that is subject to minimization). One of such simplified problems has been studied and solved in [25].
2.2. Existence, uniqueness, and regularity of the solution To formulate the conditions for well-posedness of the Stokes system (1) we follow the article [28]. We define L20 (∂ Ω ) =
⃗f ∈ L2 (∂ Ω ) :
∂Ω
⃗f · n⃗ = 0 .
1/2
Also we denote by Hσ (Ω ) the closure of Cσ∞ (Ω ) in H 1/2 (Ω ), where Cσ∞ (Ω ) = {φ ∈ C ∞ (Ω ) : ∇ · φ = 0}. The following result is a special case of Theorem 2.1 in [28]. Proposition 1. Let Ω be a bounded Lipschitz domain in Rd and let g⃗ ∈ L20 (∂ Ω ). Then system (1) admits a solution (V⃗ , p) ∈ 1/2
d
Hσ (Ω ) ∩ C ∞ (Ω )
× C ∞ (Ω ) and there exists a constant CΩ > 0 such that
∥V⃗ ∥H 1/2 (Ω ) ≤ CΩ ∥⃗g ∥L2 (∂ Ω ) .
(5)
Note that the estimate (5) ensures the uniqueness of the velocity flow V⃗ (x) in the given space. In what follows, the dimension of the phase space is set to d = 2. An example of a domain Ω and an initial condition (in the context of the already mentioned application of optimal control) is illustrated in Fig. 1. If the boundary condition g⃗ (u) is linear with respect to the control vector u: g⃗ (u) = u1 g⃗ (1) + u2 g⃗ (2) + · · · + uM g⃗ (M ) , then according to the linear homogeneity of all equations in (1) w.r.t. (V⃗ , p), a corresponding unique velocity flow V⃗ (x; u) is also linear w.r.t. u: V⃗ (x; u) = u1 V⃗ (1) (x) + u2 V⃗ (2) (x) + · · · + uM V⃗ (M ) (x), where V⃗ (i) (x) corresponds to the system (1) with boundary condition V⃗ ∂ Ω = g⃗ (i) . In this case for any control vector u such
that ess supt ∈[0,T ] |u(t )| ≤ Cu , the corresponding velocity flow V⃗ [t ] = V⃗ (x; u(t )) is uniformly bounded, i.e.,
∥V⃗ [t ]∥C ∞ (Ω ) ≤ cΩ Cu ,
0≤t ≤T
with a constant cΩ > 0 independent of Cu . Based on this fact it is reasonable to assume that u(t ) ∈ L∞ [0, T ].
278
N.R. Bayramov, J.K. Kraus / Journal of Computational and Applied Mathematics 280 (2015) 275–293
Now, let us come back to the variational problem (3). In order to construct a weak solution we follow [29, Chapter 11]. For convenience we denote by ∥ · ∥0 and ∥ · ∥1 the L2 - and H 1 -norm, respectively. We introduce the spaces
L2 0, T ; V = v : (0, T ) → V | v is measurable and
T
∥v(t )∥21 dt < ∞ ,
0
C [0, T ]; L2 (Ω ) = v : [0, T ] → L2 (Ω ) | v is measurable and max ∥v(t )∥0 dt < ∞ . 0
0≤t ≤T
Similarly one can define the space L (0, T ; V ) where V is the dual space of V (i.e., the set of all linear continuous functionals from V to R). An appropriate bilinear form for our problem is ′
2
b(c , ϕ) = ε∇ c − V⃗ c , ∇ϕ
L2 (Ω )
′
.
Note that b(c , ϕ) depends on t (because of V⃗ = V⃗ [t ]), but for simplicity we assume by default that b(c (t ), ϕ) corresponds to the velocity field V⃗ [t ]. Since ∥V⃗ [t ]∥C ∞ (Ω ) is bounded uniformly for t ∈ [0, T ], there exist two constants α > 0 and λ ≥ 0 (independent of t) such that b(v, v) + λ∥v∥20 ≥ α∥v∥21 ,
∀v ∈ V .
(6)
This means that the bilinear form b(·, ·) is continuous and weakly coercive on V × V . In order to study the well-posedness of the variational problem (3), it is rewritten in a slightly more general form with a nonzero right-hand side:
∂
∂t
⟨c (t ), ϕ⟩L2 (Ω ) + b(c (t ), ϕ) = ⟨f (t ), ϕ⟩V ′ ,V
(7)
c (0) = c0 .
The next theorem is an analog of Theorem 11.1.1 in [29]. Theorem 2. Assume that f ∈ L2 (0, T ; V ′ ), c0 ∈ L2 (Ω ) and the bilinear form b(·, ·) is continuous and weakly coercive with
λ = 0 (see condition (6)). ∂c ∂t
Then problem (7) has a unique solution c ∈ L2 (0, T ; V ) ∩ C 0 [0, T ]; L2 (Ω ) . Moreover, there exists a pointwise derivative
∈ L2 (0, T ; V ′ ) and the energy estimate t t 1 2 2 2 ∥f (τ )∥2V ′ dτ ∥c (t )∥0 + α ∥c (τ )∥1 dτ ≤ ∥c0 ∥0 + α 0 0
(8)
holds for each t ∈ [0, T ]. Proof. The estimate (8) is proven in the same way as Theorem 11.1.1 in [29], cf. [29, Remark 11.1.1].
Remark 3. In the general case of weak coercivity with λ > 0 one can use the change of variable cλ (x, t ) = e−λt c (x, t ). Then the new variable cλ (x, t ) satisfies the conditions of Theorem 2, which shows that e−2λt ∥c (t )∥20 + α
t
e−2λτ ∥c (τ )∥21 dτ ≤ ∥c0 ∥20 +
0
1
t
α
e−2λτ ∥f (τ )∥2V ′ dτ
0
and hence
∥c (t )∥ + α 2 0
t
e
2λ(t −τ )
∥c (τ )∥ dτ ≤ e 2 1
2λt
c0 20
∥ ∥ +
0
1
α
t
e2λ(t −τ ) ∥f (τ )∥2V ′ dτ .
0
t
∥c (τ )∥21 dτ is given by t t 1 2 2 2λt 2 ∥c (t )∥0 + α ∥c (τ )∥1 dτ ≤ e ∥c0 ∥0 + e2λ(t −τ ) ∥f (τ )∥2V ′ dτ . α 0 0
As a result, an upper bound for ∥c (t )∥20 + α
0
3. Numerical solution of convection–diffusion equation 3.1. The SUPG method Let Th = {Te }e be a non-degenerate shape-regular subdivision of Ω . As usual, we denote he to be the diameter of element Te and h the maximum diameter of any element in Th . Further let p ≥ 1 be a positive integer and define a conforming finite element space Vh = {vh ∈ H 1 (Ω ) : ∀Te ∈ Th (Ω ), vh |Te ∈ Pp (Te )},
N.R. Bayramov, J.K. Kraus / Journal of Computational and Applied Mathematics 280 (2015) 275–293
279
where Pp (Te ) is the set of polynomials (defined on Te ) of total degree less than or equal to p. Moreover, let (·, ·)Te denote the L2 -inner-product over Te ∈ Th and (·, ·)Ω = Te ∈Th (·, ·)Te . Now we move to a fully discretized formulation of the variational problem (3) using backward Euler time discretization. For ch , vh ∈ Vh , define the bilinear form a(ch , vh ) = △t (ε∇ ch − V⃗ ch , ∇vh )Ω + (ch , vh )Ω .
(9)
we seek an approximation ch (tk ) at time tk = k△t choosing V⃗ = V⃗ [tk ] in (9). We obtain Then for k = 1, . . . , m and △t = the following problem: Find ch (tk ) ∈ Vh such that T , m
a(ch (tk ), vh ) = ch (tk−1 ), vh Ω ∀vh ∈ Vh , (ch (0), vh )Ω = (c0 , vh )Ω ∀vh ∈ Vh .
(10)
It is well known that in case of dominant convection the standard Galerkin finite element formulation (10) becomes unstable and its solution shows spurious oscillations in the whole domain. One possible remedy is to use a streamline-upwind Petrov–Galerkin (SUPG) method, which is designed for convectiondominated problems. We briefly present one particular description of this method here, which can be found for example in [20]. The SUPG method adds a consistent diffusion term Te ∈Th τe △t ‘‘residual’’, V⃗ · ∇vh T in streamline direction to the left e hand side of the first equation of (10). Here {τe }e is a set of parameters depending on the mesh cells Te ∈ Th and the residual is the difference of the left and right hand side of the backward Euler scheme, that is, ‘‘residual’’ = ch (tk ) − ch (tk−1 ) − △t ∇ · ε∇ ch (tk ) − V⃗ ch (tk ) .
Hence, we deduce the following SUPG method from (10): For k = 1, . . . , m find ch (tk ) ∈ Vh such that
τe △t ch (tk ) − ch (tk−1 ), V⃗ · ∇vh Te a(ch (tk ), vh ) + Te ∈Th 2 ⃗ ⃗ τ + e (△t ) −∇ · ε∇ ch (tk ) − V ch (tk ) , V · ∇vh T = ch (tk−1 ), vh Ω , e Te ∈Th (ch (0), vh )Ω = (c0 , vh )Ω
(11)
for all vh ∈ Vh . A crucial question in the application of the SUPG stabilization is the choice of the parameters {τe }e . Among different options presented in [20] we have selected the simplest one (which does not require the estimation of any constants from inverse estimates) for convenience, i.e.,
τeCod =
h2e 4ε△t + 2he △t ∥V⃗ ∥2 + h2e
.
Here ∥V⃗ ∥2 is the Euclidean norm of the vector V⃗ ∈ R2 and he is the size of the mesh cell Ωe . As recommended in [20], the parameter he is chosen as the size of the mesh cell in the direction of convection. Note that the coupling of the SUPG method with implicit time stepping schemes leads to a stable discretization, regardless of the length of the time step, however spurious oscillations may be expected for very small time steps, cf. [20]. 3.2. EAFE scheme In this subsection we recall the edge-averaged finite element (EAFE) scheme that has been introduced in [26] and generalized in [27]. The target of the present work is a time-dependent (transient) problem of the form (3). In the following presentation of the EAFE method we mainly follow the article [26]. We introduce the following notation thereby assuming that the triangulation Th is shape regular: For a given Te ∈ Th denote by
• • • •
qj (1 ≤ j ≤ 3) the vertices of Te ; Eij (or simply E) the edge connecting two vertices qi and qj ; δE φ = φ(qi ) − φ(qj ) for any continuous function φ on E = Eij ; τE = δE x = qi − qj the directional vector of E.
We define the broken Sobolev space H 1 (Th ) = {v ∈ L2 (Ω ) | ∀Te ∈ Th : v ∈ H 1 (Te )} with corresponding norm and semi-norm
∥v∥2H 1 (T ) =
h
Te ∈Th
∥v∥2H 1 (Te ) ,
|v|2H 1 (T ) =
h
Te ∈Th
|v|2H 1 (Te ) .
280
N.R. Bayramov, J.K. Kraus / Journal of Computational and Applied Mathematics 280 (2015) 275–293
Further, for any element Te ∈ Th let the stiffness coefficients (aije ) be defined through the relation T
ε∇ ch · ∇vh dx =
aije ch (qi )vh (qj ) ∀ch , vh ∈ Vh , T
i ,j
Te
where the summation is only over pairs (i, j) with qi , qj ∈ Te . Note that since aiie = − T
ε∇ ch · ∇vh dx = −
T
j̸=i
aije , there holds the identity
aije (ch (qi ) − ch (qj ))(vh (qi ) − vh (qj )), T
i
Te
and thus, by summing over all elements, we obtain
ε∇ ch · ∇vh dx =
Ω
ωETe δE ch δE vh ,
(12)
Te ∈Th E ⊂Te
where ωEe = −aije and E = Eij with i < j. T
T
Consider now any simplex Te ∈ Th and any edge E ∈ Te and let us define a function ψE e on E via T
∂ψETe 1 =− (ε −1 V⃗ · τE ). ∂τE |τE | Te
Next, let HE e (V⃗ ) denote the corresponding harmonic average of e−ψE , that is, T
HE e (V⃗ ) = T
1 −1 Te . eψE ds |τE | E
We set Vh ⊂ V = H 1 (Ω ) to be the usual piecewise linear finite element space and approximate the flux j(c ) ≡ ∇ c − ε −1 V⃗ c by a constant vector jTe (c ) on each Te . Then, based on Lemma 3.1 in [26], we have jTe (c ) · τE ≈ HE e (V⃗ )δE (eψE c ) T
and hence, using (12), for any vh ∈ Vh we get
ε jTe (c ) · ∇vh dx = Te
ωETe (jTe (c ) · τE )δE vh ≈
E ⊂T e
ωETe HETe (V⃗ )δE (eψE c )δE vh .
E ⊂T e
Due to the continuity of V⃗ (x) in Ω we have that HE e ≡ HE . T
Next, define ωE = bh (ch , vh ) =
T e ⊃E
ω
Te E ,
and the approximating bilinear form bh (·, ·) on Vh × Vh by
ωE HE (V⃗ )δE (eψE ch )δE vh .
E ∈Th
Then the following finite element method can be formulated: For all t ∈ (0, T ) find ch (t ) ∈ Vh , such that
∂ ∂t
ch (t ), vh
ch (0) = c0,h
Ω
+ bh (ch (t ), vh ) = 0 ∀vh ∈ Vh , ∀vh ∈ Vh .
(13)
Now we recall another result from [26] regarding the difference between the continuous and the approximating bilinear forms (in terms of our notation), cf. Lemma 6.1 in [26]. Proposition 4. For an arbitrary function c ∈ C (Ω ), let cI ∈ Vh denote its nodal interpolant on the mesh Th . Further, assume that j(c ) ∈ H 1 (Th ). Then the following estimate holds for all vh ∈ Vh :
|b(c , vh ) − bh (cI , vh )| ≤ Cb h|j(c )|H 1 (Th ) ∥vh ∥1 .
(14)
With the estimate (14) at hand, it is easy to prove the well-posedness of problem (13). ∂c
Proposition 5. The problem (13) has a unique solution ch ∈ L2 (0, T ; Vh ) ∩ C 0 ([0, T ]; L2 (Ω )). In addition, ∂ th ∈ L2 (0, T ; Vh′ ) and the following energy estimate holds:
∥ch (t )∥20 + αˆ
t
0
∥ch (τ )∥21 dτ ≤ e2λt ∥c0,h ∥20 .
(15)
N.R. Bayramov, J.K. Kraus / Journal of Computational and Applied Mathematics 280 (2015) 275–293
281
Proof. The statement of the proposition is a particular case of Theorem 2 and the subsequent remark applied to the discrete problem (13) with zero right-hand side. It only remains to show that the bilinear form bh (·, ·) is continuous and weakly coercive on the discrete space Vh . From the proof of Lemma 6.2 in [26] we obtain for all ch , vh ∈ Vh
(16) |b(ch , vh ) − bh (ch , vh )| ≤ Cb h ε + ∥V⃗ ∥W 1,∞ (Th ) ∥ch ∥1 · ∥vh ∥1 , ⃗ ∞ where ∥V⃗ ∥W 1,∞ (Th ) = maxTe ∈Th ∥V⃗ ∥L∞ (Te ) , i ∥∇ Vi ∥L (Te ) . ⃗ Taking into account the regularity properties of V (x), cf. Proposition 1, we deduce that ∥V⃗ ∥W 1,∞ (Th ) is bounded indepen-
dently of h. Next, note that the bilinear form b(·, ·) is continuous and weakly coercive on Vh × Vh for any space Vh ⊂ V with the coercivity constant α > 0. Then (16) implies that for h ≤ h0 , where h0 is sufficiently small, the discrete bilinear form bh (·, ·) is also continuous and weakly coercive on Vh × Vh . The coercivity constant αˆ < α for bh (·, ·) can be chosen arbitrarily close to α if the mesh size is small enough. The following two lemmata will be useful for proving an error estimate later. Lemma 6. Assume that the initial condition c (0) = c0 of the variational problem (7) satisfies c0 ∈ C (Ω ) ∩ L2 (Ω ) and, moreover, c0 is a pointwise non-negative function. Then one can choose a consistent initial condition c0,h ∈ Vh for the discrete problem (13) such that
∥c0,h − cI (0)∥0 ≤ C0 h.
(17)
Proof. Without loss of generality we assume that c0 is a non-zero function. Defining c0,h
c0 (x) dx = cI (0) Ω , c (x, 0) dx Ω I
we obtain the consistency property Ω c0,h (x) dx = Ω c0 (x) dx. Denoting c¯ = Ω (cI (x, 0) − c0,h (x)) dx, from Cauchy–Schwarz inequality we get the estimate
|¯c | ≤ |Ω |1/2 · ∥c0 − cI (0)∥0 .
(18)
Moreover, since c0 ∈ Vh ⊂ V , the interpolation error c0 − cI (0) satisfies the inverse inequality
∥c0 − cI (0)∥0 ≤ C˜ |c0 |1 h.
(19)
Using inequalities (18) and (19), we can estimate
c (x, 0) − c0 (x) dx Ω I ∥c0,h − cI (0)∥0 = ∥cI (0)∥0 c (x, 0)dx Ω I ∥cI (0)∥0 ∥cI (0)∥0 = |¯c | ≤ C˜ |Ω |1/2 |c0 |1 h, c ( x , 0 ) dx c (x, 0)dx Ω I Ω I therefore obtaining (17) with C0 = C˜ |Ω |1/2 |c0 |1
∥cI (0)∥0
Ω cI (x,0) dx
.
In the next lemma let ct , ρt denote the weak time derivatives of the functions (see Definition 15 in the Appendix). Lemma 7. Assume that the solution of the variational problem (7) additionally satisfies the regularity conditions ct ∈ L2 (0, T ; V ),
c (·, t ), ct (·, t ) ∈ C (Ω ) for a.e. t ∈ (0, T ).
Then for the time derivative of the interpolation error ρ(t ) = c (t ) − cI (t ) the following estimate holds:
∥ρt (t )∥0 ≤ Cρ |ct (t )|1 h for a.e. t ∈ (0, T ).
(20)
Proof. Since c , ct ∈ L1 (0, T ; V ), one can use Theorem 19 from the Appendix to obtain c (t2 ) − c (t1 ) =
t2
ct (t ) dt
in V
(21)
t1
t
for a.e. t1 , t2 ∈ (0, T ), where t 2 ct (t ) dt is the Bochner integral. 1 We need to switch in (21) to the nodal interpolants cI (t ), (ct )I (t ) ∈ Vh . Consider an arbitrary approximating sequence of simple functions (see Definition 11 in the Appendix) n→∞
rn (t ) −−−→ ct (t ) in V
282
N.R. Bayramov, J.K. Kraus / Journal of Computational and Applied Mathematics 280 (2015) 275–293
for a.e. t ∈ (0, T ) such that limn→∞ T
ct (t ) dt = lim
0
∥rn (t ) − ct (t )∥V dt = 0. Then by definition
T
n→∞
0
T
rn (t ) dt
in V .
0
From Theorem 4.4.20 in [30] we obtain the boundedness of the nodal interpolant (ct )I (t ) (in the norm ∥ · ∥V ). Hence for the corresponding sequence of interpolants (rn )I (t ) ∈ Vh we have n→∞
(rn )I (t ) −−−→ (ct )I (t ) in Vh T for a.e. t ∈ (0, T ) and limn→∞ 0 ∥(rn )I (t ) − (ct )I (t )∥V dt = 0. This yields T T T rn (t ) dt (rn )I (t ) dt = lim (ct )I (t ) dt = lim 0
=
n→∞
0
n→∞
T
rn (t ) dt
lim
n→∞
0
=
I
0
I
0 T
ct (t ) dt
I
.
Note that such an equality holds on any interval (t1 , t2 ) ⊂ (0, T ). As a result, for a.e. t1 , t2 ∈ (0, T ) cI (t2 ) − cI (t1 ) =
t2
(ct )I (t ) dt in Vh .
(22)
t1
Again, due to the boundedness of the (nodal) interpolation operator (on the space V ), we have cI , (ct )I ∈ L2 (0, T ; V ), and, in particular, cI , (ct )I ∈ L1 (0, T ; V ). Using equality (22) together with Theorem 19, it follows that there exists a weak derivative (cI )t ∈ L1 (0, T ; V ) and (cI )t = (ct )I . Further, from (ct )I ∈ L2 (0, T ; V ), we conclude that (cI )t ∈ L2 (0, T ; V ). Finally, for ρt = ct − (cI )t = ct − (ct )I ∈ L2 (0, T ; V ) we can apply the inverse estimate
∥ρt (t )∥0 ≤ Cρ |ct (t )|1 h, see, e.g., Theorem 4.4.20 in [30].
Theorem 8. Assume that the solution of the variational problem (7) in addition to Lemmata 6–7 satisfies j(c ) ∈ L2 (0, T ; H 1 (Th )). Then for the approximation error (of the EAFE scheme) eh (t ) = ch (t ) − cI (t ) the estimate
∥eh (t )∥20 + α1
t
∥eh (τ )∥21 dτ
1/2
≤ Ceλt h
0
holds for each t ∈ [0, T ]. Proof. Note that according to Theorem 19 from the Appendix, the weak derivative ct ∈ L2 (0, T ; V ) coincides with a (strong) pointwise time derivative of c ∈ L2 (0, T ; V ). Hence for a.e. t ∈ (0, T )
∂ c (t ), ϕ Ω = ct (t ), ϕ Ω , ∂t
∀ϕ ∈ Vh .
The same property holds for the time derivative of cI (t ), see proof of Lemma 7. 1. First, we derive an equation for the error eh (t ): for any ϕ ∈ Vh
∂ ∂ ∂ ∂ eh (t ), ϕ Ω = ch (t ), ϕ Ω − cI (t ), ϕ Ω = −bh (ch (t ), ϕ) − cI (t ), ϕ Ω ∂t ∂t ∂t ∂t ∂ ∂ c (t ), ϕ Ω − cI (t ), ϕ Ω = −bh (ch (t ), ϕ) + b(c (t ), ϕ) + ∂t ∂t = −bh (ch (t ), ϕ) + b(c (t ), ϕ) + (ρt (t ), ϕ)Ω = −bh (eh (t ), ϕ) + [b(c (t ), ϕ) − bh (cI (t ), ϕ)] + (ρt (t ), ϕ)Ω = −b(eh (t ), ϕ) + [b(eh (t ), ϕ) − bh (eh (t ), ϕ)] + [b(c (t ), ϕ) − bh (cI (t ), ϕ)] + (ρt (t ), ϕ)Ω , where ρ(t ) again denotes the interpolation error c (t ) − cI (t ). Introduce a functional f ∈ L2 (0, T ; Vh′ ) defined for a.e. t ∈ (0, T ) by
⟨f (t ), ϕ⟩Vh′ ,Vh = [b(eh (t ), ϕ) − bh (eh (t ), ϕ)] + [b(c (t ), ϕ) − bh (cI (t ), ϕ)] + (ρt (t ), ϕ)Ω . In this case the error eh (t ) satisfies the system (7) with a phase space Vh and initial condition e0 = c0,h − cI (0).
N.R. Bayramov, J.K. Kraus / Journal of Computational and Applied Mathematics 280 (2015) 275–293
283
2. The regularity properties of j(c ) and j(eh ) allow us to apply the estimates from Proposition 4 for the following terms:
|b(c (t ), ϕ) − bh (cI (t ), ϕ)| ≤ Cb |j(c (t ))|H 1 (Th ) h, ∥ϕ∥1 |b(eh (t ), ϕ) − bh (eh (t ), ϕ)| ≤ Cb |j(eh (t ))|H 1 (Th ) h. ∥ϕ∥1 The third term can be estimated using Lemma 7:
|(ρt (t ), ϕ)Ω | ≤ ∥ρt (t )∥0 ≤ Cρ |ct (t )|1 h. ∥ϕ∥1 Then the conditions of the theorem guarantee that
|j(c (t ))|H 1 (Th ) ,
|ct (t )|1 ∈ L2 (0, T ).
From (the proof of) Lemma 6.2 in [26] it follows that
|j(eh (t ))|H 1 (Th ) ≤ Cj ε + ∥V⃗ ∥W 1,∞ (Th ) ∥eh (t )∥1 . As in Proposition 5 one can show that ∥V⃗ ∥W 1,∞ (Th ) is bounded independently of h and thus
|b(eh (t ), ϕ) − bh (eh (t ), ϕ)| ≤ C˜ b ∥eh (t )∥1 h. ∥ϕ∥1 As a result, we have
∥f (t )∥Vh′ = sup
⟨f (t ), ϕ⟩Vh′ ,Vh ∥ϕ∥1
ϕ∈Vh \{0}
≤ (Cf (t ) + C˜ b ∥eh (t )∥1 )h
with Cf (t ) ∈ L2 (0, T ). 3. Finally, analogously to the energy estimate (8) for the solution of (7) we have:
∥eh (t )∥20 + α
t
0
t 1 ∥eh (τ )∥21 dτ ≤ e2λt ∥e0 ∥20 + ∥f (τ )∥2V ′ dτ , h α 0
where e0 = c0,h − cI (0). Since ∥e0 ∥0 = O(h) by Lemma 6, and
∥f (τ )∥2V ′ ≤ 2 Cf2 (τ ) + C˜ b2 ∥eh (τ )∥21 h2 , h
we obtain the required estimate for eh (t ) with
α1 = α − e2λT · 2C˜ b2 /α h2 > α/2 > 0 for sufficiently small h.
3.3. Numerical schemes for more general equations In this section we briefly consider more general convection–diffusion equations for the concentration of particles c (x, t ) that are not mass-conserving. The velocity flow V⃗ (x) in the subsequent formulations remains the same. The first generalization has the following form:
∂ c + ∇ · V⃗ c − ε∇ c + γ c = f , ∂ t c (x, 0) = c0 (x), c = 0, ε∇ c − V⃗ c · n⃗ = 0,
(x, t ) ∈ Ω × (0, T ) x∈Ω (x, t ) ∈ Γ D × (0, T ) (x, t ) ∈ ΓR × (0, T ).
(23)
Here γ ≥ 0 and the domain boundary is split into Dirichlet and Robin parts: ∂ Ω = Γ D ∪ Γ R . Denote HD1 (Ω ) = {v ∈ H 1 (Ω ) : v|ΓD = 0}, where v|ΓD = 0 is interpreted as
v · χΓD = 0 in L2 (∂ Ω ) for the usual characteristic function χ . Using the space V = HD1 (Ω ), we arrive at the former variational problem (7) with analogous regularity properties. In the same manner as earlier one can define a piecewise polynomial finite element space
284
N.R. Bayramov, J.K. Kraus / Journal of Computational and Applied Mathematics 280 (2015) 275–293
Fig. 2. Domain Ω , its subdomain ω and parts of its boundary.
Vh ⊂ HD1 (Ω ). The SUPG and EAFE schemes will remain the same except that the bilinear forms b(c , v) and bh (ch , vh ) on V × V and Vh × Vh will now include additionally a mass term (γ c , v)Ω and (γ ch , vh )Ω , respectively. The second generalization of problem (2), which has been considered in [27], reads:
∂ c + ∇ · V⃗ c − ε∇ c + γ c = f , ∂ t c (x, 0) = c0 (x), c = 0, ε∇ c − V⃗ c · n⃗ = q, ε∇ c · n⃗ = 0,
(x, t ) ∈ Ω × (0, T ) x∈Ω (x, t ) ∈ Γ D × (0, T ) (x, t ) ∈ ΓNin × (0, T ) (x, t ) ∈ ΓNout × (0, T ).
(24)
Here q ∈ L2 (ΓNin × (0, T )) and the Neumann boundary is divided into two parts: ΓN = ΓNin ∪ ΓNout , where
ΓNin = {x ∈ ΓN : n⃗(x) · V⃗ (x) > 0},
ΓNout = {x ∈ ΓN : n⃗(x) · V⃗ (x) ≤ 0}.
In this case the following bilinear form on V × V (with V = HD1 (Ω )) is obtained
⃗)c , v b˜ (c , v) = b(c , v) + (γ c , v)Ω + (V⃗ · n
L2 (ΓNout )
and the weak formulation is given by
∂ ⟨c (t ), v⟩L2 (Ω ) + b˜ (c (t ), v) = ⟨f (t ), v⟩V ′ ,V + ⟨q(t ), v⟩L2 (Γ in ) ∀v ∈ V , N ∂t with initial condition c (0) = c0 . ⃗)c , v L2 (Γ out ) satisfies an additional condition of the form Note that if the boundary term (V⃗ · n N
(V⃗ · n⃗)(c − cI ), v
= O(h),
L2 (ΓNout ×(0,T ))
the error estimate provided by Theorem 8 holds for problems (23) and (24), too. This can be guaranteed for example if c ∈ L2 (0, T ; H 1 (ΓNout )). Moreover, due to the presence of the mass term (γ c , v)Ω , the approximation properties of the bilinear form bh (ch , vh ) for problem (23) (and thus also for the bilinear form bh (ch , vh )) are changed. Following Lemma 6.1 in [26], and imposing the additional regularity condition γ c ∈ L2 (0, T ; W 1,r (Th )) with r > 2, the estimate (14) converts to
1/2 |b(c , vh ) − bh (cI , vh )| ≤ Cb h |j(c )|2H 1 (T ) + |γ c |2W 1,r (T ) ∥vh ∥1 , h
where |v|W 1,r (Th ) =
Te ∈Th
|v|rW 1,r (T
e)
1/r
h
with |v|rW 1,r (T ) = e
Te
|∇v(x)|r dx.
4. Numerical results The numerical results in this section were obtained with the two schemes, EAFE and SUPG, discussed in this paper. The problem setting corresponds to a slightly simplified model as compared to the (real life) industrial application. We consider here a boundary-value problem for the flow-field in a rectangular domain Ω as shown in Fig. 2. The specification of the
N.R. Bayramov, J.K. Kraus / Journal of Computational and Applied Mathematics 280 (2015) 275–293
285
Stokes problem is given in (25). The reason to consider this simple setting is that in this situation the exact solution V⃗ , p of problem (1) is determined explicitly in terms of the problem parameters. Let
∇ · ν∇ V⃗ = ∇ p in Ω , ∇ · V⃗ = 0 in Ω , L2 · (u1 + u2 ) ⟨n, V⃗ ⟩|∂ , = membrane L1 ⟨n, V⃗ ⟩|∂ wall = 0, x2 x2 ⃗ ⟨ n , V ⟩| = − 6 · · 1 − · u1 , ∂ inlet1 L2 L2 ⃗ ⟩|∂ inlet = −6 · x2 · 1 − x2 · u2 , ⟨ n , V 2 L2 L2 ⟨τ , V⃗ ⟩|∂ = 0. Then the solution V⃗ (x), p(x) of the Stokes problem (1) is given by x1 x2 x2 x2 x2 V1 = −6(u1 + u2 ) · 1− + 6u1 1− , L1
L2
L2
L2
(25)
L2
x 2 x 3 L 2 2 2 V2 = −(u1 + u2 ) 1 − 3 +2 · · I{0
p = p0 − 12 ·
u1 · ν L2
·
x1 L2
+6·
L2
(26)
L1
x2 (u1 + u2 ) · ν x21 x2 · 2 + · 1− , L1
L2
L2
L2
where p0 is a constant. It is evident that we do not need to fix the length T of the entire time period in the numerical experiments. For the time discretization in the EAFE scheme we use the backward Euler method since it is unconditionally stable. Thus we do not need to relate the time step △t to the mesh size h. For both schemes, EAFE and SUPG, the resulting linear systems of equations were solved with the sparse direct solver UMFPACK integrated in the software package FreeFem++. 4.1. Model examples for testing the EAFE method Consider a generalization of the problem (2) with nonzero right-hand side f and non-homogeneous boundary conditions. For any given model solution cˆ ∈ H 1 (0, T ; V ) one can calculate explicitly the corresponding right-hand side functional f ∈ L2 (0, T ; V ′ ) as in (7). We can then compute the numerical solution ch (x, t ) for a given initial condition c0 (x) using the EAFE method and calculate the error ∥eh (t )∥ in different norms. The following model solutions were considered: cˆ (1) (x, t ) = x1 · (1 + sin(x21 + x22 ))e−t , 2 cˆ (2) (x, t ) = x2 · (1 + sin(x21 + x22 ))e−t /2 .
The parameters of the system are chosen as follows:
• • • •
L1 = 5, L2 = 4 — sizes of the rectangular domain Ω , ε ∈ {10−2 , 10−1 , 1, 10} — values of diffusion coefficient, △t = 10−2 — time step, (u1 , u2 ) = (3, 7) — control function is constant.
For clarity we have computed the ratio log(∥eh (t )∥0,Ω /h), which is bounded above by C +λt (according to Theorem 8). For comparison, Figs. 3–6 show also the growth of log(∥eh (t )∥1,Ω /h) and log(∥eh (t )∥∞,Ω /h) with respect to t. These numerical tests were run only for the mesh size h ≈ 10−3 since for smaller h > 0 the behavior of log(∥eh (t )∥/h) changes only slightly. We enhance the plots with the estimated theoretical maximum growth rate λ > 0 for the values of log(∥eh (t )∥0,Ω /h). Note that the error estimate provided in Theorem 8 covers the test problems. Estimates for ∥eh (t )∥1,Ω or ∥eh (t )∥∞,Ω , but under stronger assumptions, have also been provided in [29, Proposition 11.1.1, Remark 11.2.4]. However, it can observed on the Figs. 3–6 that the error in H 1 - and L∞ -norms behaves similarly in the particular case of exact solutions cˆ (1) (x, t ), cˆ (2) (x, t ). In addition, the growth rate of log(∥eh (t )∥/h) at any time moment does not exceed our theoretical upper bound λ, thus the behavior of the errors complies with the theory. Finally, the tests show that the errors increase for smaller diffusion coefficients ε , which means that in the convection-dominated regime one should choose a sufficiently small mesh size h.
286
N.R. Bayramov, J.K. Kraus / Journal of Computational and Applied Mathematics 280 (2015) 275–293
a
b
Fig. 3. Evolution of the ratios log(∥eh (t )∥0,Ω /h), log(∥eh (t )∥1,Ω /h) and log(∥eh (t )∥∞,Ω /h) for ε = 10−2 and an exact solution (a) cˆ (1) (x, t ), (b) cˆ (2) (x, t ).
a
b
Fig. 4. Evolution of the ratios log(∥eh (t )∥0,Ω /h), log(∥eh (t )∥1,Ω /h) and log(∥eh (t )∥∞,Ω /h) for ε = 10−1 and an exact solution (a) cˆ (1) (x, t ), (b) cˆ (2) (x, t ).
a
b
Fig. 5. Evolution of the ratios log(∥eh (t )∥0,Ω /h), log(∥eh (t )∥1,Ω /h) and log(∥eh (t )∥∞,Ω /h) for ε = 1 and an exact solution (a) cˆ (1) (x, t ), (b) cˆ (2) (x, t ).
4.2. Comparison of the two schemes The SUPG and the EAFE approximations were compared for the transient convection–diffusion problem that is driven by the velocity field V⃗ (x) from (26). The method and problem parameters were chosen as follows:
N.R. Bayramov, J.K. Kraus / Journal of Computational and Applied Mathematics 280 (2015) 275–293
a
287
b
Fig. 6. Evolution of the ratios log(∥eh (t )∥0,Ω /h), log(∥eh (t )∥1,Ω /h) and log(∥eh (t )∥∞,Ω /h) for ε = 10 and an exact solution (a) cˆ (1) (x, t ), (b) cˆ (2) (x, t ).
Fig. 7. c (x, t ) for t = 0.05, u = (2.1, 8.5) obtained by (a) EAFE, (b) SUPG method.
• L1 = 5, L2 = 0.05 — sizes of the rectangular domain Ω , • ε = 4 · 10−6 — diffusion coefficient, • △t = 5 · 10−3 — time step. Since the considered problem arises from technical applications, one of which is described and studied in [25], the control parameter u is selected such that after a certain period of time (not larger then T ) the concentration is accumulated mostly within a subdomain ω. We were interested only in the part of the mesh Th that corresponds to a truncated domain
Ωθ = [0, θ · L1 ] × [0, L2 ],
θ = 0.2,
which covers the subdomain ω and can be regarded as the support of c (·, t ) at any time moment t ∈ [0, T ]. The concentration c (x, t ) was evaluated only on the truncated domain as it is physically negligible on the rest of the domain. A uniform triangulation Th of the domain Ωθ into triangles was used with nx = 1000 and ny = 500 nodes along the x-axis and y-axis, respectively. For the SUPG method the finite element space Vh of piecewise linear functions was used (higher order approximations, such as for p = 2, have not resulted in the disappearance of small oscillations). In the subsequent plots of the solutions (illustrated by isolines) it is clearly seen that the EAFE and SUPG schemes produce very similar stable solutions ch (tk ) before a boundary layer is formed (driven by the velocity field V⃗ for a properly selected control parameter u(t )). In this situation spurious oscillations of the discrete solution at the boundary layer can be observed for the SUPG method. On the other hand, the monotone EAFE scheme is free from these defects. This additional stability and also the ease of implementation, in particular, the fact that the EAFE scheme does not require any stabilization parameters, certainly makes this method a good candidate for solving transient convection–diffusion problems, as they were considered in the present work. See Figs. 7–11.
288
N.R. Bayramov, J.K. Kraus / Journal of Computational and Applied Mathematics 280 (2015) 275–293
Fig. 8. c (x, t ) for t = 0.1, u = (2.1, 8.5) obtained by (a) EAFE, (b) SUPG method.
Fig. 9. c (x, t ) for t = 0.4, u = (2.1, 8.5) obtained by (a) EAFE, (b) SUPG method.
Fig. 10. c (x, t ) for t = 0.8, u = (1.6, 8.5) with switching moment T1 = 0.5, obtained by (a) EAFE, (b) SUPG method.
5. Conclusions We have compared two stable finite element schemes (SUPG and EAFE) for solving transient convection-dominated convection–diffusion problems. In general, both schemes provide similar results (for isotropic diffusion). However, for certain choices of the control parameters the SUPG solution exhibits spurious oscillations at boundary layers. This does not happen for the monotone EAFE scheme (namely, for its time-dependent generalization). This explains the choice of the numerical method made in [25] for the numerical solution of optimal control problems arising in the asymmetric flow
N.R. Bayramov, J.K. Kraus / Journal of Computational and Applied Mathematics 280 (2015) 275–293
a
289
b
Fig. 11. c (x, t ) for t = 1.0, u = (1.6, 8.5) with switching moment T1 = 0.5, obtained by (a) EAFE, (b) SUPG method.
field-flow fractionation (AF4) process. To justify the applicability of the method, a new error estimate has been derived for the EAFE scheme in the transient case. The design and analysis of efficient iterative methods for solving the arising systems of linear algebraic equations is subject of ongoing and future research. Acknowledgments The authors are grateful to Prof. Oleg Iliev, who initiated the present work motivated by an industrial application and to Prof. Ludmil Zikatanov for his valuable comments on the EAFE scheme. Thanks also go to Tigran Nagapetyan for providing the analytical expression for the solution of the Stokes system (25). This work has been supported by the Austrian Science Fund (FWF), Grant P22989-N18. Appendix In this appendix, we summarize some results about the integration and differentiation of functions f : (0, T ) → X , where X is a separable Banach space with norm ∥ · ∥, cf. [31,32]. Definition 9. A simple function f : (0, T ) → X is a function of the form f =
N
cj χEj
j =1
where E1 , . . . , EN are Lebesgue measurable subsets of (0, T ), c1 , . . . , cN ∈ X and
χE (t ) =
1, 0,
if t ∈ E if t ̸∈ E
denotes the characteristic function of E. Definition 10. A function f : (0, T ) → X is (strongly) measurable, if there is a sequence {fn } of simple functions such that fn (t ) → f (t ) strongly in X (i.e. in norm) for t a.e. in (0, T ). Definition 11. A strongly measurable function f : (0, T ) → X is (Bochner) integrable, if there is a sequence of simple functions such that fn (t ) → f (t ) pointwise a.e. in (0, T ) and T
∥f − fn ∥ dt = 0.
lim
n→∞
0
The integral of f is defined by T
T
fn dt ,
f dt = lim 0
n→∞
0
where the limit exists strongly in X .
290
N.R. Bayramov, J.K. Kraus / Journal of Computational and Applied Mathematics 280 (2015) 275–293
The value of the Bochner integral of f is independent of the sequence {fn } of approximating simple functions (which can be shown analogously to the case of real-valued functions), and
0
T
T
f dt ≤
∥f ∥ dt . 0
Definition 12. For 1 ≤ p < ∞ the space Lp (0, T ; X ) consists of all strongly measurable functions f : (0, T ) → X such that T
∥f ∥p dt < ∞ 0
and it is equipped with the norm T
∥f ∥Lp (0,T ;X ) =
∥f ∥p dt
1/p
.
0
The following result (see [32]) characterizes (Bochner) integrable functions as elements of L1 (0, T ; X ). Theorem 13 (Bochner, 1933). A function f : (0, T ) → X is Bochner integrable if and only if it is strongly measurable and T
∥f ∥ dt < ∞. 0
To prove the main theorem of this appendix, we need to introduce several auxiliary lemmata. Lemma 14. Suppose that fn : (0, T ) → X is Bochner integrable for each n ∈ N, fn (t ) → f (s) as n → ∞ strongly in X for t a.e. in (0, T ), and there is an integrable function g : (0, T ) → R such that either
∥fn (t )∥ ≤ g (t ) for t a.e. in (0, T ) and every n ∈ N, or the sequence {∥fn (t )∥}n∈N is uniformly integrable. Then f : (0, T ) → X is Bochner integrable and T
T
0
T
f dt ,
fn dt →
∥fn − f ∥ dt → 0 as n → ∞.
0
0
Proof. The proof is the same as for real-valued functions.
Let L1loc (0, T ; X ) denote the space of measurable functions f : (0, T ) → X that are integrable on every compactly supported interval (a, b), i.e., 0 < a < b < T . Also, let Cc∞ (0, T ) denote the space of smooth, real-valued functions φ : (0, T ) → R with compact support in (0, T ). Definition 15. A function f ∈ L1loc (0, T ; X ) is weakly differentiable with weak derivative ft = g ∈ L1loc (0, T ; X ) if T
φ f dt = − ′
0
T
φ g dt for every φ ∈ Cc∞ (0, T ). 0
Lemma 16. Suppose that f ∈ L1 (0, T ; X ), then f (t ) = lim
h→0
1 h
t +h
f (s) ds t
for t pointwise a.e. in (0, T ). Proof. Since X is a separable space, let {cn ∈ X : n ∈ N} be a dense subset of X . Then by the Lebesgue differentiation theorem for real-valued functions
∥f (t ) − cn ∥ = lim
h→0
1 h
t +h
∥f (s) − cn ∥ ds t
for every n ∈ N and t pointwise a.e. in (0, T ). Thus, for all such t ∈ (0, T ) and every n ∈ N, we have lim sup h→0
1 h
t +h
∥f (s) − f (t )∥ ds ≤ lim sup t
h→0
1 h
t +h
∥f (s) − cn ∥ + ∥f (t ) − cn ∥ ds = 2∥f (t ) − cn ∥. t
N.R. Bayramov, J.K. Kraus / Journal of Computational and Applied Mathematics 280 (2015) 275–293
291
Since this holds for every element cn of the dense subset, it follows that 1
lim sup
h
h→0
t +h
∥f (s) − f (t )∥ ds = 0. t
Therefore
1 lim sup h
h→0
t +h
1 f (s) ds − f (t ) ≤ lim sup h→0
t
which proves the result.
h
t +h
∥f (s) − f (t )∥ ds = 0, t
Lemma 17. Suppose that f : (0, T ) → X is locally integrable and T
φ f dt = 0 for every φ ∈ Cc∞ (0, T ). 0
Then f = 0 pointwise a.e. on (0, T ). Proof. Choose a sequence of test functions 0 ≤ φn ≤ 1 whose supports are contained inside a fixed compact subset of (0, T ) such that φn → χ(t ,t +h) pointwise, where χ(t ,t +h) is the characteristic function of the interval (t , t + h) ⊂ (0, T ). If f ∈ L1loc (0, T ; X ), then we can apply the analog of the dominated convergence theorem, Lemma 14, to obtain t +h
f (s) ds = lim
n→∞
t
But
T
φn (s)f (s) ds. 0
T 0
φ f ds = 0 for every φ ∈ Cc∞ (0, T ), therefore t +h f (s) ds = 0 t
for every (t , t + h) ⊂ (0, T ). It then follows from Lemma 16 that f = 0 pointwise a.e. in (0, T ).
Lemma 18. Suppose that f : (0, T ) → X is weakly differentiable and f ′ = 0. Then f is equivalent to a constant function. Proof. The condition that the weak derivative f ′ is zero means that (see Definition 15) T
f φ ′ dt = 0
for all φ ∈ Cc∞ (0, T ).
(27)
0
Choose a fixed test function η ∈ Cc∞ (0, T ) whose integral is equal to one, and represent an arbitrary test function φ ∈ Cc∞ (0, T ) as
φ = K η + ψ ′, where K ∈ R and ψ ∈ Cc∞ (0, T ) are given by T
φ dt ,
K =
ψ(t ) =
0
t
φ(s) − K η(s) ds.
0
If we define T
ηf dt ∈ X ,
c= 0
then (27) implies that T
(f − c )φ dt = 0 for all φ ∈ Cc∞ (0, T ), 0
and Lemma 17 implies that f = c pointwise a.e. on (0, T ).
Theorem 19. Suppose that f ∈ L1 (0, T ; X ). Then f is weakly differentiable with integrable derivative ft = g ∈ L1 (0, T ; X ) if and only if f (t ) = c0 +
t
g (s) ds
(28)
0
pointwise a.e. in (0, T ). In that case, f is differentiable pointwise a.e. and its pointwise (strong) derivative coincides with its weak derivative.
292
N.R. Bayramov, J.K. Kraus / Journal of Computational and Applied Mathematics 280 (2015) 275–293
Proof. If f is given by (28), then f (t + h) − f (t ) h
=
1
t +h
h
g (s) ds, t
and Lemma 16 implies that the strong derivative of f exists pointwise a.e. and is equal to g. We also have that
f ( t + h) − f ( t ) 1 ≤ h
h
t +h
∥g (s)∥ ds.
(29)
t
Extending f by zero to a function f : R → X , and using Fubini’s theorem, we get
1 f (t + h) − f (t ) dt ≤ h
R
h
≤
R
t +h t h
1 h
R
∥g (s)∥ ds dt
∥g (s + t )∥ ds dt ≤ 0
1
h
h
0
∥g (s + t )∥ dt ds =
∥g (t )∥ dt .
(30)
R
R
Choose φ ∈ Cc∞ (0, T ). Due to the analog of the dominated convergence theorem, Lemma 14, T
φ ′ (t )f (t ) dt = lim
T φ(t + h) − φ(t )
h→0
0
h
0
f (t ) dt .
Further we can reorder the terms in the integration (since φ is zero on [0, h1 ] ∪ [T − h2 , T ]): lim
T φ(t + h) − φ(t )
h→0
h
0
f (t ) dt = − lim
T
h→0
f (t ) − f (t − h) φ(t ) dt . h
0
The estimates (29) and (30) imply that the class of functions
f (t +h)−f (t ) h
apply again Lemma 14:
f (t ) − f (t − h) dt = − φ(t )
T
− lim
h→0
h
0
is uniformly integrable, which allows us to |h|
T
φ(t )g (t ) dt . 0
This finally means that g is the weak derivative of f : T
φ ′ (t )f (t ) dt = − 0
T
φ(t )g (t ) dt . 0
Conversely, if ft = g ∈ L1 (0, T ; X ) in the sense of weak derivatives, let
f (t ) =
t
g (s) ds. 0
ft = g, so the weak derivative (f − f )t is zero. Lemma 18 then implies that f − f Then the previous argument implies that is constant pointwise a.e., which gives (28). References [1] S. Flotron, J. Rappaz, Conservation schemes for convection–diffusion equations with Robins boundary conditions. MATHICSE Technical Report Nr. 11.2012, École Polytechnique Fédérale de Lausanne, 2012, March. [2] H.-G. Roos, M. Stynes, L. Tobiska, Robust Numerical Methods for Singularly Perturbed Differential Equations, second ed., in: Springer Series in Computational Mathematics, vol. 24, Springer-Verlag, Berlin, 2008. [3] C. Johnson, Numerical Solution of Partial Differential Equations by the Finite Element Method, Dover Publications, Inc., Mineola, NY, 2009, Reprint of the 1987 edition. [4] A.A. Samarskii, Theory of Difference Schemes, Nauka, Moscow, 1977. [5] I. Christie, D.F. Griffiths, A.R. Mitchell, O.C. Zienkiewicz, Finite element methods for second order differential equations with significant first derivatives, Internat. J. Numer. Methods Engrg. 10 (1976) 1389–1396. [6] T.J.R. Hughes, A.N. Brooks, A multidimensional upwind scheme with no crosswind diffusion. Finite element methods for convection dominated flows, in: AMD, 34, Amer. Soc. Mech. Engrs. (ASME), New York, 1979, pp. 19–35. [7] A.N. Brooks, T.J.R. Hughes, Streamline upwind/Petrov–Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier–Stokes equations, Comput. Methods Appl. Mech. Engrg. 32 (1–3) (1982) 199–259. [8] V. John, J.M. Maubach, L. Tobiska, Nonconforming streamline-diffusion-finite-element-methods for convection–diffusion problems, Numer. Math. 78 (2) (1997) 165–188. [9] V. John, J. Novo, Error analysis of the SUPG finite element discretization of evolutionary convection–diffusion–reaction equations, SIAM J. Numer. Anal. 49 (3) (2011) 1149–1176. [10] B. Semper, Numerical crosswind smear in the streamline diffusion method, Comput. Methods Appl. Mech. Engrg. 113 (1–2) (1994) 99–108. [11] V. John, P. Knobloch, On the performance of SOLD methods for convection–diffusion problems with interior layers, Int. J. Comput. Sci. Math. 1 (2–4) (2007) 245–258.
N.R. Bayramov, J.K. Kraus / Journal of Computational and Applied Mathematics 280 (2015) 275–293
293
[12] V. John, P. Knobloch, On spurious oscillations at layers diminishing (SOLD) methods for convection–diffusion equations: part II—analysis for P1 and Q1 finite elements, Comput. Methods Appl. Mech. Engrg. 197 (2008) 1997–2014. [13] V. John, P. Knobloch, On spurious oscillations at layers diminishing (SOLD) methods for convection–diffusion equations: part I—a review, Comput. Methods Appl. Mech. Engrg. 196 (2007) 2197–2215. [14] B. Ayuso, D.L. Marini, Discontinuous Galerkin methods for advection-diffusion–reaction problems, SIAM J. Numer. Anal. 47 (2) (2009) 1391–1420. [15] J. Proft, B.M. Rivière, Discontinuous Galerkin methods for convection–diffusion equations for varying and vanishing diffusivity, Int. J. Numer. Anal. Model. 6 (4) (2009) 533–561. [16] M. Braack, E. Burman, Local projection stabilization for the Oseen problem and its interpretation as a variational multiscale method, SIAM J. Numer. Anal. 43 (2006) 2544–2566. [17] G. Matthies, P. Skrzypacz, L. Tobiska, A unified convergence analysis for local projection stabilisations applied to the Oseen problem, M2AN Math. Model. Numer. Anal. 41 (2007) 713–742. [18] R. Löhner, K. Morgan, J. Peraire, M. Vahdati, Finite element flux-corrected transport (FEM-FCT) for the Euler and Navier–Stokes equations, Internat. J. Numer. Methods Fluids 7 (1987) 1093–1109. [19] D. Kuzmin, M. Möller, S. Turek, High-resolution FEM-FCT schemes for multidimensional conservation laws, Comput. Methods Appl. Mech. Engrg. 193 (2004) 4915–4946. [20] V. John, E. Schmeyer, Finite element methods for time-dependent convection–diffusion–reaction equations with small diffusion, Comput. Methods Appl. Mech. Engrg. 198 (2008) 475–494. [21] M. Augustin, A. Caiazzo, A. Fiebach, J. Fuhrmann, V. John, A. Linke, R. Umla, An assessment of discretizations for convection-dominated convection–diffusion equations, Comput. Methods Appl. Mech. Engrg. 200 (2011) 3395–3409. [22] F. Brezzi, L.D. Marini, P. Pietra, Numerical simulation of semiconductor devices, in: Proceedings of the Eighth International Conference on Computing Methods in Applied Sciences and Engineering (Versailles, 1987), vol. 75, 1989, pp. 493–514. [23] F. Brezzi, L.D. Marini, P. Pietra, Two-dimensional exponential fitting and applications to drift-diffusion models, SIAM J. Numer. Anal. 26 (1989) 1342–1355. [24] M.E. Schimpf, K.D. Caldwell, J.C. Giddings, The FFF Handbook, Wiley, New York, 2000. [25] N. Bayramov, T. Nagapetyan, R. Pinnau, Fast optimal control of asymmetric flow field flow fractionation processes, in: Proceedings of SIAM Conference on Control and Its Applications, San Diego, California (2013), pp. 207–213. [26] J. Xu, L.T. Zikatanov, A monotone finite element scheme for convection–diffusion equations, Math. Comp. 68 (228) (1999) 1429–1446. [27] R.D. Lazarov, L.T. Zikatanov, An exponential fitting scheme for general convection–diffusion equations on tetrahedral meshes, Comput. Appl. Math. 1 (92) (2005) 60–69. [28] A. Russo, G. Starita, On the existence of steady-state solutions to the Navier–Stokes system for large fluxes, Ann. Sc. Norm. Super. Pisa Cl. Sci. 7 (2008) 171–180. [29] A. Quarteroni, A. Valli, Numerical Approximation of Partial Differential Equations, in: Springer Series in Computational Mathematics, vol. 23, SpringerVerlag, 1994. [30] S.C. Brenner, L.R. Scott, The Mathematical Theory of Finite Element Methods, third ed., Springer, 2008. [31] S. Bochner, Integration von Funktionen deren Werte die Elemente eines Vektorraumes sind, Fund. Math. 20 (1933) 262–276. [32] K. Yosida, Theory of Difference Schemes, Springer-Verlag, Berlin, 1980.
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
On the stable solution of transient convection–diffusion equations N.R. Bayramov ∗ , J.K. Kraus Johann Radon Institute for Computational and Applied Mathematics, Austrian Academy of Sciences, Altenberger Straße 69, 4040 Linz, Austria
highlights • • • • •
We consider transient convection–diffusion equation with dominating convection. The model has an application in optimal control of the asymmetric flow field-flow fractionation process. For a stable discretization we propose a monotone edge-averaged finite element (EAFE) scheme. EAFE is generalized to a time-dependent case and a new error estimate is proved. We present numerical results for comparison with a popular SUPG method.
article
info
Article history: Received 5 September 2013 Received in revised form 30 November 2014 Keywords: Transient convection–diffusion equation Monotone finite element scheme Dominating convection Error estimates for time-dependent problem
abstract A transient convection–diffusion equation is considered, which particularly arises in optimization problems for the asymmetric flow field-flow fractionation (AF4) process. A timedependent generalization of the monotone edge-averaged finite element (EAFE) scheme is used to obtain a stable discretization in the convection-dominated regime. New error estimates are proved for this scheme and a comparison with the popular SUPG method is presented. Numerical experiments demonstrate that the EAFE method is more suitable for problems where boundary layers are formed. © 2014 Elsevier B.V. All rights reserved.
1. Introduction The dynamics of the concentration of small particles ruled by an external velocity flow field plays an important role in various applications in natural sciences and industry. Typically such problems are described mathematically by convection–diffusion equations coupled with Stokes or Navier–Stokes equations, see e.g., [1,2]. It is well known that in the convection-dominated regime standard Galerkin approximations suffer from a loss of monotonicity, which may result in spurious oscillations and instabilities [3,2]. These difficulties have already been observed in context of finite difference approximations in the early 1960s. Pioneering works on the analysis of stable numerical schemes for one-dimensional convection-dominated convection–diffusion equations has been conducted by Samarskii [4] in context of finite difference and by Christie et al. [5] in context of finite element discretizations. The Streamline-Upwind Petrov–Galerkin (SUPG) or Streamline-Diffusion (SD) finite element method, introduced by Hughes and Brooks [6], see also [7–9], can be viewed as a multidimensional generalization of such optimal one-dimensional upwind schemes. By adding a consistent diffusion term in streamline direction, this method provides a stable discretization
∗
Corresponding author. E-mail addresses: [email protected] (N.R. Bayramov), [email protected] (J.K. Kraus).
http://dx.doi.org/10.1016/j.cam.2014.12.001 0377-0427/© 2014 Elsevier B.V. All rights reserved.
276
N.R. Bayramov, J.K. Kraus / Journal of Computational and Applied Mathematics 280 (2015) 275–293
tool for convection dominated transport equations. However, SUPG approximations may still develop boundary or internal layers in the crosswind direction (orthogonal to the flow direction) resulting in numerical crosswind layers [10]. Several so-called spurious oscillations at layers diminishing (SOLD) methods have been proposed since the mid 1980s with the aim to remove this drawback [11,12], for an overview see also [13]. Other approaches to a stable numerical solution of (evolutionary) convection–diffusion–reaction equations include discontinuous Galerkin (DG) finite element methods, see, e.g., [14,15] and the references therein, local projection stabilization (LPS) schemes [16,17], which are based on a scale separation and local projections into a larger-scale space, and finite element flux corrected transport schemes (FEM-FCT) [18,19]. Unlike in SUPG, SOLD, DG, or LPS schemes, FEM-FCT achieves the stabilization by modifying the system matrix and the right hand side vector on the algebraic level. For an overview and performance comparison of the above mentioned approaches see [20,21]. There are many applications, for example in semi-conductor device modeling [22,23] or in field-flow fractionation processes [24], where it is desirable that the numerical scheme is monotone, which means that it preserves the maximum principle. Here we consider a transient convection–diffusion problem arising in the optimal control of the asymmetric flow field-flow fractionation (AF4) process, see [25]. The simulation of the AF4 process plays an important role in a number of different applications in medicine, biology, and chemistry. We study the edge-averaged finite element (EAFE) scheme, which has been proposed in [26] and generalized to problems with tensor coefficients in [27], for transient convection–diffusion problems with a main focus on the aforementioned application. The main contribution of this work concerns the application of this monotone scheme to time-dependent problems, the derivation of a new error estimate (for the transient case), and an experimental comparison with the SUPG method. The remainder of the paper is organized as follows. In Section 2 we introduce a mathematical description of the problem. Two stable discretizations (SUPG and EAFE schemes) of the state equation are discussed in Section 3, where we also present a new error estimate for the EAFE scheme and comment on some more general convection–diffusion problems. Section 4 is devoted to numerical experiments. 2. Mathematical model 2.1. Problem statement We consider a flow that satisfies the following Stokes equation:
△V⃗ − ∇ p = 0, ∇ · V⃗ = 0, ⃗ V = g⃗ (u(t )),
in Ω , in Ω , on ∂ Ω .
(1)
Here V⃗ is the velocity field, p is the pressure, g⃗ (u) is a certain prescribed boundary velocity function (depending on a vector parameter u) and Ω is a bounded polygonal domain in Rd . One can see that, determined by the Dirichlet boundary data g⃗ (u(t )), the velocity flow V⃗ (·) = V⃗ (·; u(t )) is a function of a time-dependent control parameter u = u(t ). We study a mass-conserving convection–diffusion equation for the concentration c (x, t ) of particles or a substance driven by the velocity flow V⃗ (x; u(t )) satisfying the system (1), that is,
∂ ∂ t c + ∇ · V⃗ c − ε∇ c = 0, c (x, 0) = c0 (x), ε∇ c − V⃗ c · n⃗ = 0,
(x, t ) ∈ Ω × (0, T ) x∈Ω (x, t ) ∈ ∂ Ω × (0, T )
(2)
⃗ is the outward unit vector normal to the boundary ∂ Ω . where ε is a diffusion coefficient and n According to the boundary condition (last equation in (2)) an appropriate choice of the space for a variational formulation of problem (2) is V = H 1 (Ω ). In other words we are looking for a function c (·, t ) ∈ V which satisfies the system (2) at any time moment t ∈ (0, T ). Multiplying the first equation in (2) by an arbitrary test function ϕ ∈ V , integrating over Ω , using the divergence theorem, and taking into account the incompressibility of the flow V⃗ (i.e., ∇ · V⃗ = 0) yields 0 = ⟨ ∂∂t c , ϕ⟩L2 (Ω ) + ⟨∇ · (V⃗ c − ε∇ c ), ϕ⟩L2 (Ω )
=
∂ ⟨c , ϕ⟩L2 (Ω ) ∂t
+ ⟨ϕ, (V⃗ c − ε∇ c ) · n⃗⟩L2 (∂ Ω ) − ⟨V⃗ c − ε∇ c , ∇ϕ⟩L2 (Ω ) .
Further, by using the boundary condition we obtain
∂ ⟨c (t ), ϕ⟩L2 (Ω ) ∂t
c (0) = c0 ,
+ ε∇ c (t ) − V⃗ c (t ), ∇ϕ L2 (Ω ) = 0 ∀ϕ ∈ V ,
(3)
where the first equation in (3) should be understood in the sense of distributions in (0, T ), that is, for a.e. t ∈ (0, T ), and the initial condition is regarded as an equivalence of two elements of L2 (Ω ).
N.R. Bayramov, J.K. Kraus / Journal of Computational and Applied Mathematics 280 (2015) 275–293
277
Fig. 1. Domain Ω of the system and support of initial condition (as considered in applied optimization problem).
The need to consider the introduced state equation arises, for instance, in the problem of optimal control of a steady accumulation of particles within a subdomain ω ⊂ Ω . A simple form of this abstract problem is given by
c (x, T ) dx −→ min u(t )∈U Ω \ω c (x, T ) dx ≤ κi , i = 1, . . . , M
(4)
ωi
where {κi }M i=1 are given constants, U is a control set and ω = 1≤i≤M ωi . Note that problem (4) in general is rather difficult to solve because of the presence of both, control and state constraints. It can be modified in order to obtain a simpler optimal control problem without state constraints (by adding a penalty term in the functional that is subject to minimization). One of such simplified problems has been studied and solved in [25].
2.2. Existence, uniqueness, and regularity of the solution To formulate the conditions for well-posedness of the Stokes system (1) we follow the article [28]. We define L20 (∂ Ω ) =
⃗f ∈ L2 (∂ Ω ) :
∂Ω
⃗f · n⃗ = 0 .
1/2
Also we denote by Hσ (Ω ) the closure of Cσ∞ (Ω ) in H 1/2 (Ω ), where Cσ∞ (Ω ) = {φ ∈ C ∞ (Ω ) : ∇ · φ = 0}. The following result is a special case of Theorem 2.1 in [28]. Proposition 1. Let Ω be a bounded Lipschitz domain in Rd and let g⃗ ∈ L20 (∂ Ω ). Then system (1) admits a solution (V⃗ , p) ∈ 1/2
d
Hσ (Ω ) ∩ C ∞ (Ω )
× C ∞ (Ω ) and there exists a constant CΩ > 0 such that
∥V⃗ ∥H 1/2 (Ω ) ≤ CΩ ∥⃗g ∥L2 (∂ Ω ) .
(5)
Note that the estimate (5) ensures the uniqueness of the velocity flow V⃗ (x) in the given space. In what follows, the dimension of the phase space is set to d = 2. An example of a domain Ω and an initial condition (in the context of the already mentioned application of optimal control) is illustrated in Fig. 1. If the boundary condition g⃗ (u) is linear with respect to the control vector u: g⃗ (u) = u1 g⃗ (1) + u2 g⃗ (2) + · · · + uM g⃗ (M ) , then according to the linear homogeneity of all equations in (1) w.r.t. (V⃗ , p), a corresponding unique velocity flow V⃗ (x; u) is also linear w.r.t. u: V⃗ (x; u) = u1 V⃗ (1) (x) + u2 V⃗ (2) (x) + · · · + uM V⃗ (M ) (x), where V⃗ (i) (x) corresponds to the system (1) with boundary condition V⃗ ∂ Ω = g⃗ (i) . In this case for any control vector u such
that ess supt ∈[0,T ] |u(t )| ≤ Cu , the corresponding velocity flow V⃗ [t ] = V⃗ (x; u(t )) is uniformly bounded, i.e.,
∥V⃗ [t ]∥C ∞ (Ω ) ≤ cΩ Cu ,
0≤t ≤T
with a constant cΩ > 0 independent of Cu . Based on this fact it is reasonable to assume that u(t ) ∈ L∞ [0, T ].
278
N.R. Bayramov, J.K. Kraus / Journal of Computational and Applied Mathematics 280 (2015) 275–293
Now, let us come back to the variational problem (3). In order to construct a weak solution we follow [29, Chapter 11]. For convenience we denote by ∥ · ∥0 and ∥ · ∥1 the L2 - and H 1 -norm, respectively. We introduce the spaces
L2 0, T ; V = v : (0, T ) → V | v is measurable and
T
∥v(t )∥21 dt < ∞ ,
0
C [0, T ]; L2 (Ω ) = v : [0, T ] → L2 (Ω ) | v is measurable and max ∥v(t )∥0 dt < ∞ . 0
0≤t ≤T
Similarly one can define the space L (0, T ; V ) where V is the dual space of V (i.e., the set of all linear continuous functionals from V to R). An appropriate bilinear form for our problem is ′
2
b(c , ϕ) = ε∇ c − V⃗ c , ∇ϕ
L2 (Ω )
′
.
Note that b(c , ϕ) depends on t (because of V⃗ = V⃗ [t ]), but for simplicity we assume by default that b(c (t ), ϕ) corresponds to the velocity field V⃗ [t ]. Since ∥V⃗ [t ]∥C ∞ (Ω ) is bounded uniformly for t ∈ [0, T ], there exist two constants α > 0 and λ ≥ 0 (independent of t) such that b(v, v) + λ∥v∥20 ≥ α∥v∥21 ,
∀v ∈ V .
(6)
This means that the bilinear form b(·, ·) is continuous and weakly coercive on V × V . In order to study the well-posedness of the variational problem (3), it is rewritten in a slightly more general form with a nonzero right-hand side:
∂
∂t
⟨c (t ), ϕ⟩L2 (Ω ) + b(c (t ), ϕ) = ⟨f (t ), ϕ⟩V ′ ,V
(7)
c (0) = c0 .
The next theorem is an analog of Theorem 11.1.1 in [29]. Theorem 2. Assume that f ∈ L2 (0, T ; V ′ ), c0 ∈ L2 (Ω ) and the bilinear form b(·, ·) is continuous and weakly coercive with
λ = 0 (see condition (6)). ∂c ∂t
Then problem (7) has a unique solution c ∈ L2 (0, T ; V ) ∩ C 0 [0, T ]; L2 (Ω ) . Moreover, there exists a pointwise derivative
∈ L2 (0, T ; V ′ ) and the energy estimate t t 1 2 2 2 ∥f (τ )∥2V ′ dτ ∥c (t )∥0 + α ∥c (τ )∥1 dτ ≤ ∥c0 ∥0 + α 0 0
(8)
holds for each t ∈ [0, T ]. Proof. The estimate (8) is proven in the same way as Theorem 11.1.1 in [29], cf. [29, Remark 11.1.1].
Remark 3. In the general case of weak coercivity with λ > 0 one can use the change of variable cλ (x, t ) = e−λt c (x, t ). Then the new variable cλ (x, t ) satisfies the conditions of Theorem 2, which shows that e−2λt ∥c (t )∥20 + α
t
e−2λτ ∥c (τ )∥21 dτ ≤ ∥c0 ∥20 +
0
1
t
α
e−2λτ ∥f (τ )∥2V ′ dτ
0
and hence
∥c (t )∥ + α 2 0
t
e
2λ(t −τ )
∥c (τ )∥ dτ ≤ e 2 1
2λt
c0 20
∥ ∥ +
0
1
α
t
e2λ(t −τ ) ∥f (τ )∥2V ′ dτ .
0
t
∥c (τ )∥21 dτ is given by t t 1 2 2 2λt 2 ∥c (t )∥0 + α ∥c (τ )∥1 dτ ≤ e ∥c0 ∥0 + e2λ(t −τ ) ∥f (τ )∥2V ′ dτ . α 0 0
As a result, an upper bound for ∥c (t )∥20 + α
0
3. Numerical solution of convection–diffusion equation 3.1. The SUPG method Let Th = {Te }e be a non-degenerate shape-regular subdivision of Ω . As usual, we denote he to be the diameter of element Te and h the maximum diameter of any element in Th . Further let p ≥ 1 be a positive integer and define a conforming finite element space Vh = {vh ∈ H 1 (Ω ) : ∀Te ∈ Th (Ω ), vh |Te ∈ Pp (Te )},
N.R. Bayramov, J.K. Kraus / Journal of Computational and Applied Mathematics 280 (2015) 275–293
279
where Pp (Te ) is the set of polynomials (defined on Te ) of total degree less than or equal to p. Moreover, let (·, ·)Te denote the L2 -inner-product over Te ∈ Th and (·, ·)Ω = Te ∈Th (·, ·)Te . Now we move to a fully discretized formulation of the variational problem (3) using backward Euler time discretization. For ch , vh ∈ Vh , define the bilinear form a(ch , vh ) = △t (ε∇ ch − V⃗ ch , ∇vh )Ω + (ch , vh )Ω .
(9)
we seek an approximation ch (tk ) at time tk = k△t choosing V⃗ = V⃗ [tk ] in (9). We obtain Then for k = 1, . . . , m and △t = the following problem: Find ch (tk ) ∈ Vh such that T , m
a(ch (tk ), vh ) = ch (tk−1 ), vh Ω ∀vh ∈ Vh , (ch (0), vh )Ω = (c0 , vh )Ω ∀vh ∈ Vh .
(10)
It is well known that in case of dominant convection the standard Galerkin finite element formulation (10) becomes unstable and its solution shows spurious oscillations in the whole domain. One possible remedy is to use a streamline-upwind Petrov–Galerkin (SUPG) method, which is designed for convectiondominated problems. We briefly present one particular description of this method here, which can be found for example in [20]. The SUPG method adds a consistent diffusion term Te ∈Th τe △t ‘‘residual’’, V⃗ · ∇vh T in streamline direction to the left e hand side of the first equation of (10). Here {τe }e is a set of parameters depending on the mesh cells Te ∈ Th and the residual is the difference of the left and right hand side of the backward Euler scheme, that is, ‘‘residual’’ = ch (tk ) − ch (tk−1 ) − △t ∇ · ε∇ ch (tk ) − V⃗ ch (tk ) .
Hence, we deduce the following SUPG method from (10): For k = 1, . . . , m find ch (tk ) ∈ Vh such that
τe △t ch (tk ) − ch (tk−1 ), V⃗ · ∇vh Te a(ch (tk ), vh ) + Te ∈Th 2 ⃗ ⃗ τ + e (△t ) −∇ · ε∇ ch (tk ) − V ch (tk ) , V · ∇vh T = ch (tk−1 ), vh Ω , e Te ∈Th (ch (0), vh )Ω = (c0 , vh )Ω
(11)
for all vh ∈ Vh . A crucial question in the application of the SUPG stabilization is the choice of the parameters {τe }e . Among different options presented in [20] we have selected the simplest one (which does not require the estimation of any constants from inverse estimates) for convenience, i.e.,
τeCod =
h2e 4ε△t + 2he △t ∥V⃗ ∥2 + h2e
.
Here ∥V⃗ ∥2 is the Euclidean norm of the vector V⃗ ∈ R2 and he is the size of the mesh cell Ωe . As recommended in [20], the parameter he is chosen as the size of the mesh cell in the direction of convection. Note that the coupling of the SUPG method with implicit time stepping schemes leads to a stable discretization, regardless of the length of the time step, however spurious oscillations may be expected for very small time steps, cf. [20]. 3.2. EAFE scheme In this subsection we recall the edge-averaged finite element (EAFE) scheme that has been introduced in [26] and generalized in [27]. The target of the present work is a time-dependent (transient) problem of the form (3). In the following presentation of the EAFE method we mainly follow the article [26]. We introduce the following notation thereby assuming that the triangulation Th is shape regular: For a given Te ∈ Th denote by
• • • •
qj (1 ≤ j ≤ 3) the vertices of Te ; Eij (or simply E) the edge connecting two vertices qi and qj ; δE φ = φ(qi ) − φ(qj ) for any continuous function φ on E = Eij ; τE = δE x = qi − qj the directional vector of E.
We define the broken Sobolev space H 1 (Th ) = {v ∈ L2 (Ω ) | ∀Te ∈ Th : v ∈ H 1 (Te )} with corresponding norm and semi-norm
∥v∥2H 1 (T ) =
h
Te ∈Th
∥v∥2H 1 (Te ) ,
|v|2H 1 (T ) =
h
Te ∈Th
|v|2H 1 (Te ) .
280
N.R. Bayramov, J.K. Kraus / Journal of Computational and Applied Mathematics 280 (2015) 275–293
Further, for any element Te ∈ Th let the stiffness coefficients (aije ) be defined through the relation T
ε∇ ch · ∇vh dx =
aije ch (qi )vh (qj ) ∀ch , vh ∈ Vh , T
i ,j
Te
where the summation is only over pairs (i, j) with qi , qj ∈ Te . Note that since aiie = − T
ε∇ ch · ∇vh dx = −
T
j̸=i
aije , there holds the identity
aije (ch (qi ) − ch (qj ))(vh (qi ) − vh (qj )), T
i
Te
and thus, by summing over all elements, we obtain
ε∇ ch · ∇vh dx =
Ω
ωETe δE ch δE vh ,
(12)
Te ∈Th E ⊂Te
where ωEe = −aije and E = Eij with i < j. T
T
Consider now any simplex Te ∈ Th and any edge E ∈ Te and let us define a function ψE e on E via T
∂ψETe 1 =− (ε −1 V⃗ · τE ). ∂τE |τE | Te
Next, let HE e (V⃗ ) denote the corresponding harmonic average of e−ψE , that is, T
HE e (V⃗ ) = T
1 −1 Te . eψE ds |τE | E
We set Vh ⊂ V = H 1 (Ω ) to be the usual piecewise linear finite element space and approximate the flux j(c ) ≡ ∇ c − ε −1 V⃗ c by a constant vector jTe (c ) on each Te . Then, based on Lemma 3.1 in [26], we have jTe (c ) · τE ≈ HE e (V⃗ )δE (eψE c ) T
and hence, using (12), for any vh ∈ Vh we get
ε jTe (c ) · ∇vh dx = Te
ωETe (jTe (c ) · τE )δE vh ≈
E ⊂T e
ωETe HETe (V⃗ )δE (eψE c )δE vh .
E ⊂T e
Due to the continuity of V⃗ (x) in Ω we have that HE e ≡ HE . T
Next, define ωE = bh (ch , vh ) =
T e ⊃E
ω
Te E ,
and the approximating bilinear form bh (·, ·) on Vh × Vh by
ωE HE (V⃗ )δE (eψE ch )δE vh .
E ∈Th
Then the following finite element method can be formulated: For all t ∈ (0, T ) find ch (t ) ∈ Vh , such that
∂ ∂t
ch (t ), vh
ch (0) = c0,h
Ω
+ bh (ch (t ), vh ) = 0 ∀vh ∈ Vh , ∀vh ∈ Vh .
(13)
Now we recall another result from [26] regarding the difference between the continuous and the approximating bilinear forms (in terms of our notation), cf. Lemma 6.1 in [26]. Proposition 4. For an arbitrary function c ∈ C (Ω ), let cI ∈ Vh denote its nodal interpolant on the mesh Th . Further, assume that j(c ) ∈ H 1 (Th ). Then the following estimate holds for all vh ∈ Vh :
|b(c , vh ) − bh (cI , vh )| ≤ Cb h|j(c )|H 1 (Th ) ∥vh ∥1 .
(14)
With the estimate (14) at hand, it is easy to prove the well-posedness of problem (13). ∂c
Proposition 5. The problem (13) has a unique solution ch ∈ L2 (0, T ; Vh ) ∩ C 0 ([0, T ]; L2 (Ω )). In addition, ∂ th ∈ L2 (0, T ; Vh′ ) and the following energy estimate holds:
∥ch (t )∥20 + αˆ
t
0
∥ch (τ )∥21 dτ ≤ e2λt ∥c0,h ∥20 .
(15)
N.R. Bayramov, J.K. Kraus / Journal of Computational and Applied Mathematics 280 (2015) 275–293
281
Proof. The statement of the proposition is a particular case of Theorem 2 and the subsequent remark applied to the discrete problem (13) with zero right-hand side. It only remains to show that the bilinear form bh (·, ·) is continuous and weakly coercive on the discrete space Vh . From the proof of Lemma 6.2 in [26] we obtain for all ch , vh ∈ Vh
(16) |b(ch , vh ) − bh (ch , vh )| ≤ Cb h ε + ∥V⃗ ∥W 1,∞ (Th ) ∥ch ∥1 · ∥vh ∥1 , ⃗ ∞ where ∥V⃗ ∥W 1,∞ (Th ) = maxTe ∈Th ∥V⃗ ∥L∞ (Te ) , i ∥∇ Vi ∥L (Te ) . ⃗ Taking into account the regularity properties of V (x), cf. Proposition 1, we deduce that ∥V⃗ ∥W 1,∞ (Th ) is bounded indepen-
dently of h. Next, note that the bilinear form b(·, ·) is continuous and weakly coercive on Vh × Vh for any space Vh ⊂ V with the coercivity constant α > 0. Then (16) implies that for h ≤ h0 , where h0 is sufficiently small, the discrete bilinear form bh (·, ·) is also continuous and weakly coercive on Vh × Vh . The coercivity constant αˆ < α for bh (·, ·) can be chosen arbitrarily close to α if the mesh size is small enough. The following two lemmata will be useful for proving an error estimate later. Lemma 6. Assume that the initial condition c (0) = c0 of the variational problem (7) satisfies c0 ∈ C (Ω ) ∩ L2 (Ω ) and, moreover, c0 is a pointwise non-negative function. Then one can choose a consistent initial condition c0,h ∈ Vh for the discrete problem (13) such that
∥c0,h − cI (0)∥0 ≤ C0 h.
(17)
Proof. Without loss of generality we assume that c0 is a non-zero function. Defining c0,h
c0 (x) dx = cI (0) Ω , c (x, 0) dx Ω I
we obtain the consistency property Ω c0,h (x) dx = Ω c0 (x) dx. Denoting c¯ = Ω (cI (x, 0) − c0,h (x)) dx, from Cauchy–Schwarz inequality we get the estimate
|¯c | ≤ |Ω |1/2 · ∥c0 − cI (0)∥0 .
(18)
Moreover, since c0 ∈ Vh ⊂ V , the interpolation error c0 − cI (0) satisfies the inverse inequality
∥c0 − cI (0)∥0 ≤ C˜ |c0 |1 h.
(19)
Using inequalities (18) and (19), we can estimate
c (x, 0) − c0 (x) dx Ω I ∥c0,h − cI (0)∥0 = ∥cI (0)∥0 c (x, 0)dx Ω I ∥cI (0)∥0 ∥cI (0)∥0 = |¯c | ≤ C˜ |Ω |1/2 |c0 |1 h, c ( x , 0 ) dx c (x, 0)dx Ω I Ω I therefore obtaining (17) with C0 = C˜ |Ω |1/2 |c0 |1
∥cI (0)∥0
Ω cI (x,0) dx
.
In the next lemma let ct , ρt denote the weak time derivatives of the functions (see Definition 15 in the Appendix). Lemma 7. Assume that the solution of the variational problem (7) additionally satisfies the regularity conditions ct ∈ L2 (0, T ; V ),
c (·, t ), ct (·, t ) ∈ C (Ω ) for a.e. t ∈ (0, T ).
Then for the time derivative of the interpolation error ρ(t ) = c (t ) − cI (t ) the following estimate holds:
∥ρt (t )∥0 ≤ Cρ |ct (t )|1 h for a.e. t ∈ (0, T ).
(20)
Proof. Since c , ct ∈ L1 (0, T ; V ), one can use Theorem 19 from the Appendix to obtain c (t2 ) − c (t1 ) =
t2
ct (t ) dt
in V
(21)
t1
t
for a.e. t1 , t2 ∈ (0, T ), where t 2 ct (t ) dt is the Bochner integral. 1 We need to switch in (21) to the nodal interpolants cI (t ), (ct )I (t ) ∈ Vh . Consider an arbitrary approximating sequence of simple functions (see Definition 11 in the Appendix) n→∞
rn (t ) −−−→ ct (t ) in V
282
N.R. Bayramov, J.K. Kraus / Journal of Computational and Applied Mathematics 280 (2015) 275–293
for a.e. t ∈ (0, T ) such that limn→∞ T
ct (t ) dt = lim
0
∥rn (t ) − ct (t )∥V dt = 0. Then by definition
T
n→∞
0
T
rn (t ) dt
in V .
0
From Theorem 4.4.20 in [30] we obtain the boundedness of the nodal interpolant (ct )I (t ) (in the norm ∥ · ∥V ). Hence for the corresponding sequence of interpolants (rn )I (t ) ∈ Vh we have n→∞
(rn )I (t ) −−−→ (ct )I (t ) in Vh T for a.e. t ∈ (0, T ) and limn→∞ 0 ∥(rn )I (t ) − (ct )I (t )∥V dt = 0. This yields T T T rn (t ) dt (rn )I (t ) dt = lim (ct )I (t ) dt = lim 0
=
n→∞
0
n→∞
T
rn (t ) dt
lim
n→∞
0
=
I
0
I
0 T
ct (t ) dt
I
.
Note that such an equality holds on any interval (t1 , t2 ) ⊂ (0, T ). As a result, for a.e. t1 , t2 ∈ (0, T ) cI (t2 ) − cI (t1 ) =
t2
(ct )I (t ) dt in Vh .
(22)
t1
Again, due to the boundedness of the (nodal) interpolation operator (on the space V ), we have cI , (ct )I ∈ L2 (0, T ; V ), and, in particular, cI , (ct )I ∈ L1 (0, T ; V ). Using equality (22) together with Theorem 19, it follows that there exists a weak derivative (cI )t ∈ L1 (0, T ; V ) and (cI )t = (ct )I . Further, from (ct )I ∈ L2 (0, T ; V ), we conclude that (cI )t ∈ L2 (0, T ; V ). Finally, for ρt = ct − (cI )t = ct − (ct )I ∈ L2 (0, T ; V ) we can apply the inverse estimate
∥ρt (t )∥0 ≤ Cρ |ct (t )|1 h, see, e.g., Theorem 4.4.20 in [30].
Theorem 8. Assume that the solution of the variational problem (7) in addition to Lemmata 6–7 satisfies j(c ) ∈ L2 (0, T ; H 1 (Th )). Then for the approximation error (of the EAFE scheme) eh (t ) = ch (t ) − cI (t ) the estimate
∥eh (t )∥20 + α1
t
∥eh (τ )∥21 dτ
1/2
≤ Ceλt h
0
holds for each t ∈ [0, T ]. Proof. Note that according to Theorem 19 from the Appendix, the weak derivative ct ∈ L2 (0, T ; V ) coincides with a (strong) pointwise time derivative of c ∈ L2 (0, T ; V ). Hence for a.e. t ∈ (0, T )
∂ c (t ), ϕ Ω = ct (t ), ϕ Ω , ∂t
∀ϕ ∈ Vh .
The same property holds for the time derivative of cI (t ), see proof of Lemma 7. 1. First, we derive an equation for the error eh (t ): for any ϕ ∈ Vh
∂ ∂ ∂ ∂ eh (t ), ϕ Ω = ch (t ), ϕ Ω − cI (t ), ϕ Ω = −bh (ch (t ), ϕ) − cI (t ), ϕ Ω ∂t ∂t ∂t ∂t ∂ ∂ c (t ), ϕ Ω − cI (t ), ϕ Ω = −bh (ch (t ), ϕ) + b(c (t ), ϕ) + ∂t ∂t = −bh (ch (t ), ϕ) + b(c (t ), ϕ) + (ρt (t ), ϕ)Ω = −bh (eh (t ), ϕ) + [b(c (t ), ϕ) − bh (cI (t ), ϕ)] + (ρt (t ), ϕ)Ω = −b(eh (t ), ϕ) + [b(eh (t ), ϕ) − bh (eh (t ), ϕ)] + [b(c (t ), ϕ) − bh (cI (t ), ϕ)] + (ρt (t ), ϕ)Ω , where ρ(t ) again denotes the interpolation error c (t ) − cI (t ). Introduce a functional f ∈ L2 (0, T ; Vh′ ) defined for a.e. t ∈ (0, T ) by
⟨f (t ), ϕ⟩Vh′ ,Vh = [b(eh (t ), ϕ) − bh (eh (t ), ϕ)] + [b(c (t ), ϕ) − bh (cI (t ), ϕ)] + (ρt (t ), ϕ)Ω . In this case the error eh (t ) satisfies the system (7) with a phase space Vh and initial condition e0 = c0,h − cI (0).
N.R. Bayramov, J.K. Kraus / Journal of Computational and Applied Mathematics 280 (2015) 275–293
283
2. The regularity properties of j(c ) and j(eh ) allow us to apply the estimates from Proposition 4 for the following terms:
|b(c (t ), ϕ) − bh (cI (t ), ϕ)| ≤ Cb |j(c (t ))|H 1 (Th ) h, ∥ϕ∥1 |b(eh (t ), ϕ) − bh (eh (t ), ϕ)| ≤ Cb |j(eh (t ))|H 1 (Th ) h. ∥ϕ∥1 The third term can be estimated using Lemma 7:
|(ρt (t ), ϕ)Ω | ≤ ∥ρt (t )∥0 ≤ Cρ |ct (t )|1 h. ∥ϕ∥1 Then the conditions of the theorem guarantee that
|j(c (t ))|H 1 (Th ) ,
|ct (t )|1 ∈ L2 (0, T ).
From (the proof of) Lemma 6.2 in [26] it follows that
|j(eh (t ))|H 1 (Th ) ≤ Cj ε + ∥V⃗ ∥W 1,∞ (Th ) ∥eh (t )∥1 . As in Proposition 5 one can show that ∥V⃗ ∥W 1,∞ (Th ) is bounded independently of h and thus
|b(eh (t ), ϕ) − bh (eh (t ), ϕ)| ≤ C˜ b ∥eh (t )∥1 h. ∥ϕ∥1 As a result, we have
∥f (t )∥Vh′ = sup
⟨f (t ), ϕ⟩Vh′ ,Vh ∥ϕ∥1
ϕ∈Vh \{0}
≤ (Cf (t ) + C˜ b ∥eh (t )∥1 )h
with Cf (t ) ∈ L2 (0, T ). 3. Finally, analogously to the energy estimate (8) for the solution of (7) we have:
∥eh (t )∥20 + α
t
0
t 1 ∥eh (τ )∥21 dτ ≤ e2λt ∥e0 ∥20 + ∥f (τ )∥2V ′ dτ , h α 0
where e0 = c0,h − cI (0). Since ∥e0 ∥0 = O(h) by Lemma 6, and
∥f (τ )∥2V ′ ≤ 2 Cf2 (τ ) + C˜ b2 ∥eh (τ )∥21 h2 , h
we obtain the required estimate for eh (t ) with
α1 = α − e2λT · 2C˜ b2 /α h2 > α/2 > 0 for sufficiently small h.
3.3. Numerical schemes for more general equations In this section we briefly consider more general convection–diffusion equations for the concentration of particles c (x, t ) that are not mass-conserving. The velocity flow V⃗ (x) in the subsequent formulations remains the same. The first generalization has the following form:
∂ c + ∇ · V⃗ c − ε∇ c + γ c = f , ∂ t c (x, 0) = c0 (x), c = 0, ε∇ c − V⃗ c · n⃗ = 0,
(x, t ) ∈ Ω × (0, T ) x∈Ω (x, t ) ∈ Γ D × (0, T ) (x, t ) ∈ ΓR × (0, T ).
(23)
Here γ ≥ 0 and the domain boundary is split into Dirichlet and Robin parts: ∂ Ω = Γ D ∪ Γ R . Denote HD1 (Ω ) = {v ∈ H 1 (Ω ) : v|ΓD = 0}, where v|ΓD = 0 is interpreted as
v · χΓD = 0 in L2 (∂ Ω ) for the usual characteristic function χ . Using the space V = HD1 (Ω ), we arrive at the former variational problem (7) with analogous regularity properties. In the same manner as earlier one can define a piecewise polynomial finite element space
284
N.R. Bayramov, J.K. Kraus / Journal of Computational and Applied Mathematics 280 (2015) 275–293
Fig. 2. Domain Ω , its subdomain ω and parts of its boundary.
Vh ⊂ HD1 (Ω ). The SUPG and EAFE schemes will remain the same except that the bilinear forms b(c , v) and bh (ch , vh ) on V × V and Vh × Vh will now include additionally a mass term (γ c , v)Ω and (γ ch , vh )Ω , respectively. The second generalization of problem (2), which has been considered in [27], reads:
∂ c + ∇ · V⃗ c − ε∇ c + γ c = f , ∂ t c (x, 0) = c0 (x), c = 0, ε∇ c − V⃗ c · n⃗ = q, ε∇ c · n⃗ = 0,
(x, t ) ∈ Ω × (0, T ) x∈Ω (x, t ) ∈ Γ D × (0, T ) (x, t ) ∈ ΓNin × (0, T ) (x, t ) ∈ ΓNout × (0, T ).
(24)
Here q ∈ L2 (ΓNin × (0, T )) and the Neumann boundary is divided into two parts: ΓN = ΓNin ∪ ΓNout , where
ΓNin = {x ∈ ΓN : n⃗(x) · V⃗ (x) > 0},
ΓNout = {x ∈ ΓN : n⃗(x) · V⃗ (x) ≤ 0}.
In this case the following bilinear form on V × V (with V = HD1 (Ω )) is obtained
⃗)c , v b˜ (c , v) = b(c , v) + (γ c , v)Ω + (V⃗ · n
L2 (ΓNout )
and the weak formulation is given by
∂ ⟨c (t ), v⟩L2 (Ω ) + b˜ (c (t ), v) = ⟨f (t ), v⟩V ′ ,V + ⟨q(t ), v⟩L2 (Γ in ) ∀v ∈ V , N ∂t with initial condition c (0) = c0 . ⃗)c , v L2 (Γ out ) satisfies an additional condition of the form Note that if the boundary term (V⃗ · n N
(V⃗ · n⃗)(c − cI ), v
= O(h),
L2 (ΓNout ×(0,T ))
the error estimate provided by Theorem 8 holds for problems (23) and (24), too. This can be guaranteed for example if c ∈ L2 (0, T ; H 1 (ΓNout )). Moreover, due to the presence of the mass term (γ c , v)Ω , the approximation properties of the bilinear form bh (ch , vh ) for problem (23) (and thus also for the bilinear form bh (ch , vh )) are changed. Following Lemma 6.1 in [26], and imposing the additional regularity condition γ c ∈ L2 (0, T ; W 1,r (Th )) with r > 2, the estimate (14) converts to
1/2 |b(c , vh ) − bh (cI , vh )| ≤ Cb h |j(c )|2H 1 (T ) + |γ c |2W 1,r (T ) ∥vh ∥1 , h
where |v|W 1,r (Th ) =
Te ∈Th
|v|rW 1,r (T
e)
1/r
h
with |v|rW 1,r (T ) = e
Te
|∇v(x)|r dx.
4. Numerical results The numerical results in this section were obtained with the two schemes, EAFE and SUPG, discussed in this paper. The problem setting corresponds to a slightly simplified model as compared to the (real life) industrial application. We consider here a boundary-value problem for the flow-field in a rectangular domain Ω as shown in Fig. 2. The specification of the
N.R. Bayramov, J.K. Kraus / Journal of Computational and Applied Mathematics 280 (2015) 275–293
285
Stokes problem is given in (25). The reason to consider this simple setting is that in this situation the exact solution V⃗ , p of problem (1) is determined explicitly in terms of the problem parameters. Let
∇ · ν∇ V⃗ = ∇ p in Ω , ∇ · V⃗ = 0 in Ω , L2 · (u1 + u2 ) ⟨n, V⃗ ⟩|∂ , = membrane L1 ⟨n, V⃗ ⟩|∂ wall = 0, x2 x2 ⃗ ⟨ n , V ⟩| = − 6 · · 1 − · u1 , ∂ inlet1 L2 L2 ⃗ ⟩|∂ inlet = −6 · x2 · 1 − x2 · u2 , ⟨ n , V 2 L2 L2 ⟨τ , V⃗ ⟩|∂ = 0. Then the solution V⃗ (x), p(x) of the Stokes problem (1) is given by x1 x2 x2 x2 x2 V1 = −6(u1 + u2 ) · 1− + 6u1 1− , L1
L2
L2
L2
(25)
L2
x 2 x 3 L 2 2 2 V2 = −(u1 + u2 ) 1 − 3 +2 · · I{0
p = p0 − 12 ·
u1 · ν L2
·
x1 L2
+6·
L2
(26)
L1
x2 (u1 + u2 ) · ν x21 x2 · 2 + · 1− , L1
L2
L2
L2
where p0 is a constant. It is evident that we do not need to fix the length T of the entire time period in the numerical experiments. For the time discretization in the EAFE scheme we use the backward Euler method since it is unconditionally stable. Thus we do not need to relate the time step △t to the mesh size h. For both schemes, EAFE and SUPG, the resulting linear systems of equations were solved with the sparse direct solver UMFPACK integrated in the software package FreeFem++. 4.1. Model examples for testing the EAFE method Consider a generalization of the problem (2) with nonzero right-hand side f and non-homogeneous boundary conditions. For any given model solution cˆ ∈ H 1 (0, T ; V ) one can calculate explicitly the corresponding right-hand side functional f ∈ L2 (0, T ; V ′ ) as in (7). We can then compute the numerical solution ch (x, t ) for a given initial condition c0 (x) using the EAFE method and calculate the error ∥eh (t )∥ in different norms. The following model solutions were considered: cˆ (1) (x, t ) = x1 · (1 + sin(x21 + x22 ))e−t , 2 cˆ (2) (x, t ) = x2 · (1 + sin(x21 + x22 ))e−t /2 .
The parameters of the system are chosen as follows:
• • • •
L1 = 5, L2 = 4 — sizes of the rectangular domain Ω , ε ∈ {10−2 , 10−1 , 1, 10} — values of diffusion coefficient, △t = 10−2 — time step, (u1 , u2 ) = (3, 7) — control function is constant.
For clarity we have computed the ratio log(∥eh (t )∥0,Ω /h), which is bounded above by C +λt (according to Theorem 8). For comparison, Figs. 3–6 show also the growth of log(∥eh (t )∥1,Ω /h) and log(∥eh (t )∥∞,Ω /h) with respect to t. These numerical tests were run only for the mesh size h ≈ 10−3 since for smaller h > 0 the behavior of log(∥eh (t )∥/h) changes only slightly. We enhance the plots with the estimated theoretical maximum growth rate λ > 0 for the values of log(∥eh (t )∥0,Ω /h). Note that the error estimate provided in Theorem 8 covers the test problems. Estimates for ∥eh (t )∥1,Ω or ∥eh (t )∥∞,Ω , but under stronger assumptions, have also been provided in [29, Proposition 11.1.1, Remark 11.2.4]. However, it can observed on the Figs. 3–6 that the error in H 1 - and L∞ -norms behaves similarly in the particular case of exact solutions cˆ (1) (x, t ), cˆ (2) (x, t ). In addition, the growth rate of log(∥eh (t )∥/h) at any time moment does not exceed our theoretical upper bound λ, thus the behavior of the errors complies with the theory. Finally, the tests show that the errors increase for smaller diffusion coefficients ε , which means that in the convection-dominated regime one should choose a sufficiently small mesh size h.
286
N.R. Bayramov, J.K. Kraus / Journal of Computational and Applied Mathematics 280 (2015) 275–293
a
b
Fig. 3. Evolution of the ratios log(∥eh (t )∥0,Ω /h), log(∥eh (t )∥1,Ω /h) and log(∥eh (t )∥∞,Ω /h) for ε = 10−2 and an exact solution (a) cˆ (1) (x, t ), (b) cˆ (2) (x, t ).
a
b
Fig. 4. Evolution of the ratios log(∥eh (t )∥0,Ω /h), log(∥eh (t )∥1,Ω /h) and log(∥eh (t )∥∞,Ω /h) for ε = 10−1 and an exact solution (a) cˆ (1) (x, t ), (b) cˆ (2) (x, t ).
a
b
Fig. 5. Evolution of the ratios log(∥eh (t )∥0,Ω /h), log(∥eh (t )∥1,Ω /h) and log(∥eh (t )∥∞,Ω /h) for ε = 1 and an exact solution (a) cˆ (1) (x, t ), (b) cˆ (2) (x, t ).
4.2. Comparison of the two schemes The SUPG and the EAFE approximations were compared for the transient convection–diffusion problem that is driven by the velocity field V⃗ (x) from (26). The method and problem parameters were chosen as follows:
N.R. Bayramov, J.K. Kraus / Journal of Computational and Applied Mathematics 280 (2015) 275–293
a
287
b
Fig. 6. Evolution of the ratios log(∥eh (t )∥0,Ω /h), log(∥eh (t )∥1,Ω /h) and log(∥eh (t )∥∞,Ω /h) for ε = 10 and an exact solution (a) cˆ (1) (x, t ), (b) cˆ (2) (x, t ).
Fig. 7. c (x, t ) for t = 0.05, u = (2.1, 8.5) obtained by (a) EAFE, (b) SUPG method.
• L1 = 5, L2 = 0.05 — sizes of the rectangular domain Ω , • ε = 4 · 10−6 — diffusion coefficient, • △t = 5 · 10−3 — time step. Since the considered problem arises from technical applications, one of which is described and studied in [25], the control parameter u is selected such that after a certain period of time (not larger then T ) the concentration is accumulated mostly within a subdomain ω. We were interested only in the part of the mesh Th that corresponds to a truncated domain
Ωθ = [0, θ · L1 ] × [0, L2 ],
θ = 0.2,
which covers the subdomain ω and can be regarded as the support of c (·, t ) at any time moment t ∈ [0, T ]. The concentration c (x, t ) was evaluated only on the truncated domain as it is physically negligible on the rest of the domain. A uniform triangulation Th of the domain Ωθ into triangles was used with nx = 1000 and ny = 500 nodes along the x-axis and y-axis, respectively. For the SUPG method the finite element space Vh of piecewise linear functions was used (higher order approximations, such as for p = 2, have not resulted in the disappearance of small oscillations). In the subsequent plots of the solutions (illustrated by isolines) it is clearly seen that the EAFE and SUPG schemes produce very similar stable solutions ch (tk ) before a boundary layer is formed (driven by the velocity field V⃗ for a properly selected control parameter u(t )). In this situation spurious oscillations of the discrete solution at the boundary layer can be observed for the SUPG method. On the other hand, the monotone EAFE scheme is free from these defects. This additional stability and also the ease of implementation, in particular, the fact that the EAFE scheme does not require any stabilization parameters, certainly makes this method a good candidate for solving transient convection–diffusion problems, as they were considered in the present work. See Figs. 7–11.
288
N.R. Bayramov, J.K. Kraus / Journal of Computational and Applied Mathematics 280 (2015) 275–293
Fig. 8. c (x, t ) for t = 0.1, u = (2.1, 8.5) obtained by (a) EAFE, (b) SUPG method.
Fig. 9. c (x, t ) for t = 0.4, u = (2.1, 8.5) obtained by (a) EAFE, (b) SUPG method.
Fig. 10. c (x, t ) for t = 0.8, u = (1.6, 8.5) with switching moment T1 = 0.5, obtained by (a) EAFE, (b) SUPG method.
5. Conclusions We have compared two stable finite element schemes (SUPG and EAFE) for solving transient convection-dominated convection–diffusion problems. In general, both schemes provide similar results (for isotropic diffusion). However, for certain choices of the control parameters the SUPG solution exhibits spurious oscillations at boundary layers. This does not happen for the monotone EAFE scheme (namely, for its time-dependent generalization). This explains the choice of the numerical method made in [25] for the numerical solution of optimal control problems arising in the asymmetric flow
N.R. Bayramov, J.K. Kraus / Journal of Computational and Applied Mathematics 280 (2015) 275–293
a
289
b
Fig. 11. c (x, t ) for t = 1.0, u = (1.6, 8.5) with switching moment T1 = 0.5, obtained by (a) EAFE, (b) SUPG method.
field-flow fractionation (AF4) process. To justify the applicability of the method, a new error estimate has been derived for the EAFE scheme in the transient case. The design and analysis of efficient iterative methods for solving the arising systems of linear algebraic equations is subject of ongoing and future research. Acknowledgments The authors are grateful to Prof. Oleg Iliev, who initiated the present work motivated by an industrial application and to Prof. Ludmil Zikatanov for his valuable comments on the EAFE scheme. Thanks also go to Tigran Nagapetyan for providing the analytical expression for the solution of the Stokes system (25). This work has been supported by the Austrian Science Fund (FWF), Grant P22989-N18. Appendix In this appendix, we summarize some results about the integration and differentiation of functions f : (0, T ) → X , where X is a separable Banach space with norm ∥ · ∥, cf. [31,32]. Definition 9. A simple function f : (0, T ) → X is a function of the form f =
N
cj χEj
j =1
where E1 , . . . , EN are Lebesgue measurable subsets of (0, T ), c1 , . . . , cN ∈ X and
χE (t ) =
1, 0,
if t ∈ E if t ̸∈ E
denotes the characteristic function of E. Definition 10. A function f : (0, T ) → X is (strongly) measurable, if there is a sequence {fn } of simple functions such that fn (t ) → f (t ) strongly in X (i.e. in norm) for t a.e. in (0, T ). Definition 11. A strongly measurable function f : (0, T ) → X is (Bochner) integrable, if there is a sequence of simple functions such that fn (t ) → f (t ) pointwise a.e. in (0, T ) and T
∥f − fn ∥ dt = 0.
lim
n→∞
0
The integral of f is defined by T
T
fn dt ,
f dt = lim 0
n→∞
0
where the limit exists strongly in X .
290
N.R. Bayramov, J.K. Kraus / Journal of Computational and Applied Mathematics 280 (2015) 275–293
The value of the Bochner integral of f is independent of the sequence {fn } of approximating simple functions (which can be shown analogously to the case of real-valued functions), and
0
T
T
f dt ≤
∥f ∥ dt . 0
Definition 12. For 1 ≤ p < ∞ the space Lp (0, T ; X ) consists of all strongly measurable functions f : (0, T ) → X such that T
∥f ∥p dt < ∞ 0
and it is equipped with the norm T
∥f ∥Lp (0,T ;X ) =
∥f ∥p dt
1/p
.
0
The following result (see [32]) characterizes (Bochner) integrable functions as elements of L1 (0, T ; X ). Theorem 13 (Bochner, 1933). A function f : (0, T ) → X is Bochner integrable if and only if it is strongly measurable and T
∥f ∥ dt < ∞. 0
To prove the main theorem of this appendix, we need to introduce several auxiliary lemmata. Lemma 14. Suppose that fn : (0, T ) → X is Bochner integrable for each n ∈ N, fn (t ) → f (s) as n → ∞ strongly in X for t a.e. in (0, T ), and there is an integrable function g : (0, T ) → R such that either
∥fn (t )∥ ≤ g (t ) for t a.e. in (0, T ) and every n ∈ N, or the sequence {∥fn (t )∥}n∈N is uniformly integrable. Then f : (0, T ) → X is Bochner integrable and T
T
0
T
f dt ,
fn dt →
∥fn − f ∥ dt → 0 as n → ∞.
0
0
Proof. The proof is the same as for real-valued functions.
Let L1loc (0, T ; X ) denote the space of measurable functions f : (0, T ) → X that are integrable on every compactly supported interval (a, b), i.e., 0 < a < b < T . Also, let Cc∞ (0, T ) denote the space of smooth, real-valued functions φ : (0, T ) → R with compact support in (0, T ). Definition 15. A function f ∈ L1loc (0, T ; X ) is weakly differentiable with weak derivative ft = g ∈ L1loc (0, T ; X ) if T
φ f dt = − ′
0
T
φ g dt for every φ ∈ Cc∞ (0, T ). 0
Lemma 16. Suppose that f ∈ L1 (0, T ; X ), then f (t ) = lim
h→0
1 h
t +h
f (s) ds t
for t pointwise a.e. in (0, T ). Proof. Since X is a separable space, let {cn ∈ X : n ∈ N} be a dense subset of X . Then by the Lebesgue differentiation theorem for real-valued functions
∥f (t ) − cn ∥ = lim
h→0
1 h
t +h
∥f (s) − cn ∥ ds t
for every n ∈ N and t pointwise a.e. in (0, T ). Thus, for all such t ∈ (0, T ) and every n ∈ N, we have lim sup h→0
1 h
t +h
∥f (s) − f (t )∥ ds ≤ lim sup t
h→0
1 h
t +h
∥f (s) − cn ∥ + ∥f (t ) − cn ∥ ds = 2∥f (t ) − cn ∥. t
N.R. Bayramov, J.K. Kraus / Journal of Computational and Applied Mathematics 280 (2015) 275–293
291
Since this holds for every element cn of the dense subset, it follows that 1
lim sup
h
h→0
t +h
∥f (s) − f (t )∥ ds = 0. t
Therefore
1 lim sup h
h→0
t +h
1 f (s) ds − f (t ) ≤ lim sup h→0
t
which proves the result.
h
t +h
∥f (s) − f (t )∥ ds = 0, t
Lemma 17. Suppose that f : (0, T ) → X is locally integrable and T
φ f dt = 0 for every φ ∈ Cc∞ (0, T ). 0
Then f = 0 pointwise a.e. on (0, T ). Proof. Choose a sequence of test functions 0 ≤ φn ≤ 1 whose supports are contained inside a fixed compact subset of (0, T ) such that φn → χ(t ,t +h) pointwise, where χ(t ,t +h) is the characteristic function of the interval (t , t + h) ⊂ (0, T ). If f ∈ L1loc (0, T ; X ), then we can apply the analog of the dominated convergence theorem, Lemma 14, to obtain t +h
f (s) ds = lim
n→∞
t
But
T
φn (s)f (s) ds. 0
T 0
φ f ds = 0 for every φ ∈ Cc∞ (0, T ), therefore t +h f (s) ds = 0 t
for every (t , t + h) ⊂ (0, T ). It then follows from Lemma 16 that f = 0 pointwise a.e. in (0, T ).
Lemma 18. Suppose that f : (0, T ) → X is weakly differentiable and f ′ = 0. Then f is equivalent to a constant function. Proof. The condition that the weak derivative f ′ is zero means that (see Definition 15) T
f φ ′ dt = 0
for all φ ∈ Cc∞ (0, T ).
(27)
0
Choose a fixed test function η ∈ Cc∞ (0, T ) whose integral is equal to one, and represent an arbitrary test function φ ∈ Cc∞ (0, T ) as
φ = K η + ψ ′, where K ∈ R and ψ ∈ Cc∞ (0, T ) are given by T
φ dt ,
K =
ψ(t ) =
0
t
φ(s) − K η(s) ds.
0
If we define T
ηf dt ∈ X ,
c= 0
then (27) implies that T
(f − c )φ dt = 0 for all φ ∈ Cc∞ (0, T ), 0
and Lemma 17 implies that f = c pointwise a.e. on (0, T ).
Theorem 19. Suppose that f ∈ L1 (0, T ; X ). Then f is weakly differentiable with integrable derivative ft = g ∈ L1 (0, T ; X ) if and only if f (t ) = c0 +
t
g (s) ds
(28)
0
pointwise a.e. in (0, T ). In that case, f is differentiable pointwise a.e. and its pointwise (strong) derivative coincides with its weak derivative.
292
N.R. Bayramov, J.K. Kraus / Journal of Computational and Applied Mathematics 280 (2015) 275–293
Proof. If f is given by (28), then f (t + h) − f (t ) h
=
1
t +h
h
g (s) ds, t
and Lemma 16 implies that the strong derivative of f exists pointwise a.e. and is equal to g. We also have that
f ( t + h) − f ( t ) 1 ≤ h
h
t +h
∥g (s)∥ ds.
(29)
t
Extending f by zero to a function f : R → X , and using Fubini’s theorem, we get
1 f (t + h) − f (t ) dt ≤ h
R
h
≤
R
t +h t h
1 h
R
∥g (s)∥ ds dt
∥g (s + t )∥ ds dt ≤ 0
1
h
h
0
∥g (s + t )∥ dt ds =
∥g (t )∥ dt .
(30)
R
R
Choose φ ∈ Cc∞ (0, T ). Due to the analog of the dominated convergence theorem, Lemma 14, T
φ ′ (t )f (t ) dt = lim
T φ(t + h) − φ(t )
h→0
0
h
0
f (t ) dt .
Further we can reorder the terms in the integration (since φ is zero on [0, h1 ] ∪ [T − h2 , T ]): lim
T φ(t + h) − φ(t )
h→0
h
0
f (t ) dt = − lim
T
h→0
f (t ) − f (t − h) φ(t ) dt . h
0
The estimates (29) and (30) imply that the class of functions
f (t +h)−f (t ) h
apply again Lemma 14:
f (t ) − f (t − h) dt = − φ(t )
T
− lim
h→0
h
0
is uniformly integrable, which allows us to |h|
T
φ(t )g (t ) dt . 0
This finally means that g is the weak derivative of f : T
φ ′ (t )f (t ) dt = − 0
T
φ(t )g (t ) dt . 0
Conversely, if ft = g ∈ L1 (0, T ; X ) in the sense of weak derivatives, let
f (t ) =
t
g (s) ds. 0
ft = g, so the weak derivative (f − f )t is zero. Lemma 18 then implies that f − f Then the previous argument implies that is constant pointwise a.e., which gives (28). References [1] S. Flotron, J. Rappaz, Conservation schemes for convection–diffusion equations with Robins boundary conditions. MATHICSE Technical Report Nr. 11.2012, École Polytechnique Fédérale de Lausanne, 2012, March. [2] H.-G. Roos, M. Stynes, L. Tobiska, Robust Numerical Methods for Singularly Perturbed Differential Equations, second ed., in: Springer Series in Computational Mathematics, vol. 24, Springer-Verlag, Berlin, 2008. [3] C. Johnson, Numerical Solution of Partial Differential Equations by the Finite Element Method, Dover Publications, Inc., Mineola, NY, 2009, Reprint of the 1987 edition. [4] A.A. Samarskii, Theory of Difference Schemes, Nauka, Moscow, 1977. [5] I. Christie, D.F. Griffiths, A.R. Mitchell, O.C. Zienkiewicz, Finite element methods for second order differential equations with significant first derivatives, Internat. J. Numer. Methods Engrg. 10 (1976) 1389–1396. [6] T.J.R. Hughes, A.N. Brooks, A multidimensional upwind scheme with no crosswind diffusion. Finite element methods for convection dominated flows, in: AMD, 34, Amer. Soc. Mech. Engrs. (ASME), New York, 1979, pp. 19–35. [7] A.N. Brooks, T.J.R. Hughes, Streamline upwind/Petrov–Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier–Stokes equations, Comput. Methods Appl. Mech. Engrg. 32 (1–3) (1982) 199–259. [8] V. John, J.M. Maubach, L. Tobiska, Nonconforming streamline-diffusion-finite-element-methods for convection–diffusion problems, Numer. Math. 78 (2) (1997) 165–188. [9] V. John, J. Novo, Error analysis of the SUPG finite element discretization of evolutionary convection–diffusion–reaction equations, SIAM J. Numer. Anal. 49 (3) (2011) 1149–1176. [10] B. Semper, Numerical crosswind smear in the streamline diffusion method, Comput. Methods Appl. Mech. Engrg. 113 (1–2) (1994) 99–108. [11] V. John, P. Knobloch, On the performance of SOLD methods for convection–diffusion problems with interior layers, Int. J. Comput. Sci. Math. 1 (2–4) (2007) 245–258.
N.R. Bayramov, J.K. Kraus / Journal of Computational and Applied Mathematics 280 (2015) 275–293
293
[12] V. John, P. Knobloch, On spurious oscillations at layers diminishing (SOLD) methods for convection–diffusion equations: part II—analysis for P1 and Q1 finite elements, Comput. Methods Appl. Mech. Engrg. 197 (2008) 1997–2014. [13] V. John, P. Knobloch, On spurious oscillations at layers diminishing (SOLD) methods for convection–diffusion equations: part I—a review, Comput. Methods Appl. Mech. Engrg. 196 (2007) 2197–2215. [14] B. Ayuso, D.L. Marini, Discontinuous Galerkin methods for advection-diffusion–reaction problems, SIAM J. Numer. Anal. 47 (2) (2009) 1391–1420. [15] J. Proft, B.M. Rivière, Discontinuous Galerkin methods for convection–diffusion equations for varying and vanishing diffusivity, Int. J. Numer. Anal. Model. 6 (4) (2009) 533–561. [16] M. Braack, E. Burman, Local projection stabilization for the Oseen problem and its interpretation as a variational multiscale method, SIAM J. Numer. Anal. 43 (2006) 2544–2566. [17] G. Matthies, P. Skrzypacz, L. Tobiska, A unified convergence analysis for local projection stabilisations applied to the Oseen problem, M2AN Math. Model. Numer. Anal. 41 (2007) 713–742. [18] R. Löhner, K. Morgan, J. Peraire, M. Vahdati, Finite element flux-corrected transport (FEM-FCT) for the Euler and Navier–Stokes equations, Internat. J. Numer. Methods Fluids 7 (1987) 1093–1109. [19] D. Kuzmin, M. Möller, S. Turek, High-resolution FEM-FCT schemes for multidimensional conservation laws, Comput. Methods Appl. Mech. Engrg. 193 (2004) 4915–4946. [20] V. John, E. Schmeyer, Finite element methods for time-dependent convection–diffusion–reaction equations with small diffusion, Comput. Methods Appl. Mech. Engrg. 198 (2008) 475–494. [21] M. Augustin, A. Caiazzo, A. Fiebach, J. Fuhrmann, V. John, A. Linke, R. Umla, An assessment of discretizations for convection-dominated convection–diffusion equations, Comput. Methods Appl. Mech. Engrg. 200 (2011) 3395–3409. [22] F. Brezzi, L.D. Marini, P. Pietra, Numerical simulation of semiconductor devices, in: Proceedings of the Eighth International Conference on Computing Methods in Applied Sciences and Engineering (Versailles, 1987), vol. 75, 1989, pp. 493–514. [23] F. Brezzi, L.D. Marini, P. Pietra, Two-dimensional exponential fitting and applications to drift-diffusion models, SIAM J. Numer. Anal. 26 (1989) 1342–1355. [24] M.E. Schimpf, K.D. Caldwell, J.C. Giddings, The FFF Handbook, Wiley, New York, 2000. [25] N. Bayramov, T. Nagapetyan, R. Pinnau, Fast optimal control of asymmetric flow field flow fractionation processes, in: Proceedings of SIAM Conference on Control and Its Applications, San Diego, California (2013), pp. 207–213. [26] J. Xu, L.T. Zikatanov, A monotone finite element scheme for convection–diffusion equations, Math. Comp. 68 (228) (1999) 1429–1446. [27] R.D. Lazarov, L.T. Zikatanov, An exponential fitting scheme for general convection–diffusion equations on tetrahedral meshes, Comput. Appl. Math. 1 (92) (2005) 60–69. [28] A. Russo, G. Starita, On the existence of steady-state solutions to the Navier–Stokes system for large fluxes, Ann. Sc. Norm. Super. Pisa Cl. Sci. 7 (2008) 171–180. [29] A. Quarteroni, A. Valli, Numerical Approximation of Partial Differential Equations, in: Springer Series in Computational Mathematics, vol. 23, SpringerVerlag, 1994. [30] S.C. Brenner, L.R. Scott, The Mathematical Theory of Finite Element Methods, third ed., Springer, 2008. [31] S. Bochner, Integration von Funktionen deren Werte die Elemente eines Vektorraumes sind, Fund. Math. 20 (1933) 262–276. [32] K. Yosida, Theory of Difference Schemes, Springer-Verlag, Berlin, 1980.