Published in Comput. Meth. Appl. Mech. Engng., Vol. 195, pp. 3926–3946, 2006
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.
Considerable effort has been spent in recent years to derive finite element methods (FEM)  for the solution of the advection-diffusion-reaction equation. The physical behaviour of this equation is varied. In the propagation regime, originated by large values of the productive term (the Helmholtz equation), solutions are exponentially modulated sinusoidal functions. Typical problems found here in the numerical solution are those of phase, amplitude and pollution errors . In the exponential regime, solutions are of the form of real exponential functions where absorptive (dissipative) or productive source terms are possible. Numerical schemes here find difficulties to approximate the sharp gradients appearing in the neighborhood of boundary and internal layers due to high Peclet and/or Damköhler numbers.
Stabilized methods to tackle above problems have been based on different kinds of streamline-upwind-Petrov-Galerkin (SUPG) [2-7], Galerkin Least Square [5-11], Subgrid Scale (SGS) [5,6,12,13] and Residual Free Bubbles  finite element methods. While a single stabilization parameter suffices to yield stabilized (and even nodally exact results) for the one-dimensional (1D) advection-diffusion and the diffusion-reaction equations (Vol. 3 in  and [8-15]), this is not the case for the diffusion-advection-reaction equation. Here, in general, two stabilization parameters are needed in order to ensure a stabilized solution for all range of physical parameters and boundary conditions [4,10,14]. As reported in [12,13] both SUPG, GLS and SGS methods with a single stabilization parameter (and indeed the standard Galerkin scheme) fail to obtain an stabilized solution for some specific boundary conditions in the exponential regime with negative (absorption) terms when there is a negative streamwise gradient of the solution.
Apart from the many practical applications of the advection-diffusion-reaction problem, the interest to derive stabilized numerical schemes for this equation extends to the solution of the transient advection-diffusion equation, as the discretized form of this equation in time can be interpreted as an analogous steady-state advection-diffusion-absorption problem. This analogy has been analyzed in [16,17].
In this paper we present a reliable and accurate stabilized formulation for the advection-diffusion-reaction equation in the exponential regime with absorption (destructive) reaction terms using a single stabilization parameter.
The stabilized formulation is based on the standard Galerkin FEM solution of the modified governing differential equations derived via the Finite Calculus (FIC) method . The FIC equations are obtained by expressing the balance of fluxes in a domain of finite size. This introduces additional stabilizing terms in the differential equations of the infinitessimal theory which are a function of the balance domain dimensions. These dimensions are termed characteristic length parameters and play a key role in the stabilization process.
The modified governing equations lead to stabilized numerical schemes using whatever numerical method. It is interesting that many of the stabilized FEM can be recovered using the FIC formulation. The FIC method has been successfully applied to problems of convection-diffusion [18-24], diffusion-absorption , incompressible fluid flow [24-30] and standard and incompressible solid mechanics [31-33] using finite element and meshless procedures.
The Galerkin FEM-FIC formulation here proposed introduces naturally an additional non linear dissipation term into the discretized equations which is governed by a single stabilization parameter. In the absence of the absorption term the formulation simplifies to the standard Petrov-Galerkin approach for the advection-diffusion problem For the diffusion-absorption case the diffusion-type stabilization term is identical to that recently obtained by Felippa and Oñate using a variational FIC approach and modified equation methods . The limit value of the stabilization parameter ensuring a stable numerical solution is derived. More accurate numerical results can be obtained using an expression of the stabilization parameter which is a function of the sign of the solution and that of its first and second derivatives. This introduces a non-linearity in the solution scheme and a simple iterative algorithm leading to stabilized and accurate results in just two iterations is presented. The merit of the new formulation is that it yields stabilized and accurate numerical solutions for the advection-diffusion-absorption equation for all range of the physical parameters and boundary conditions using a single stabilization parameter. As previously mentioned, this property was only attained by previous similar stabilized formulations using two stabilization parameters. We present a collection of 1D examples showing the efficiency and accuracy of the new FEM-FIC formulation for different values of the advective and absorptive terms.
In the last part of the paper the FIC formulation is extended to multi-dimensional advection-diffusion problems. The key ingredients of the formulation are given and some examples of the good performance of the new stabilized method for the solution of 2D problems are presented.
The governing equation for the 1D stationary advection-diffusion-absorption problem can be written in the FIC formulation as
In above equations is the state variable, is the problem domain, is the domain length, is the velocity field, is the diffusion, is the absorption, dissipation or destruction source parameter, is a constant source term, and are the prescribed values of the total flux and the unknown function at the Neumann and Dirichlet boundaries and , respectively and is a characteristic length which plays the role of a stabilization parameter. In the 1D problem and consist of four combinations at the end points of the problem domain.
Eqs.(1) and (2) are obtained by expressing the balance of fluxes in an arbitrary 1D domain of finite size within the problem domain and at the Neumann boundary, respectively. The variations of the transported variables within the balance domain are approximated by Taylor series expansions retaining one order higher terms than in the infinitessimal theory . The underlined terms in Eqs.(1) and (2) emanate from these higher order expansions and they are essential to derive stabilized numerical schemes using whatever discretization method. Note that as the characteristic length parameter tends to zero the FIC differential equations gradually recover the standard infinitessimal form giving in the limit (for ) in and .
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 . At the discrete level 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 as a proportion of a typical element dimension (i.e. the element length for 1D problems) .
Successful applications of the FIC method to a variety of problems in computational mechanics can be found in [18-31].
We will construct a standard finite element discretization of the 1D analysis domain length with index ranging from 1 to the number of elements . The state variable is approximated by over the analysis domain. The approximated variable is interpolated within each element with nodes in the standard manner, i.e.
where are the element shape functions and are nodal values of the approximate function . Substituting Eq.(5a) into Eqs.(1) and (2) gives
where and and are the residuals of the approximate solution in the problem domain and on the Neumann and Dirichlet boundaries and , respectively.
The weighted residual form of Eqs.(6)–(8) is written as
where and are test functions satisfying on . In Eq.(9a) and in the following and denotes the value of at the boundary . Using Eqs.(6) and (7) into Eq.(9a) we have
Assuming smooth enough solutions and integrating by parts the term involving in the first integral gives for
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 to be discontinuous. Also in Eq.(10) has been assumed to be constant within each element, (i.e. within ).
The weak form is obtained by integrating by parts the advective and diffusive terms within in the first integral of Eq.(10). This gives
Observation of the stabilization terms in the last integral of Eq.(11) shows that they introduce an additional diffusion (of value ) as well as an additional advection (of value ). 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 from the terms depending on in as
We define now
where is a stabilization parameter acting as a proportion of the actual physical diffusion .
Eqs. (12) and (13) allow to rewrite Eq.(11) as
Wee see clearly that the last term of Eq.(14) introduces within each element an additional diffusion of value .
Substituting expression (5b) into (14) and choosing a Galerkin method with within each element gives the discrete system of FE equations written in the standard matrix form as
where is the vector of nodal unknowns and the element contributions to matrix K and vector are
The amount of balancing diffusion in Eq.(16) clearly depends on the (non linear) stabilization parameter . The element and critical values of are deduced in the next section for linear two node elements.
We note that the integral of the term in Eq.(17) vanishes after asssembly when and are uniform over a patch of linear elements.
The matrix and the vector for two node linear elements are (for constant values of , and )
In Eqs.(18) and (19) index denotes element values.
Let us consider a mesh of just two elements 1 and 2 of equal size with nodes 1, 2 and 3.
For the particular case of the following difference equation is found from Eq.(18)
where and are the element Peclet number and a velocity independent dimensionless element number, respectively. The Damköler number can be obtained as .
We will solve now for the value of for Dirichlet boundary conditions at nodes 1 and 3.
The two following cases are considered.
Application of the stencil of Eq.(20) gives the following value for
A physically meaningful solution for any values of and must necessarily satisfy . Clearly this requires
The stencil of Eq.(20) gives
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 is proposed (see Appendix A):
where denotes values computed at the element mid point.
It is clear from Eq.(26) that the computation of the stabilization parameter requires the knowledge of the sign of the numerical solution and that of the first and second derivatives of 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.
The following constant value of over the mesh yields a stabilized solution (for a uniform distribution of and )
where is the critical stabilization parameter. The proof of Eq.(28) is given in the Appendix B.
The use of for all elements ensures a stabilized numerical solution for all ranges of and although it can yield slight overdiffusive results in some cases.
The use of 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)
The following two steps iterative scheme is proposed in order to obtain a stabilized and accurate solution.
Step 1. Compute a first stabilized solution using the critical value given by Eq.(28). This ensures a stabilized, although sometimes slightly overdiffusive solution.
Step 2.1. Compute the signs of the first and second derivatives of within each element. The second derivative field is obtained as follows
where denotes averaged values of the first derivative at node of element . At the boundary nodes the constant value of the derivative of within the element is taken in Eq.(30); i.e. .
Step 2.2. Compute the enhanced stabilized solution using the element value of 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.
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
Constants and depend on the boundary conditions of the problem. For the Dirichlet problem ( at and at ) we obtain
Let us now write the stencil of Eq.(20) for three arbitrary nodes , , and a uniform mesh with and as
where is the nodal value of the stabilization parameter or the th node. Note that in Eq.(33) has assumed to be constant within the tributary domain of node .
The optimal nodal value of the stabilization parameter is found from Eq.(33) as
where denotes the exact analytical value of as obtained from Eq.(31a) by
Substituting the exact nodal values of Eq.(35) into Eq.(34) gives the analytical expression of as
Eqs.(34) and (36) indicate that depends on the coordinate of node . The advective () and absorptive () limits of are given below.
In the examples presented in the next section we have found more convenient to compute by Eq.(34) with the numerical values of () computed from Eq.(35).
If the absorption terms are zero, , and from Eq.(36)
For , , , and Eq.(36) gives
Note that both the advective and absorptive limits of are independent of the nodal coordinates.
Remark. The optimal stabilization parameter for the advective limit (Eq.(37)) coincides with the classical expression obtained by many authors (see Vol. 3 in ).
The expression of Eq.(38) coincides with that obtained by Felippa and Oñate  using the FIC variational approach and the modified equation method.
The expression of the element characteristic length parameter is deduced from Eq.(13) as
where is defined in Eq.(25).
It can be readily deduced from Eq.(39) that vanishes in regions where is small. In high gradient zones for (advective limit), whereas for (absorptive limit). The sign of in the latter case depends on the sign of where and are defined in Eq.(27).
The expression of is needed in Eqs.(17) and (19) for the case of .
The stabilized formulation described above can be derived using either or as the stabilization parameter. The choice of instead of has the advantage that is defined as a positive number (see Eq.(25)), whereas is a distance that can be either positive or negative, depending on the values of the absorptive parameter and of the sign ratio (see previous section). The change of sign of within the mesh introduces a difficulty in the iterative solution process. This difficulty is avoided by selecting as the stabilization parameter, which is the option chosen in this work.
Indeed, the characteristic length parameter is still needed in the expression of (see Eqs.(17) and (19)). This does not pose any problem as the value (and sign) of can be effectively computed in the second step from Eq.(39) using the stabilized solution obtained in the first step.
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 and and different values of and .
Tables 1–11 show the numerical results obtained at the nodes for the standard Galerkin method () and using the element () and critical () values of . Recall that the results for and 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 yielding nodaly exact results.
We note that the numbers in the tables for the values of refer to element values for and and to nodal values for the case of .
Figures 1–6 show the distribution of along the domain for (Galerkin), and 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 (Table 9 and Figure 5) small oscillations around the exact distribution of are obtained for . 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 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 and and and , respectively. Both cases are solved for and . The nodal results are shown in Tables 12 and 13, while the distributions of are plotted in Figures 7 and 8.
|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 , with and . 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 stabilization procedure presented can be easily extended to multi-dimensional problems. Only the main ingredients of the formulation are presented here. For a more detailed description of the multi-dimensional formulation see .
Consider the general steady-state advection-diffusion-reaction equations written for the sourceless case as
where for 2D problems
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
where is the original differential equation as defined in Eq.(40).
The Dirichlet and boundary conditions of the FIC formulation are
where is the normal vector to the boundary where the normal flux is prescribed. As usual index denotes the prescribed values.
In Eqs.(42) and (44) is the characteristic vector of the FIC formulation which components play the role of stabilization parameters. The underlined terms in Eqs.(42) and (44) introduce the necessary stability in the numerical solution.
As shown in  the components of the characteristic length h allow to obtain stabilized finite element solutions in the presence of strong gradients of near the boundaries (boundary layers) and within the analysis domain (internal layers). If vector h is taken to be parallel to the velocity u the standard SUPG method is recovered . Keeping the more general form of h allows to reproduce the best features of the so called shock capturing or transverse dissipation schemes.
Oñate et al.  have recently shown that a general form of the characteristic lengths can be found by writting the FIC equations along the principal curvature directions of the solution denoted hereafter as and . The resulting FIC-FEM formulation is equivalent in this case (for linear elements) to adding a non linear stabilizing diffusion matrix to the standard infinitessimal equation. The stabilizing diffusion matrix depends on the velocities along the principal curvature directions of the solution which must be found iteratively. The iterative scheme can be simplified by assuming that the main principal curvature direction at each integration point within an element are coincident with the gradient vector direction at that point. For details see . This technique is extended here accounting for the absorption term.
Following the ideas described in  the FIC equations can be written as
where the stabilization matrix is given by
where is the angle that the main principal curvature direction of the solution forms with the global axis. The different terms in Eq.(47–48) are
where , are the velocity components along each of the principal curvature directions and and are the maximum values of the projection of the element sides along the principal directions and , respectively.
The values of and 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 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 with a zero value of the the stabilization matrix (i.e. ). 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 to be parallel to the gradient vector for simplicity. The direction is chosen orthogonal to .
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 .
The analysis domain in the first two 2D examples presented is a square of size 8 units. A relatively coarse mesh of four node bi-linear square elements is used in all cases.
The boundary conditions for the first and second examples are at the boundary and at and zero normal flux at and . This reproduces the condition of the 1D examples solved in previous sections. The first example is analized for , and giving , and 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 where is prescribed to a value of . The numerical results obtained with the standard (non stabilized) Galerkin solution (i.e. ) 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 .
|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 , and giving , and . 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 .
|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 at and and at and . The square domain has a side of 10 units in this case. A mesh of four node elements of unit length is used. The physical parameters are , and giving , and .
The correct solution for this absorption dominated problem has two boundary layers next to the boundaries at and . 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 at and and at and . The physical parameters are , and giving , and .
The correct solution has two boundary layers next to the boundaries at and . The Galerkin solution gives oscillatory results next to these boundaries and also next to the boundary at (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 similarly as for the 1D case. Details of the iterative scheme and examples of applications are given in .
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.
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).
 Zienkiewicz O.C. and Taylor R.L. The Finite Element Method. Volume 3. 5th Edition, Butterworth-Heinemann, 2001.
 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.
 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.
 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.
 Codina R. Comparison of some finite element methods for solving the diffusion-convection-reaction equation. Comput. Methods Appl. Mech. Engrg., 156, 186–210, 1998.
 Codina R. On stabilized finite element methods for linear systems of convection-diffusion-reaction equations. Comput. Meth. Appl. Mech. Engrg., 188, 61–82, 2000.
 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.
 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.
 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.
 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.
 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.
 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.
 Hauke G. A simple subgrid scale stabilized method for the advection-diffusion reaction equation. Comput. Methods Appl. Mech. Engrg. 2002; 191:2925–2947.
 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.
 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.
 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.
 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.
 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.
 Oñate E. Possibilities of finite calculus in computational mechanics. Int. J. Num. Meth. Engng., 60, 255–281, 2004.
 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.
 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.
 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.
 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.
 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.
 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.
 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.
 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.
 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.
 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.
 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.
 Oñate E., Perazzo F. and Miquel J. A finite point method for elasticity problems. Computers and Structures, 79, 22–25, 2151–2163, 2001.
 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.
 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.
 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.
Let us consider the modified form of the convection-diffusion reaction equation with the stabilization term having the form of a diffusion of value i.e.
where is the stabilization parameter. For simplicity the source term has been taken equal to zero in Eq.(A.1).
From Eq.(A.1) we can obtain
Let us estimate now the value of over a patch of finite elements. From the observation of the finite element matrices (Eq.(18)) we can write
where is the maximum value of the approximate finite element function over a patch of elements of equal size and , and 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
The discretized equation for the two linear element stencil of Eq.(20) can be written as
The characteristic polinomial is
and the characteristic roots are
A non oscilatory (stable) numerical solution is found if . This is fulfilled if
The inequality (B.4) implies
The second inequality (B.5) implies
Hence the stability limit of is
where is the critical stabilization parameter.
Indeed a stabilized solution will be obtained if .
The condition for the onset of instability for the original (non-stabilized) problem is obtained by making 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
Hence, the stability condition is found as