Adaptive Finite Element Method for Steady Convection-Diffusion Equation

This paper examines the numerical solution of the convection-diffusion equation in 2-D. The solution of this equation possesses singularities in the form of boundary or interior layers due to non-smooth boundary conditions. To overcome such singularities arising from these critical regions, the adaptive finite element method is employed. This scheme is based on the streamline diffusion method combined with Neumann-type posteriori estimator. The effectiveness of this approach is illustrated by different examples with several numerical experiments.


Introduction
This paper deals with the scalar convection-diffusion equation.This equation describes the transport of scalar quantity, e.g., temperature or concentration.We are interested in the convection dominated case.In this case, the solution of a convection-diffusion equation frequently has boundary or interior layers.It is well known that the standard Galerkin finite element discretization on uniform grids produces inaccurate oscillatory solutions to convection diffusion problems.Therefore several stabilized finite element methods have been developed, e.g., the streamline-upwind Petrov-Galerkin (SUPG) method [1] & [2] or streamline-diffusion finite element method (SDFEM) [3] is designed to overcome these problems by introducing a small amount of artificial diffusion in the direction of streamlines.The numerical solution obtained from the SDFEM has the desirable property that the accuracy in regions where the exact solution is smooth will not be degraded as a result of discontinuities and layers in the exact solution [4] & [5].However, the numerical solution obtained from the SDFEM can be oscillatory in regions where there are layers.One common technique to increase the accuracy of the finite element solution in these critical regions is through local grid refinement, the so-called h-method.The question is how to identify those regions and how to obtain a good balance between the refined and unrefined regions such that the overall accuracy is optimal.
Another related problem is to obtain reliable estimates of the accuracy of the computed numerical solution.A priori estimate are often insufficient and can't be used to estimate the exact error.Therefore, it is natural to acquire a posteriori error estimators to pinpoint where the error is large and, at the same time, properly bound the exact error on the whole domain.The error estimator should be local and should yield reliable upper and lower bounds for the true error in a user-specified norm.Global upper bounds are sufficient to obtain a numerical solution with accuracy below a prescribed tolerance.Local lower bounds are necessary to ensure that the grid is correctly refined so that one obtains a numerical solution with a prescribed tolerance using a nearly minimal number of grid-points.
For two-dimensional problems, several estimators have been shown to be asymptotically exact when used on uniform meshes provided the solution of the problem is smooth enough [6]- [8].Estimators based on computing residuals, so-called residual-type estimators, and estimators based on solving a local Dirichlet problem, so-called Dirichlet-type estimators, were introduced in [9].Estimators based on solving a local Neumann problem, socalled Neumann-type estimators, were first given in [10].These estimators have been studied by many researchers in [11]- [16].The Zienkiewicz-Zhu (ZZ) type of estimators based on recovery of gradient and Hessian are also well developed, see [17] & [18], and articles cited therein.
In this paper we introduce and analyze from theoretical and experimental points of view an adaptive scheme to efficiently solve the convection-diffusion equation.This scheme is based on the streamline-diffusion finite element method (SDFEM) introduced in [3] combined with an error estimator similar to the one developed in [14].We prove global upper and local lower error estimates in the energy norm, with constants which only depend on the shape-regularity of the mesh and the polynomial degree of the finite element approximating space.We perform several numerical experiments to show the effectiveness of our approach to capture boundary and inner layers sharply and without significant oscillations.
The paper is organized as follows.In Section 2 we recall the convection-diffusion problem under consideration and the Streamline Diffusion Finite Element Method.In Section 3 we define a posteriori error estimator with the energy norm of the finite element approximation error.Finally, in Section 4, we introduce the adaptive scheme and report the results of the numerical tests.

Linear Convection-Diffusion Equation
We consider the following steady linear convection-diffusion equation . We are interested in the convection dominated case and assume that (A.1) 0 L norm and the 1 H semi-norm (also called Energy Norm) are defined as Ω , respectively.We shall denote the above norm and semi-norm by the following convention if no subscript index is given then we assume an ordinary 2 L norm, 0, .Ω , and if no subscript index is given then we shall assume it is the whole of Ω .
To define weak form of Equation (2.1), we need two classes of functions: the trial functions : The standard variational formulation of Equation (2.1) is given by: Find where be a decomposition of Ω into triangles.We need to make the following geometrical assumptions on the family of triangulations h ℑ 1) Admissibility: whenever 1 K and 2 which means for any 0 h > and for any : , for triangular elements, where ( ) P K is the space of polynomials of degree not greater than 1 on K.In the case of convection-dominated problem, the standard Galerkin approximation of Equation (2.6) may produce unphysical behavior, oscillation, if the mesh is too coarse in critical regions.To circumvent these difficulties, stability of the discretization has to be increased by introducing artificial diffusion along streamlines.The Streamline-Diffusion Finite Element Method (SDFEM) [1]- [3] stabilizes a convection-dominated problem by adding weighted residuals to the standard Galerkin finite element method for hyperbolic equations which combines good stability with high order accuracy, convergence results are available (see [19]).
The SDFEM yields the following discrete problem obtained: where In Equation (2.11), a constant K δ must be chosen for every element K. Let the mesh Peclet number be de- fined by, From the analysis of the SDFEM, the following choice of K δ are optimal; see [20]: where k h is a measure of the element length in the direction of the convection flow b.For other parameter choice, see [21]- [24].

A Posteriori Error Estimator
In this topic, we introduce the analysis of a Neumann-type error estimator proposed in [14] which is an extension of the work [25].In their work, they modify the well-known Bank and Weiser estimator [10] and using the idea of Ainsworth & Oden in [26], they solve a local (element) Poisson problem over a suitably chosen (higher order) approximation space with data from interior residuals and flux jumps along element edges.
We now introduce some definitions and notations that will be needed for the error estimates.
We denote by ( ) the set of all element edges and the subsets relating to internal, Dirichlet and Neumann edges respectively as

{ }
, : We denote For the lowest order 1 P approximations over a triangular element subdivision, 2 0 K u ∇ = , so that the interior residual of element K is given by ( ) and the internal residual is approximated by ( ) where 0 K  is the to be the outward normals with respect to the edge E from element K and S respectively, we have [ ]

,
, The approximation space is denoted by consisting of edge and interior bubble functions respectively: where each member of the space is a quadratic (or biquadratic) edge bubble function E Ψ that is nonzero on edge E of element K, but non zero valued on all other edges of K. K B is the space spanned by interior cubic (or biquadratic) bubbles K φ i.e., { } : 0 1 , 0 on where each function is associated with an element K, and is zero on all edges of K, nonzero on the interior of K, and 1 The upshot is that the local problems are always well posed and that for each triangular element a 4 × 4 system of equations must be solved to compute K e .For an element h K ∈ ℑ , the local error estimate is the energy norm of K e given by where ( ) In the following, we use the short-hand notation L S .The Kay and Silvester's a posteriori error estimation can be read as following: Theorem 1.If the variational Equation (2.6) solved with a grid of linear triangular elements, and if the triangle aspect ratio condition is satisfied with β Ω , then, the estimator K η computed via Equation (3.9) satisfies the (global) upper bound property ( ) where C is independent of  and h and K h is the length of the longest edge of element K. Proof.See the details in [14].Theorem 2. If the variational Equation (2.6) with 1 b ∞ = is solved via either the Galerkin formulation or the SD formulation Equation (2.11), using a grid of linear triangular elements, and if the triangle aspect ratio condition is satisfied, then the estimator K η computed via Equation (3.9) is a local lower bound for ( ) where c is independent of  , and K ω represents the patch of four elements that have at least one boundary edge E from the set ( ) K  .Proof.See the details in [14].

Numerical Experiments
In this section we report three series of numerical experiments with the Streamline Diffusion stabilization method described in Section (2) and an h-adaptive mesh-refinement strategy based on the error estimator analyzed in Section (3).In all the experiments we have used piecewise linear finite elements (i.e., k P polynomial degree 1 k = ) and we have taken as geometric domain the unit square ( ) ( ) with different boundary conditions.We have considered varying values of the coefficients  , 1 b , and 2 b of the convection-diffusion equation.
The adaptive procedure consists in solving Equation (2.11) on a sequence of meshes up to finally attain a solution with an estimated error within a prescribed tolerance.To attain this purpose, we initiate the process with a quasi-uniform mesh and, at each step, a new mesh better adapted to the solution of Equation (2.6) must be created.This is done by computing the local error estimators K η for all K in the "old" mesh h ℑ , and refining those elements K * with { } * max : is a prescribed parameter.In all our experiments we have chosen 0.5 θ = . For other marking strategies, we refer to [27].The implementation used in this paper is derived from iFEM [28].This software package is the successor of AFEM@MATLAB [29], which contains an advanced refinement tool. .Figure 1 shows the successfully refined meshes created in the adaptive process for , as well as the corresponding computed solution.Figure 2 shows the error curves for the exact and estimated errors.  .The results show that the estimated error is well bounded as described in [22].
Example 2 (Interior layers) We consider Equation (2.1) with Figure 5 and Figure 6 clearly show that the adaptive method has successfully refined the correct elements using a greater concentration of elements in the interior layer.Figure 7 and Figure 8 show the estimated and exact error curves decrease monotonically for        We do not include error curves because no analytical solution is known in this case.
Figure 9 and Figure 10 show some of the successively refined meshes created in the adaptive process for  3), it is hard to fully resolve the internal layers and the numerical solution display a small oscillatory pattern in the internal layer.

Conclusions
An adaptive finite element scheme for the convection-diffusion equation has been introduced and analyzed.This scheme is based on the Streamline Diffusion Finite element method combined with a Neumann-type error estimator.
Several numerical experiments are reported.For in Example (3), the numerical solution displays small oscillatory pattern in the internal layer which requires a high computing cost to produce an accurate internal layer.In general, it is  quite evident that our error estimator provides an effective refinement indicator even in the presence of internal layers.
of all element vertices (that do not lie on the Dirichlet boundary D Γ ).Let E  be the set of vertices of h E ∈  , and for h K ∈ ℑ , h E ∈  and h X ∈  we define the local "patches" of elements as For any edge E of an element h K ∈ ℑ , we define the flux jump as , function on the inter-element edge E and [ ] E v measures the jump of v across E , that is, for

Example 1 (
Exponential boundary layer) The first test problem contains an exponential boundary layer.This problem corresponds to the case of

Figure 2 .
Figure 2.Estimated and exact error curves using f and boundary condition are determined from the exact solution:

Example 3 (
Interior and boundary layers) We consider Equation (2

Figure 4 .
Figure 4.Estimated and exact error curves using

Figure 7 .
Figure 7.Estimated and exact error curves using

Figure 8 .
Figure 8.Estimated and exact error curves using

3 )
Discontinuities at ( ) 0,1 causes u to have an internal layer of width ( ) the left and 1 u = to the right, as well as a boundary layer along the outflow boundary.
them show the effectiveness of this scheme to capture boundary and inner layers very sharply and without significant oscillations.But in the case of