Stabilized formulation for the advection-diffusion-absorption equation using finite calculus and linear finite elements



==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) <span id='citeF-1'></span>[[#cite-1|[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 <span id='citeF-2'></span>[[#cite-2|[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)  <span id='citeF-2'></span>[[#cite-2|[2]]], Galerkin Least Square <span id='citeF-5'></span>[[#cite-5|[5]]], Subgrid Scale (SGS) <span id='citeF-5'></span>[[#cite-5|[5]]] and Residual Free Bubbles <span id='citeF-14'></span>[[#cite-14|[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 <span id='citeF-1'></span>[[#cite-1|[1]]] and <span id='citeF-8'></span>[[#cite-8|[8]]]), 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 <span id='citeF-4'></span>[[#cite-4|[4]]]. As reported in <span id='citeF-12'></span>[[#cite-12|[12]]] 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 <span id='citeF-16'></span>[[#cite-16|[16]]].

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 <span id='citeF-20'></span>[[#cite-20|[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 <span id='citeF-20'></span>[[#cite-20|[20]]], diffusion-absorption <span id='citeF-15'></span>[[#cite-15|[15]]], incompressible fluid flow <span id='citeF-26'></span>[[#cite-26|[26]]] and standard and incompressible solid mechanics <span id='citeF-31'></span>[[#cite-31|[31]]] 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 <span id='citeF-15'></span>[[#cite-15|[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

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>r - \underline{{h\over 2} {dr\over dx}}{h\over 2} {dr\over dx}=0\quad \hbox{in } x \in (0,L) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>-u\phi + k {d\phi \over dx} +q^p - \underline{{h\over 2} r}{h\over 2} r=0\quad \hbox{on }\Gamma _q  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\phi -\phi ^p =0 \quad \hbox{on }\Gamma _\phi  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
|}

where

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>r =:-u {d\phi \over dx} +{d\over dx} \left(k{d\phi \over dx}\right)- s\phi + Q </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
|}

In above equations <math display="inline">\phi </math> is the state variable, <math display="inline">x \in [0,L]</math> is the problem domain, <math display="inline">L</math> is the domain length, <math display="inline">u</math> is the velocity field, <math display="inline">k\ge 0</math> is the diffusion, <math display="inline">s\ge 0</math> is the absorption, dissipation or destruction source parameter, <math display="inline">Q</math> is a constant source term, <math display="inline">q^p</math> and <math display="inline">\phi ^p</math> are the prescribed values of the total flux and the unknown function at the Neumann and Dirichlet boundaries <math display="inline">\Gamma _q</math> and <math display="inline">\Gamma _\phi </math>, respectively and <math display="inline">h</math> is a ''characteristic length'' which plays the role of a stabilization parameter. In the 1D problem <math display="inline">\Gamma _\phi </math> and <math display="inline">\Gamma _q</math> 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 <span id='citeF-20'></span>[[#cite-20|[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 <math display="inline">h</math> tends to zero the FIC differential equations gradually recover the standard infinitessimal form giving in the limit (for <math display="inline">h=0</math>) <math display="inline">r=0</math> in <math display="inline">x\in (0,L)</math> and <math display="inline">-u \phi + k{d\phi \over dx}+ q^p=0\quad \hbox{on }\Gamma _q</math>.

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 <math display="inline">h</math>. At the discrete level <math display="inline">h</math> 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 <math display="inline">h</math> as a proportion of a typical element dimension (i.e. the element length for 1D problems) <span id='citeF-20'></span>[[#cite-20|20]].

Successful applications of the FIC method to a variety of problems in computational mechanics can be found in <span id='citeF-20'></span>[[#cite-20|[20]]].

==3 FINITE ELEMENT FORMULATION==

We will construct a standard finite element discretization <math display="inline">\left\{l^e\right\}</math> of the 1D analysis domain length <math display="inline">L</math> with index <math display="inline">e</math> ranging from 1 to the number of elements <math display="inline">N</math> <span id='citeF-1'></span>[[#cite-1|[1]]]. The state variable <math display="inline">\phi </math> is approximated by <math display="inline">\bar \phi </math> over the analysis domain. The approximated variable <math display="inline">\bar \phi </math> is interpolated within each element with <math display="inline">n</math> nodes in the standard manner, i.e.

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\phi \simeq \bar \phi \quad \hbox{for}\quad x \in [0,L]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" |  (5a)
|}

with

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\bar \phi =\sum \limits _{i=1}^n N_i \phi _i</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" |  (5b)
|}

where <math display="inline">N_i</math> are the element shape functions and <math display="inline">\phi _i</math> are nodal values of the approximate function <math display="inline">\bar \phi </math>. Substituting Eq.(5a) into Eqs.(1) and (2) gives

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\bar r - {h\over 2} {d\bar r\over dx} =r_\Omega \quad \hbox{in } x\in (0,L) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>-u\bar \phi + k {d\bar \phi \over dx} +q^p - {h\over 2} \bar r=r_q\quad \hbox{on }\Gamma _q  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\bar \phi -\phi ^p =r_\phi \quad \hbox{on }\Gamma _\phi  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
|}

where <math display="inline">\bar r =r(\bar \phi )</math> and <math display="inline">r_\Omega , r_q</math> and <math display="inline">r_\phi </math> are the residuals of the approximate solution in the problem domain and on the Neumann and Dirichlet boundaries <math display="inline">\Gamma _q</math> and <math display="inline">\Gamma _\phi </math>, respectively.

The weighted residual form of Eqs.(6)&#8211;(8) is written as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int _L W_i r_\Omega dx +   [\hat W_i r_q]_{\Gamma _q} =0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" |  (9a)
|}

where <math display="inline">W_i(x)</math> and <math display="inline">\hat W_i</math> are test functions satisfying <math display="inline">W_i =\hat W_i =0</math> on <math display="inline">\Gamma _\phi </math>. In Eq.(9a) and in the following <math display="inline">\int _L (\cdot ) dx :=  \int _0^L (\cdot ) dx</math> and <math display="inline">[g]_{\Gamma _q}</math> denotes the value of <math display="inline">g</math> at the boundary <math display="inline">\Gamma _q</math>. Using Eqs.(6) and (7) into Eq.(9a) we have

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\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</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" |  (9b)
|}

Assuming smooth enough solutions and integrating by parts the term involving <math display="inline">h</math> in the first integral gives for <math display="inline">\hat W_i=-W_i</math>

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int _L W_i \bar r d x - \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</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (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 <math display="inline">\bar r</math> to be discontinuous. Also in  Eq.(10) <math display="inline">h</math> has been assumed to be constant within each element, (i.e. <math display="inline">{dh\over dx}=0</math> within <math display="inline">l^e</math>).

The  weak form is obtained by integrating by parts the advective and diffusive terms within <math display="inline">\bar r</math> in the first integral of Eq.(10). This gives

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\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 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
|}

Observation of the stabilization terms in the last integral of Eq.(11) shows that they introduce an additional diffusion (of value <math display="inline">{uh\over 2}</math>) as well as an additional advection (of value <math display="inline">{sh\over 2}</math>). 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 <math display="inline">Q</math> from the terms depending on <math display="inline">\bar \phi </math> in <math display="inline">\bar r</math> as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{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=</math>
|-
| style="text-align: center;" | <math> = -{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 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
|}

We define now

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{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 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (13)
|}

where <math display="inline">\beta (\bar \phi )</math> is a stabilization parameter acting as a proportion of the actual physical diffusion <math display="inline">k</math>.

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

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\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 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
|}

Wee see clearly that the last term of Eq.(14) introduces within each element an additional diffusion of value <math display="inline">\beta k</math>.

Substituting expression (5b) into (14) and choosing a Galerkin method with <math display="inline">W_i =N_i</math> within each element gives the discrete system of FE equations written in the standard matrix form as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\textbf K}\bar{\boldsymbol \phi } ={\textbf f} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
|}

where <math display="inline">\bar{\boldsymbol \phi }</math> is the vector of nodal unknowns and the element contributions to matrix <math display="inline">\textbf  K</math> and vector <math display="inline">\textbf  f</math> are

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>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</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (16)
|-
| style="text-align: center;" | <math> f_i^e\!\! = \!\! \int _{l^e} \left(N_i + {h\over 2} {dW_i\over dx} \right)Qdx - (N_i q^p)_{\Gamma _q} </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (17)
|}
|}

The amount of balancing diffusion in Eq.(16) clearly depends on the (non linear) stabilization parameter <math display="inline">\beta </math>. The element and critical  values of <math display="inline">\beta </math> are deduced in the next section for linear two node elements.

We note that the integral of the term <math display="inline">{h\over 2} {dW_i\over dx}Q</math> in Eq.(17) vanishes after asssembly when <math display="inline">h</math> and <math display="inline">Q</math> are uniform over a patch of linear elements.

==4 COMPUTATION OF THE STABILIZATION PARAMETER FOR LINEAR ELEMENTS==

The matrix <math display="inline">{\textbf K}^e</math> and the vector <math display="inline">{\textbf f}^e</math> for two node linear elements are (for constant values of <math display="inline">u,k</math>, <math display="inline">s</math> and <math display="inline">Q</math>)

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\textbf 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] </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (18)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\textbf 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} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (19)
|}

In Eqs.(18) and (19) index <math display="inline">e</math> denotes element values.

===4.1 <span id='lb-4.1'></span>Computation of the element stabilization parameter β<sup>e</sup>===

Let us consider a mesh of just two elements 1 and 2 of equal size <math display="inline">l^e =l</math>  with nodes 1, 2 and 3.

For the particular case  of <math display="inline">Q=0</math> the following difference equation is found from Eq.(18)

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>-\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  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (20)
|}

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

We will solve now for the value of <math display="inline">\phi _2</math> for Dirichlet boundary conditions at nodes 1 and 3.

The two following cases are considered.

===Case A. ϕ₁=0, ϕ₃= ϕ₃<sup>p</sup>===

Application of the stencil of Eq.(20) gives the following value for <math display="inline">\phi _2</math>

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\phi _2 = {\phi _3^p \over (2+\beta ^1+\beta ^2 + {2w\over 3})} \left(\beta ^2 +1 -\gamma - {w\over 6} \right) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (21)
|}

A physically meaningful solution for any values of <math display="inline">\gamma </math> and <math display="inline">w</math> must necessarily satisfy <math display="inline">\phi _2\ge 0</math>. Clearly this requires

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\beta ^2\ge \gamma + {w\over 6} -1 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (22)
|}

===Case B. ϕ₁=ϕ<sup>p</sup>₁, ϕ₃= 0===

The stencil of Eq.(20) gives

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\phi _2 = { \phi _1^p \over (2+\beta ^1+\beta ^2 + {2w\over 3})} \left(1 + \beta ^1  +\gamma - {w\over 6} \right) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (23)
|}

Now <math display="inline">\phi _2 \ge 0</math> if

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\beta ^1 \ge {w\over 6} - \gamma - 1  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (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 <math display="inline">\beta ^e</math> is proposed (see Appendix A):

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\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 \le{0} \end{array} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (25)
|}

with

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\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]} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (26)
|}

and

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>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 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (27)
|}

where <math display="inline">\bar{(\cdot )}^e</math> denotes values computed at the element mid point.

It is clear from Eq.(26) that the computation of the stabilization parameter <math display="inline">\beta ^e</math> requires the knowledge of the sign of the numerical solution <math display="inline">\bar \phi </math> and that of the  first and second derivatives of <math display="inline">\bar \phi </math> 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 <math display="inline">\beta ^e</math> over the mesh yields a stabilized solution (for a uniform distribution of <math display="inline">\gamma </math> and <math display="inline">w</math>)

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\displaystyle{\beta ^e \ge \beta _c = {w\over 6} +\vert \gamma \vert -1} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (28)
|}

where <math display="inline">\beta _c </math> is the ''critical stabilization parameter''. The proof of Eq.(28) is given in the Appendix B.

The use of <math display="inline">\beta ^e =\beta _c</math> for all elements ''ensures a stabilized numerical solution'' for all ranges of <math display="inline">\gamma </math> and <math display="inline">w</math> although it can yield slight overdiffusive results in some cases.

The use  of <math display="inline">\beta ^e =\beta _c </math> 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)

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{w\over 6} +\vert \gamma \vert -1>0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (29)
|}

===4.3 Iterative solution scheme===

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


'''Step 1.'''
* Compute a first stabilized solution <math display="inline">\bar{\boldsymbol \phi }^1</math> using the critical value <math display="inline">\beta ^e = \beta _c</math> 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 <math display="inline">\bar{\boldsymbol \phi }^1</math> within each element. The second derivative field is obtained as follows

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\left({d^2 \bar \phi ^1 \over dx^2}\right)^e= {1\over l^e} \left[\left({d\hat \phi ^1 \over dx}\right)^e_2 - \left({d\hat \phi ^1 \over dx}\right)^e_1\right] </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (30)
|}

where <math display="inline">({\hat \cdot })_i^e</math> denotes averaged values of the first derivative at node <math display="inline">i</math> of element <math display="inline">e</math>. At the boundary nodes the constant value of the derivative of <math display="inline">\bar \phi </math> within the element is taken in Eq.(30); i.e. <math display="inline">(\hat{\cdot })_i^e = \left({d\bar \phi \over dx}\right)^{(e)} = {\bar \phi _2 - \bar \phi _1 \over l^e}</math>.

* ''Step 2.2''. Compute the enhanced stabilized solution <math display="inline">{\boldsymbol \phi }^2</math> using the element value of <math display="inline">\beta ^e</math> 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

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\phi = A e^{{ax\over l}} + Be^{{bx\over l}}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" |  (31a)
|}

where

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>a=\gamma + \delta \qquad , \qquad b=\gamma - \delta \qquad , \qquad \delta = (\gamma ^2 + w)^{1/2}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" |  (31b)
|}

Constants <math display="inline">A</math> and <math display="inline">B</math> depend on the boundary conditions of the problem. For the Dirichlet problem (<math display="inline">\phi = \phi ^p_1</math> at <math display="inline">x=0</math> and <math display="inline">\phi = \phi ^p_N</math> at <math display="inline">x=L</math>) we obtain

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>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}} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (32)
|}

Let us now write the stencil of Eq.(20) for three arbitrary nodes <math display="inline">i-1</math>, <math display="inline">i</math>, <math display="inline">i+1</math> and a uniform mesh with <math display="inline">l^e =l</math> and <math display="inline">\beta ^{i-1}=\beta ^i=\beta _i</math> as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>-\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 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (33)
|}

where <math display="inline">\beta _i</math> is the nodal value of the stabilization parameter or the <math display="inline">i</math>th node. Note that in Eq.(33) <math display="inline">\beta _i</math> has assumed to be constant within the tributary domain of node <math display="inline">i</math>.

The optimal nodal value of the stabilization parameter <math display="inline">\beta _i^{opt}</math> is found  from Eq.(33) as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\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 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (34)
|}

where <math display="inline">\phi _{i}^{ex}</math> denotes the exact analytical value of <math display="inline">\phi _i</math> as obtained from Eq.(31a) by

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\phi _{i}^{ex}=Ae^{{ax_i\over l}} + Be^{{bx_i\over l}} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (35)
|}

Substituting the exact nodal values of Eq.(35) into  Eq.(34)  gives the analytical expression of <math display="inline">\beta _i^{opt}</math> as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\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)\over \sinh \left(\displaystyle{a\over 2}\right)}\right)^2 e^{{b-a\over l}x_i}}\right]  + {w\over 6} \left[1+  { \displaystyle{3\over 2}  \left(1+ \displaystyle{A\over B}e^{{a-b\over l}x_i}\right)\over  \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 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (36)
|}

Eqs.(34) and (36) indicate that <math display="inline">\beta _i^{opt}</math> depends on the coordinate of node <math display="inline">i</math>. The advective (<math display="inline">w=0</math>) and absorptive (<math display="inline">\gamma =0</math>) limits of <math display="inline">\beta _i^{opt}</math> are given below.

In the examples presented  in the next section we have found more convenient to compute <math display="inline">\beta _i^{opt}</math> by Eq.(34) with the numerical values of <math display="inline">\phi _i^{ex}</math> (<math display="inline">i=1,2\cdots N+1</math>) computed from Eq.(35).

===4.5 <span id='lb-4.5'></span>Advective limit (w=0)===

If the absorption terms are zero, <math display="inline">w=0</math>, <math display="inline">a=2\gamma , b=0</math> and  from Eq.(36)

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\beta ^{opt}_i = \gamma \coth \gamma - 1  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (37)
|}

===4.6 <span id='lb-4.6'></span>Absorptive limit (γ=0)===

For <math display="inline">\gamma =0</math>, <math display="inline">a=\delta </math>, <math display="inline">b=-\delta </math>, <math display="inline">\delta = \sqrt w</math> and Eq.(36) gives

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\beta ^{opt}_i = {w\over 6} + {w\over 4\sinh ^2 \left({\sqrt w\over 2}\right)} -1  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (38)
|}

Note that both the advective and absorptive limits of <math display="inline">\beta _i^{opt}</math> 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 <span id='citeF-1'></span>[[#cite-1|[1]]]).<br />The expression of Eq.(38) coincides with that obtained by Felippa and Oñate <span id='citeF-15'></span>[[#cite-15|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 <math display="inline">h^e</math> is deduced from Eq.(13) as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>h^e={2\beta ^e k \left({d\bar \phi \over dx}\right)^e\over \left(u{d\bar \phi \over dx} - k {d^2\bar \phi \over dx^2} + s\bar \phi \right)^e} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (39)
|}

where <math display="inline">\beta ^e</math> is defined in Eq.(25).

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

The expression of <math display="inline">h^e</math> is needed in Eqs.(17) and (19) for the case of <math display="inline">Q\not =0</math>.

==6 <span id='lb-6'></span>ADVANTAGES OF USING β VERSUS h AS THE STABILIZATION PARAMETER==

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

Indeed, the characteristic length parameter <math display="inline">h</math> is still needed in the expression of <math display="inline">{f}^e</math> (see Eqs.(17) and (19)). This does  not pose any problem as the value (and sign) of <math display="inline">h</math> 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 <math display="inline">\bar \phi _1 =8</math> and <math display="inline">\bar \phi _9=3</math> and different values of <math display="inline">\gamma </math> and <math display="inline">w</math>.

Tables 1&#8211;8 show the numerical results obtained at the nodes for the standard  Galerkin method (<math display="inline">\beta =0</math>) and using the element (<math display="inline">\beta ^e</math>) and critical (<math display="inline">\beta _c</math>) values of <math display="inline">\beta </math>. Recall that the results for <math display="inline">\beta _c</math> and <math display="inline">\beta ^e</math> 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 <math display="inline">\beta _i^{opt}</math> yielding nodaly exact results.

We  note that the numbers in the tables for the values of <math display="inline">\beta </math> refer to ''element values'' for <math display="inline">\beta ^e</math> and <math display="inline">\beta _c</math> and to ''nodal values'' for the case of <math display="inline">\beta ^{opt}_i</math>.

Figures 1&#8211;6 show the distribution of <math display="inline">\bar \phi </math> along the domain for <math display="inline">\beta =0</math> (Galerkin), <math display="inline">\beta ^e</math> and <math display="inline">\beta _c</math>  for some of the cases studied, as well as the exact solution.

We note that for the case <math display="inline">\gamma =10, \beta =4</math> (Table 9 and Figure 5) small oscillations around the exact distribution of <math display="inline">\phi </math> are obtained for <math display="inline">\beta =\beta ^e</math>. 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 <math display="inline">\phi </math> 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 <math display="inline">\bar \phi _1=0</math> and <math display="inline">\bar \phi _8=1</math> and <math display="inline">\bar \phi _1=1</math> and <math display="inline">\bar \phi _8=0</math>, respectively. Both cases are solved for <math display="inline">\gamma =10</math> and <math display="inline">w=200</math>. The nodal results  are shown in Tables 12 and 13, while the distributions of <math display="inline">\bar \phi </math> are plotted in Figures 7 and 8.

In references <span id='citeF-12'></span>[[#cite-12|[12]]] it was shown that the SUPG, GLS and SGS methods with a single stabilization parameter fail to provide a stabilized solution for the case <math display="inline">\bar \phi _1=1</math>, <math display="inline">\bar \phi _8=0</math> with <math display="inline">\gamma =10</math> and <math display="inline">w=200</math>. 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.

<ol>

<li>The Galerkin solution (<math display="inline">\beta =0</math>) is unstable for the cases where the instability condition (Eq.(29))  is violated, as expected.  </li>

<li>The solution obtained with the critical value <math display="inline">\beta _c</math> is always stable, although it yields slightly overdiffusive results in some cases (see Tables 2, 6, 9 and 10).  </li>

<li>The results obtained with <math display="inline">\beta ^e</math>  are less diffusive than those obtained with <math display="inline">\beta _c</math> and quite accurate for all ranges of <math display="inline">\gamma </math> and <math display="inline">w</math>.  </li>

<li>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 <span id='citeF-12'></span>[[#cite-12|[12]]].  </li>

</ol>

==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 <span id='citeF-34'></span>[[#cite-34|[34]]].

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

<span id="eq-40"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>r(\phi ): =- {u}^T {\boldsymbol \nabla } \phi + {\boldsymbol \nabla }^T {D}{\boldsymbol \nabla }\phi - s\phi =0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (40)
|}

where for 2D problems

<span id="eq-41"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{u}=[u,v]^T\quad ,\quad {\boldsymbol \nabla }=\left[{\partial  \over \partial x},{\partial  \over \partial y}\right]^T\quad ,\quad {D} =k \left[\begin{matrix}1 &0\\ 0&1\\\end{matrix}\right] </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (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.([[#eq-40|[40]]]) is written

<span id="eq-42"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>r - \underline{{1\over 2} {h}^T {\boldsymbol \nabla }r}{1\over 2} {h}^T {\boldsymbol \nabla }r=0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (42)
|}

where <math display="inline">r</math> is the original differential equation as defined in Eq.([[#eq-40|[40]]]).

The Dirichlet and boundary conditions of the FIC formulation are

<span id="eq-43"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\phi - \phi ^p =0 \quad \hbox{on}\quad \Gamma _\phi   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (43)
|}

<span id="eq-44"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>- {u}^T {n}\phi + {n}^T {D} {\boldsymbol \nabla } \phi + q^p - \underline{{1\over 2} {h}^T {n}r}{1\over 2} {h}^T {n}r=0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (44)
|}

where <math display="inline">n</math> is the normal vector to the boundary where the normal flux is prescribed. As usual index <math display="inline">p</math> denotes the prescribed values.

In Eqs.([[#eq-42|42]]) and ([[#eq-44|44]]) <math display="inline">{h}=[h_x,h_y]^T</math> is the characteristic vector of the FIC formulation which components play the role of stabilization parameters. The underlined terms in Eqs.([[#eq-42|42]]) and ([[#eq-44|44]]) introduce the necessary stability in the numerical solution.

As shown in <span id='citeF-23'></span>[[#cite-23|23]] the components of the characteristic length '''h''' allow to obtain stabilized finite element solutions in the presence of strong gradients of <math display="inline">\phi </math> 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 <span id='citeF-20'></span>[[#cite-20|20]]. 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.'' <span id='citeF-23'></span>[[#cite-23|23]] 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 <math display="inline">\xi </math> and <math display="inline">\eta </math>. 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 <span id='citeF-23'></span>[[#cite-23|23]]. This technique is  extended here accounting for the absorption term.

Following the ideas described in <span id='citeF-23'></span>[[#cite-23|23]] the FIC equations can be written as

<span id="eq-45"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>- {u}^T {\boldsymbol \nabla } \phi + {\boldsymbol \nabla }^T ({D}+\bar {D}) {\boldsymbol \nabla }\phi - s\phi =0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (45)
|}

where the stabilization matrix <math display="inline">\bar {D}</math> is given by

<span id="eq-46"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\bar {D} = {T}^T \bar {D}' {T} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (46)
|}

<span id="eq-47"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\bar {D}' =\left[\begin{matrix}k_\xi{+} k_\xi ^s &\vdots &0 \\  \cdots & \cdots & \cdots \\  0& \vdots & k_\eta + k_\eta ^s\\\end{matrix}  \right] </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (47)
|}

and

<span id="eq-48"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{T}=\left[\begin{matrix}\cos \alpha & \sin \alpha \\ -\sin \alpha &\cos \alpha \\\end{matrix}\right]  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (48)
|}

where <math display="inline">\alpha </math> is the angle that the main principal curvature direction of the solution <math display="inline">\xi </math> forms with the global <math display="inline">x</math> axis. The different terms in Eq.(47&#8211;48) are

<span id="eq-49"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\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} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (49)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>k^s_\xi ={sl_\xi ^2\over 6}~,~ k^s_\eta ={sl_\eta ^2\over 6}</math>
|}
|}

where <math display="inline">u_\xi </math>, <math display="inline">u_\eta </math> are the velocity components along each of the principal curvature directions and <math display="inline">l_\xi </math> and <math display="inline">l_\eta </math> are the maximum values of the projection of the element sides along the principal directions <math display="inline">\xi </math> and <math display="inline">\eta </math>, respectively.

The values of <math display="inline"> k^s_\xi </math> and <math display="inline">k^s_\eta </math> in Eq.([[#eq-49|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 <span id='citeF-34'></span>[[#cite-34|34]].

The derivation of the finite element equations follows the standard Galerkin  procedure applied to Eq.([[#eq-43|43]]&#8211;[[#eq-45|45]]). For details see <span id='citeF-23'></span>[[#cite-23|23]].

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

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

<br/>

'''Step 2'''. Compute  the principal curvature directions of the solution at the element center. As mentioned above we have chosen the main principal direction <math display="inline">\xi </math> to be parallel to the gradient vector <math display="inline">{\boldsymbol \nabla }\phi </math> for simplicity. The direction <math display="inline">\eta </math> is chosen orthogonal to <math display="inline">\xi </math>.

Compute <math display="inline"> \bar {D}'</math>, <math display="inline">{T}</math> and <math display="inline">\bar{D}</math> from Eqs.([[#eq-46|46]]&#8211;[[#eq-48|48]]). Solve for <math display="inline">{}^2\bar{\boldsymbol \phi }</math>.<br/>

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 <span id='citeF-34'></span>[[#cite-34|34]].

===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 <math display="inline">8\times 8</math> four node bi-linear square elements is used in all cases.

The boundary conditions for the first and second examples are <math display="inline">\phi ^p =8</math> at the boundary <math display="inline">x=0</math> and <math display="inline">\phi =3</math> at <math display="inline">x=8</math> and zero normal flux at <math display="inline">y=0</math> and <math display="inline">y=8</math>. This reproduces the condition of the 1D examples solved in previous sections. The first example is analized for <math display="inline">{u} = [2,0]^T</math>, <math display="inline">k=1</math> and <math display="inline">s=1</math> giving <math display="inline">w=20</math>, <math display="inline">\gamma _x=1</math> and <math display="inline">\gamma _y=0</math> 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 <math display="inline">x=8</math> where <math display="inline">\phi </math> is prescribed to a value of <math display="inline">3</math>. The numerical results obtained with the standard (non stabilized) Galerkin solution (i.e. <math display="inline">\bar {D}=0</math>) 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 <math display="inline">\beta = \beta _c</math>.

The second example is similar to the first one with <math display="inline">{u} =[20,0]^T</math>, <math display="inline">k=1</math> and <math display="inline">s=20</math> giving <math display="inline">w=20</math>, <math display="inline">\gamma _x =10</math> and <math display="inline">\gamma _y =0</math>. 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 <math display="inline">\beta =\beta _c</math>.

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

The correct solution for this absorption dominated problem has two boundary layers next to the boundaries at <math display="inline">x=0</math> and <math display="inline">y=8</math>. 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.

The boundary conditions for the last example are <math display="inline">\phi =50</math> at <math display="inline">x=0</math> and <math display="inline">y =10</math> and <math display="inline">\phi =100</math> at <math display="inline">y=0</math> and <math display="inline">x=10</math>. The physical parameters are <math display="inline">{u} =[0,5]^T</math>, <math display="inline">k=1</math> and <math display="inline">s=12</math> giving <math display="inline">w=12</math>, <math display="inline">\gamma _x =0</math>  and <math display="inline">\gamma _y =2.5</math>.

The correct solution has two boundary layers next to the boundaries at <math display="inline">x=10</math> and <math display="inline">y=10</math>. The Galerkin solution gives oscillatory results next to these boundaries and also next to the boundary at <math display="inline">x=0</math> (Figure 12). All these oscillations disappear with the stabilized FIC formulation here proposed.

More accurate results for all these problems can be found by using the non linear form of the element  absorptive stabilization parameter <math display="inline">\beta ^e</math> similarly as for the 1D case. Details of the iterative scheme and examples of applications are given in <span id='citeF-34'></span>[[#cite-34|34]].

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

===BIBLIOGRAPHY===

<div id="cite-1"></div>
'''[[#citeF-1|[1]]]'''  Zienkiewicz O.C.  and  Taylor R.L. ''The Finite Element Method''. Volume 3. 5th Edition, Butterworth-Heinemann, 2001.

<div id="cite-2"></div>
'''[[#citeF-2|[2]]]''' Tezduyar T.E. and Park Y.J. Discontinuity-capturing finite element formulations for nonlinear convection-diffusion-reaction equations. ''Comput. Methods Appl. Mech. Engrg.'', '''59''', 307&#8211;325, 1986.

<div id="cite-3"></div>
'''[3]''' Tezduyar T.E., Park Y.J. and Deans H.A. Finite element procedures for time-dependent convection-diffusion-reaction systems. ''Int. J. Num. Meth. Fluids'', '''7''', 1013&#8211;1033, 1987.

<div id="cite-4"></div>
'''[[#citeF-4|[4]]]'''  Idelsohn S., Nigro N. Storti M. and Buscaglia G. A Petrov-Galerkin formulation for advection-reaction-diffusion problems. ''Comput. Methods Appl. Mech. Engrg.'', '''136''', 27&#8211;46, 1996.

<div id="cite-5"></div>
'''[[#citeF-5|[5]]]'''  Codina R. Comparison of some finite element methods for solving the diffusion-convection-reaction equation. ''Comput. Methods Appl. Mech. Engrg.'', '''156''', 186&#8211;210, 1998.

<div id="cite-6"></div>
'''[6]'''  Codina R. On stabilized finite element methods for linear systems of convection-diffusion-reaction equations. ''Comput. Meth. Appl. Mech. Engrg.'', '''188''', 61&#8211;82, 2000.

<div id="cite-7"></div>
'''[7]'''  Burman E. and Ern A. Nonlinear diffusion and discrete maxium principle for stabilized Galerkin approximations of the convection-diffusion-reaction equation. ''Comput. Meth. Appl. Mech. Engrg.'', '''191''', 3833&#8211;3855, 2002.

<div id="cite-8"></div>
'''[[#citeF-8|[8]]]'''  Harari I. and  Hughes T.J.R. Finite element methods for the Helmholtz equation in an exterior domain: model problems. ''Comput. Meth. Appl. Mech. Engrg.'', '''87''', 59&#8211;96, 1991.

<div id="cite-9"></div>
'''[9]'''   Harari I. and  Hughes T.J.R. Galerkin/least squares finite element methods for the reduced wave equation with non-reflecting boundary conditions in unbounded domains. ''Comput. Meth. Appl. Mech. Engrg.'', '''98''', 411&#8211;454, 1992.

<div id="cite-10"></div>
'''[10]'''   Harari I. and  Hughes T.J.R. Stabilized finite element method for steady advection-diffusion with production. ''Comput. Meth. Appl. Mech. Engrg.'', '''115''', 165&#8211;191, 1994.

<div id="cite-11"></div>
'''[11]'''  Harari I., Grosh K., Hughes T.J.R., Malhotra M.,  Pinsky P.M., Stewart J.R. and  Thompson L.L. Recent development in finite element methods for structural acoustics. ''Archives of Computational Mechanics in Engineering'', '''3, 2-3''', 131&#8211;309, 1996.

<div id="cite-12"></div>
'''[[#citeF-12|[12]]]'''  Hauke G. and Garcia-Olivares A. Variational subgrid scale formulations for the advection-diffusion-reaction equation. ''Comput. Methods Appl. Mech. Engrg.'' 2000; '''190''':6847&#8211;6865.

<div id="cite-13"></div>
'''[13]''' Hauke  G. A simple subgrid scale stabilized method for the advection-diffusion reaction equation. ''Comput. Methods Appl. Mech. Engrg.'' 2002; '''191''':2925&#8211;2947.

<div id="cite-14"></div>
'''[[#citeF-14|[14]]]'''  Brezzi F., Hauke G., Marin L.D. and Sangalli S. Link-cutting bubbles for the stabilization of convection-diffusion-reaction problems. ''Mathematical Models and Methods in Applied Sciences'', World Scientific Publishing Company, 2002.

<div id="cite-15"></div>
'''[[#citeF-15|[15]]]'''  Felippa C.A. and  Oñate E. Nodally exact Ritz discretization of 1D diffusion-absorption and Helmholtz equations by variational FIC and modified equation methods. Research Report No. PI 237, CIMNE, Barcelona 2004. Submitted to ''Comput. Mechanics''.

<div id="cite-16"></div>
'''[[#citeF-16|[16]]]'''  Idelsohn S.R.,  Heinrich J.C. and Oñate E. Petrov-Galerkin methods for the transient advective-diffusive equation with sharp gradients. ''Int. J. Num. Meth. Engng.'', '''39''', 1455&#8211;1473, 1996.

<div id="cite-17"></div>
'''[17]'''  Harari I. Stability of semidiscrete advection-diffusion in transient computation. Proceedings of ''6th World Congress on Computational Mechanics'', Beijing, Sept. 2004, Z.H. Yao, M.W. Yuan and W.X. Zhong (Eds.), Tsinghua Univ. Press-Springer.

<div id="cite-18"></div>
'''[18]''' Oñate E. Derivation of stabilized equations for  advective-diffusive transport and fluid flow problems.  ''Comput. Methods Appl. Mech. Engrg.'', '''151''' (1-2), 233&#8211;267, 1998.

<div id="cite-19"></div>
'''[19]'''  Oñate E. Possibilities of finite calculus in computational mechanics. ''Int. J. Num. Meth. Engng.'', '''60''', 255&#8211;281, 2004.

<div id="cite-20"></div>
'''[[#citeF-20|[20]]]''' Oñate E., García J. and Idelsohn S.R. Computation of the stabilization  parameter for the finite element solution of advective-diffusive problems. ''Int. J. Num. Meth. Fluids'', '''25''', 1385&#8211;1407, 1997.

<div id="cite-21"></div>
'''[21]''' Oñate E., Manzan M. Stabilization techniques for finite element analysis of convection diffusion problems. In ''Comput. Anal. of Heat Transfer'',  WIT Press, G. Comini and B. Sunden (Eds.), 2000.

<div id="cite-22"></div>
'''[22]''' Oñate E., Manzan M. A general procedure for deriving stabilized space-time finite element methods for advective-diffusive problems. ''Int. J. Num. Meth. Fluids'', '''31''', 203&#8211;221, 1999.

<div id="cite-23"></div>
'''[[#citeF-23|[23]]]''' Oñate E., Zárate F. and Idelsohn S.R. Finite element formulation for convection-diffusion problems with sharp gradients using finite calculus. ''Comput. Meth. Appl. Mech. Engng.'' Submitted Nov. 2004.

<div id="cite-24"></div>
'''[24]'''  Oñate E., Idelsohn S.R., Zienkiewicz O.C. and Taylor R.L. A finite point method in computational mechanics. Applications to convective transport and fluid flow. ''Int. J. Num. Meth. Engng.'', '''39''', 3839&#8211;3866, 1996.

<div id="cite-25"></div>
'''[25]'''  Oñate E., Idelsohn S.R. and Sacco C. A finite point method for incompressible flow problems. ''Computing and Visualization in Science'', '''17''', 127&#8211;154, 2000.

<div id="cite-26"></div>
'''[[#citeF-26|[26]]]'''  Oñate E. A stabilized finite element method for incompressible viscous flows using a finite increment calculus formulation. ''Comput. Methods Appl. Mech. Engrg.'', '''182''', (1&#8211;2), 355&#8211;370, 2000.

<div id="cite-27"></div>
'''[27]'''  Oñate E.,  García J. A finite element method for fluid-structure interaction with surface waves using a finite calculus formulation. in ''Comput. Methods Appl. Mech. Engrg.'', '''191''' (6-7), 635-660, 2001.

<div id="cite-28"></div>
'''[28]'''  Idelsohn S.R., Oñate E. and Del Pin F. A Lagrangian meshless finite element method applied to fluid-structure interaction problems. ''Computers and Structures'', '''81''', 655&#8211;671, 2003.

<div id="cite-29"></div>
'''[29]'''   Oñate E., Idelsohn S.R., Del Pin F. and  Aubry R. The particle finite element method. An overview (PFEM). ''Int. J. Comput. Methods'', '''1''' ('''2'''), 267&#8211;307, 2004.

<div id="cite-30"></div>
'''[30]'''   Oñate E., García J. and Idelsohn S.R. Ship Hydrodynamics. ''Encyclopedia of Computational Mechanics'', T. Hughes, R. de Borst and E. Stein (Eds.), Vol. 3, Chapter 18, 579&#8211;607, J. Wiley, 2004.

<div id="cite-31"></div>
'''[[#citeF-31|[31]]]'''  Oñate E., Perazzo F. and Miquel J. A finite point method for elasticity problems. ''Computers and Structures'', '''79''', 22&#8211;25, 2151&#8211;2163, 2001.

<div id="cite-32"></div>
'''[32]'''  Oñate E., Taylor R.L., Zienkiewicz O.C. and  Rojek J. A residual correction method based on finite calculus. ''Engineering Computations'', '''20''' (5/6), 629&#8211;658, 2003.

<div id="cite-33"></div>
'''[33]''' Oñate E., Rojek J., Taylor R.L. and Zienkiewicz O.C. Finite calculus formulation for analysis of incompressible solids using linear triangles and tetrahedra. ''Int. J. Num. Meth. Engng.'', '''59''' (11), 1473&#8211;1500, 2004.

<div id="cite-34"></div>
'''[[#citeF-34|[34]]]''' Oñate E., Miquel, J. and Zárate, F. Stabilized FIC/FEM formulation for multi-dimensional advection-diffusion-reaction problems. Research Report No. 269, CIMNE,  Barcelona, May 2005.

==APPENDIX A==

Let us consider the modified form of the convection-diffusion reaction equation with the stabilization term having the form of a diffusion of value <math display="inline">\beta k</math> i.e.

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>-u \phi' + k(1+\beta )\phi'' -s\phi=0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (50)
|}

where <math display="inline">\beta </math> is the stabilization parameter. For simplicity the source term <math display="inline">Q</math> has been taken equal to zero in Eq.(A.1).

From Eq.(A.1) we can obtain

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\beta = {s\phi \over k\phi'' } + {u\phi' \over k\phi'' } -1 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (51)
|}

Let us estimate now the value of <math display="inline">\beta </math> over a patch of finite elements. From the observation of the finite element matrices (Eq.(18)) we can write

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\phi \simeq \bar \phi =S_0 O\left({\bar \phi _{max}\over 6}\right)\quad ,\quad \phi' \simeq \bar \phi' =S_1 O\left({\bar \phi _{max}\over 2l}\right)\quad ,\quad \phi'' \simeq \bar \phi'' =S_2 O\left({\bar \phi _{max}\over l^2}\right) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (52)
|}

where <math display="inline">\phi _{max}</math> is the maximum value of the approximate finite element function <math display="inline">\bar \phi </math> over a patch of elements of equal size <math display="inline">l</math> and <math display="inline">S_0</math>, <math display="inline">S_1</math> and <math display="inline">S_2</math> denote the signs of the solution and of its first and second derivatives, respectively.

Substitution of Eq.(A.3) into (A.2) gives the elemental value of the stabilization parameter as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\beta ^e \simeq \left({S_0\over S_2} \right){w\over 6} + \left({S_1\over S_2} \right)\gamma -1 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (53)
|}

==APPENDIX B==

===STABILITY ANALYSIS OF THE DISCRETIZED FIC ADVECTION-DIFFUSION-ABSORPTION EQUATION FOR LINEAR ELEMENTS===

The discretized equation for the two linear element stencil of Eq.(20) can be written as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left({w\over 6} - 1 -\beta ^2-\gamma \right)\phi _3 + 2\left({w\over 3}+1+{\beta ^1+\beta ^2\over 2}\right)\phi _2 + \left({w\over 6} - 1 -\beta ^1+\gamma \right)\phi _1=0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" |  (B.1)
|}