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)  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 . 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  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  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 . 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 , 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 . 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

 $r-{\underline {{h \over 2}{dr \over dx}}}{h \over 2}{dr \over dx}=0\quad {\hbox{in }}x\in (0,L)$
(1)
 $-u\phi +k{d\phi \over dx}+q^{p}-{\underline {{h \over 2}r}}{h \over 2}r=0\quad {\hbox{on }}\Gamma _{q}$
(2)
 $\phi -\phi ^{p}=0\quad {\hbox{on }}\Gamma _{\phi }$
(3)

where

 $r=:-u{d\phi \over dx}+{d \over dx}\left(k{d\phi \over dx}\right)-s\phi +Q$
(4)

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 . 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) .

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}$ . 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.

 $\phi \simeq {\bar {\phi }}\quad {\hbox{for}}\quad x\in [0,L]$
(5a)

with

 ${\bar {\phi }}=\sum \limits _{i=1}^{n}N_{i}\phi _{i}$
(5b)

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

 ${\bar {r}}-{h \over 2}{d{\bar {r}} \over dx}=r_{\Omega }\quad {\hbox{in }}x\in (0,L)$
(6)
 $-u{\bar {\phi }}+k{d{\bar {\phi }} \over dx}+q^{p}-{h \over 2}{\bar {r}}=r_{q}\quad {\hbox{on }}\Gamma _{q}$
(7)
 ${\bar {\phi }}-\phi ^{p}=r_{\phi }\quad {\hbox{on }}\Gamma _{\phi }$
(8)

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

 $\int _{L}W_{i}r_{\Omega }dx+[{\hat {W}}_{i}r_{q}]_{\Gamma _{q}}=0$
(9a)

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

 $\int _{L}W_{i}\left({\bar {r}}-{h \over 2}{d{\bar {r}} \over dx}\right)dx+\left[{\hat {W}}_{i}\left(-u{\bar {\phi }}+k{d{\bar {\phi }} \over dx}+q_{p}-{h \over 2}{\bar {r}}\right)\right]_{\Gamma _{q}}=0$
(9b)

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}}$

 $\int _{L}W_{i}{\bar {r}}dx-\left[W_{i}\left(-u{\bar {\phi }}+k{d{\bar {\phi }} \over dx}+q^{p}\right)\right]_{\Gamma _{q}}+\sum \limits _{e}\int _{l^{e}}{h \over 2}{dW_{i} \over dx}{\bar {r}}dx=0$
(10)

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

 $\int _{L}\left[u{dW_{i} \over dx}{\bar {\phi }}-{dW_{i} \over dx}k{d{\bar {\phi }} \over dx}-W_{i}s{\bar {\phi }}+W_{i}Q\right]dx-[W_{i}q^{p}]_{\Gamma _{q}}+\sum \limits _{e}\int _{l^{e}}{h \over 2}{dW_{i} \over dx}{\bar {r}}dx=0$
(11)

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

 ${h{\bar {r}} \over 2}{dW_{i} \over dx}={h \over 2}\left(-u{d{\bar {\phi }} \over dx}+k{d^{2}{\bar {\phi }} \over dx^{2}}-s{\bar {\phi }}\right){dW_{i} \over dx}+{h \over 2}{dW_{i} \over dx}Q=$ $=-{h \over 2}\left[u-\left(k{d^{2}{\bar {\phi }} \over dx^{2}}-s{\bar {\phi }}\right)\left({d{\bar {\phi }} \over dx}\right)^{-1}\right]{dW_{i} \over dx}{d{\bar {\phi }} \over dx}+{h \over 2}{dW_{i} \over dx}Q$
(12)

We define now

 ${h \over 2}\left[u-\left(k{d^{2}{\bar {\phi }} \over dx^{2}}-s{\bar {\phi }}\right)\left({d{\bar {\phi }} \over dx}\right)^{-1}\right]=\beta ({\bar {\phi }})k$
(13)

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

 $\int _{L}\left[u{dW_{i} \over dx}{\bar {\phi }}-{dW_{i} \over dx}k{d{\bar {\phi }} \over dx}-W_{i}s{\bar {\phi }}+W_{i}Q\right]dx-[W_{i}q^{p}]_{\Gamma _{q}}-\sum \limits _{e}\int _{l^{e}}\left(\beta k{dW_{i} \over dx}{d{\bar {\phi }} \over dx}-{h \over 2}{dW_{i} \over dx}Q\right)dx=0$
(14)

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

 ${\boldsymbol {K}}{\bar {\boldsymbol {\phi }}}={\boldsymbol {f}}$
(15)

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

 $K_{ij}^{e}\!\!=\!\!\int _{l^{e}}\left(-u{dN_{i} \over dx}N_{j}+k(1+\beta ){dN_{i} \over dx}{dN_{j} \over dx}+sN_{i}N_{j}\right)dx$ (16) $f_{i}^{e}\!\!=\!\!\int _{l^{e}}\left(N_{i}+{h \over 2}{dW_{i} \over dx}\right)Qdx-(N_{i}q^{p})_{\Gamma _{q}}$ (17)

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}$)

 ${\boldsymbol {K}}^{e}=-{u \over 2}\left[{\begin{matrix}-1&-1\\1&1\\\end{matrix}}\right]+{k \over l^{e}}(1+\beta ^{e})\left[{\begin{matrix}1&-1\\-1&1\\\end{matrix}}\right]+{sl^{e} \over 6}\left[{\begin{matrix}2&1\\1&2\\\end{matrix}}\right]$
(18)
 ${\boldsymbol {f}}^{e}={Ql^{e} \over 2}\left\{{\begin{matrix}1-\displaystyle {h^{e} \over 2}\\1+\displaystyle {h^{e} \over 2}\\\end{matrix}}\right\}\qquad +\quad {\hbox{boundary term}}$
(19)

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)

 $-\gamma \left(\phi _{1}-\phi _{3}\right)-(1+\beta ^{1})\phi _{1}+(2+\beta ^{1}+\beta ^{2})\phi _{2}-(1+\beta ^{2})\phi _{3}+{w \over 6}(\phi _{1}+4\phi _{2}+\phi _{3})=0$
(20)

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}}$

 $\phi _{2}={\phi _{3}^{p} \over (2+\beta ^{1}+\beta ^{2}+{2w \over 3})}\left(\beta ^{2}+1-\gamma -{w \over 6}\right)$
(21)

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

 $\beta ^{2}\geq \gamma +{w \over 6}-1$
(22)

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

The stencil of Eq.(20) gives

 $\phi _{2}={\phi _{1}^{p} \over (2+\beta ^{1}+\beta ^{2}+{2w \over 3})}\left(1+\beta ^{1}+\gamma -{w \over 6}\right)$
(23)

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

 $\beta ^{1}\geq {w \over 6}-\gamma -1$
(24)

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):

 ${\begin{array}{ll}\beta ^{e}={\bar {\beta }}^{e}\qquad {\hbox{if }}\quad {\bar {\beta }}^{e}>0\\\beta ^{e}=0\qquad {\hbox{if }}\quad {\bar {\beta }}^{e}\leq {0}\end{array}}$
(25)

with

 $\displaystyle {{\bar {\beta }}^{e}=\left[\left({S_{0} \over S_{2}}\right){w \over 6}+\left({S_{1} \over S_{2}}\right)\gamma -1\right]}$
(26)

and

 $S_{0}={\hbox{sign }}({\bar {\phi }})^{e},\quad S_{1}={\hbox{sign }}\left({d{\bar {\phi }} \over dx}\right)^{e},\quad S_{2}={\hbox{sign }}\left({d^{2}{\bar {\phi }} \over dx^{2}}\right)^{e}$
(27)

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}$)

 $\displaystyle {\beta ^{e}\geq \beta _{c}={w \over 6}+\vert \gamma \vert -1}$
(28)

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)

 ${w \over 6}+\vert \gamma \vert -1>0$
(29)

### 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

 $\left({d^{2}{\bar {\phi }}^{1} \over dx^{2}}\right)^{e}={1 \over l^{e}}\left[\left({d{\hat {\phi }}^{1} \over dx}\right)_{2}^{e}-\left({d{\hat {\phi }}^{1} \over dx}\right)_{1}^{e}\right]$
(30)

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

 $\phi =Ae^{ax \over l}+Be^{bx \over l}$
(31a)

where

 $a=\gamma +\delta \qquad ,\qquad b=\gamma -\delta \qquad ,\qquad \delta =(\gamma ^{2}+w)^{1/2}$
(31b)

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

 $A={\phi _{1}^{p}e^{bN}-\phi _{N}^{p} \over e^{bN}-e^{aN}}\quad ,\quad B={\phi _{N}^{p}-\phi _{1}^{p}e^{aN} \over e^{bN}-e^{aN}}$
(32)

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

 $-\gamma (\phi _{i-1}-\phi _{i+1})+(1+\beta _{i})(-\phi _{i-1}+2\phi _{i}-\phi _{i+1})+{w \over 6}(\phi _{i-1}+4\phi _{i}+\phi _{i+1})=0$
(33)

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

 $\beta _{i}^{opt}=\gamma {\phi _{i+1}^{ex}-\phi _{i-1}^{ex} \over \phi _{i+1}^{ex}-2\phi _{i}^{ex}+\phi _{i-1}^{ex}}+{w \over 6}\left(1+{6\phi _{i}^{ex} \over \phi _{i+1}^{ex}-2\phi _{i}^{ex}+\phi _{i-1}^{ex}}\right)-1$
(34)

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

 $\phi _{i}^{ex}=Ae^{ax_{i} \over l}+Be^{bx_{i} \over l}$
(35)

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

 $\beta _{i}^{opt}=\gamma \left[{\coth \displaystyle {a \over 2}+{B \over 2A}\displaystyle {\sinh b \over \sinh ^{2}\left(\displaystyle {a \over 2}\right)}e^{{b-a \over l}x_{i}} \over 1+{B \over A}\left(\displaystyle {\sinh \left(\displaystyle {b \over 2}\right)/\sinh \left(\displaystyle {a \over 2}\right)}\right)^{2}e^{{b-a \over l}x_{i}}}\right]+{w \over 6}\left[1+{3 \over 2}\left(1+\displaystyle {A \over B}e^{{a-b \over l}x_{i}}\right)/\sinh ^{2}\left(\displaystyle {b \over 2}\right)+\displaystyle {A \over B}e^{{a-b \over l}x_{i}}\sinh ^{2}\left(\displaystyle {a \over 2}\right)}\right]-1$
(36)

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).

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

 $\beta _{i}^{opt}=\gamma \coth \gamma -1$
(37)

### 4.6 Absorptive limit (γ=0)

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

 $\beta _{i}^{opt}={w \over 6}+{w \over 4\sinh ^{2}\left({{\sqrt {w}} \over 2}\right)}-1$
(38)

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 ).

The expression of Eq.(38) coincides with that obtained by Felippa and Oñate  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

 $h^{e}={2\beta ^{e}k\left({d{\bar {\phi }} \over dx}\right)^{e}/\left(u{d{\bar {\phi }} \over dx}-k{d^{2}{\bar {\phi }} \over dx^{2}}+s{\bar {\phi }}\right)^{e}}$
(39)

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.

 $\gamma =0.1$, $\omega =20$ ${\bar {\phi }}(\beta =0$) ${\bar {\phi }}(\beta ^{e}$) ${\bar {\phi }}(\beta _{c}$) $\phi$(exact) $\beta ^{e}$ $\beta _{c}$ $\beta _{i}^{opt}$ 1 8,00E+00 8,00E+00 8,00E+00 8,00E+00 2,233334 2,433334 – 2 -1,19E+00 2,69E-10 7,92E-02 1,01E-01 2,233334 2,433334 2,49 3 1,78E-01 9,07E-21 7,84E-04 1,27E-03 2,233334 2,433334 2,49 4 -2,69E-02 3,05E-31 7,76E-06 1,60E-05 2,233334 2,433334 2,49 5 6,06E-03 1,38E-41 7,69E-08 2,36E-07 2,433334 2,433334 2,51633 6 -1,35E-02 1,08E-31 7,61E-10 3,30E-06 2,433334 2,433334 2,64614 7 7,93E-02 3,27E-21 7,54E-12 3,20E-04 2,433334 2,433334 2,64624 8 -4,88E-01 9,90E-08 9,90E-08 3,10E-02 2,433334 2,433334 2,64624 9 3,00E+00 3,00E+00 3,00E+00 3,00E+00 – – –

 $\gamma =1$, $\omega =5$ ${\bar {\phi }}(\beta =0$) ${\bar {\phi }}(\beta ^{e}$) ${\bar {\phi }}(\beta _{c}$) $\phi$(exact) $\beta ^{e}$ $\beta _{c}$ $\beta _{i}^{opt}$ 1 8,00E+00 8,00E+00 8,00E+00 8,00E+00 0 0,8333 2 1,69E+00 1,69E+00 2,29E+00 1,88E+00 0 0,8333 0,2235 3 3,59E-01 3,59E-01 6,53E-01 4,41E-01 0 0,8333 0,2235 4 7,57E-02 7,59E-02 1,87E-01 1,03E-01 0 0,8333 0,2235 5 1,77E-02 1,61E-02 5,33E-02 2,43E-02 0 0,8333 0,2248 6 -6,97E-03 3,42E-03 1,52E-02 5,79E-03 0 0,8333 0,3641 7 6,93E-02 6,46E-04 4,35E-03 4,36E-03 0,8333 0,8333 1,038 8 -4,54E-01 1,85E-04 1,23E-03 9,56E-02 0,8333 0,8333 1,0681 9 3,00E+00 3,00E+00 3,00E+00 3,00E+00 – – –

 $\gamma =1$, $\omega =20$ ${\bar {\phi }}(\beta =0$) ${\bar {\phi }}(\beta ^{e}$) ${\bar {\phi }}(\beta _{c}$) $\phi$(exact) $\beta ^{e}$ $\beta _{c}$ $\beta _{i}^{opt}$ 1 8,00E+00 8,00E+00 8,00E+00 8,00E+00 1,3333 3,3333 – 2 -7,09E-01 2,96E-10 7,27E-01 2,22E-01 1,3333 3,3333 1,8645 3 6,32E-02 1,10E-20 6,61E-02 6,19E-03 1,3333 3,3333 1,8645 4 -7,18E-03 4,08E-31 6,01E-03 1,72E-04 1,3333 3,3333 1,8645 5 7,74E-03 1,02E-32 5,46E-04 4,78E-06 1,3333 3,3333 1,866 6 -3,27E-02 9,18E-32 4,97E-05 2,93E-07 3,3333 3,3333 3,2664 7 1,47E-01 2,75E-21 4,52E-06 4,25E-05 3,3333 3,3333 3,4167 8 -6,65E-01 9,09E-11 1,32E-06 1,13E-02 3,3333 3,3333 3,4167 9 3,00E+00 3,00E+00 3,00E+00 3,00E+00 – – –

 $\gamma =1$, $\omega =120$ ${\bar {\phi }}(\beta =0$) ${\bar {\phi }}(\beta ^{e}$) ${\bar {\phi }}(\beta _{c}$) $\phi$(exact) $\beta ^{e}$ $\beta _{c}$ $\beta _{i}^{opt}$ 1 8,00E+00 8,00E+00 8,00E+00 8,00E+00 18 20 2 -1,86E+00 0,00E+00 1,31E-01 3,66E-04 18 20 18,0054 3 4,34E-01 0,00E+00 2,15E-03 1,68E-08 18 20 18,0054 4 -1,04E-01 0,00E+00 3,52E-05 7,66E-13 18 20 18,0054 5 3,69E-02 0,00E+00 5,78E-07 7,09E-17 18 20 18,0072 6 -5,73E-02 0,00E+00 9,47E-09 6,97E-15 20 20 20,00075 7 2,02E-01 0,00E+00 1,55E-10 1,13E-09 20 20 20,00075 8 -7,76E-01 0,00E+00 2,55E-12 1,84E-04 20 20 20,0075 9 3,00E+00 3,00E+00 3,00E+00 3,00E+00 – – –

 $\gamma =2$, $\omega =0.1$ ${\bar {\phi }}(\beta =0$) ${\bar {\phi }}(\beta ^{e}$) ${\bar {\phi }}(\beta _{c}$) $\phi$(exact) $\beta ^{e}$ $\beta _{c}$ $\beta _{i}^{opt}$ 1 8,00E+00 8,00E+00 8,00E+00 8,00E+00 0,00E+00 1,02E+00 2 7,81E+00 7,80E+00 7,80E+00 7,80E+00 0,00E+00 1,02E+00 5,17E-05 3 7,61E+00 7,61E+00 7,61E+00 7,61E+00 0,00E+00 1,02E+00 5,03E-05 4 7,44E+00 7,43E+00 7,43E+00 7,43E+00 0,00E+00 1,02E+00 -3,15E-05 5 7,20E+00 7,24E+00 7,25E+00 7,24E+00 0,00E+00 1,02E+00 -4,75E-03 6 7,20E+00 7,07E+00 7,07E+00 7,07E+00 0,00E+00 1,02E+00 -3,66E-01 7 6,50E+00 6,89E+00 6,90E+00 6,89E+00 9,83E-01 1,02E+00 1,17E+00 8 7,91E+00 6,75E+00 6,73E+00 6,66E+00 9,83E-01 1,02E+00 1,09E+00 9 3,00E+00 3,00E+00 3,00E+00 3,00E+00 – – –

 $\gamma =2$, $\omega =2$ ${\bar {\phi }}(\beta =0$) ${\bar {\phi }}(\beta ^{e}$) ${\bar {\phi }}(\beta _{c}$) $\phi$(exact) $\beta ^{e}$ $\beta _{c}$ $\beta _{i}^{opt}$ 1 8,00E+00 8,00E+00 8,00E+00 8,00E+00 0 0,01667 – 2 5,09855 5,099631 5,333333 5,10362871 0 1,33333333 0,01900914 3 3,253624 3,250922 3,555556 3,25587825 0 1,33333333 0,01900914 4 2,063041 2,071956 2,37037 2,07709922 0 1,33333333 0,0190093 5 1,349646 1,321954 1,580247 1,32509295 0 1,33333333 0,01903119 6 0,7519671 0,8390264 1,053498 0,84535221 0 1,33333333 0,02196073 7 0,819374 0,5463428 0,702332 0,53967226 0 1,33333333 0,32747624 8 -0,5445008 0,3121959 0,4710313 0,3765327 0 1,33333333 1,36940125 9 3,00E+00 3,00E+00 3,00E+00 3,00E+00 1,33333333 1,33333333 1,36940125

 $\gamma =10$, $\omega =0.1$ ${\bar {\phi }}(\beta =0$) ${\bar {\phi }}(\beta ^{e}$) ${\bar {\phi }}(\beta _{c}$) $\phi$(exact) $\beta ^{e}$ $\beta _{c}$ $\beta _{i}^{opt}$ 1 8,00E+00 8,00E+00 8,00E+00 8,00E+00 0 9,016667 – 2 1,05E+01 7,97E+00 7,96E+00 7,96E+00 0 9,016667 2,10E-06 3 7,34E+00 7,92E+00 7,92E+00 7,92E+00 0 9,016667 2,10E-06 4 1,11E+01 7,90E+00 7,88E+00 7,88E+00 0 9,016667 2,10E-06 5 6,39E+00 7,84E+00 7,84E+00 7,84E+00 0 9,016667 2,10E-06 6 1,21E+01 7,82E+00 7,80E+00 7,80E+00 0 9,016667 2,10E-06 7 5,01E+00 7,75E+00 7,76E+00 7,76E+00 8,983333 9,016667 -4,44E-04 8 1,36E+01 7,72E+00 7,73E+00 7,72E+00 8,983333 9,016667 9,02E+00 9 3,00E+00 3,00E+00 3,00E+00 3,00E+00 – – –

 $\gamma =10$, $\omega =1$ ${\bar {\phi }}(\beta =0$) ${\bar {\phi }}(\beta ^{e}$) ${\bar {\phi }}(\beta _{c}$) $\phi$(exact) $\beta ^{e}$ $\beta _{c}$ $\beta _{i}^{opt}$ 1 8,00E+00 8,00E+00 8,00E+00 8,00E+00 0 9,16667 – 2 8,65E+00 7,69E+00 7,62E+00 7,61E+00 0 9,16667 0,00022111 3 6,94E+00 7,22E+00 7,26E+00 7,24E+00 0 9,16667 0,00022111 4 8,20E+00 6,98E+00 6,91E+00 6,89E+00 0 9,16667 0,00022111 5 5,81E+00 6,50E+00 6,58E+00 6,55E+00 0 9,16667 0,00022111 6 8,00E+00 6,36E+00 6,27E+00 6,23E+00 0 9,16667 0,00022111 7 4,54E+00 5,83E+00 5,97E+00 5,93E+00 8,833333 9,16667 0,00021822 8 8,14E+00 5,59E+00 5,69E+00 5,64E+00 8,833333 9,16667 9,22133708 9 3,00E+00 3,00E+00 3,00E+00 3,00E+00 – – –

 $\gamma =10$, $\omega =4$ $\phi (\beta =0$) $\phi (\beta ^{e}$) $\phi (\beta _{c}$) $\phi$(exact) $\beta ^{e}$ $\beta _{c}$ $\beta _{i}^{opt}$ 1 8,00E+00 8,00E+00 8,00E+00 8,00E+00 0 9,16667 – 2 6,206664 6,510888 6,666698 6,562702696 0 9,666666667 0,004136478 3 5,555403 5,408537 5,555607 5,383633335 0 9,666666667 0,004136478 4 3,952791 4,348897 4,629694 4,416398125 0 9,666666667 0,004136478 5 4,030291 3,682072 3,858096 3,622938486 0 9,666666667 0,004136478 6 2,27974 2,871269 3,215095 2,972033521 0 9,666666667 0,004136478 7 3,207678 2,549878 2,679258 2,43807155 0 9,666666667 0,004136709 8 0,8884292 1,838311 2,232884 2,000042344 0 9,666666667 9,137861599 9 3,00E+00 3,00E+00 3,00E+00 3,00E+00 9,666666667 – –

 $\gamma =10$, $\omega =20$ ${\bar {\phi }}(\beta =0$) ${\bar {\phi }}(\beta ^{e}$) ${\bar {\phi }}(\beta _{c}$) $\phi$(exact) $\beta ^{e}$ $\beta _{c}$ $\beta _{i}^{opt}$ 1 8,00E+00 8,00E+00 8,00E+00 8,00E+00 0 12,333334 2 2,94E+00 3,06E+00 4,00E+00 3,08E+00 0 12,333334 0,17281 3 1,32E+00 1,17E+00 2,00E+00 1,19E+00 0 12,333334 0,17281 4 1,80E-01 4,47E-01 1,00E+00 4,57E-01 0 12,333334 0,17281 5 5,99E-01 1,72E-01 5,00E-01 1,76E-01 0 12,333334 0,17281 6 -6,33E-01 6,46E-02 2,50E-01 6,77E-02 0 12,333334 0,17281 7 1,16E+00 2,64E-02 1,25E-01 2,61E-02 0 12,333334 0,17281 8 -1,83E+00 7,31E-03 6,25E-02 1,00E-02 12,333334 12,333334 12,2935 9 3,00E+00 3,00E+00 3,00E+00 3,00E+00 – – –
 $\gamma =10$, $\omega =120$ ${\bar {\phi }}(\beta =0$) ${\bar {\phi }}(\beta ^{e}$) ${\bar {\phi }}(\beta _{c}$) $\phi$(exact) $\beta ^{e}$ $\beta _{c}$ $\beta _{i}^{opt}$ 1 8,00E+00 8,00E+00 8,00E+00 8,00E+00 9 29 – 2 -9,18E-01 0,00E+00 1,14E+00 6,37E-02 9 29 9,8109 3 1,12E-01 0,00E+00 1,63E-01 5,08E-04 9 29 9,8109 4 -3,24E-02 0,00E+00 2,33E-02 4,05E-06 9 29 9,8109 5 5,67E-02 0,00E+00 3,33E-03 3,22E-08 9 29 9,8109 6 -1,50E-01 0,00E+00 4,76E-04 2,56E-10 9 29 9,8109 7 4,08E-01 0,00E+00 6,80E-05 9,14E-12 9 29 12,9409 8 -1,11E+00 0,00E+00 9,71E-06 1,31E-12 29 29 29 9 3,00E+00 3,00E+00 3,00E+00 3,00E+00 – – –

We note that the numbers in the tables for the values of ${\textstyle \beta }$ refer to element values for ${\textstyle \beta ^{e}}$ and ${\textstyle \beta _{c}}$ and to nodal values for the case of ${\textstyle \beta _{i}^{opt}}$.

Figures 1–6 show the distribution of ${\textstyle {\bar {\phi }}}$ along the domain for ${\textstyle \beta =0}$ (Galerkin), ${\textstyle \beta ^{e}}$ and ${\textstyle \beta _{c}}$ for some of the cases studied, as well as the exact solution. Figure 1: $\phi _{1}=8,\phi _{9}=3,\gamma =1$ and $\omega =5$. FIC results for a mesh of 8 linear elements obtained for $\beta =0$ (Galerkin), $\beta ^{e}$ and $\beta _{c}$. Comparison with the analytical solution (see Table 2) Figure 2: $\phi _{1}=8,\phi _{9}=3,\gamma =1$ and $\omega =20$. FIC results for a mesh of 8 linear elements obtained for $\beta =0$ (Galerkin), $\beta ^{e}$ and $\beta _{c}$. Comparison with the analytical solution (see Table 3) Figure 3: $\phi _{1}=8,\phi _{9}=3,\gamma =1$ and $\omega =120$. FIC results for a mesh of 8 linear elements obtained for $\beta =0$ (Galerkin), $\beta ^{e}$ and $\beta _{c}$. Comparison with the analytical solution (see Table 4) Figure 4: $\phi _{1}=8,\phi _{9}=3,\gamma =2$ and $\omega =2$. FIC results for a mesh of 8 linear elements obtained for $\beta =0$ (Galerkin), $\beta ^{e}$ and $\beta _{c}$. Comparison with the analytical solution (see Table 6) Figure 5: $\phi _{1}=8,\phi _{9}=3,\gamma =10$ and $\omega =4$. FIC results for a mesh of 8 linear elements obtained for $\beta =0$ (Galerkin), $\beta ^{e}$ and $\beta _{c}$. Comparison with the analytical solution (see Table 9) Figure 6: $\phi _{1}=8,\phi _{9}=3,\gamma =10$ and $\omega =20$. FIC results for a mesh of 8 linear elements obtained for $\beta =0$ (Galerkin), $\beta ^{e}$ and $\beta _{c}$. Comparison with the analytical solution (see Table 10)

We note that for the case ${\textstyle \gamma =10,\beta =4}$ (Table 9 and Figure 5) small oscillations around the exact distribution of ${\textstyle \phi }$ are obtained for ${\textstyle \beta =\beta ^{e}}$. This is due to the fact that in this case the stabilizing diffusion has been introduced in the last element of the mesh only. This suffices to stabilize the solution in the vecinity of the boundary layer at the right end of the domain. However it does not suffices to capture accurately the smooth distribution of ${\textstyle \phi }$ in the rest of the domain (the errors osciallate between 9% for node 8 to a maximum of 4% for the rest of the nodes). The oscillations are very much reduced if an additional diffusion of similar value is introduced in the next element in the mesh.

The first set of examples is complemented by two similar problems with Dirichlet conditions ${\textstyle {\bar {\phi }}_{1}=0}$ and ${\textstyle {\bar {\phi }}_{8}=1}$ and ${\textstyle {\bar {\phi }}_{1}=1}$ and ${\textstyle {\bar {\phi }}_{8}=0}$, respectively. Both cases are solved for ${\textstyle \gamma =10}$ and ${\textstyle w=200}$. The nodal results are shown in Tables 12 and 13, while the distributions of ${\textstyle {\bar {\phi }}}$ are plotted in Figures 7 and 8.

 $\gamma =10$, $\omega =200$ ${\bar {\phi }}(\beta =0$) ${\bar {\phi }}(\beta ^{e}$) ${\bar {\phi }}(\beta _{c}$) $\phi$(exact) $\beta ^{e}$ $\beta _{c}$ $\beta _{i}^{opt}$ 1 1 1 1 1 22,3333333 42,3333333 – 2 -0,1745558 3,7037E-11 0,09090909 0,00066183 22,3333333 42,3333333 22,4526286 3 0,03046973 1,3717E-21 0,00826446 4,3801E-07 22,3333333 42,3333333 22,4526286 4 -0,00531866 5,0805E-32 0,00075131 2,8989E-10 22,3333333 42,3333333 22,4526286 5 0,00092839 1,8817E-42 6,8301E-05 1,9186E-13 22,3333333 42,3333333 22,4526286 6 -0,00016203 6,9692E-53 6,2092E-06 1,2698E-16 22,3333333 42,3333333 22,4526286 7 2,8194E-05 2,5812E-63 5,6447E-07 8,4035E-20 22,3333333 42,3333333 22,4526286 8 -4,6527E-06 9,5599E-74 5,1316E-08 5,5617E-23 22,3333333 42,3333333 22,4526199 9 0 0 0 0 – – –

 $\gamma =10$, $\omega =200$ ${\bar {\phi }}(\beta =0$) ${\bar {\phi }}(\beta ^{e}$) ${\bar {\phi }}(\beta _{c}$) $\phi$(exact) $\beta ^{e}$ $\beta _{c}$ $\beta _{i}^{opt}$ 1 0 0 0 0 42,3333333 42,3333333 – 2 -0,00040908 2,3464E-81 2,3464E-81 8,7898E-84 42,3333333 42,3333333 42,3333333 3 0,00130776 7,7431E-70 7,7431E-70 6,4435E-72 42,3333333 42,3333333 42,3333333 4 -0,0039649 2,5552E-58 2,5552E-58 4,7236E-60 42,3333333 42,3333333 42,3333333 5 0,01198527 8,4323E-47 8,4323E-47 3,4627E-48 42,3333333 42,3333333 42,3333333 6 -0,03622341 2,7826E-35 2,7826E-35 2,5384E-36 42,3333333 42,3333333 42,3333333 7 0,1094779 9,1827E-24 9,1827E-24 1,8608E-20 42,3333333 42,3333333 42,3333333 8 -0,3308744 3,0303E-12 3,0303E-12 1,3641E-12 42,3333333 42,3333333 42,3333333 9 1 1 1 1 – – – Figure 7: $\phi _{1}=8,\phi _{9}=3,\gamma =10$ and $\omega =200$. FIC results for a mesh of 8 linear elements obtained for $\beta =0$ (Galerkin), $\beta ^{e}$ and $\beta _{c}$. Comparison with the analytical solution (see Table 12) Figure 8: $\phi _{1}=8,\phi _{9}=3,\gamma =1$ and $\omega =20$. FIC results for a mesh of 8 linear elements obtained for $\beta =0$ (Galerkin), $\beta ^{e}$ and $\beta _{c}$. Comparison with the analytical solution (see Table 13)

In references [12,13] it was shown that the SUPG, GLS and SGS methods with a single stabilization parameter fail to provide a stabilized solution for the case ${\textstyle {\bar {\phi }}_{1}=1}$, ${\textstyle {\bar {\phi }}_{8}=0}$ with ${\textstyle \gamma =10}$ and ${\textstyle w=200}$. The FIC results, however, are physically correct as shown in Table 13 and Figure 8.

The following conclusions are drawn from the results obtained for the 1D examples.

1. The Galerkin solution (${\textstyle \beta =0}$) is unstable for the cases where the instability condition (Eq.(29)) is violated, as expected.
2. The solution obtained with the critical value ${\textstyle \beta _{c}}$ is always stable, although it yields slightly overdiffusive results in some cases (see Tables 2, 6, 9 and 10).
3. The results obtained with ${\textstyle \beta ^{e}}$ are less diffusive than those obtained with ${\textstyle \beta _{c}}$ and quite accurate for all ranges of ${\textstyle \gamma }$ and ${\textstyle w}$.
4. The FIC algorithm provides a stabilized solution for Dirichlet problems when there is a negative streamwise gradient of the solution. This is an advantage versus SUPG, GLS and SGS methods using a single stabilization parameter which fail for these type of problems [12,13].

## 8 EXTENSION TO MULTI-DIMENSIONAL PROBLEMS

The stabilization procedure presented can be easily extended to multi-dimensional problems. Only the main ingredients of the formulation are presented here. For a more detailed description of the multi-dimensional formulation see .

Consider the general steady-state advection-diffusion-reaction equations written for the sourceless case as

 $r(\phi ):=-{\boldsymbol {u}}^{T}{\boldsymbol {\nabla }}\phi +{\boldsymbol {\nabla }}^{T}{\boldsymbol {D}}{\boldsymbol {\nabla }}\phi -s\phi =0$
(40)

where for 2D problems

 ${\boldsymbol {u}}=[u,v]^{T}\quad ,\quad {\boldsymbol {\nabla }}=\left[{\partial \over \partial x},{\partial \over \partial y}\right]^{T}\quad ,\quad {\boldsymbol {D}}=k\left[{\begin{matrix}1&0\\0&1\\\end{matrix}}\right]$
(41)

are respectively the velocity vector, the gradient vector and the diffusivity matrix, respectively. For simplicity we have assumed an isotropic diffusion matrix.

The FIC form of Eq.(40) is written

 $r-{\underline {{1 \over 2}{\boldsymbol {h}}^{T}{\boldsymbol {\nabla }}r}}{1 \over 2}{\boldsymbol {h}}^{T}{\boldsymbol {\nabla }}r=0$
(42)

where ${\textstyle r}$ is the original differential equation as defined in Eq.(40).

The Dirichlet and boundary conditions of the FIC formulation are

 $\phi -\phi ^{p}=0\quad {\hbox{on}}\quad \Gamma _{\phi }$
(43)

 $-{\boldsymbol {u}}^{T}{\boldsymbol {n}}\phi +{\boldsymbol {n}}^{T}{\boldsymbol {D}}{\boldsymbol {\nabla }}\phi +q^{p}-{\underline {{1 \over 2}{\boldsymbol {h}}^{T}{\boldsymbol {n}}r}}{1 \over 2}{\boldsymbol {h}}^{T}{\boldsymbol {n}}r=0$
(44)

where ${\textstyle n}$ is the normal vector to the boundary where the normal flux is prescribed. As usual index ${\textstyle p}$ denotes the prescribed values.

In Eqs.(42) and (44) ${\textstyle {\boldsymbol {h}}=[h_{x},h_{y}]^{T}}$ is the characteristic vector of the FIC formulation which components play the role of stabilization parameters. The underlined terms in Eqs.(42) and (44) introduce the necessary stability in the numerical solution.

As shown in  the components of the characteristic length h allow to obtain stabilized finite element solutions in the presence of strong gradients of ${\textstyle \phi }$ near the boundaries (boundary layers) and within the analysis domain (internal layers). If vector h is taken to be parallel to the velocity u the standard SUPG method is recovered . Keeping the more general form of h allows to reproduce the best features of the so called shock capturing or transverse dissipation schemes.

Oñate et al.  have recently shown that a general form of the characteristic lengths can be found by writting the FIC equations along the principal curvature directions of the solution denoted hereafter as ${\textstyle \xi }$ and ${\textstyle \eta }$. The resulting FIC-FEM formulation is equivalent in this case (for linear elements) to adding a non linear stabilizing diffusion matrix to the standard infinitessimal equation. The stabilizing diffusion matrix depends on the velocities along the principal curvature directions of the solution which must be found iteratively. The iterative scheme can be simplified by assuming that the main principal curvature direction at each integration point within an element are coincident with the gradient vector direction at that point. For details see . This technique is extended here accounting for the absorption term.

Following the ideas described in  the FIC equations can be written as

 $-{\boldsymbol {u}}^{T}{\boldsymbol {\nabla }}\phi +{\boldsymbol {\nabla }}^{T}({\boldsymbol {D}}+{\bar {\boldsymbol {D}}}){\boldsymbol {\nabla }}\phi -s\phi =0$
(45)

where the stabilization matrix ${\textstyle {\bar {\boldsymbol {D}}}}$ is given by

 ${\bar {\boldsymbol {D}}}={\boldsymbol {T}}^{T}{\bar {\boldsymbol {D}}}'{\boldsymbol {T}}$
(46)

 ${\bar {\boldsymbol {D}}}'=\left[{\begin{matrix}k_{\xi }{+}k_{\xi }^{s}&\vdots &0\\\cdots &\cdots &\cdots \\0&\vdots &k_{\eta }+k_{\eta }^{s}\\\end{matrix}}\right]$
(47)

and

 ${\boldsymbol {T}}=\left[{\begin{matrix}\cos \alpha &\sin \alpha \\-\sin \alpha &\cos \alpha \\\end{matrix}}\right]$
(48)

where ${\textstyle \alpha }$ is the angle that the main principal curvature direction of the solution ${\textstyle \xi }$ forms with the global ${\textstyle x}$ axis. The different terms in Eq.(47–48) are

 ${\begin{array}{l}\displaystyle k_{\xi }=\alpha _{\xi }{u_{\xi }l_{\xi } \over 2}~,~k_{\eta }=\alpha _{\eta }{u_{\eta }l_{\eta } \over 2}\\\displaystyle \alpha _{\xi }=\coth \gamma _{\xi }-{1 \over \gamma _{\xi }}~,~\alpha _{\eta }=\coth \gamma _{\eta }-{1 \over \gamma _{\eta }}~,~\gamma _{\xi }={u_{\xi }l_{\xi } \over 2k}~,~\gamma _{\eta }={u_{\eta }l_{\eta } \over 2k}\end{array}}$
(49)
 $k_{\xi }^{s}={sl_{\xi }^{2} \over 6}~,~k_{\eta }^{s}={sl_{\eta }^{2} \over 6}$

where ${\textstyle u_{\xi }}$, ${\textstyle u_{\eta }}$ are the velocity components along each of the principal curvature directions and ${\textstyle l_{\xi }}$ and ${\textstyle l_{\eta }}$ are the maximum values of the projection of the element sides along the principal directions ${\textstyle \xi }$ and ${\textstyle \eta }$, respectively.

The values of ${\textstyle k_{\xi }^{s}}$ and ${\textstyle k_{\eta }^{s}}$ in Eq.(49) correspond to the critical stabilization parameters for the absorptive limit, as defined for the 1D problem. A non linear form of these parameters depending on the sign of the solution derivatives along the principal curvature directions can be obtained [23,34].

The derivation of the finite element equations follows the standard Galerkin procedure applied to Eq.(4345). For details see [23,34].

The stabilized numerical results presented in the next section have been found by the following two steps algorithm.

Step 1. (Galerkin step). Solve the problem for ${\textstyle {}^{1}{\bar {\boldsymbol {\phi }}}}$ with a zero value of the the stabilization matrix (i.e. ${\textstyle {\bar {\boldsymbol {D}}}={\boldsymbol {0}}}$). Unstable results characterized by undershoots and overshoots in the numerical solution will be found for high values of the advective and the absorption terms.

Step 2. Compute the principal curvature directions of the solution at the element center. As mentioned above we have chosen the main principal direction ${\textstyle \xi }$ to be parallel to the gradient vector ${\textstyle {\boldsymbol {\nabla }}\phi }$ for simplicity. The direction ${\textstyle \eta }$ is chosen orthogonal to ${\textstyle \xi }$.

Compute ${\textstyle {\bar {\boldsymbol {D}}}'}$, ${\textstyle {\boldsymbol {T}}}$ and ${\textstyle {\bar {\boldsymbol {D}}}}$ from Eqs.(4648). Solve for ${\textstyle {}^{2}{\bar {\boldsymbol {\phi }}}}$.

Excellent results are obtained with this simple two steps algorithm for problems with boundary layers orthogonal to the velocity vector as shown in the examples presented in the next section. A refined version of the algorithm applicable to problems with internal layers is described in .

### 8.1 2D Examples

The analysis domain in the first two 2D examples presented is a square of size 8 units. A relatively coarse mesh of ${\textstyle 8\times 8}$ four node bi-linear square elements is used in all cases.

The boundary conditions for the first and second examples are ${\textstyle \phi ^{p}=8}$ at the boundary ${\textstyle x=0}$ and ${\textstyle \phi =3}$ at ${\textstyle x=8}$ and zero normal flux at ${\textstyle y=0}$ and ${\textstyle y=8}$. This reproduces the condition of the 1D examples solved in previous sections. The first example is analized for ${\textstyle {\boldsymbol {u}}=[2,0]^{T}}$, ${\textstyle k=1}$ and ${\textstyle s=1}$ giving ${\textstyle w=20}$, ${\textstyle \gamma _{x}=1}$ and ${\textstyle \gamma _{y}=0}$ which corresponds to the 1D example analyzed in Table 3. The correct solution for this absorptive dominant case is a boundary layer in the vecinity of the side at ${\textstyle x=8}$ where ${\textstyle \phi }$ is prescribed to a value of ${\textstyle 3}$. The numerical results obtained with the standard (non stabilized) Galerkin solution (i.e. ${\textstyle {\bar {\boldsymbol {D}}}=0}$) are oscillatory as expected. The stabilized FIC formulation elliminates the oscillations and yields the correct physical solution as shown in Figure 9. Note that the 2D solution coincides precisely with the 1D results of Table 3 and Figure 2 for ${\textstyle \beta =\beta _{c}}$. Figure 9: 2D advection-conduction-absorption problem over a square domain of size equal to 8 units. $\phi ^{p}=8$ at $x=0$, $\phi ^{p}=3$ at $x=8$, $q_{n}=0$ at $y=0$ and $y=8$. ${\boldsymbol {u}}=[2,0]^{T}$, $k=1$, $s=20$, $w=20$, $\gamma _{x}=1$ and $\gamma _{y}=0$. Galerkin and FIC solutions obtained with a mesh of $8\times 8$ four node square elements.

The second example is similar to the first one with ${\textstyle {\boldsymbol {u}}=[20,0]^{T}}$, ${\textstyle k=1}$ and ${\textstyle s=20}$ giving ${\textstyle w=20}$, ${\textstyle \gamma _{x}=10}$ and ${\textstyle \gamma _{y}=0}$. This problem corresponds to the 1D case of Table 10. The Galerkin solution is again oscillatory, whrereas the FIC results are physically sound (see Figure 10). Once more the 2D FIC results coincide with the 1D values shown in Table 10 and Figure 6 for ${\textstyle \beta =\beta _{c}}$. Figure 10: 2D advection-conduction-absorption problem over a square domain of size equal to 8 units. $\phi ^{p}=8$ at $x=0$, $\phi ^{p}=3$ at $x=8$, $q_{n}=0$ at $y=0$ and $y=8$. ${\boldsymbol {u}}=[20,0]^{T}$, $k=1$, $s=20$, $w=20$, $\gamma _{x}=10$ and $\gamma _{y}=0$. Galerkin and FIC solutions obtained with a mesh of $8\times 8$ four node square elements.

The boundary conditions for the third 2D example are ${\textstyle \phi =10}$ at ${\textstyle x=0}$ and ${\textstyle y=10}$ and ${\textstyle \phi =0}$ at ${\textstyle y=0}$ and ${\textstyle x=10}$. The square domain has a side of 10 units in this case. A mesh of ${\textstyle 10\times 10}$ four node elements of unit length is used. The physical parameters are ${\textstyle {\boldsymbol {u}}=[0,1\times 10^{-3}]^{T}}$, ${\textstyle k=1}$ and ${\textstyle s=12}$ giving ${\textstyle w=12}$, ${\textstyle \gamma _{x}=0}$ and ${\textstyle \gamma _{y}=5\times 10^{-4}}$.

The correct solution for this absorption dominated problem has two boundary layers next to the boundaries at ${\textstyle x=0}$ and ${\textstyle y=8}$. The Galerkin solution is oscillatory in these zones as expected (Figure 11) whereas the FIC results capture the boundary layer within the elements adjacent to the boundary. Figure 11: 2D advection-conduction-absorption problem over a square domain of size equal to 10 units. $\phi ^{p}=8$ at $x=0$, $\phi ^{p}=3$ at $x=10$, $q_{n}=0$ at $y=0$ and $y=10$. ${\boldsymbol {u}}=[0,1\times 10^{-3}]^{T}$, $k=1$, $s=12$, $w=12$, $\gamma _{x}=0$ and $\gamma _{y}=5\times 10^{-4}$. Galerkin and FIC solutions obtained with a mesh of $10\times 10$ four node square elements.

The boundary conditions for the last example are ${\textstyle \phi =50}$ at ${\textstyle x=0}$ and ${\textstyle y=10}$ and ${\textstyle \phi =100}$ at ${\textstyle y=0}$ and ${\textstyle x=10}$. The physical parameters are ${\textstyle {\boldsymbol {u}}=[0,5]^{T}}$, ${\textstyle k=1}$ and ${\textstyle s=12}$ giving ${\textstyle w=12}$, ${\textstyle \gamma _{x}=0}$ and ${\textstyle \gamma _{y}=2.5}$.

The correct solution has two boundary layers next to the boundaries at ${\textstyle x=10}$ and ${\textstyle y=10}$. The Galerkin solution gives oscillatory results next to these boundaries and also next to the boundary at ${\textstyle x=0}$ (Figure 12). All these oscillations disappear with the stabilized FIC formulation here proposed. Figure 12: 2D advection-conduction-absorption problem over a square domain of size equal to 10 units. $\phi ^{p}=8$ at $x=0$, $\phi ^{p}=3$ at $x=10$, $q_{n}=0$ at $y=0$ and $y=10$. ${\boldsymbol {u}}=[0,5]^{T}$, $k=1$, $s=12$, $w=12$, $\gamma _{x}=0$ and $\gamma _{y}=2.5$. Galerkin and FIC solutions obtained with a mesh of $10\times 10$ four node square elements.

More accurate results for all these problems can be found by using the non linear form of the element absorptive stabilization parameter ${\textstyle \beta ^{e}}$ similarly as for the 1D case. Details of the iterative scheme and examples of applications are given in .

## 9 CONCLUSIONS

The FEM-FIC formulation presented allows to obtain a stabilized solution for the advection-diffusion-absorption equation. For the 1D problem the formulation is equivalent to adding a non-linear diffusion term to the standard discretized equations, together with a consistent modification of the constant source term for each element. The stabilizing diffusion term is governed by a single stabilization parameter. The use of the constant critical value of the stabilization parameter for all the elements in the mesh allows to obtain a stabilized numerical solution in a single step. A more accurate (less diffusive) solution can be obtained using the two-step iterative scheme proposed. This obviously increases the cost of the solution and for many practical cases the use of the critical stabilization parameter may be more attractive. The 1D examples presented demonstrate that the method provides stabilized and very accurate results for all range of the physical parameters and boundary conditions using a quite coarse mesh.

The equivalence of the FIC method with a stabilizing diffusion term extends naturally to multidimensional problems, where the advantages of the proposed approach should prove to be very advantageous as indicated by the excellent results obtained in the 2D examples solved.

## ACKNOWLEDGEMENTS

The numerical solutions for the 2D problems were obtained by Dr. F. Zarate from CIMNE whose contribution is gratefully acknowledged.

The authors also thank Profs. C. Felippa and S.R. Idelsohn for many useful discussions.

This work has been sponsored by the Spanish Ministry of Educación y Ciencia (Plan Nacional, Project number: DPI2001-2193).

### Document information Published on 01/01/2006

DOI: 10.1016/j.cma.2005.07.020

### Document Score 0

Times cited: 20
Views 50
Recommendations 0