The paper studies finite element methods for the simulation of time–dependent convection–diffusion–reaction equations with small diffusion: the SUPG method, a SOLD method and two types of FEM–FCT methods. The methods are assessed, in particular with respect to the size of the spurious oscillations in the computed solutions, at a 3D example with nonhomogeneous Dirichlet boundary conditions and homogeneous Neumann boundary conditions

John, Volker

Schmeyer, Ellen

Numerical studies of stabilized finite element methods for solving scalar time--dependent convection--diffusion--reaction equations with small diffusion are presented in this paper. These studies include the Streamline--Upwind Petrov--Galerkin (SUPG) method with different parameters, various Spurious Oscillations at Layers Diminishing (SOLD) methods, a Local Projection Stabilization (LPS) scheme based on enrichment and two Finite Element Method Flux--Corrected--Transport (FEM--FCT) methods. The focus of the evaluation of the numerical results is on the reduction of spurious oscillations

John, V.

Schmeyer, E.

Repository: Freie Universität Berlin (FU), Math Department (fu_mi_publications)

Finite element methods for time-dependent convection-diffusion-reaction equations with small diffusion

Universaar

On finite element methods for 3D time–dependent convection–diffusion–reaction equations with small diffusion

Acronym

Universita¨t des SaarlandesUN IV E R S IT A SSAR A V I E NS ISFachrichtung 6.1 – MathematikPreprint Nr. 219On finite element methods for 3Dtime–dependentconvection–diffusion–reaction equations withsmall diffusionVolker John and Ellen SchmeyerSaarbru¨cken 2008Fachrichtung 6.1 – Mathematik Preprint No. 219Universita¨t des Saarlandes submitted: September 8, 2008On finite element methods for 3Dtime–dependentconvection–diffusion–reaction equations withsmall diffusionVolker JohnSaarland UniversityDepartment of MathematicsP.O. Box 15 11 5066041 Saarbru¨ckenGermanyjohn@math.uni-sb.deEllen SchmeyerSaarland UniversityDepartment of MathematicsP.O. Box 15 11 5066041 Saarbru¨ckenGermanyschmeyer@math.uni-sb.deEdited byFR 6.1 – MathematikUniversita¨t des SaarlandesPostfach 15 11 5066041 Saarbru¨ckenGermanyFax: + 49 681 302 4443e-Mail: preprint@math.uni-sb.deWWW: http://www.math.uni-sb.de/AbstractThe paper studies finite element methods for the simulation oftime–dependent convection–diffusion–reaction equations with smalldiffusion: the SUPG method, a SOLD method and two types of FEM–FCT methods. The methods are assessed, in particular with respectto the size of the spurious oscillations in the computed solutions, at a3D example with nonhomogeneous Dirichlet boundary conditions andhomogeneous Neumann boundary conditions.1 IntroductionThe simulation of various applications requires the numerical solution oftime–dependent convection–diffusion–reaction equations. A typical exam-ple are processes which involve a chemical reaction in a flow field [5]. Sucha reaction can be modeled with a coupled system of time–dependent non-linear convection–diffusion–reaction equations for the concentrations of thereactants and the products.Typically, the solution of these equations possesses layers. A numericalmethod for the simulation of these equations, whose results can be consideredto be useful, should meet two requirements:• sharp layers (with respect to the used mesh size) should be computed,• spurious oscillations in the solution must not occur.The second requirement means in particular that the computed solution doesnot have negative values if, for instance, the behavior of concentrations issimulated. A number of finite element methods have been developed for thesimulation of convection–diffusion–reaction equations with small diffusion.One of the most popular ones is the streamline upwind Petrov–Galerkin(SUPG) method from [2, 1]. This method leads to solutions with sharplayers, however also with sometimes considerable spurious oscillations. Toreduce these oscillations, a number of so–called Spurious Oscillations at Lay-ers Diminishing (SOLD) schemes have been proposed, see the reviews [3, 4].SOLD schemes add additional, in general nonlinear, stabilization terms tothe SUPG method. A completely different finite element approach of treatingequations with small diffusion are Finite Element Method Flux–Corrected–Transport (FEM–FCT) schemes [10, 8]. These methods do not modify thebilinear form but they manipulate the matrix and the right hand side of aGalerkin finite element method.A first comparison of finite element methods for time–dependent convection–diffusion–reaction equations was presented in [6]. The numerical examples of1[6] studied problems in 2D with homogeneous Dirichlet boundary conditions.The present paper extends the studies of [6] to 3D problems with inhomo-geneous Dirichlet and Neumann boundary conditions. This is a realisticsituation in applications.2 Finite element methods time–dependent convection–diffusion–reaction equationsWe consider a linear time–dependent convection–diffusion–reaction equationut − ε∆u+ b · ∇u+ cu = 0 in (0, T ]× Ω, (1)where ε > 0 is the diffusion coefficient, b ∈ L∞(0, T ; (W 1,∞(Ω))) is theconvection field, c ∈ L∞(0, T ;L∞(Ω)) is the non–negative reaction coefficient,T > 0 is the final time and Ω ⊂ R3 is a bounded domain. This equation hasto be equipped with an initial condition u0 = u(0,x) and with appropriateboundary conditions.In the numerical studies, (1) will be discretized in time with the Crank-Nicolson scheme using equidistant time steps ∆t. This leads in the discretetime tk to the equationuk + 0.5∆t (−ε∆uk + b · ∇uk + cuk) (2)= uk−1 − 0.5∆t (−ε∆uk−1 + b · ∇uk−1 + cuk−1) + 0.5∆tfk−1 + 0.5∆tfk.Equation (2) can be considered as a steady–state convection–diffusion–reactionequation, with the diffusion, convection and reaction, respectively, given byD = 0.5∆tε, C = 0.5∆tb, R = 1 + 0.5∆tc.The Galerkin finite element method reads as follows: Find uhk ∈ Vhanssuchthat(uhk, vh) + 0.5∆t((ε∇uhk,∇vh) + (b · ∇uhk + cuhk, vh))= (uhk−1, vh)− 0.5∆t((ε∇uhk−1,∇vh) + (b · ∇uhk−1 + cuhk−1, vh))(3)+0.5∆t(fk−1, vh) + 0.5∆t(fk, vh)for all V htest, where V hansand V htestare appropriate finite element spaces.The SUPG method adds a consistent diffusion term in streamline directionto the left hand side of (3)∑K∈T hτK(Rh(uhk),C · ∇vh)K,2where {τK} is a set of parameters depending on the mesh cells {K}. Theresidual Rh(uhk) is defined by the difference of the left hand side and the righthand side of (2). Different proposals for the choice of the parameters {τK}can be found in the literature. In the numerical studies of [6], the choicefrom [7]τK = min{hK∆t‖b‖2,11 + 0.5∆tc,2h2K∆tε},where ‖ · ‖2 denotes the Euclidean norm of a vector, has been proven tobe the best one. However, it is well known that numerical solutions whichare computed with the SUPG method often possess non–negligible spuriousoscillations at the layers.SOLD methods are trying to reduce the spurious oscillations of SUPG meth-ods by adding another stabilization term to the SUPG method. This sta-bilization term is in general nonlinear. There are several classes of SOLDmethods, see [3, 4]. It was found in the numerical studies of [6] that the bestresults among the SOLD methods were obtained with a method adding ananisotropic diffusion term(ε˜Cos∇uhk,∇vh) with Cos =I −C⊗C‖C‖22if C 6= 0,0 else,and the parameterε˜|K = max{0, Cdiam(K)|Rh(uhk)|2‖∇uhk‖2−D}, (4)where diam(K) is the diameter of a mesh cell K. This type of parameter wasproposed in [7] and modified to the form (4) in [3]. The SOLD parameter(4) contains a free parameter C which has to be chosen by the user.The last approach which will be studied in our numerical tests are FEM–FCT schemes. They start with the algebraic equation which corresponds tothe Galerkin finite element method (3)(MC + 0.5∆tA)uk = (MC − 0.5∆tA)uk−1 + 0.5∆tfk−1 + 0.5∆tfk, (5)where (MC)ij = (mij) = (ϕj, ϕi) is the consistent mass matrix. The firstgoal of FEM–FCT schemes consists in manipulating (5) such that a stablebut low order scheme is represented. To this end, define L = A+D withD = (dij), dij = −max{0, aij, aji} = for i 6= j, dii = −N∑j=1,j 6=idij ,3and ML = diag(mi) with mi =∑Nj=1mij , where N is the number of degreesof freedom. ML is called lumped mass matrix. The low order scheme reads(ML + 0.5∆tL)uk = (ML − 0.5∆tL)uk−1 + 0.5∆tfk−1 + 0.5∆tfk. (6)The second goal of FEM–FCT schemes consists in the modification of theright hand side of (6) such that diffusion is removed in (6) where it is notneeded but spurious oscillations are still suppressed(ML + 0.5∆tL)uk = (ML − 0.5∆tL)uk−1 + 0.5∆tfk−1 + 0.5∆tfk+f ∗(uk, uk−1). (7)The computation of the anti–diffusive flux vector f ∗(uk, uk−1) is somewhatinvolved and we refer to [8, 9, 10, 6] for details. Its computation relies on apredictor step which uses an explicit and stable low order scheme. Thus, astability issue arises in FEM–FCT schemes which leads to the CFL–like con-dition ∆t < 2mini mi/lii. This condition was fulfilled in the numerical testspresented in Section 3. We will consider a nonlinear approach for computingf ∗(uk, uk−1) [10, 9] and a linear approach [8].3 Numerical studiesWe consider a situation which has some typical features of a chemical re-action in application. First, the domain is three dimensional, Ω = (0, 1)3.There is an inlet at (0, 5/8, 5/8)×(0, 6/8, 6/8) and an outlet at (1, 3/8, 4/8)×(1, 4/8, 5/8). The convection is given by b = (1,−1/4,−1/8)T , which corre-sponds to the vector pointing from the center of the inlet to the center of theoutlet. Thus, the convection will not be aligned to the mesh. The diffusionis given by ε = 10−6 and the reaction byc(x) ={1 if ‖x− g‖2 ≤ 0.10 else,where g is the line through the center of the inlet and the center of the outlet.That means, a reaction takes place only where the solution (concentration)is expected to be transported. The inflow boundary condition isuin =sin(pit/2) if t ∈ [0, 1]1 if t ∈ (1, 2]sin(pi(t− 1)/2) if t ∈ (2, 3].At the outflow, homogeneous Neumann boundary conditions are prescribed.The right hand side was set to be f = 0 in Ω for all times and the final time4in our numerical studies was T = 3. The orders of magnitudes for diffusion,convection, reaction and concentration correspond to the situation of [5].Results will be presented for the P1 finite element on a tetrahedral meshand the Q1 finite element on a hexahedral mesh. The number of degrees offreedom on both meshes is 35 937, including Dirichlet nodes. The diameterof the mesh cells is about 0.054 for the hexahedral mesh and between 0.054and 0.076 for the tetrahedral mesh. The Crank–Nicolson scheme was appliedwith ∆t = 0.001.From the construction of the problem, it is expected that the solution istransported from the inflow to the outflow with a little smearing due to thediffusion. It should take values in [0, 1]. The size of the spurious oscillationsin the numerical schemes will be illustrated with uhmin(t), Fig. 1, the size ofthe undershoots. The undershoots are particularly dangerous in applicationssince they represent non–physical situations, like negative concentrations.Fig. 2 shows the distribution of the undershoots with uhmin(t) 6 0.01 for theSUPG method at t = 2. Cut planes of the solutions at t = 2 are given inFigs. 3 – 7. These cut planes contain the centers of the inlet and the outletand they are parallel to the z–axis. Note, some wiggles which can be seen inthe contour lines, might be due to the rather coarse meshes. For illustratingthe spurious oscillations, a color bar is given for each cut plane.0 0.5 1 1.5 2 2.5 3−0.6−0.5−0.4−0.3−0.2−0.100.1  SUPGFEM−FCT−linFEM−FCT−nonlinKLR C=0.2KLR C=0.40 0.5 1 1.5 2 2.5 3−0.6−0.5−0.4−0.3−0.2−0.100.1  SUPGFEM−FCT−linFEM−FCT−nonlinKLR C=0.2KLR C=0.4Figure 1: Minimal value of finite element solutions uhmin(t), left Q1, right P1.The numerical results show the large amount of spurious oscillations in thesolutions computed with the SUPG method. Fig. 2 demonstrates that thesolutions are globally polluted with spurious oscillations. The oscillationswere considerably reduced and localized (not shown here) with the SOLDmethod KLR02. Increasing the constant leads to a decrease of the spuriousoscillations, Fig 1. From the numerical studies of [3, 4] it is known that anincrease of the constant results to somewhat more smearing of the solutions.However, this is rather tolerable in applications than spurious oscillations.5Figure 2: Distribution of negative oscillations uhmin(t) 6 0.01 for the SUPGmethod at t = 2, left Q1, right P1.Figure 3: Cut of the solution, SUPG method at t = 2, left Q1, right P1.The solutions obtained with the FEM–FCT methods are almost free of spuri-ous oscillations. These schemes gave the best results in the numerical studies.Computing times for the methods are given in the Table 1. For solvingthe equations in the nonlinear schemes, the same fixed point iteration asdescribed in [4, 6] was used. The iterations were stopped if the Euclideannorm of the residual was less than 10−8. It can be observed that the non-linear schemes are considerably more expensive than the linear methods. InKLR02, the computing times increase with increasing size of the user–chosenparameter. All observations correspond to the results obtained in [6] for 2Dproblems.4 Summary and ConclusionsThe paper studied several finite element method for solving time–dependentconvection–diffusion–reaction equations in a 3D domain with inhomogeneous6Figure 4: Cut of the solution, SOLD method (4), C = 0.2 at t = 2, left Q1,right P1.Figure 5: Cut of the solution, SOLD method (4), C = 0.4 at t = 2, left Q1,right P1.Dirichlet and homogeneous Neumann boundary conditions. The SUPGmethodled to globally polluted solutions with large spurious oscillations. These os-cillations were reduced considerably with a SOLD method, however on theexpense of much larger computing times. FEM–FCT methods led to almostoscillation–free solutions. From the aspects of solution quality and computingtime, the linear FEM–FCT scheme seems to be, among the studied methods,the most appropriate method to be used in applications.Acknowledgment. The research of E. Schmeyer was supported by theDeutsche Forschungsgemeinschaft (DFG) by grant No. Jo 329/8–17Figure 6: Cut of the solution linear FEM–FCT method at t = 2, left Q1,right P1.Figure 7: Cut of the solution, nonlinear FEM–FCT method at t = 2, left Q1,right P1.References[1] A.N. Brooks and T.J.R. Hughes. Streamline upwind/Petrov-Galerkinformulations for convection dominated flows with particular emphasison the incompressible Navier–Stokes equations. Comput. Methods Appl.Mech. Engrg., 32:199 – 259, 1982.[2] T.J.R. Hughes and A.N. Brooks. A multidimensional upwind schemewith no crosswind diffusion. In T.J.R. Hughes, editor, Finite ElementMethods for Convection Dominated Flows, AMD vol.34, pages 19 – 35.ASME, New York, 1979.[3] V. John and P. Knobloch. A comparison of spurious oscillations at layersdiminishing (sold) methods for convection–diffusion equations: Part I –a review. Comput. Methods Appl. Mech. Engrg., 196:2197 – 2215, 2007.8method Q1 P1SUPG 5989 9473SOLD (4) C = 0.2 24832 25050SOLD (4) C = 0.4 33688 30932FEM FCT linear 5920 6509FEM FCT nonlinear 9768 10398Table 1: Computing times in seconds.[4] V. John and P. Knobloch. A comparison of spurious oscillations at layersdiminishing (sold) methods for convection–diffusion equations: Part II– analysis for P1 and Q1 finite elements. Comput. Methods Appl. Mech.Engrg., 197:1997 – 2014, 2008.[5] V. John, M. Roland, T. Mitkova, K. Sundmacher, L. Tobiska, andA. Voigt. Simulations of population balance systems with one inter-nal coordinate using finite element methods. Chem. Engrg. Sci., 2008.in press.[6] V. John and E. Schmeyer. Stabilized finite element methods for time–dependent convection–diffusion–reaction equations. Comput. MethodsAppl. Mech. Engrg., 2008. in press.[7] T. Knopp, G. Lube, and G. Rapin. Stabilized finite element methodswith shock capturing for advection–diffusion problems. Comput. Meth-ods Appl. Mech. Engrg., 191:2997 – 3013, 2002.[8] D. Kuzmin. Explicit and implicit FEM–FCT algorithms with flux lin-earization. Ergebnisberichte Angew. Math. 358, University of Dort-mund, 2008.[9] D. Kuzmin and M. Mo¨ller. Algebraic flux correction I. Scalar conserva-tion laws. In R. Lo¨hner D. Kuzmin and S. Turek, editors, Flux–CorrectedTransport: Principles, Algorithms and Applications, pages 155 – 206.Springer, 2005.[10] D. Kuzmin, M. Mo¨ller, and S. Turek. High–resolution FEM–FCTschemes for multidimensional conservation laws. Comput. Methods Appl.Mech. Engrg., 193:4915 – 4946, 2004.9

The paper studies finite element methods for the simulation of time-dependent convection-diffusion-reaction equations with small diffusion: the SUPG method, a SOLD method and two types of FEM-FCT methods. The methods are assessed, in particular with respect to the size of the spurious oscillations in the computed solutions, at a 3D example with nonhomogeneous Dirichlet boundary conditions and homogeneous Neumann boundary conditions

On finite element methods for 3D time-dependent convection-diffusion-reaction equations with small diffusion

Scientific documents from the Saarland University,

Universita¨t des Saarlandes
U
N I
V E R S IT A S
S
A
R A V I E N
S I
S
Fachrichtung 6.1 – Mathematik
Preprint Nr. 219
On finite element methods for 3D
time–dependent
convection–diffusion–reaction equations with
small diffusion
Volker John and Ellen Schmeyer
Saarbru¨cken 2008
Fachrichtung 6.1 – Mathematik Preprint No. 219
Universita¨t des Saarlandes submitted: September 8, 2008
On finite element methods for 3D
time–dependent
convection–diffusion–reaction equations with
small diffusion
Volker John
Saarland University
Department of Mathematics
P.O. Box 15 11 50
66041 Saarbru¨cken
Germany
john@math.uni-sb.de
Ellen Schmeyer
Saarland University
Department of Mathematics
P.O. Box 15 11 50
66041 Saarbru¨cken
Germany
schmeyer@math.uni-sb.de
Edited by
FR 6.1 – Mathematik
Universita¨t des Saarlandes
Postfach 15 11 50
66041 Saarbru¨cken
Germany
Fax: + 49 681 302 4443
e-Mail: preprint@math.uni-sb.de
WWW: http://www.math.uni-sb.de/
Abstract
The paper studies finite element methods for the simulation of
time–dependent convection–diffusion–reaction equations with small
diffusion: the SUPG method, a SOLD method and two types of FEM–
FCT methods. The methods are assessed, in particular with respect
to the size of the spurious oscillations in the computed solutions, at a
3D example with nonhomogeneous Dirichlet boundary conditions and
homogeneous Neumann boundary conditions.
1 Introduction
The simulation of various applications requires the numerical solution of
time–dependent convection–diffusion–reaction equations. A typical exam-
ple are processes which involve a chemical reaction in a flow field [5]. Such
a reaction can be modeled with a coupled system of time–dependent non-
linear convection–diffusion–reaction equations for the concentrations of the
reactants and the products.
Typically, the solution of these equations possesses layers. A numerical
method for the simulation of these equations, whose results can be considered
to be useful, should meet two requirements:
• sharp layers (with respect to the used mesh size) should be computed,
• spurious oscillations in the solution must not occur.
The second requirement means in particular that the computed solution does
not have negative values if, for instance, the behavior of concentrations is
simulated. A number of finite element methods have been developed for the
simulation of convection–diffusion–reaction equations with small diffusion.
One of the most popular ones is the streamline upwind Petrov–Galerkin
(SUPG) method from [2, 1]. This method leads to solutions with sharp
layers, however also with sometimes considerable spurious oscillations. To
reduce these oscillations, a number of so–called Spurious Oscillations at Lay-
ers Diminishing (SOLD) schemes have been proposed, see the reviews [3, 4].
SOLD schemes add additional, in general nonlinear, stabilization terms to
the SUPG method. A completely different finite element approach of treating
equations with small diffusion are Finite Element Method Flux–Corrected–
Transport (FEM–FCT) schemes [10, 8]. These methods do not modify the
bilinear form but they manipulate the matrix and the right hand side of a
Galerkin finite element method.
A first comparison of finite element methods for time–dependent convection–
diffusion–reaction equations was presented in [6]. The numerical examples of
1
[6] studied problems in 2D with homogeneous Dirichlet boundary conditions.
The present paper extends the studies of [6] to 3D problems with inhomo-
geneous Dirichlet and Neumann boundary conditions. This is a realistic
situation in applications.
2 Finite element methods time–dependent convection–
diffusion–reaction equations
We consider a linear time–dependent convection–diffusion–reaction equation
ut − ε∆u+ b · ∇u+ cu = 0 in (0, T ]× Ω, (1)
where ε > 0 is the diffusion coefficient, b ∈ L∞(0, T ; (W 1,∞(Ω))) is the
convection field, c ∈ L∞(0, T ;L∞(Ω)) is the non–negative reaction coefficient,
T > 0 is the final time and Ω ⊂ R3 is a bounded domain. This equation has
to be equipped with an initial condition u0 = u(0,x) and with appropriate
boundary conditions.
In the numerical studies, (1) will be discretized in time with the Crank-
Nicolson scheme using equidistant time steps ∆t. This leads in the discrete
time tk to the equation
uk + 0.5∆t (−ε∆uk + b · ∇uk + cuk) (2)
= uk−1 − 0.5∆t (−ε∆uk−1 + b · ∇uk−1 + cuk−1) + 0.5∆tfk−1 + 0.5∆tfk.
Equation (2) can be considered as a steady–state convection–diffusion–reaction
equation, with the diffusion, convection and reaction, respectively, given by
D = 0.5∆tε, C = 0.5∆tb, R = 1 + 0.5∆tc.
The Galerkin finite element method reads as follows: Find uhk ∈ V
h
ans
such
that
(uhk, v
h) + 0.5∆t
(
(ε∇uhk,∇v
h) + (b · ∇uhk + cu
h
k, v
h)
)
= (uhk−1, v
h)− 0.5∆t
(
(ε∇uhk−1,∇v
h) + (b · ∇uhk−1 + cu
h
k−1, v
h)
)
(3)
+0.5∆t(fk−1, v
h) + 0.5∆t(fk, v
h)
for all V h
test
, where V h
ans
and V h
test
are appropriate finite element spaces.
The SUPG method adds a consistent diffusion term in streamline direction
to the left hand side of (3)
∑
K∈T h
τK
(
Rh(uhk),C · ∇v
h
)
K
,
2
where {τK} is a set of parameters depending on the mesh cells {K}. The
residual Rh(uhk) is defined by the difference of the left hand side and the right
hand side of (2). Different proposals for the choice of the parameters {τK}
can be found in the literature. In the numerical studies of [6], the choice
from [7]
τK = min
{
hK
∆t‖b‖2
,
1
1 + 0.5∆tc
,
2h2K
∆tε
}
,
where ‖ · ‖2 denotes the Euclidean norm of a vector, has been proven to
be the best one. However, it is well known that numerical solutions which
are computed with the SUPG method often possess non–negligible spurious
oscillations at the layers.
SOLD methods are trying to reduce the spurious oscillations of SUPG meth-
ods by adding another stabilization term to the SUPG method. This sta-
bilization term is in general nonlinear. There are several classes of SOLD
methods, see [3, 4]. It was found in the numerical studies of [6] that the best
results among the SOLD methods were obtained with a method adding an
anisotropic diffusion term
(ε˜Cos∇u
h
k,∇v
h) with Cos =


I −
C⊗C
‖C‖2
2
if C 6= 0,
0 else,
and the parameter
ε˜|K = max
{
0, C
diam(K)|Rh(uhk)|
2‖∇uhk‖2
−D
}
, (4)
where diam(K) is the diameter of a mesh cell K. This type of parameter was
proposed in [7] and modified to the form (4) in [3]. The SOLD parameter
(4) contains a free parameter C which has to be chosen by the user.
The last approach which will be studied in our numerical tests are FEM–
FCT schemes. They start with the algebraic equation which corresponds to
the Galerkin finite element method (3)
(MC + 0.5∆tA)uk = (MC − 0.5∆tA)uk−1 + 0.5∆tfk−1 + 0.5∆tfk, (5)
where (MC)ij = (mij) = (ϕj, ϕi) is the consistent mass matrix. The first
goal of FEM–FCT schemes consists in manipulating (5) such that a stable
but low order scheme is represented. To this end, define L = A+D with
D = (dij), dij = −max{0, aij, aji} = for i 6= j, dii = −
N∑
j=1,j 6=i
dij ,
3
and ML = diag(mi) with mi =
∑N
j=1mij , where N is the number of degrees
of freedom. ML is called lumped mass matrix. The low order scheme reads
(ML + 0.5∆tL)uk = (ML − 0.5∆tL)uk−1 + 0.5∆tfk−1 + 0.5∆tfk. (6)
The second goal of FEM–FCT schemes consists in the modification of the
right hand side of (6) such that diffusion is removed in (6) where it is not
needed but spurious oscillations are still suppressed
(ML + 0.5∆tL)uk = (ML − 0.5∆tL)uk−1 + 0.5∆tfk−1 + 0.5∆tfk
+f ∗(uk, uk−1). (7)
The computation of the anti–diffusive flux vector f ∗(uk, uk−1) is somewhat
involved and we refer to [8, 9, 10, 6] for details. Its computation relies on a
predictor step which uses an explicit and stable low order scheme. Thus, a
stability issue arises in FEM–FCT schemes which leads to the CFL–like con-
dition ∆t < 2mini mi/lii. This condition was fulfilled in the numerical tests
presented in Section 3. We will consider a nonlinear approach for computing
f ∗(uk, uk−1) [10, 9] and a linear approach [8].
3 Numerical studies
We consider a situation which has some typical features of a chemical re-
action in application. First, the domain is three dimensional, Ω = (0, 1)3.
There is an inlet at (0, 5/8, 5/8)×(0, 6/8, 6/8) and an outlet at (1, 3/8, 4/8)×
(1, 4/8, 5/8). The convection is given by b = (1,−1/4,−1/8)T , which corre-
sponds to the vector pointing from the center of the inlet to the center of the
outlet. Thus, the convection will not be aligned to the mesh. The diffusion
is given by ε = 10−6 and the reaction by
c(x) =
{
1 if ‖x− g‖2 ≤ 0.1
0 else,
where g is the line through the center of the inlet and the center of the outlet.
That means, a reaction takes place only where the solution (concentration)
is expected to be transported. The inflow boundary condition is
uin =


sin(pit/2) if t ∈ [0, 1]
1 if t ∈ (1, 2]
sin(pi(t− 1)/2) if t ∈ (2, 3].
At the outflow, homogeneous Neumann boundary conditions are prescribed.
The right hand side was set to be f = 0 in Ω for all times and the final time
4
in our numerical studies was T = 3. The orders of magnitudes for diffusion,
convection, reaction and concentration correspond to the situation of [5].
Results will be presented for the P1 finite element on a tetrahedral mesh
and the Q1 finite element on a hexahedral mesh. The number of degrees of
freedom on both meshes is 35 937, including Dirichlet nodes. The diameter
of the mesh cells is about 0.054 for the hexahedral mesh and between 0.054
and 0.076 for the tetrahedral mesh. The Crank–Nicolson scheme was applied
with ∆t = 0.001.
From the construction of the problem, it is expected that the solution is
transported from the inflow to the outflow with a little smearing due to the
diffusion. It should take values in [0, 1]. The size of the spurious oscillations
in the numerical schemes will be illustrated with uh
min
(t), Fig. 1, the size of
the undershoots. The undershoots are particularly dangerous in applications
since they represent non–physical situations, like negative concentrations.
Fig. 2 shows the distribution of the undershoots with uh
min
(t) 6 0.01 for the
SUPG method at t = 2. Cut planes of the solutions at t = 2 are given in
Figs. 3 – 7. These cut planes contain the centers of the inlet and the outlet
and they are parallel to the z–axis. Note, some wiggles which can be seen in
the contour lines, might be due to the rather coarse meshes. For illustrating
the spurious oscillations, a color bar is given for each cut plane.
0 0.5 1 1.5 2 2.5 3
−0.6
−0.5
−0.4
−0.3
−0.2
−0.1
0
0.1
 
 
SUPG
FEM−FCT−lin
FEM−FCT−nonlin
KLR C=0.2
KLR C=0.4
0 0.5 1 1.5 2 2.5 3
−0.6
−0.5
−0.4
−0.3
−0.2
−0.1
0
0.1
 
 
SUPG
FEM−FCT−lin
FEM−FCT−nonlin
KLR C=0.2
KLR C=0.4
Figure 1: Minimal value of finite element solutions uh
min
(t), left Q1, right P1.
The numerical results show the large amount of spurious oscillations in the
solutions computed with the SUPG method. Fig. 2 demonstrates that the
solutions are globally polluted with spurious oscillations. The oscillations
were considerably reduced and localized (not shown here) with the SOLD
method KLR02. Increasing the constant leads to a decrease of the spurious
oscillations, Fig 1. From the numerical studies of [3, 4] it is known that an
increase of the constant results to somewhat more smearing of the solutions.
However, this is rather tolerable in applications than spurious oscillations.
5
Figure 2: Distribution of negative oscillations uh
min
(t) 6 0.01 for the SUPG
method at t = 2, left Q1, right P1.
Figure 3: Cut of the solution, SUPG method at t = 2, left Q1, right P1.
The solutions obtained with the FEM–FCT methods are almost free of spuri-
ous oscillations. These schemes gave the best results in the numerical studies.
Computing times for the methods are given in the Table 1. For solving
the equations in the nonlinear schemes, the same fixed point iteration as
described in [4, 6] was used. The iterations were stopped if the Euclidean
norm of the residual was less than 10−8. It can be observed that the non-
linear schemes are considerably more expensive than the linear methods. In
KLR02, the computing times increase with increasing size of the user–chosen
parameter. All observations correspond to the results obtained in [6] for 2D
problems.
4 Summary and Conclusions
The paper studied several finite element method for solving time–dependent
convection–diffusion–reaction equations in a 3D domain with inhomogeneous
6
Figure 4: Cut of the solution, SOLD method (4), C = 0.2 at t = 2, left Q1,
right P1.
Figure 5: Cut of the solution, SOLD method (4), C = 0.4 at t = 2, left Q1,
right P1.
Dirichlet and homogeneous Neumann boundary conditions. The SUPGmethod
led to globally polluted solutions with large spurious oscillations. These os-
cillations were reduced considerably with a SOLD method, however on the
expense of much larger computing times. FEM–FCT methods led to almost
oscillation–free solutions. From the aspects of solution quality and computing
time, the linear FEM–FCT scheme seems to be, among the studied methods,
the most appropriate method to be used in applications.
Acknowledgment. The research of E. Schmeyer was supported by the
Deutsche Forschungsgemeinschaft (DFG) by grant No. Jo 329/8–1
7
Figure 6: Cut of the solution linear FEM–FCT method at t = 2, left Q1,
right P1.
Figure 7: Cut of the solution, nonlinear FEM–FCT method at t = 2, left Q1,
right P1.
References
[1] A.N. Brooks and 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:199 – 259, 1982.
[2] T.J.R. Hughes and A.N. Brooks. A multidimensional upwind scheme
with no crosswind diffusion. In T.J.R. Hughes, editor, Finite Element
Methods for Convection Dominated Flows, AMD vol.34, pages 19 – 35.
ASME, New York, 1979.
[3] V. John and P. Knobloch. A comparison of spurious oscillations at layers
diminishing (sold) methods for convection–diffusion equations: Part I –
a review. Comput. Methods Appl. Mech. Engrg., 196:2197 – 2215, 2007.
8
method Q1 P1
SUPG 5989 9473
SOLD (4) C = 0.2 24832 25050
SOLD (4) C = 0.4 33688 30932
FEM FCT linear 5920 6509
FEM FCT nonlinear 9768 10398
Table 1: Computing times in seconds.
[4] V. John and P. Knobloch. A comparison of 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:1997 – 2014, 2008.
[5] V. John, M. Roland, T. Mitkova, K. Sundmacher, L. Tobiska, and
A. Voigt. Simulations of population balance systems with one inter-
nal coordinate using finite element methods. Chem. Engrg. Sci., 2008.
in press.
[6] V. John and E. Schmeyer. Stabilized finite element methods for time–
dependent convection–diffusion–reaction equations. Comput. Methods
Appl. Mech. Engrg., 2008. in press.
[7] T. Knopp, G. Lube, and G. Rapin. Stabilized finite element methods
with shock capturing for advection–diffusion problems. Comput. Meth-
ods Appl. Mech. Engrg., 191:2997 – 3013, 2002.
[8] D. Kuzmin. Explicit and implicit FEM–FCT algorithms with flux lin-
earization. Ergebnisberichte Angew. Math. 358, University of Dort-
mund, 2008.
[9] D. Kuzmin and M. Mo¨ller. Algebraic flux correction I. Scalar conserva-
tion laws. In R. Lo¨hner D. Kuzmin and S. Turek, editors, Flux–Corrected
Transport: Principles, Algorithms and Applications, pages 155 – 206.
Springer, 2005.
[10] D. Kuzmin, M. Mo¨ller, and S. Turek. High–resolution FEM–FCT
schemes for multidimensional conservation laws. Comput. Methods Appl.
Mech. Engrg., 193:4915 – 4946, 2004.
9


http://scidok.sulb.uni-saarland.de/volltexte/2013/4749/pdf/preprint_219_08.pdf

On finite element methods for 3D time–dependent convection–diffusion–reaction equations with small diffusion

Abstract

Similar works

Full text

Available Versions

Repository: Freie Universität Berlin (FU), Math Department (fu_mi_publications)

Universaar

Acronym

Repository: Freie Universität Berlin (FU), Math Department (fu_mi_publications)

Scientific documents from the Saarland University,