Latest revision as of 11:15, 12 November 2018

Published in Comput. Meth. Appl. Mech. Engng., Vol. 195, pp. 3926–3946, 2006
doi: 10.1016/j.cma.2005.07.020

Abstract

A stabilized finite element method (FEM) for the steady state advection-diffusion-absorption equation is presented. The stabilized formulation is based on the modified governing differential equations derived via the Finite Calculus (FIC) method. The basis of the method is detailed for the 1D problem. It is shown that the stabilization terms act as a non linear additional diffusion governed by a single stabilization parameter. A critical constant value of this parameter ensuring a stabilized solution using two node linear elements is given. More accurate numerical results can be obtained by using a simple expression of the non linear stabilization parameter depending on the signs of the numerical solution and of its derivatives. A straightforward two steps algorithm yielding a stable and accurate solution for all the range of physical parameters and boundary conditions is described. The extension to the multi-dimensional case is briefly described. Numerical results for 1D and 2D problems are presented showing the efficiency and accuracy of the new stabilized formulation.

1 INTRODUCTION

Considerable effort has been spent in recent years to derive finite element methods (FEM) [1] for the solution of the advection-diffusion-reaction equation. The physical behaviour of this equation is varied. In the propagation regime, originated by large values of the productive term (the Helmholtz equation), solutions are exponentially modulated sinusoidal functions. Typical problems found here in the numerical solution are those of phase, amplitude and pollution errors [2]. In the exponential regime, solutions are of the form of real exponential functions where absorptive (dissipative) or productive source terms are possible. Numerical schemes here find difficulties to approximate the sharp gradients appearing in the neighborhood of boundary and internal layers due to high Peclet and/or Damköhler numbers.

Stabilized methods to tackle above problems have been based on different kinds of streamline-upwind-Petrov-Galerkin (SUPG) [2-7], Galerkin Least Square [5-11], Subgrid Scale (SGS) [5,6,12,13] and Residual Free Bubbles [14] finite element methods. While a single stabilization parameter suffices to yield stabilized (and even nodally exact results) for the one-dimensional (1D) advection-diffusion and the diffusion-reaction equations (Vol. 3 in [1] and [8-15]), this is not the case for the diffusion-advection-reaction equation. Here, in general, two stabilization parameters are needed in order to ensure a stabilized solution for all range of physical parameters and boundary conditions [4,10,14]. As reported in [12,13] both SUPG, GLS and SGS methods with a single stabilization parameter (and indeed the standard Galerkin scheme) fail to obtain an stabilized solution for some specific boundary conditions in the exponential regime with negative (absorption) terms when there is a negative streamwise gradient of the solution.

Apart from the many practical applications of the advection-diffusion-reaction problem, the interest to derive stabilized numerical schemes for this equation extends to the solution of the transient advection-diffusion equation, as the discretized form of this equation in time can be interpreted as an analogous steady-state advection-diffusion-absorption problem. This analogy has been analyzed in [16,17].

In this paper we present a reliable and accurate stabilized formulation for the advection-diffusion-reaction equation in the exponential regime with absorption (destructive) reaction terms using a single stabilization parameter.

The stabilized formulation is based on the standard Galerkin FEM solution of the modified governing differential equations derived via the Finite Calculus (FIC) method [20]. The FIC equations are obtained by expressing the balance of fluxes in a domain of finite size. This introduces additional stabilizing terms in the differential equations of the infinitessimal theory which are a function of the balance domain dimensions. These dimensions are termed characteristic length parameters and play a key role in the stabilization process.

The modified governing equations lead to stabilized numerical schemes using whatever numerical method. It is interesting that many of the stabilized FEM can be recovered using the FIC formulation. The FIC method has been successfully applied to problems of convection-diffusion [18-24], diffusion-absorption [15], incompressible fluid flow [24-30] and standard and incompressible solid mechanics [31-33] using finite element and meshless procedures.

The Galerkin FEM-FIC formulation here proposed introduces naturally an additional non linear dissipation term into the discretized equations which is governed by a single stabilization parameter. In the absence of the absorption term the formulation simplifies to the standard Petrov-Galerkin approach for the advection-diffusion problem For the diffusion-absorption case the diffusion-type stabilization term is identical to that recently obtained by Felippa and Oñate using a variational FIC approach and modified equation methods [15]. The limit value of the stabilization parameter ensuring a stable numerical solution is derived. More accurate numerical results can be obtained using an expression of the stabilization parameter which is a function of the sign of the solution and that of its first and second derivatives. This introduces a non-linearity in the solution scheme and a simple iterative algorithm leading to stabilized and accurate results in just two iterations is presented. The merit of the new formulation is that it yields stabilized and accurate numerical solutions for the advection-diffusion-absorption equation for all range of the physical parameters and boundary conditions using a single stabilization parameter. As previously mentioned, this property was only attained by previous similar stabilized formulations using two stabilization parameters. We present a collection of 1D examples showing the efficiency and accuracy of the new FEM-FIC formulation for different values of the advective and absorptive terms.

In the last part of the paper the FIC formulation is extended to multi-dimensional advection-diffusion problems. The key ingredients of the formulation are given and some examples of the good performance of the new stabilized method for the solution of 2D problems are presented.

2 FIC FORMULATION OF THE 1D STATIONARY ADVECTION-DIFFUSION-ABSORPTION EQUATION

The governing equation for the 1D stationary advection-diffusion-absorption problem can be written in the FIC formulation as

where

In above equations ${\textstyle \phi }$ is the state variable, ${\textstyle x\in [0,L]}$ is the problem domain, ${\textstyle L}$ is the domain length, ${\textstyle u}$ is the velocity field, ${\textstyle k\geq 0}$ is the diffusion, ${\textstyle s\geq 0}$ is the absorption, dissipation or destruction source parameter, ${\textstyle Q}$ is a constant source term, ${\textstyle q^{p}}$ and ${\textstyle \phi ^{p}}$ are the prescribed values of the total flux and the unknown function at the Neumann and Dirichlet boundaries ${\textstyle \Gamma _{q}}$ and ${\textstyle \Gamma _{\phi }}$, respectively and ${\textstyle h}$ is a characteristic length which plays the role of a stabilization parameter. In the 1D problem ${\textstyle \Gamma _{\phi }}$ and ${\textstyle \Gamma _{q}}$ consist of four combinations at the end points of the problem domain.

Eqs.(1) and (2) are obtained by expressing the balance of fluxes in an arbitrary 1D domain of finite size within the problem domain and at the Neumann boundary, respectively. The variations of the transported variables within the balance domain are approximated by Taylor series expansions retaining one order higher terms than in the infinitessimal theory [20]. The underlined terms in Eqs.(1) and (2) emanate from these higher order expansions and they are essential to derive stabilized numerical schemes using whatever discretization method. Note that as the characteristic length parameter ${\textstyle h}$ tends to zero the FIC differential equations gradually recover the standard infinitessimal form giving in the limit (for ${\textstyle h=0}$) ${\textstyle r=0}$ in ${\textstyle x\in (0,L)}$ and ${\textstyle -u\phi +k{d\phi \over dx}+q^{p}=0\quad {\hbox{on }}\Gamma _{q}}$.

Similarly as in all stabilized methods, the stability and accuracy of the numerical solution depends on the values of the stabilization parameter, i.e. of the characteristic length ${\textstyle h}$. At the discrete level ${\textstyle h}$ can be interpreted as a distance related to the (macroscopic) domain within which the space derivatives are computed. In the finite element practice it is usual to express ${\textstyle h}$ as a proportion of a typical element dimension (i.e. the element length for 1D problems) [18].

Successful applications of the FIC method to a variety of problems in computational mechanics can be found in [18-31].

3 FINITE ELEMENT FORMULATION

We will construct a standard finite element discretization ${\textstyle \left\{l^{e}\right\}}$ of the 1D analysis domain length ${\textstyle L}$ with index ${\textstyle e}$ ranging from 1 to the number of elements ${\textstyle N}$ [1]. The state variable ${\textstyle \phi }$ is approximated by ${\textstyle {\bar {\phi }}}$ over the analysis domain. The approximated variable ${\textstyle {\bar {\phi }}}$ is interpolated within each element with ${\textstyle n}$ nodes in the standard manner, i.e.

with

where ${\textstyle N_{i}}$ are the element shape functions and ${\textstyle \phi _{i}}$ are nodal values of the approximate function ${\textstyle {\bar {\phi }}}$. Substituting Eq.(5a) into Eqs.(1) and (2) gives

where ${\textstyle {\bar {r}}=r({\bar {\phi }})}$ and ${\textstyle r_{\Omega },r_{q}}$ and ${\textstyle r_{\phi }}$ are the residuals of the approximate solution in the problem domain and on the Neumann and Dirichlet boundaries ${\textstyle \Gamma _{q}}$ and ${\textstyle \Gamma _{\phi }}$, respectively.

The weighted residual form of Eqs.(6)–(8) is written as

where ${\textstyle W_{i}(x)}$ and ${\textstyle {\hat {W}}_{i}}$ are test functions satisfying ${\textstyle W_{i}={\hat {W}}_{i}=0}$ on ${\textstyle \Gamma _{\phi }}$. In Eq.(9a) and in the following ${\textstyle \int _{L}(\cdot )dx:=\int _{0}^{L}(\cdot )dx}$ and ${\textstyle [g]_{\Gamma _{q}}}$ denotes the value of ${\textstyle g}$ at the boundary ${\textstyle \Gamma _{q}}$. Using Eqs.(6) and (7) into Eq.(9a) we have

Assuming smooth enough solutions and integrating by parts the term involving ${\textstyle h}$ in the first integral gives for ${\textstyle {\hat {W}}_{i}=-W_{i}}$

Note that the residual term has dissappeared from the boundary integral. This is due to the consistency of the FIC governing equations in the domain and on the Neumann boundary (Eqs.(1) and (2)).

The third term in Eq.(10) is computed as the sum of the integrals over the element interiors, therefore allowing for the space derivatives of ${\textstyle {\bar {r}}}$ to be discontinuous. Also in Eq.(10) ${\textstyle h}$ has been assumed to be constant within each element, (i.e. ${\textstyle {dh \over dx}=0}$ within ${\textstyle l^{e}}$).

The weak form is obtained by integrating by parts the advective and diffusive terms within ${\textstyle {\bar {r}}}$ in the first integral of Eq.(10). This gives

Observation of the stabilization terms in the last integral of Eq.(11) shows that they introduce an additional diffusion (of value ${\textstyle {uh \over 2}}$) as well as an additional advection (of value ${\textstyle {sh \over 2}}$). The advective -type stabilization term is quite inconvenient for computational purposes as it leads to an antisymmetric stabilization matrix. In order to overcome this problem we will reformulate the stabilization terms in Eq.(11) with the aim of obtaining a single stabilizing term having the form of a diffusion.

With that purpose in mind we rewrite the integrand of the last term of Eq.(11) by splitting the constant source ${\textstyle Q}$ from the terms depending on ${\textstyle {\bar {\phi }}}$ in ${\textstyle {\bar {r}}}$ as

We define now

where ${\textstyle \beta ({\bar {\phi }})}$ is a stabilization parameter acting as a proportion of the actual physical diffusion ${\textstyle k}$.

Eqs. (12) and (13) allow to rewrite Eq.(11) as

Wee see clearly that the last term of Eq.(14) introduces within each element an additional diffusion of value ${\textstyle \beta k}$.

Substituting expression (5b) into (14) and choosing a Galerkin method with ${\textstyle W_{i}=N_{i}}$ within each element gives the discrete system of FE equations written in the standard matrix form as

where ${\textstyle {\bar {\boldsymbol {\phi }}}}$ is the vector of nodal unknowns and the element contributions to matrix K and vector ${\textstyle {\boldsymbol {f}}}$ are

The amount of balancing diffusion in Eq.(16) clearly depends on the (non linear) stabilization parameter ${\textstyle \beta }$. The element and critical values of ${\textstyle \beta }$ are deduced in the next section for linear two node elements.

We note that the integral of the term ${\textstyle {h \over 2}{dW_{i} \over dx}Q}$ in Eq.(17) vanishes after asssembly when ${\textstyle h}$ and ${\textstyle Q}$ are uniform over a patch of linear elements.

4 COMPUTATION OF THE STABILIZATION PARAMETER FOR LINEAR ELEMENTS

The matrix ${\textstyle {\boldsymbol {K}}^{e}}$ and the vector ${\textstyle {\boldsymbol {f}}^{e}}$ for two node linear elements are (for constant values of ${\textstyle u,k}$, ${\textstyle s}$ and ${\textstyle Q}$)

In Eqs.(18) and (19) index ${\textstyle e}$ denotes element values.

4.1 Computation of the element stabilization parameter β^e

Let us consider a mesh of just two elements 1 and 2 of equal size ${\textstyle l^{e}=l}$ with nodes 1, 2 and 3.

For the particular case of ${\textstyle Q=0}$ the following difference equation is found from Eq.(18)

where ${\textstyle \gamma =\displaystyle {ul \over 2k}}$ and ${\textstyle w=\displaystyle {sl^{2} \over k}}$ are the element Peclet number and a velocity independent dimensionless element number, respectively. The Damköler number ${\textstyle \sigma }$ can be obtained as ${\textstyle \sigma ={w \over 2\gamma }}$.

We will solve now for the value of ${\textstyle \phi _{2}}$ for Dirichlet boundary conditions at nodes 1 and 3.

The two following cases are considered.

Case A. ϕ₁=0, ϕ₃= ϕ₃^p

Application of the stencil of Eq.(20) gives the following value for ${\textstyle \phi _{2}}$

A physically meaningful solution for any values of ${\textstyle \gamma }$ and ${\textstyle w}$ must necessarily satisfy ${\textstyle \phi _{2}\geq 0}$. Clearly this requires

Case B. ϕ₁=ϕ^p₁, ϕ₃= 0

The stencil of Eq.(20) gives

Now ${\textstyle \phi _{2}\geq 0}$ if

From the results of above patch tests and the inspection of the modified governing equations assuming that the stabilization terms have the form of a diffusion, the following expression for the element stabilization parameter ${\textstyle \beta ^{e}}$ is proposed (see Appendix A):

with

and

where ${\textstyle {\bar {(\cdot )}}^{e}}$ denotes values computed at the element mid point.

It is clear from Eq.(26) that the computation of the stabilization parameter ${\textstyle \beta ^{e}}$ requires the knowledge of the sign of the numerical solution ${\textstyle {\bar {\phi }}}$ and that of the first and second derivatives of ${\textstyle {\bar {\phi }}}$ within each element. This necessarily leads to an iterative scheme. A simple algorithm which provides an stabilized and accurate solution in just two steps is presented in Section 4.3.

4.2 Critical stabilization parameter and unstability conditions

The following constant value of ${\textstyle \beta ^{e}}$ over the mesh yields a stabilized solution (for a uniform distribution of ${\textstyle \gamma }$ and ${\textstyle w}$)

where ${\textstyle \beta _{c}}$ is the critical stabilization parameter. The proof of Eq.(28) is given in the Appendix B.

The use of ${\textstyle \beta ^{e}=\beta _{c}}$ for all elements ensures a stabilized numerical solution for all ranges of ${\textstyle \gamma }$ and ${\textstyle w}$ although it can yield slight overdiffusive results in some cases.

The use of ${\textstyle \beta ^{e}=\beta _{c}}$ is also an excellent choice to obtain a first stabilized solution in the iterative scheme described in the next section.

The following condition is an indicator of a unstable solution (see the Appendix B)

4.3 Iterative solution scheme

The following two steps iterative scheme is proposed in order to obtain a stabilized and accurate solution.

step

Step 1. Compute a first stabilized solution ${\textstyle {\bar {\boldsymbol {\phi }}}^{1}}$ using the critical value ${\textstyle \beta ^{e}=\beta _{c}}$ given by Eq.(28). This ensures a stabilized, although sometimes slightly overdiffusive solution.

Step 2.

Step 2.1. Compute the signs of the first and second derivatives of ${\textstyle {\bar {\boldsymbol {\phi }}}^{1}}$ within each element. The second derivative field is obtained as follows

where ${\textstyle ({\hat {\cdot }})_{i}^{e}}$ denotes averaged values of the first derivative at node ${\textstyle i}$ of element ${\textstyle e}$. At the boundary nodes the constant value of the derivative of ${\textstyle {\bar {\phi }}}$ within the element is taken in Eq.(30); i.e. ${\textstyle ({\hat {\cdot }})_{i}^{e}=\left({d{\bar {\phi }} \over dx}\right)^{(e)}={{\bar {\phi }}_{2}-{\bar {\phi }}_{1} \over l^{e}}}$.

Step 2.2. Compute the enhanced stabilized solution ${\textstyle {\boldsymbol {\phi }}^{2}}$ using the element value of ${\textstyle \beta ^{e}}$ given by Eqs.(25) and (26).

In all the examples solved the above two steps have sufficied to obtain a converged stabilized solution. The reason of this is that the signs of the first and second derivative fields typically do not change any further after the second step over the relevant regions of the analysis domain (i.e. over regions adjacent to high gradient zones). We note that the stabilized solution is also remarkably accurate as shown in the examples presented in a next section.

4.4 Optimal value of the stabilization parameter

The optimal value of the stabilization parameter is defined as such that it guarantees a nodally exact solution. For 1D problems this parameter can be readily found making use of the exact analytical solution. For the zero source case this is written as

where

Constants ${\textstyle A}$ and ${\textstyle B}$ depend on the boundary conditions of the problem. For the Dirichlet problem (${\textstyle \phi =\phi _{1}^{p}}$ at ${\textstyle x=0}$ and ${\textstyle \phi =\phi _{N}^{p}}$ at ${\textstyle x=L}$) we obtain

Let us now write the stencil of Eq.(20) for three arbitrary nodes ${\textstyle i-1}$, ${\textstyle i}$, ${\textstyle i+1}$ and a uniform mesh with ${\textstyle l^{e}=l}$ and ${\textstyle \beta ^{i-1}=\beta ^{i}=\beta _{i}}$ as

where ${\textstyle \beta _{i}}$ is the nodal value of the stabilization parameter or the ${\textstyle i}$th node. Note that in Eq.(33) ${\textstyle \beta _{i}}$ has assumed to be constant within the tributary domain of node ${\textstyle i}$.

The optimal nodal value of the stabilization parameter ${\textstyle \beta _{i}^{opt}}$ is found from Eq.(33) as

where ${\textstyle \phi _{i}^{ex}}$ denotes the exact analytical value of ${\textstyle \phi _{i}}$ as obtained from Eq.(31a) by

Substituting the exact nodal values of Eq.(35) into Eq.(34) gives the analytical expression of ${\textstyle \beta _{i}^{opt}}$ as

Eqs.(34) and (36) indicate that ${\textstyle \beta _{i}^{opt}}$ depends on the coordinate of node ${\textstyle i}$. The advective (${\textstyle w=0}$) and absorptive (${\textstyle \gamma =0}$) limits of ${\textstyle \beta _{i}^{opt}}$ are given below.

In the examples presented in the next section we have found more convenient to compute ${\textstyle \beta _{i}^{opt}}$ by Eq.(34) with the numerical values of ${\textstyle \phi _{i}^{ex}}$ (${\textstyle i=1,2\cdots N+1}$) computed from Eq.(35).

4.5 Advective limit (w=0)

If the absorption terms are zero, ${\textstyle w=0}$, ${\textstyle a=2\gamma ,b=0}$ and from Eq.(36)

4.6 Absorptive limit (γ=0)

For ${\textstyle \gamma =0}$, ${\textstyle a=\delta }$, ${\textstyle b=-\delta }$, ${\textstyle \delta ={\sqrt {w}}}$ and Eq.(36) gives

Note that both the advective and absorptive limits of ${\textstyle \beta _{i}^{opt}}$ are independent of the nodal coordinates.

Remark. The optimal stabilization parameter for the advective limit (Eq.(37)) coincides with the classical expression obtained by many authors (see Vol. 3 in [1]).

The expression of Eq.(38) coincides with that obtained by Felippa and Oñate [15] using the FIC variational approach and the modified equation method.

5 VALUE OF THE ELEMENT CHARACTERISTIC LENGTH PARAMETER

The expression of the element characteristic length parameter ${\textstyle h^{e}}$ is deduced from Eq.(13) as

where ${\textstyle \beta ^{e}}$ is defined in Eq.(25).

It can be readily deduced from Eq.(39) that ${\textstyle h^{e}}$ vanishes in regions where ${\textstyle {d{\bar {\phi }} \over dx}}$ is small. In high gradient zones ${\textstyle h^{e}\to l^{e}}$ for ${\textstyle u\to \infty }$ (advective limit), whereas ${\textstyle h^{e}\to \pm {2 \over 3}l^{e}}$ for ${\textstyle s\to \infty }$ (absorptive limit). The sign of ${\textstyle h^{e}}$ in the latter case depends on the sign of ${\textstyle \displaystyle {S_{1} \over S_{0}}}$ where ${\textstyle S_{1}}$ and ${\textstyle S_{0}}$ are defined in Eq.(27).

The expression of ${\textstyle h^{e}}$ is needed in Eqs.(17) and (19) for the case of ${\textstyle Q\not =0}$.

6 ADVANTAGES OF USING β VERSUS h AS THE STABILIZATION PARAMETER

The stabilized formulation described above can be derived using either ${\textstyle \beta }$ or ${\textstyle h}$ as the stabilization parameter. The choice of ${\textstyle \beta }$ instead of ${\textstyle h}$ has the advantage that ${\textstyle \beta }$ is defined as a positive number (see Eq.(25)), whereas ${\textstyle h}$ is a distance that can be either positive or negative, depending on the values of the absorptive parameter ${\textstyle s}$ and of the sign ratio ${\textstyle \displaystyle {S_{1} \over S_{0}}}$ (see previous section). The change of sign of ${\textstyle h}$ within the mesh introduces a difficulty in the iterative solution process. This difficulty is avoided by selecting ${\textstyle \beta }$ as the stabilization parameter, which is the option chosen in this work.

Indeed, the characteristic length parameter ${\textstyle h}$ is still needed in the expression of ${\textstyle {\boldsymbol {f}}^{e}}$ (see Eqs.(17) and (19)). This does not pose any problem as the value (and sign) of ${\textstyle h}$ can be effectively computed in the second step from Eq.(39) using the stabilized solution obtained in the first step.

7 1D NUMERICAL EXAMPLES

The efficiency and accuracy of the iterative scheme described in the previous section is proved in a number of examples of application.

All examples presented in this section are solved in a 1D domain of discretized with eight two-node linear elements of unit size. The first set of examples is solved with the same Dirichlet conditions ${\textstyle {\bar {\phi }}_{1}=8}$ and ${\textstyle {\bar {\phi }}_{9}=3}$ and different values of ${\textstyle \gamma }$ and ${\textstyle w}$.

Tables 1–11 show the numerical results obtained at the nodes for the standard Galerkin method (${\textstyle \beta =0}$) and using the element (${\textstyle \beta ^{e}}$) and critical (${\textstyle \beta _{c}}$) values of ${\textstyle \beta }$. Recall that the results for ${\textstyle \beta _{c}}$ and ${\textstyle \beta ^{e}}$ correspond to those obtained respectively in the first and second step of the iterative algorithm described in Section 4.3. The exact analytical nodal solution is also presented as well as the optimal nodal value ${\textstyle \beta _{i}^{opt}}$ yielding nodaly exact results.