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

Latest revision as of 12:15, 12 November 2018

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

Abstract

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

1 INTRODUCTION

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

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

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

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

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

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

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

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

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

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

(1)

(2)

(3)

where

(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 [20]. The underlined terms in Eqs.(1) and (2) emanate from these higher order expansions and they are essential to derive stabilized numerical schemes using whatever discretization method. Note that as the characteristic length parameter ${\textstyle h}$ tends to zero the FIC differential equations gradually recover the standard infinitessimal form giving in the limit (for ${\textstyle h=0}$ ) ${\textstyle r=0}$ in ${\textstyle x\in (0,L)}$ and ${\textstyle -u\phi +k{d\phi \over dx}+q^{p}=0\quad {\hbox{on }}\Gamma _{q}}$ .

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

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

3 FINITE ELEMENT FORMULATION

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

(5a)

with

(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

(6)

(7)

(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

(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

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

(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 \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

(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

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

(18)

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

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

(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

(22)

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

The stencil of Eq.(20) gives

(23)

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

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

(25)

with

(26)

and

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

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

(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

(31a)

where

(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

(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

(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

(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+{\displaystyle {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).

4.5 Advective limit (w=0)

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

(37)

4.6 Absorptive limit (γ=0)

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

(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 [1]).

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

5 VALUE OF THE ELEMENT CHARACTERISTIC LENGTH PARAMETER

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

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

Table. 1 $\phi _{1}=8,\phi _{9}=3,\gamma =0.1$ and $\omega =20$
$\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	–	–	–

Table. 2 $\phi _{1}=8,\phi _{9}=3,\gamma =1$ and $\omega =5$ (see Figure 1)
$\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	–	–	–

Table. 3 $\phi _{1}=8,\phi _{9}=3,\gamma =1$ and $\omega =20$ (see Figure 2)
$\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	–	–	–

Table. 4 $\phi _{1}=8,\phi _{9}=3,\gamma =1$ and $\omega =120$ (see Figure 3)
$\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	–	–	–

Table. 5 $\phi _{1}=8,\phi _{9}=3,\gamma =2$ and $\omega =0.1$
$\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	–	–	–

Table. 6 $\phi _{1}=8,\phi _{9}=3,\gamma =2$ and $\omega =2$ (see Figure 4)
$\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

Table. 7 $\phi _{1}=8,\phi _{9}=3,\gamma =10$ and $\omega =1$
$\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	–	–	–

Table. 8 $\phi _{1}=8,\phi _{9}=3,\gamma =10$ and $\omega =1$
$\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	–	–	–

Table. 9 $\phi _{1}=8,\phi _{9}=3,\gamma =10$ and $\omega =4$ (see Figure 5)
$\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	–	–

Table. 10 $\phi _{1}=8,\phi _{9}=3,\gamma =10$ and $\omega =20$ (see Figure 6)
$\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	–	–	–

Table. 11 $\phi _{1}=8,\phi _{9}=3,\gamma =10$ and $\omega =120$
$\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:

and

. FIC results for a mesh of 8 linear elements obtained for

(Galerkin),

and

. Comparison with the analytical solution (see Table 2)

Figure 2:

and

. FIC results for a mesh of 8 linear elements obtained for

(Galerkin),

and

. Comparison with the analytical solution (see Table 3)

Figure 3:

and

. FIC results for a mesh of 8 linear elements obtained for

(Galerkin),

and

. Comparison with the analytical solution (see Table 4)

Figure 4:

and

. FIC results for a mesh of 8 linear elements obtained for

(Galerkin),

and

. Comparison with the analytical solution (see Table 6)

Figure 5:

and

. FIC results for a mesh of 8 linear elements obtained for

(Galerkin),

and

. Comparison with the analytical solution (see Table 9)

Figure 6:

and

. FIC results for a mesh of 8 linear elements obtained for

(Galerkin),

and

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

Table. 12 $\phi _{1}=0,\phi _{9}=1,\gamma =10$ and $\omega =200$ (see Figure 7)
$\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	–	–	–

Table. 13 $\phi _{1}=0,\phi _{9}=0,\gamma =10$ and $\omega =200$ (see Figure 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	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:

and

. FIC results for a mesh of 8 linear elements obtained for

(Galerkin),

and

. Comparison with the analytical solution (see Table 12)

Figure 8:

and

. FIC results for a mesh of 8 linear elements obtained for

(Galerkin),

and

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

The Galerkin solution ( ${\textstyle \beta =0}$ ) is unstable for the cases where the instability condition (Eq.(29)) is violated, as expected.
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).
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}$ .
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 [34].

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

(40)

where for 2D problems

(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

(42)

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

The Dirichlet and boundary conditions of the FIC formulation are

(43)

(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 [23] 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 [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. [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 ${\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 [23]. This technique is extended here accounting for the absorption term.

Following the ideas described in [23] the FIC equations can be written as

(45)

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

(46)

(47)

and

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

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.(43–45). 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.(46–48). 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 [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 ${\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.

at

,

at

,

at

and

.

,

and

. Galerkin and FIC solutions obtained with a mesh of

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.

at

,

at

,

at

and

.

,

and

. Galerkin and FIC solutions obtained with a mesh of

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.

at

,

at

,

at

and

.

,

and

. Galerkin and FIC solutions obtained with a mesh of

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.

at

,

at

,

at

and

.

,

and

. Galerkin and FIC solutions obtained with a mesh of

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

REFERENCES

[1] Zienkiewicz O.C. and Taylor R.L. The Finite Element Method. Volume 3. 5th Edition, Butterworth-Heinemann, 2001.

[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–325, 1986.

[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–1033, 1987.

[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–46, 1996.

[5] Codina R. Comparison of some finite element methods for solving the diffusion-convection-reaction equation. Comput. Methods Appl. Mech. Engrg., 156, 186–210, 1998.

[6] Codina R. On stabilized finite element methods for linear systems of convection-diffusion-reaction equations. Comput. Meth. Appl. Mech. Engrg., 188, 61–82, 2000.

[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–3855, 2002.

[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–96, 1991.

[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–454, 1992.

[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–191, 1994.

[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–309, 1996.

[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–6865.

[13] Hauke G. A simple subgrid scale stabilized method for the advection-diffusion reaction equation. Comput. Methods Appl. Mech. Engrg. 2002; 191:2925–2947.

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

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

[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–1473, 1996.

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

[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–267, 1998.

[19] Oñate E. Possibilities of finite calculus in computational mechanics. Int. J. Num. Meth. Engng., 60, 255–281, 2004.

[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–1407, 1997.

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

[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–221, 1999.

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

[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–3866, 1996.

[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–154, 2000.

[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–2), 355–370, 2000.

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

[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–671, 2003.

[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–307, 2004.

[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–607, J. Wiley, 2004.

[31] Oñate E., Perazzo F. and Miquel J. A finite point method for elasticity problems. Computers and Structures, 79, 22–25, 2151–2163, 2001.

[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–658, 2003.

[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–1500, 2004.

[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 ${\textstyle \beta k}$ i.e.

(A.1)

where ${\textstyle \beta }$ is the stabilization parameter. For simplicity the source term ${\textstyle Q}$ has been taken equal to zero in Eq.(A.1).

From Eq.(A.1) we can obtain

(A.2)

Let us estimate now the value of ${\textstyle \beta }$ over a patch of finite elements. From the observation of the finite element matrices (Eq.(18)) we can write

(A.3)

where ${\textstyle \phi _{max}}$ is the maximum value of the approximate finite element function ${\textstyle {\bar {\phi }}}$ over a patch of elements of equal size ${\textstyle l}$ and ${\textstyle S_{0}}$ , ${\textstyle S_{1}}$ and ${\textstyle S_{2}}$ 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

(A.4)

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

(B.1)

The characteristic polinomial is

(B.2)

and the characteristic roots are

l_{1,2}={-\left({w \over 3}+1+{\beta ^{1}+\beta ^{2} \over 2}\right)\pm {\sqrt {\left({w \over 3}+1+{\beta ^{1}+\beta ^{2} \over 2}\right)^{2}-\left({w \over 6}-1-\beta ^{2}-\gamma \right)\left({w \over 6}-1-\beta ^{1}+\gamma \right)}}/\left({w \over 6}-1-\beta ^{2}-\gamma \right)}

(B.3)

A non oscilatory (stable) numerical solution is found if ${\textstyle l_{1,2}\geq 0}$ . This is fulfilled if

(B.4)

(B.5)

The inequality (B.4) implies

(B.6)

or

(B.7)

The second inequality (B.5) implies

(B.8)

or

(B.9)

Hence the stability limit of ${\textstyle \beta ^{e}}$ is

(B.10)

where ${\textstyle \beta _{c}}$ is the critical stabilization parameter.

Indeed a stabilized solution will be obtained if ${\textstyle \beta ^{e}\geq \beta _{c}}$ .

The condition for the onset of instability for the original (non-stabilized) problem is obtained by making ${\textstyle \beta =0}$ in Eq.(B.6) and (B.8) and reversing the sing of the inequalities, i.e. an oscilatory solution will be found if one of the two conditions is fullfiled

(B.11)

Hence, the stability condition is found as

(B.12)

Latest revision as of 12:15, 12 November 2018

Abstract

1 INTRODUCTION

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

3 FINITE ELEMENT FORMULATION

4 COMPUTATION OF THE STABILIZATION PARAMETER FOR LINEAR ELEMENTS

4.1 Computation of the element stabilization parameter β^e

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

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

4.2 Critical stabilization parameter and unstability conditions

4.3 Iterative solution scheme

4.4 Optimal value of the stabilization parameter

4.5 Advective limit (w=0)

4.6 Absorptive limit (γ=0)

5 VALUE OF THE ELEMENT CHARACTERISTIC LENGTH PARAMETER

6 ADVANTAGES OF USING β VERSUS h AS THE STABILIZATION PARAMETER

7 1D NUMERICAL EXAMPLES

8 EXTENSION TO MULTI-DIMENSIONAL PROBLEMS

8.1 2D Examples

9 CONCLUSIONS

ACKNOWLEDGEMENTS

REFERENCES

APPENDIX A

APPENDIX B. STABILITY ANALYSIS OF THE DISCRETIZED FIC ADVECTION-DIFFUSION-ABSORPTION EQUATION FOR LINEAR ELEMENTS

Document information

Document Score

Share this document

claim authorship

@@ Line 1: / Line 1: @@
-<!-- metadata commented in wiki content
+Published in ''Comput. Meth. Appl. Mech. Engng.'', Vol. 195, pp. 3926–3946, 2006<br />
-==STABILIZED FORMULATION FOR THE ADVECTION-DIFFUSION-ABSORPTION EQUATION USING FINITE CALCULUS AND LINEAR FINITE ELEMENTS==
+doi: 10.1016/j.cma.2005.07.020
-'''Eugenio Oñate<math>^1</math>, Juan Miquel<math>^1</math>  and Guillermo Hauke<math>^2</math>'''
-{|
-|-
-|<math>^1</math> International Center for Numerical Methods in Engineering (CIMNE)
-|-
-| Universidad Politécnica de Cataluña
-|-
-| Edificio C1, Gran Capitán s/n, 08034, Barcelona, Spain
-|-
-| e-mail: [mailto:onate@cimne.upc.edu onate@cimne.upc.edu]
-|}
-{|
-|-
-|<math>^2</math> Departamento de Mecánica de Fluidos
-|-
-| Centro Politécnico Superior
-|-
-| Maria Luna 3, 50015 Zaragoza, Spain
-|}
--->
 ==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.
-<br/><br/>
 ==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.
+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.
+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-7]]], Galerkin Least Square <span id='citeF-5'></span>[[#cite-5|[5-11]]], Subgrid Scale (SGS) <span id='citeF-5'></span>[[#cite-5|[5,6,12,13]]] 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-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 <span id='citeF-4'></span>[[#cite-4|[4,10,14]]]. As reported in <span id='citeF-12'></span>[[#cite-12|[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 <span id='citeF-16'></span>[[#cite-16|16]].
+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,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 <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 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 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|[18-24]]], diffusion-absorption <span id='citeF-15'></span>[[#cite-15|[15]]], incompressible fluid flow <span id='citeF-26'></span>[[#cite-26|[24-30]]] and standard and incompressible solid mechanics <span id='citeF-31'></span>[[#cite-31|[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 <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.
+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.
@@ Line 97: / Line 73: @@
 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>.
+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]].
+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|[18]]].
-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]].
+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|[18-31]]].
 ==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.
+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;"
@@ Line 468: / Line 444: @@
 step
-'''Step step.'''step
+'''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.
-* 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.step.''step
+''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
-* 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;"
@@ Line 490: / Line 464: @@
 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>.
-* 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).
+''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.
@@ Line 616: / Line 590: @@
 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]]).
+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]]]).
-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.
+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==
@@ Line 652: / Line 626: @@
 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.
+Tables 1&#8211;11 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 ): =- {\boldsymbol u}^T {\boldsymbol \nabla } \phi + {\boldsymbol \nabla }^T {\boldsymbol 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>{\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] </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} {\boldsymbol h}^T {\boldsymbol \nabla }r}{1\over 2} {\boldsymbol 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>- {\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 </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">{\boldsymbol 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>- {\boldsymbol u}^T {\boldsymbol \nabla } \phi + {\boldsymbol \nabla }^T ({\boldsymbol D}+\bar {\boldsymbol 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 {\boldsymbol 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 {\boldsymbol D} = {\boldsymbol T}^T \bar {\boldsymbol D}' {\boldsymbol 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 {\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] </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>{\boldsymbol 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 {\boldsymbol D}={\boldsymbol 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 {\boldsymbol D}'</math>, <math display="inline">{\boldsymbol T}</math> and <math display="inline">\bar{\boldsymbol 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">{\boldsymbol 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 {\boldsymbol 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">{\boldsymbol 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">{\boldsymbol 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">{\boldsymbol 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)
-|}
-The characteristic polinomial is
-{| 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)l^2 + 2\left({w\over 3}+1+{\beta ^1+\beta ^2\over 2}\right)l + \left({w\over 6} - 1 -\beta ^1+\gamma \right)=0</math>
-|}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (B.2)
-|}
-and  the characteristic roots are
-{| class="formulaSCP" style="width: 100%; text-align: left;"
-|-
-|
-{| style="text-align: left; margin:auto;width: 100%;"
-|-
-| style="text-align: center;" | <math>l_{1,2} = {-\left({w\over 3} +1+{\beta ^1+\beta ^2\over 2}\right)\pm \sqrt{\left({w\over 3} +1+{\beta ^1+\beta ^2\over 2}\right)^2 - \left({w\over 6} - 1 -\beta ^2-\gamma \right)\left({w\over 6} - 1 -\beta ^1+\gamma \right)}/\left({w\over 6} - 1 -\beta ^2-\gamma \right)}</math>
-|}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (B.3)
-|}
-A non oscilatory (stable) numerical solution is found if <math display="inline">l_{1,2} \ge 0</math>. This is fulfilled if
-{| class="formulaSCP" style="width: 100%; text-align: left;"
-|-
-|
-{| style="text-align: left; margin:auto;width: 100%;"
-|-
-| style="text-align: center;" | <math>l_1 +l_2 = - {w/3 + 1+ {\beta ^1+\beta ^2\over 2}\over w/6 -1 -\beta ^2-\gamma }>0</math>
-|}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (B.4)
-|}
-{| class="formulaSCP" style="width: 100%; text-align: left;"
-|-
-|
-{| style="text-align: left; margin:auto;width: 100%;"
-|-
-| style="text-align: center;" | <math>l_1 l_2 = {w/6 -1 -\beta ^1 +\gamma \over w/6 -1 -\beta ^2-\gamma }>0</math>
-|}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (B.5)
-|}
-The inequality (B.4) implies
-{| class="formulaSCP" style="width: 100%; text-align: left;"
-|-
-|
-{| style="text-align: left; margin:auto;width: 100%;"
-|-
-| style="text-align: center;" | <math>{w\over 6} -1 -\beta ^2-\gamma \le 0</math>
-|}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (B.6)
-|}
-or
-{| class="formulaSCP" style="width: 100%; text-align: left;"
-|-
-|
-{| style="text-align: left; margin:auto;width: 100%;"
-|-
-| style="text-align: center;" | <math>\beta ^2\ge {w\over 6} - \gamma -1</math>
-|}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (B.7)
-|}
-The second inequality (B.5) implies
-{| class="formulaSCP" style="width: 100%; text-align: left;"
-|-
-|
-{| style="text-align: left; margin:auto;width: 100%;"
-|-
-| style="text-align: center;" | <math>{w\over 6} -1 -\beta ^1+\gamma \le 0</math>
-|}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (B.8)
-|}
-or
-{| 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;" | (B.9)
-|}
-Hence the stability limit of <math display="inline">\beta ^e</math>  is
-{| class="formulaSCP" style="width: 100%; text-align: left;"
-|-
-|
-{| style="text-align: left; margin:auto;width: 100%;"
-|-
-| style="text-align: center;" | <math>\displaystyle{\beta ^e=\beta _c={w\over 6}+\vert \gamma \vert -1}</math>
-|}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (B.10)
-|}
-where <math display="inline">\beta _c</math> is the critical stabilization parameter.
-Indeed a stabilized solution will be obtained if <math display="inline">\beta ^e \ge \beta _c</math>.
-The condition for the onset of instability for the original (non-stabilized) problem is obtained by making <math display="inline">\beta =0</math> in Eq.(B.6) and (B.8) and reversing the sing of the inequalities, i.e. an ''oscilatory solution'' will be found if one of the two conditions is fullfiled
-{| 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{w\over 6}+\gamma -1\ge 0\\ \\ \displaystyle{w\over 6}-\gamma -1\ge 0 \end{array}</math>
-|}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (B.11)
-|}
-Hence, the stability condition  is found as
-{| 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\ge{0}</math>
-|}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (B.12)
-|}
@@ Line 1,968: / Line 1,409: @@
 | style="border-left: 2px solid;border-right: 2px solid;" | &#8211;
 | style="border-left: 2px solid;border-right: 2px solid;" | &#8211;
 |}
@@ Line 2,164: / Line 1,604: @@
 | style="border-left: 2px solid;border-right: 2px solid;" | &#8211;
 | style="border-left: 2px solid;border-right: 2px solid;" | &#8211;
 |}
 {|  class="floating_tableSCP wikitable" style="text-align: right; margin: 1em auto;min-width:50%;"
@@ Line 2,264: / Line 1,702: @@
 |}
+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.
+<div id='img-1'></div>
+{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
+|-
+|[[Image:Draft_Samper_267560997-test-Fig1.png|400px|ϕ₁=8, ϕ₉=3, γ=1 and ω=5. FIC results for a mesh of 8 linear elements obtained for β=0 (Galerkin), β<sup>e</sup> and β<sub>c</sub>. Comparison with the analytical solution (see Table 2)]]
+|- style="text-align: center; font-size: 75%;"
+| colspan="1" | '''Figure 1:''' <math>\phi _1=8, \phi _9=3, \gamma =1</math> and <math>\omega =5</math>. FIC results for a mesh of 8 linear elements obtained for <math>\beta =0</math> (Galerkin), <math>\beta ^e</math> and <math>\beta _c</math>. Comparison with the analytical solution (see Table 2)
+|}
+<div id='img-2'></div>
+{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
+|-
+|[[Image:Draft_Samper_267560997-test-Fig2.png|400px|ϕ₁=8, ϕ₉=3, γ=1 and ω=20. FIC results for a mesh of 8 linear elements obtained for β=0 (Galerkin), β<sup>e</sup> and β<sub>c</sub>. Comparison with the analytical solution (see Table 3)]]
+|- style="text-align: center; font-size: 75%;"
+| colspan="1" | '''Figure 2:''' <math>\phi _1=8, \phi _9=3, \gamma =1</math> and <math>\omega =20</math>. FIC results for a mesh of 8 linear elements obtained for <math>\beta =0</math> (Galerkin), <math>\beta ^e</math> and <math>\beta _c</math>. Comparison with the analytical solution (see Table 3)
+|}
+<div id='img-3'></div>
+{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
+|-
+|[[Image:Draft_Samper_267560997-test-Fig3.png|400px|ϕ₁=8, ϕ₉=3, γ=1 and ω=120. FIC results for a mesh of 8 linear elements obtained for β=0 (Galerkin), β<sup>e</sup> and β<sub>c</sub>. Comparison with the analytical solution (see Table 4)]]
+|- style="text-align: center; font-size: 75%;"
+| colspan="1" | '''Figure 3:''' <math>\phi _1=8, \phi _9=3, \gamma =1</math> and <math>\omega =120</math>. FIC results for a mesh of 8 linear elements obtained for <math>\beta =0</math> (Galerkin), <math>\beta ^e</math> and <math>\beta _c</math>. Comparison with the analytical solution (see Table 4)
+|}
+<div id='img-4'></div>
+{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
+|-
+|[[Image:Draft_Samper_267560997-test-Fig4.png|400px|ϕ₁=8, ϕ₉=3, γ=2 and ω=2. FIC results for a mesh of 8 linear elements obtained for β=0 (Galerkin), β<sup>e</sup> and β<sub>c</sub>. Comparison with the analytical solution (see Table 6)]]
+|- style="text-align: center; font-size: 75%;"
+| colspan="1" | '''Figure 4:''' <math>\phi _1=8, \phi _9=3, \gamma =2</math> and <math>\omega =2</math>. FIC results for a mesh of 8 linear elements obtained for <math>\beta =0</math> (Galerkin), <math>\beta ^e</math> and <math>\beta _c</math>. Comparison with the analytical solution (see Table 6)
+|}
+<div id='img-5'></div>
+{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
+|-
+|[[Image:Draft_Samper_267560997-test-Fig5.png|400px|ϕ₁=8, ϕ₉=3, γ=10 and ω=4. FIC results for a mesh of 8 linear elements obtained for β=0 (Galerkin), β<sup>e</sup> and β<sub>c</sub>. Comparison with the analytical solution (see Table 9)]]
+|- style="text-align: center; font-size: 75%;"
+| colspan="1" | '''Figure 5:''' <math>\phi _1=8, \phi _9=3, \gamma =10</math> and <math>\omega =4</math>. FIC results for a mesh of 8 linear elements obtained for <math>\beta =0</math> (Galerkin), <math>\beta ^e</math> and <math>\beta _c</math>. Comparison with the analytical solution (see Table 9)
+|}
+<div id='img-6'></div>
+{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
+|-
+|[[Image:Draft_Samper_267560997-test-Fig6.png|400px|ϕ₁=8, ϕ₉=3, γ=10 and ω=20. FIC results for a mesh of 8 linear elements obtained for β=0 (Galerkin), β<sup>e</sup> and β<sub>c</sub>. Comparison with the analytical solution (see Table 10)]]
+|- style="text-align: center; font-size: 75%;"
+| colspan="1" | '''Figure 6:''' <math>\phi _1=8, \phi _9=3, \gamma =10</math> and <math>\omega =20</math>. FIC results for a mesh of 8 linear elements obtained for <math>\beta =0</math> (Galerkin), <math>\beta ^e</math> and <math>\beta _c</math>. Comparison with the analytical solution (see Table 10)
+|}
+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.
@@ Line 2,458: / Line 1,955: @@
 | style="border-left: 2px solid;border-right: 2px solid;" | &#8211;
 | style="border-left: 2px solid;border-right: 2px solid;" | &#8211;
 |}
-<div id='img-1'></div>
+<div id='img-7'></div>
 {| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 |-
-|[[Image:Draft_Samper_267560997-test-Fig1.png|600px|ϕ₁=8, ϕ₉=3, γ=1 and ω=5. FIC results for a mesh of 8 linear elements obtained for β=0 (Galerkin), β<sup>e</sup> and β<sub>c</sub>. Comparison with the analytical solution (see Table 2)]]
+|[[Image:Draft_Samper_267560997-test-Fig7.png|400px|ϕ₁=8, ϕ₉=3, γ=10 and ω=200. FIC results for a mesh of 8 linear elements obtained for β=0 (Galerkin), β<sup>e</sup> and β<sub>c</sub>. Comparison with the analytical solution (see Table 12)]]
 |- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 1:''' <math>\phi _1=8, \phi _9=3, \gamma =1</math> and <math>\omega =5</math>. FIC results for a mesh of 8 linear elements obtained for <math>\beta =0</math> (Galerkin), <math>\beta ^e</math> and <math>\beta _c</math>. Comparison with the analytical solution (see Table 2)
+| colspan="1" | '''Figure 7:''' <math>\phi _1=8, \phi _9=3, \gamma =10</math> and <math>\omega =200</math>. FIC results for a mesh of 8 linear elements obtained for <math>\beta =0</math> (Galerkin), <math>\beta ^e</math> and <math>\beta _c</math>. Comparison with the analytical solution (see Table 12)
 |}
-<div id='img-2'></div>
+<div id='img-8'></div>
 {| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 |-
-|[[Image:Draft_Samper_267560997-test-Fig2.png|600px|ϕ₁=8, ϕ₉=3, γ=1 and ω=20. FIC results for a mesh of 8 linear elements obtained for β=0 (Galerkin), β<sup>e</sup> and β<sub>c</sub>. Comparison with the analytical solution (see Table 3)]]
+|[[Image:Draft_Samper_267560997-test-Fig8.png|400px|ϕ₁=8, ϕ₉=3, γ=1 and ω=20. FIC results for a mesh of 8 linear elements obtained for β=0 (Galerkin), β<sup>e</sup> and β<sub>c</sub>. Comparison with the analytical solution (see Table 13)]]
 |- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 2:''' <math>\phi _1=8, \phi _9=3, \gamma =1</math> and <math>\omega =20</math>. FIC results for a mesh of 8 linear elements obtained for <math>\beta =0</math> (Galerkin), <math>\beta ^e</math> and <math>\beta _c</math>. Comparison with the analytical solution (see Table 3)
+| colspan="1" | '''Figure 8:''' <math>\phi _1=8, \phi _9=3, \gamma =1</math> and <math>\omega =20</math>. FIC results for a mesh of 8 linear elements obtained for <math>\beta =0</math> (Galerkin), <math>\beta ^e</math> and <math>\beta _c</math>. Comparison with the analytical solution (see Table 13)
 |}
-<div id='img-3'></div>
-{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
+In references <span id='citeF-12'></span>[[#cite-12|[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 <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,13]]].  </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;"
 |-
-|[[Image:Draft_Samper_267560997-test-Fig3.png|600px|ϕ₁=8, ϕ₉=3, γ=1 and ω=120. FIC results for a mesh of 8 linear elements obtained for β=0 (Galerkin), β<sup>e</sup> and β<sub>c</sub>. Comparison with the analytical solution (see Table 4)]]
+|
-|- style="text-align: center; font-size: 75%;"
+{| style="text-align: left; margin:auto;width: 100%;"
-| colspan="1" | '''Figure 3:''' <math>\phi _1=8, \phi _9=3, \gamma =1</math> and <math>\omega =120</math>. FIC results for a mesh of 8 linear elements obtained for <math>\beta =0</math> (Galerkin), <math>\beta ^e</math> and <math>\beta _c</math>. Comparison with the analytical solution (see Table 4)
+|-
+| style="text-align: center;" | <math>r(\phi ): =- {\boldsymbol u}^T {\boldsymbol \nabla } \phi + {\boldsymbol \nabla }^T {\boldsymbol D}{\boldsymbol \nabla }\phi - s\phi =0 </math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (40)
 |}
-<div id='img-4'></div>
+where for 2D problems
-{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
+<span id="eq-41"></span>
+{| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
-|[[Image:Draft_Samper_267560997-test-Fig4.png|600px|ϕ₁=8, ϕ₉=3, γ=2 and ω=2. FIC results for a mesh of 8 linear elements obtained for β=0 (Galerkin), β<sup>e</sup> and β<sub>c</sub>. Comparison with the analytical solution (see Table 6)]]
+|
-|- style="text-align: center; font-size: 75%;"
+{| style="text-align: left; margin:auto;width: 100%;"
-| colspan="1" | '''Figure 4:''' <math>\phi _1=8, \phi _9=3, \gamma =2</math> and <math>\omega =2</math>. FIC results for a mesh of 8 linear elements obtained for <math>\beta =0</math> (Galerkin), <math>\beta ^e</math> and <math>\beta _c</math>. Comparison with the analytical solution (see Table 6)
+|-
+| style="text-align: center;" | <math>{\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] </math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (41)
 |}
-<div id='img-5'></div>
+are respectively the velocity vector, the gradient vector and the diffusivity matrix, respectively. For simplicity we have assumed an isotropic diffusion matrix.
-{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
+The FIC form of Eq.([[#eq-40|40]]) is written
+<span id="eq-42"></span>
+{| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
-|[[Image:Draft_Samper_267560997-test-Fig5.png|600px|ϕ₁=8, ϕ₉=3, γ=10 and ω=4. FIC results for a mesh of 8 linear elements obtained for β=0 (Galerkin), β<sup>e</sup> and β<sub>c</sub>. Comparison with the analytical solution (see Table 9)]]
+|
-|- style="text-align: center; font-size: 75%;"
+{| style="text-align: left; margin:auto;width: 100%;"
-| colspan="1" | '''Figure 5:''' <math>\phi _1=8, \phi _9=3, \gamma =10</math> and <math>\omega =4</math>. FIC results for a mesh of 8 linear elements obtained for <math>\beta =0</math> (Galerkin), <math>\beta ^e</math> and <math>\beta _c</math>. Comparison with the analytical solution (see Table 9)
+|-
+| style="text-align: center;" | <math>r - \underline{{1\over 2} {\boldsymbol h}^T {\boldsymbol \nabla }r}{1\over 2} {\boldsymbol h}^T {\boldsymbol \nabla }r=0 </math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (42)
 |}
-<div id='img-6'></div>
+where <math display="inline">r</math> is the original differential equation as defined in Eq.([[#eq-40|40]]).
-{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
+The Dirichlet and boundary conditions of the FIC formulation are
+<span id="eq-43"></span>
+{| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
-|[[Image:Draft_Samper_267560997-test-Fig6.png|600px|ϕ₁=8, ϕ₉=3, γ=10 and ω=20. FIC results for a mesh of 8 linear elements obtained for β=0 (Galerkin), β<sup>e</sup> and β<sub>c</sub>. Comparison with the analytical solution (see Table 10)]]
+|
-|- style="text-align: center; font-size: 75%;"
+{| style="text-align: left; margin:auto;width: 100%;"
-| colspan="1" | '''Figure 6:''' <math>\phi _1=8, \phi _9=3, \gamma =10</math> and <math>\omega =20</math>. FIC results for a mesh of 8 linear elements obtained for <math>\beta =0</math> (Galerkin), <math>\beta ^e</math> and <math>\beta _c</math>. Comparison with the analytical solution (see Table 10)
+|-
+| 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)
 |}
-<div id='img-7'></div>
+<span id="eq-44"></span>
-{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
+{| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
-|[[Image:Draft_Samper_267560997-test-Fig7.png|600px|ϕ₁=8, ϕ₉=3, γ=10 and ω=200. FIC results for a mesh of 8 linear elements obtained for β=0 (Galerkin), β<sup>e</sup> and β<sub>c</sub>. Comparison with the analytical solution (see Table 12)]]
+|
-|- style="text-align: center; font-size: 75%;"
+{| style="text-align: left; margin:auto;width: 100%;"
-| colspan="1" | '''Figure 7:''' <math>\phi _1=8, \phi _9=3, \gamma =10</math> and <math>\omega =200</math>. FIC results for a mesh of 8 linear elements obtained for <math>\beta =0</math> (Galerkin), <math>\beta ^e</math> and <math>\beta _c</math>. Comparison with the analytical solution (see Table 12)
+|-
+| style="text-align: center;" | <math>- {\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 </math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (44)
 |}
-<div id='img-8'></div>
+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.
-{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
+In Eqs.([[#eq-42|42]]) and ([[#eq-44|44]]) <math display="inline">{\boldsymbol 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;"
 |-
-|[[Image:Draft_Samper_267560997-test-Fig8.png|600px|ϕ₁=8, ϕ₉=3, γ=1 and ω=20. FIC results for a mesh of 8 linear elements obtained for β=0 (Galerkin), β<sup>e</sup> and β<sub>c</sub>. Comparison with the analytical solution (see Table 13)]]
+|
-|- style="text-align: center; font-size: 75%;"
+{| style="text-align: left; margin:auto;width: 100%;"
-| colspan="1" | '''Figure 8:''' <math>\phi _1=8, \phi _9=3, \gamma =1</math> and <math>\omega =20</math>. FIC results for a mesh of 8 linear elements obtained for <math>\beta =0</math> (Galerkin), <math>\beta ^e</math> and <math>\beta _c</math>. Comparison with the analytical solution (see Table 13)
+|-
+| style="text-align: center;" | <math>- {\boldsymbol u}^T {\boldsymbol \nabla } \phi + {\boldsymbol \nabla }^T ({\boldsymbol D}+\bar {\boldsymbol 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 {\boldsymbol 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 {\boldsymbol D} = {\boldsymbol T}^T \bar {\boldsymbol D}' {\boldsymbol 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 {\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] </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>{\boldsymbol 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|[23,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,34]]].
+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 {\boldsymbol D}={\boldsymbol 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 {\boldsymbol D}'</math>, <math display="inline">{\boldsymbol T}</math> and <math display="inline">\bar{\boldsymbol 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">{\boldsymbol 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 {\boldsymbol 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>.
 <div id='img-9'></div>
 {| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 |-
-|[[Image:Draft_Samper_267560997-test-Figure9.png|600px|2D advection-conduction-absorption problem over a square domain of size equal to 8 units. ϕ<sup>p</sup>=8 at x=0, ϕ<sup>p</sup>=3 at x=8, qₙ=0 at y=0 and y=8. u = [2,0]<sup>T</sup>, k=1, s=20, w=20, γₓ=1 and γ<sub>y</sub>=0. Galerkin and FIC solutions obtained with a mesh of 8 ×8 four node square elements.]]
+|[[File:Draft_Samper_267560997_1174_fig9.png|600px|2D advection-conduction-absorption problem over a square domain of size equal to 8 units. ϕ<sup>p</sup>=8 at x=0, ϕ<sup>p</sup>=3 at x=8, qₙ=0 at y=0 and y=8. u = [2,0]<sup>T</sup>, k=1, s=20, w=20, γₓ=1 and γ<sub>y</sub>=0. Galerkin and FIC solutions obtained with a mesh of 8 ×8 four node square elements.]]
 |- style="text-align: center; font-size: 75%;"
 | colspan="1" | '''Figure 9:''' 2D advection-conduction-absorption problem over a square domain of size equal to 8 units. <math>\phi ^p =8</math> at <math>x=0</math>, <math>\phi ^p =3</math> at <math>x=8</math>, <math>q_n =0</math> at <math>y=0</math> and <math>y=8</math>. <math>{\boldsymbol u} = [2,0]^T</math>, <math>k=1</math>, <math>s=20</math>, <math>w=20</math>, <math>\gamma _x=1</math> and <math>\gamma _y=0</math>. Galerkin and FIC solutions obtained with a mesh of <math>8 \times 8</math> four node square elements.
 |}
+The second example is similar to the first one with <math display="inline">{\boldsymbol 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>.
 <div id='img-10'></div>
 {| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 |-
-|[[Image:Draft_Samper_267560997-test-Figure_10.png|600px|2D advection-conduction-absorption problem over a square domain of size equal to 8 units. ϕ<sup>p</sup>=8 at x=0, ϕ<sup>p</sup>=3 at x=8, qₙ=0 at y=0 and y=8. u = [20,0]<sup>T</sup>, k=1, s=20, w=20, γₓ=10 and γ<sub>y</sub>=0. Galerkin and FIC solutions obtained with a mesh of 8 ×8 four node square elements.]]
+|[[File:Draft_Samper_267560997_7137_Fig10.png|600px|2D advection-conduction-absorption problem over a square domain of size equal to 8 units. ϕ<sup>p</sup>=8 at x=0, ϕ<sup>p</sup>=3 at x=8, qₙ=0 at y=0 and y=8. u = [20,0]<sup>T</sup>, k=1, s=20, w=20, γₓ=10 and γ<sub>y</sub>=0. Galerkin and FIC solutions obtained with a mesh of 8 ×8 four node square elements.]]
 |- style="text-align: center; font-size: 75%;"
 | colspan="1" | '''Figure 10:''' 2D advection-conduction-absorption problem over a square domain of size equal to 8 units. <math>\phi ^p =8</math> at <math>x=0</math>, <math>\phi ^p =3</math> at <math>x=8</math>, <math>q_n =0</math> at <math>y=0</math> and <math>y=8</math>. <math>{\boldsymbol u} = [20,0]^T</math>, <math>k=1</math>, <math>s=20</math>, <math>w=20</math>, <math>\gamma _x=10</math> and <math>\gamma _y=0</math>. Galerkin and FIC solutions obtained with a mesh of <math>8 \times 8</math> four node square elements.
 |}
+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">{\boldsymbol 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.
 <div id='img-11'></div>
 {| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 |-
-|[[Image:Draft_Samper_267560997-test-Figure_11.png|600px|2D advection-conduction-absorption problem over a square domain of size equal to 10 units. ϕ<sup>p</sup>=8 at x=0, ϕ<sup>p</sup>=3 at x=10, qₙ=0 at y=0 and y=10. u = [0,1×10⁻³]<sup>T</sup>, k=1, s=12, w=12, γₓ=0 and γ<sub>y</sub>= 5×10⁻⁴. Galerkin and FIC solutions obtained with a mesh of 10 ×10 four node square elements.]]
+|[[File:Draft_Samper_267560997_7278_Fig11.png|600px|2D advection-conduction-absorption problem over a square domain of size equal to 10 units. ϕ<sup>p</sup>=8 at x=0, ϕ<sup>p</sup>=3 at x=10, qₙ=0 at y=0 and y=10. u = [0,1×10⁻³]<sup>T</sup>, k=1, s=12, w=12, γₓ=0 and γ<sub>y</sub>= 5×10⁻⁴. Galerkin and FIC solutions obtained with a mesh of 10 ×10 four node square elements.]]
 |- style="text-align: center; font-size: 75%;"
 | colspan="1" | '''Figure 11:''' 2D advection-conduction-absorption problem over a square domain of size equal to 10 units. <math>\phi ^p =8</math> at <math>x=0</math>, <math>\phi ^p =3</math> at <math>x=10</math>, <math>q_n =0</math> at <math>y=0</math> and <math>y=10</math>. <math>{\boldsymbol u} = [0,1\times 10^{-3}]^T</math>, <math>k=1</math>, <math>s=12</math>, <math>w=12</math>, <math>\gamma _x=0</math> and <math>\gamma _y= 5\times 10^{-4}</math>. Galerkin and FIC solutions obtained with a mesh of <math>10 \times 10</math> four node square elements.
 |}
+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">{\boldsymbol 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.
 <div id='img-12'></div>
 {| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
 |-
-|[[Image:Draft_Samper_267560997-test-Figure12.png|600px|2D advection-conduction-absorption problem over a square domain of size equal to 10 units. ϕ<sup>p</sup>=8 at x=0, ϕ<sup>p</sup>=3 at x=10, qₙ=0 at y=0 and y=10. u = [0,5]<sup>T</sup>, k=1, s=12, w=12, γₓ=0 and γ<sub>y</sub>=2.5. Galerkin and FIC solutions obtained with a mesh of 10 ×10 four node square elements.]]
+|[[File:Draft_Samper_267560997_4095_Fig12.png|600px|2D advection-conduction-absorption problem over a square domain of size equal to 10 units. ϕ<sup>p</sup>=8 at x=0, ϕ<sup>p</sup>=3 at x=10, qₙ=0 at y=0 and y=10. u = [0,5]<sup>T</sup>, k=1, s=12, w=12, γₓ=0 and γ<sub>y</sub>=2.5. Galerkin and FIC solutions obtained with a mesh of 10 ×10 four node square elements.]]
 |- style="text-align: center; font-size: 75%;"
 | colspan="1" | '''Figure 12:''' 2D advection-conduction-absorption problem over a square domain of size equal to 10 units. <math>\phi ^p =8</math> at <math>x=0</math>, <math>\phi ^p =3</math> at <math>x=10</math>, <math>q_n =0</math> at <math>y=0</math> and <math>y=10</math>. <math>{\boldsymbol u} = [0,5]^T</math>, <math>k=1</math>, <math>s=12</math>, <math>w=12</math>, <math>\gamma _x=0</math> and <math>\gamma _y=2.5</math>. Galerkin and FIC solutions obtained with a mesh of <math>10 \times 10</math> 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 <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).
+==REFERENCES==
+<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>
+'''[[#citeF-3|[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>
+'''[[#citeF-6|[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>
+'''[[#citeF-7|[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>
+'''[[#citeF-9|[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>
+'''[[#citeF-10|[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>
+'''[[#citeF-11|[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>
+'''[[#citeF-13|[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>
+'''[[#citeF-17|[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>
+'''[[#citeF-10¡8|[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>
+'''[[#citeF-19|[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>
+'''[[#citeF-21|[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>
+'''[[#citeF-22|[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>
+'''[[#citeF-24|[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>
+'''[[#citeF-25|[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>
+'''[[#citeF-27|[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>
+'''[[#citeF-28|[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>
+'''[[#citeF-29|[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>
+'''[[#citeF-30|[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>
+'''[[#citeF-32|[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>
+'''[[#citeF-33|[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;" | (A.1)
+|}
+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;" | (A.2)
+|}
+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;" | (A.3)
+|}
+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;" | (A.4)
+|}
+==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)
+|}
+The characteristic polinomial is
+{| 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)l^2 + 2\left({w\over 3}+1+{\beta ^1+\beta ^2\over 2}\right)l + \left({w\over 6} - 1 -\beta ^1+\gamma \right)=0</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (B.2)
+|}
+and  the characteristic roots are
+{| class="formulaSCP" style="width: 100%; text-align: left;"
+|-
+|
+{| style="text-align: left; margin:auto;width: 100%;"
+|-
+| style="text-align: center;" | <math>l_{1,2} = {-\left({w\over 3} +1+{\beta ^1+\beta ^2\over 2}\right)\pm \sqrt{\left({w\over 3} +1+{\beta ^1+\beta ^2\over 2}\right)^2 - \left({w\over 6} - 1 -\beta ^2-\gamma \right)\left({w\over 6} - 1 -\beta ^1+\gamma \right)}/\left({w\over 6} - 1 -\beta ^2-\gamma \right)}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (B.3)
+|}
+A non oscilatory (stable) numerical solution is found if <math display="inline">l_{1,2} \ge 0</math>. This is fulfilled if
+{| class="formulaSCP" style="width: 100%; text-align: left;"
+|-
+|
+{| style="text-align: left; margin:auto;width: 100%;"
+|-
+| style="text-align: center;" | <math>l_1 +l_2 = - {w/3 + 1+ {\beta ^1+\beta ^2\over 2}\over w/6 -1 -\beta ^2-\gamma }>0</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (B.4)
+|}
+{| class="formulaSCP" style="width: 100%; text-align: left;"
+|-
+|
+{| style="text-align: left; margin:auto;width: 100%;"
+|-
+| style="text-align: center;" | <math>l_1 l_2 = {w/6 -1 -\beta ^1 +\gamma \over w/6 -1 -\beta ^2-\gamma }>0</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (B.5)
+|}
+The inequality (B.4) implies
+{| class="formulaSCP" style="width: 100%; text-align: left;"
+|-
+|
+{| style="text-align: left; margin:auto;width: 100%;"
+|-
+| style="text-align: center;" | <math>{w\over 6} -1 -\beta ^2-\gamma \le 0</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (B.6)
+|}
+or
+{| class="formulaSCP" style="width: 100%; text-align: left;"
+|-
+|
+{| style="text-align: left; margin:auto;width: 100%;"
+|-
+| style="text-align: center;" | <math>\beta ^2\ge {w\over 6} - \gamma -1</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (B.7)
+|}
+The second inequality (B.5) implies
+{| class="formulaSCP" style="width: 100%; text-align: left;"
+|-
+|
+{| style="text-align: left; margin:auto;width: 100%;"
+|-
+| style="text-align: center;" | <math>{w\over 6} -1 -\beta ^1+\gamma \le 0</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (B.8)
+|}
+or
+{| 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;" | (B.9)
+|}
+Hence the stability limit of <math display="inline">\beta ^e</math>  is
+{| class="formulaSCP" style="width: 100%; text-align: left;"
+|-
+|
+{| style="text-align: left; margin:auto;width: 100%;"
+|-
+| style="text-align: center;" | <math>\displaystyle{\beta ^e=\beta _c={w\over 6}+\vert \gamma \vert -1}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (B.10)
+|}
+where <math display="inline">\beta _c</math> is the critical stabilization parameter.
+Indeed a stabilized solution will be obtained if <math display="inline">\beta ^e \ge \beta _c</math>.
+The condition for the onset of instability for the original (non-stabilized) problem is obtained by making <math display="inline">\beta =0</math> in Eq.(B.6) and (B.8) and reversing the sing of the inequalities, i.e. an ''oscilatory solution'' will be found if one of the two conditions is fullfiled
+{| 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{w\over 6}+\gamma -1\ge 0\\ \\ \displaystyle{w\over 6}-\gamma -1\ge 0 \end{array}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (B.11)
+|}
+Hence, the stability condition  is found as
+{| 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\ge{0}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (B.12)
 |}

Latest revision as of 12:15, 12 November 2018

Abstract

1 INTRODUCTION

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

3 FINITE ELEMENT FORMULATION

4 COMPUTATION OF THE STABILIZATION PARAMETER FOR LINEAR ELEMENTS

4.1 Computation of the element stabilization parameter βe

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

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

4.2 Critical stabilization parameter and unstability conditions

4.3 Iterative solution scheme

4.4 Optimal value of the stabilization parameter

4.5 Advective limit (w=0)

4.6 Absorptive limit (γ=0)

5 VALUE OF THE ELEMENT CHARACTERISTIC LENGTH PARAMETER

6 ADVANTAGES OF USING β VERSUS h AS THE STABILIZATION PARAMETER

7 1D NUMERICAL EXAMPLES

8 EXTENSION TO MULTI-DIMENSIONAL PROBLEMS

8.1 2D Examples

9 CONCLUSIONS

ACKNOWLEDGEMENTS

REFERENCES

APPENDIX A

APPENDIX B. STABILITY ANALYSIS OF THE DISCRETIZED FIC ADVECTION-DIFFUSION-ABSORPTION EQUATION FOR LINEAR ELEMENTS

Document information

Document Score

Share this document

claim authorship

4.1 Computation of the element stabilization parameter β^e

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

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