Published in Comput. Meth. Appl. Mech. Engng. , Vol. 195, pp. 339-362, 2006
A methodology for error estimation and mesh adaptation for finite element (FE) analysis of incompressible viscous flow is presented. The error estimation method is based on the evaluation of the energy rate (the power) of the FE residuals of the momentum and incompressibility equations. The residuals are computed using recovered values of the derivatives of the velocity and pressure variable obtained via a nodal derivative recovery technique. Two mesh adaptation procedures based on: a) the equi-distribution of the residual power among the elements in the mesh and b) the equi-distribution of the density of the total residual power are presented. The stabilized form of the Navier-Stokes equations using a finite calculus (FIC) formulation is solved with a fractional step FE scheme. This allows the use of linear triangles and tetrahedra with an equal order interpolation for the velocity and pressure variables. A nodal-based approach is used for computing the residual power integrals. The examples show the ability of the error estimation and the mesh adaptation process to capture the high gradients of the solution in the vecinity of boundary layers and in zones of high vorticity of the fluid for both steady state and transient flow situations.
The estimation of the error in the finite element solution of the Navier-Stokes equations for an incompressible fluid is a difficult problem . The presence of convective terms in the momentum equations prevents the use of standard error estimators based on the energy norm of the error developed for elliptic type problems typical of solid mechanics and heat transfer analysis [2-4]. Extension of these methods to fluid mechanics problems have been reported in . A different class of mesh refinement procedures for fluid flow problems based on error indicators using approximations of the Hessian operator has been reported in [6-10]. Alternative error estimation approaches for fluid flow problems can be found in [11-16].
The error estimation method presented in this paper is based on the evaluation of a particular form of the energy rate (the power) of the finite element residuals of the momentum and mass balance (incompressibility) equations. The residuals are computed using recovered values of the derivatives of the velocity and pressure variable obtained via a nodal derivative recovery technique. Two mesh adaptation processes based on: a) the equi-distribution of the residual power among the elements in the mesh and b) the equi-distribution of the density of the total residual power are presented. The stabilized form of the Navier-Stokes equations using a finite calculus (FIC) procedure developed by the authors [17-21] is used for the FE solution of the examples presented in the paper using a fractional step scheme. This allows the use of linear triangles and tetrahedra with an equal order interpolation for the velocity and pressure variables. A nodal-based approach is used for computation of the residual power error and the total power generated by the viscous stresses and the pressure terms in the mesh. The residual power form derived from the FIC approach includes power terms emanating from the velocity field and from the gradient of the velocity field. It is shown that accounting for the terms involving the velocity gradient provides an effective error estimation and mesh adaptation procedure for CFD problems. Also the FIC form of the residual power expression can be effectively extended for error estimation and mesh adaptivity of other problems in computational mechanics .
The examples show the ability of the error estimation and the mesh adaptation process to capture the high gradients of the solution in the vecinity of boundary layers and in zones of high vorticity of the fluid for both steady state and transient flow situations.
We consider here the solution of the Navier-Stokes equations for an incompressible fluid written as
For the transient case
with when is the number of space dimensions ( for 2D flow problems).
In eq.(3) is the fluid density (here assumed to be constant), is the velocity component in the -th direction, the pressure, the body forces and the viscous stress tensor components related to the velocity gradients through the fluid viscosity by
Summation convention for repeated indexes in products and derivatives is used unless otherwise specified.
The field equations are completed by the standard boundary conditions
where and are prescribed tractions and the prescribed displacement at the Neuman and Dirichlet boundaries and , respectively.
We will choose an equal order finite element (FE) approximation for the velocity and pressure
where denotes approximate values, are the shape functions, are nodal values and is the number of nodes in the mesh . Standard three node linear triangles are used in the examples given in the paper.
It is well known that stabilized finite element schemes are needed for overcoming the numerical instabilities induced by high values of the convective terms and by the incompressibility constraint when equal order interpolation for the velocities and the presure are used. A number of stabilized FE methods are available in the literature and most of them can be casted in the following form
In above and denote the approximate finite element residuals for the momentum and mass balance equations, respectively, and are appropriate functions of the shape functions and and are stabilization parameters which depend on the element dimensions and the physical variables of the problem.
Eqs.(9) and (10) can be interpreted as extended Galerkin expressions where the last integrals over the element interiors has been added to the standard Galerkin forms of the momentum and mass balance equations in order to get a physically meanful (stabilized) numerical solution. Many popular stabilized methods (SUPG, GLS, SGS, etc.) can be written in the form of Eqs.(9) and (10) and they just differ in the choice of the weighting functions and and in the values of the stabilization parameters [1,24-30].
A somehow different class of stabilized method developed by the authors is based on a finite calculus (FIC) formulation [17-21]. Here the standard momentum and mass balance equations are modified by enforcing the balance laws in a domain of finite size. The modified governing equations contain additional terms which play the role of stabilizing terms in the numerical solution. The Galerkin FE form of the FIC equations of an incompressible fluid can also be written by integral forms analogous to Eqs.(9) and (10) and hence we will retain these equations as general expressions applicable to any stabilized method. Details on the FIC formulation used for the solution of the examples presented in the paper are given in a later section.
We will assume that a finite element solution has been found for the velocity and pressure variables using any stabilized FE scheme. The next step in the error estimation and mesh adaptivity method proposed is to compute the enhanced recovered values of the derivatives of the velocity and the pressure at the element nodes. The derivative recovery procedure used here is explained in the next section.
The recovered nodal derivative fields of the velocities and the pressure are used now to compute the finite element residuals of the discretized momentum and mass balance equations (1) and (2) and the Neumann boundary conditions (6) as:
where denote values computed using the recovered derivative fields. In Eqs.(11) can be interpreted as the ith component of an out-of-balance body force, whereas represents an spurious non-zero volumetric strain rate. Similarly, the recovered residuals are out-of-balance tractions at the Neuman boundary.
The power due to the recovered residuals can be written as
where the weights and are functions of the approximate velocity vector field and the pressure, respectively. Indeed the simplest choice is and . However, many other options are possible. In this work we have chosen
where is a characteristic length of the order of the element size. This choice is a result of the finite calculus formulation described in a next section.
Note that at the solid boundaries and at the exterior boundary on . This justifies neglecting the boundary integral in Eq.(12a) by choosing .
The weights and as defined in Eq.(13) are discontinuous for standar finite element interpolations. A continuous discrete form of these weights can be found by using the recovered values of the derivatives of and as explained in Section 4.
The integrals over the domain can be computed by summing up the individual element contributions in the usual manner. However in the examples presented in the paper the volume integrals in Eqs.(12) are computed using a nodal integration scheme as
where denotes values computed at node , is the number of internal nodes in the mesh and is the tributary area of node (Figure 1). In Eqs.(14) absolute values are used to avoid random cancelation of the residual power terms. Indeed Eqs.(14) imply nodal continuity of the weighting functions obtained via the derivative recovery process explained in the next Section. The global error for a given mesh is estimated as
where is the estimated error value for an individual node defined as
We note that the residual power expressions chosen can be interpreted as the projection of the recovered residuals of the momentum and incompressibility equations against the weighting functions and , respectively. In this sense, Eqs.(12) can be viewed as a particular class of residual projection method. A review of these procedures can be found in  and references herein.
The recovered derivative field is defined so as to satisfy the following variational form
where are the same shape functions used to approximate the velocity field and are known derivatives values directly obtained from the finite element solution.
A standard finite element approximation is chosen for the recovered derivative field; i.e.
Substituting Eq.(18) into (17) gives the following system of equations from which the recovered nodal derivative values can be obtained
where is the number of nodes in the mesh.
The element contributions to and are
We recall that in above is the (discontinuous) derivative field obtained from the finite element solution and are the recovered derivative values at node . The process is repeated for .
The full form of matrix is used for the solution of Eq.(19) using a standard Jacobian iterative scheme.
The same procedure is followed in order to recover the nodal derivatives of the pressure and the nodal second derivatives of the velocity field needed for computation of all the terms in the recovered momentum equations .
A global error parameter is defined as
where is the prescribed global error value and is the total power in the mesh defined here as
The first term in Eq.(23) is the power of the viscous stresses, whereas the second one is the power of the pressure gradient terms. The integral in Eq.(23) is computed using nodal information by an expression analogous to Eq.(15).
The mesh adaption process aims to reduce the global error (leading to a value of ) while satisfying a local error condition for each individual element. The two options considered in this work for the local error condition are presented next.
A local error parameter is defined for each element as
where is the residual power error for element and is the total number of elements in the mesh. Eq.(24) is based in the so called equi-distribution of the global error, i.e. if for all elements then the total residual power error is uniformly distributed in the mesh .
The error indicator for the element is defined as
Parameters and are found by analyzing the convergence rate of the global and local error parameters as described in . In the examples presented in the paper the following values of and have been used
where is the order of the finite element interpolation (i.e. for linear elements).
The new element size is computed in terms of the error indicator and the current element size as
An alternative local error parameter based on the equidistribution of the density of the residual power can be defined as
where and are the area of the analysis domain and the area of an individual element, respectively. Now if the density of the error (error per unit area) is uniformly distributed over the mesh.
The expression for the error indicator for the element is again given by Eq.(25). The following values of and are now chosen as deduced following the procedure described in 
The new element sizes are computed from Eq.(27).
This mesh adaptation criteria leads to meshes where the smaller elements are concentrated in zones where the higher gradients of the velocity and pressure occur, whereas much larger elemens are placed in other zones of the domain . This is adequate for capturing accurate solutions in the vecinity of boundary layers or singularity zones. The drawback of this procedure is that it tends to refine progressively (and very fast) the mesh in these zones until . This leads to regions with very small elements and hence to an extremely large number of elements in the mesh. A limiting value in the size of the smallest element in the new mesh is mandatory in these cases.
Alternatively the criterium based on the equidistribution of the global error described in the previous section leads to meshes with a “smoother” distribution of the element sizes between the regions where high gradients occur and the other zones in the domain.
A comparison of these two mesh adaptions procedures for linear structural analysis can be found in .
The mesh refinement strategy is as follows.
Step 1 Compute the nodal recovered values of the derivatives of the velocity and pressure fields using the algorithm described in Section 4.
Step 2 Compute the residual power error for each element in the mesh and the total residual power error using Eqs.(15) and (16).
Step 3 Compute the global and local (element) error parameters and .
Step 4 Compute the new element sizes using Eq.(27). Perform a smoothing of the new element sizes to avoid problems during the regeneration of the mesh.
Step 5 Regenerate a new mesh. Minimum and maximum sizes of the elements in the mesh are prescribed in order to have a better control of the mesh generation process. The mesh generation is carried out using the GiD pre/postprocessing system .
Above mesh refinement strategy is equally applicable for the two mesh adaption procedures explained in the previous sections.
For steady-state flow problems the error estimation and the mesh refinement process are performed after a steady state solution has been found. For transient problems the mesh refinement is performed after a prescribed number of time steps.
The FIC governing equations for incompressible viscous flows are obtained by applying the standard conservation laws expressing balance of momentum and mass over a control domain. Assuming that the control domain has finite dimensions and expressing the variation of mass and momentum over the domain using Taylor series expansions of one order higher than those used in the standard infinitesimal theory, the following expressions are found [17-21]:
where and were defined in Eqs.(3) and (4), respectively.
Eqs.(30) and (31) are the FIC forms of the governing differential equations for an incompressible flow. The terms underlined in Eqs.(30) and (31) introduce naturally the necessary stabilization at the discretization level. The so called characteristic length vector is defined as (for 2D problems)
where and are characteristic dimensions of the finite control domain where balance of momentum and mass conservation is enforced. The components of vector introduce the necessary stabilization along the streamline and transverse directions to the flow in the discrete problem [17-21].
The method to derive the modified differential equations (30) and (31) incorporating the stabilization terms is termed finite calculus (FIC) as a reference to the standard infinitesimal calculus where the size of the domain where balance of mass and momentum is enforced is assumed to be negligible. Note that for the standard infinitesimal form of the momentum and mass balance equations is recovered .
Eqs.(30) and (31) are completed by the following FIC boundary conditions .
Equilibrium of surface tractions at the boundary
Prescribed velocity at the boundaries
In eq.(34) and denote the tangential velocity to the boundary and its prescribed value, respectively.
Eq.(35) expresses the balance of mass on an arbitrary domain next to the boundary and and denote the velocity normal to the boundary and its prescribed value, respectively. The value of is zero on solid walls and stationary free surfaces.
Also in eqs.(34) and (35) and are the parts of the boundary of where the tangential and normal velocities are prescribed, respectively. The Dirichlet boundary is defined as .
The underlined terms in eqs.(33) and (35) introduce the necessary stabilization at the boundaries in a form consistent with that of eqs.(30) and (31). Details of the derivation of eqs.(30–35) can be found in  whereas the derivation of eq.(9) is shown in .
We will choose again an equal order finite element approximation for the velocity and pressure variables given by Eqs.(8)
Again we will assume that a finite element solution has been found for the velocity and the pressure variables. The finite element algorithm used for the examples solved in this paper is based on a fractional step scheme [19-21,30]. The next step is to obtain enhanced recovered values of the derivatives of the velocity and the pressure at the element nodes as explained in previouss section.
The recovered nodal derivative fields of the velocities and the pressure are used now to compute the FIC residuals of the momentum and mass balance equations (1) and (2):
where once more denote values computed using the recovered derivative fields.
The total power due to the FIC residuals can be written as
Substituting the expressions of and of Eqs.(36) into (37) and integrating by parts some derivative terms gives
Numerical experiments have shown that the contribution of the boundary integrals in Eqs.(38) is negligible versus the value of the integrals over the domain. Hence the contribution of the surface integrals in Eqs.(38) will not be taken into account from now onwards. An explanation for this is that at the solid boundaries in a viscous flow and hence the line integrals in Eq.(38a) vanish there. Note also that at the exterior boundary on and on and this justifies neglecting the line integral at those boundaries. In Eq.(38b) the first line integral vanishes if is made equal to the prescribed normal velocity when solving the system of discretized equations. The second integral in Eq.(38b) can be neglected by accepting that either (external boundary) or at the boundary.
In summary, the relevant residual power terms for the FIC momentum and mass balance equations are
In the examples presented in the paper the integrals in Eq.(39) are computed using a nodal integration scheme. This requires computing the recovered nodal value of the derivative of the velocity and pressure field. This is an unexpensive task as the recovered nodal values are needed to evaluate the recovered residual. The nodal computation of the residual power expressions is therefore written as
the different terms have the meaning explained for Eq.(14).
We note that in the examples presented next the characteristic length parameters have been defined as the characteristic sizes of each tributary domain for a node (see Figure 1).
|Figure 1: Tributary domain for a node and definition of characteristic lengths.|
The global and element error norm are defined by Eqs.(15) and (16), respectively. The error estimation process and the mesh adaptivity scheme follow precisely the steps explained in Section 5.
Eqs.(40) define the weighting functions in the general residual power forms of Eqs.(12) and (13) as
These expressions, obtained using arguments from the FIC formulation, can be used to compute variational residual forms of the error in order to derive effective mesh refinement procedures for other problems of computational mechanics using the general scheme presented in Section 5.3 .
We note again however that Eqs.(41) are not the only option for and and some other choices are explored in  to derive a variety of error estimation and mesh refinement procedures for convection-diffusion problems. Numerical results reported in  show that the expressions given by Eqs.(41) have a better performance.
Typical parameters used for the examples presented are: prescribed global error , maximum and minimum sizes of the new elements generated m and m.
This a transient problem. The mesh refinement process has been carried out at a time . The value of the Reynolds number for this problem is 100. Figure 2 shows the initial mesh of 2574 nodes and 5148 linear elements used.
|Figure 2: Flow past a NACA airfoil. Characteristic length =1m, . Original mesh|
The mesh adaptation procedure chosen first is that of equidistribution of the global error. Figures 3 and 4 show the mesh obtained after five refinement steps. The distribution of the velocity field for that mesh is shown in Figures 5 and 6. Table 1 shows the value of the global error and the evolution of the number of elements and nodes during the mesh refinement process. The percentage errors in the drag and lift values are also given in Table 1. We note that the solution was found using a laminar flow model. Therefore some error in the drag values is to be expected even for fine meshes. The error of 6.4% obtained in this case is quite acceptable.
|Figure 3: Mesh obtained after five refinement steps using the equidistribution of the global error.|
|Figure 4: Details of the refined mesh of Figure 3 in the vecinity of the airfoil|
|Figure 5: Velocity field after five refinement steps|
|Figure 6: Detail of the velocity field|
|Mesh||Global error (%)||Nodes||Elements||Drag error %||Lift error%|
The same problem was solved using the criterium of the equidistribution of error density. The same maximum and minimum sizes for the generated elements are used (i.e., m and m). In this case, after the third step the mesher was unable to generate a new grid from the element size distributions given by the method. Figures 7 and 8 show the mesh obtained after the second step. The failure of the mesh generator is due to the very large gradients of the element size field leading to an extremely high density of elements concentrated near the airfoil.
|Figure 7: Mesh obtained after two refinement steps using the equidistribution of the error density (min size = 0.001m)|
|Figure 8: Detail of the refined mesh shown in the previous figure|
|Mesh||Global error (%)||Nodes||Elements||Drag error %||Lift error%|
|Figure 9: Mesh obtained after four refinement steps using the equidistribution of the error density (min size = 0.002m).|
|Figure 10: Detail of the refined mesh of Figure 9.|
Table 2 shows that the criterium of equidistribution of error density leads to a smaller error for the drag and lift in fewer time steps. Note however the increase in the number of elements and nodes with respect to that of Table 1. In order to avoid this problem, the minimum element size was increased next to 0.002m. This allowed to complete the mesh refinement process and the results are shown in Table 3 and Figures 9 and 10. Note that the drag and lift errors have now slightly increased, showing the dependence of these values with the minimum element size. Again we note that these errors can be reduced using a turbulence model.
|Mesh||Global error (%)||Nodes||Elements||Drag error %||Lift error%|
Finally Figures 11 and 12 respectively show the evolution of the drag and lift errors for the three cases considered and Figure 13 shows the convergence of the global error.
|Figure 11: Flow past a NACA airfoil. Convergence of the drag error for both refinement methods|
|Figure 12: Flow past a NACA airfoil. Convergence of the lift error for both refinement methods. Numbers next to legend denote the minimum element size prescribed|
|Figure 13: Convergence of the global error for both refinement methods|
The mesh refinement process for this transient problem has been carried out at a time . Figure 14 shows the original mesh of 2135 nodes and 5365 linear elements used.
Using the equidistribution of the global error, the maximum and minimum size of the elements to be generated are 0.3m and 0.001m, respectively. Figures 15 and 16 show the mesh obtained after two refinement steps. Figures 17 and 18 show the velocity field for the third mesh. Table 4 shows the value of the global error and the evolution of the number of elements and nodes during the mesh refinement process.
|Figure 14: Flow past a cylinder. Radius=1m, . Original mesh|
|Figure 15: Mesh obtained after two refinement steps using the equidistribution of the global error|
|Figure 16: Detail of the refined mesh in the neighbourhood of the cylinder|
|Figure 17: Velocity field after two refinement steps. Equidistribution of global error|
|Figure 18: Detail of the velocity field|
|Mesh||Global error (%)||Nodes||Elements|
The same problem was solved using the criterium based on the equidistribution of the error density. Now the minimum size of the elements generated is increased to 0.003m, in order to avoid problems in the regeneration of the mesh. The mesh obtained after two steps is shown in Figures 19 and 20. Table 5 shows the evolution of the global error and the number of elements and nodes during the mesh refinement process.
Figure 21 shows the convergence of the global error for both cases.
|Figure 19: Mesh obtained after two refinement steps using the equidistribution of the error density|
|Figure 20: Details of the refined mesh of Figure 19|
|Mesh||Global error (%)||Nodes||Elements|
|Figure 21: The convergence of the global error for both cases.|
This is another well known benchmark problem in the CFD literature. The problem is solved for a Reynolds number of 1000. The mesh regeneration was performed at the time .
The initial mesh of 477 nodes and 952 linear elements is shown in Figure 22. The mesh obtained after three refinement steps and the corresponding velocity field are shown in Figure 23. The mesh adaptation criterium based on the equidistribution of the global error was used. Table 6 shows the value of the global error and the evolution of the number of elements and nodes during the mesh refinement process.
Figure 24 shows the convergence of the global error.
|Figure 22: Driven cavity flow. Side length =1m, . Original mesh|
|Figure 23: Mesh obtained after three refinement steps (equidistribution of the global error) and velocity field.|
|Mesh||Global error (%)||Nodes||Elements|
|Figure 24: Convergence of the global error in the cavity problem.|
|Figure 25: High-lift airfoil configuration. Characteristic length =0.6m, Re=100. Initial mesh|
|Figure 26: Mesh obtained after two refinement steps using equiditribution of global error ()|
|Figure 27: Mesh obtained in the fith after changing to in the fourth step|
This last example is a multiple airfoil configuration including slat, wing and flap. The problem has been solved with a laminar solver and a Reynolds number of 100. In this example the maximum size of the generated elements has been reduced to 0.05m. The initial mesh of 11703 nodes and 23480 elements is shown in Figure 25. The resulting mesh after two refinement steps is shown in Figure 26. The convergence of the global error with the mesh refinement is shown in Table 7. The algorithm based on the equidistribution of the global error has been used in this example.
|Mesh||Global error (%)||Nodes||Elements|
After the fourth step the prescribed global error was reduced from to . The result is shown in Figure 27.
The velocity field for the fifth step is shown in Figure 28.
Figures 29 and 30 show the convergence of the drag and lift error, respectively. Figure 31 shows the evolution of the global error.
|Figure 28: Velocity field obtained after five refinement steps|
|Figure 29: The convergence of the drag error|
|Figure 30: The convergence of the lift error|
|Figure 31: The convergence of the global error|
An error estimator for the finite element analysis of incompressible fluid flow problems has been proposed. The estimator is based on the weigthed form (the power) of the residuals of the stabilized momentum and mass balance equations. The particular form of the weighting functions used has been found from a finite calculus formulation of the Navier-Stokes equations. This weighted residual form derived is also applicable for error estimation and mesh adaptivity in other problems in computational mechanics .
A simple mesh refinement strategy based on the equal distribution of the total residual power among the elements in the mesh has been presented. The refinement strategy has proven to be more cost-effective than that based on the equal distribution of the density of the error also presented in the paper. The efficiency of the error estimation and the mesh adaptation process has been verified in the solution of steady state and transient flow problems using a FIC stabilized formulation for the Navier-Stokes equations and a standard fractional step scheme.
The method yields refined meshes which capture very well the features of the velocity and pressure fields in high gradient zones near the solid boundaries and within the analysis domain. Note that in the problems presented here the global error is not reduced below its limited value. The reason is that a minimun element size has been prescribed in order to facilitate the mesh generation process. Should this limit be reduced then the global error would be subsequently reduced.
Further work in this field will focus in the relationship of the residual power error values and the error in selected individual variables (velocity, pressure, etc) or other solutions outputs target (drag, lift, etc.).
The numerical results were obtained using the finite element code Tdyn provided by the company COMPASS Ingeniería y Sistemas SA . Thanks are also given to Dassault Aviation for providing the multiple airfoil geometry for the last example.
 O.C. Zienkiewicz and R.L. Taylor. The finite element method. 5th Edition, 3 Volumes, Butterworth–Heinemann, 2000.
 O.C. Zienkiewicz and J.Z. Zhu. A simple error estimator and adaptive procedure for practical engineering analysis. Int. J. Num. Meth. Eng., 24, 337-57, 1987.
 J.Z. Zhu and O.C. Zienkiewicz. Adaptive techniques in the finite element method. Comm. App. Num. Math., 4, 197–204, 1988.
 Oñate E. and G. Bugeda. A study of mesh optimality criteria in adaptive finite element analysis. Engineering Computations, Vol. 3, 307-321, 1994.
 J. Wu, Z.J. Zhu, J. Szmelter and O.C. Zienkiewicz. Error estimation and adaptivity in Navier-Stokes incompressible flows. Comput Mechancis, 6, 254–70, 1990.
 R. Löhner, K. Morgan and O.C. Zienkiewicz. An adaptive finite element procedure for compressible high speed flows. Comp. Meth. Appl. Mech. Eng., 51, 441–65, 1985.
 R. Löhner, K Morgan and O.C. Zienkiewicz. Adaptive grid refinement for the Euler and compressible Navier-Stokes equations. in Accuracy Estimates and Adaptive Refinements in Finite Element Computations (eds I. Babuska, O.C. Zienkiewicz, J. Gago y E.R. de A. Oliveira), cap. 15, pp. 281–98, Wiley, Chichester, 1986.
 J. Peraire, M. Vahdati, K. Morgan and O.C. Zienkiewicz. Adaptive remeshing for compressible flow computations. J. Comp. Phys., 72, 449–66, 1987.
 J.T. Oden, T. Strouboulis and P. Devloo. Adaptive finite element methods for high speed compressible flows. Int. J. Num. Meth. Eng., 7, 1211–28, 1987.
 R. Löhner. Adaptive remeshing for transient problems. Comput. Meth. Appl. Mech. Eng., 75, 195–214, 1989.
 J.T. Oden, W. Wu and M. Ainsworth. An a posteriori error estimate for finite element approximations of the Navier-Stokes equations, Comput. Meth. Appl. Mech. Eng., 111, 185–202, 1993.
 C. Johnson, R. Rannacher and M. Boman. Numerics and hydrodynamic stability: Towards error control in CFD. SIAM J. Numer. Anal., 32, 1058–1079, 1995.
 R.C. Almeida, R.A. Feijoo, A.C. Galeao, C. Padra and R.S. Silva. Adaptive finite element computational fluid dynamics using an anisotropic error estimator. Comput. Meth. Appl. Mech. Eng., 182, 379–400, 2000.
 J.M. Cascón, G.C. García and R. Rodríguez. A priori and a posteriori error analysis for a large-scale ocean circulation finite element model. Comput. Meth. Appl. Mech. Eng., 192, 5305–5327, 2003.
 W. Bangerth and R. Rannacher. Adaptive finite elements for differential equations. Birkhäuser Verlag, Basel, 2003.
 D.L. Darmofal, D.A. Venditti. Anisotropic grid adaptation for functional outputs: Application to two-dimensional viscous flows. Journal of Computational Physics, 187, 22–66, 2003.
 E. Oñate. Derivation of stabilized equations for advective-diffusive transport and fluid flow problems. Comput. Methods Appl. Mech. Engrg., 151:1-2, 233–267, 1998.
 E. Oñate. 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.
 E. Oñate and J. García. 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.
 E. Oñate, J. García, G. Bugeda and S.R. Idelsohn. A general stabilized formulation for incompressible fluid flow using finite calculus and the FEM. In Towards a new fluid dynamics with its challenges in aeronautics, J. Periaux, D. Chamption, O. Pironneau and Ph. Thomas (Eds.), CIMNE, Barcelona, Spain 2002.
 E. Oñate, J. García, S.R. Idelsohn and F. Del Pin. FIC formulations for finite element analysis of incompressible flows. Eulerian, ALE and Lagrangian approaches, Submitted to Comput. Methods Appl. Mech. Engng., July 2004.
 E. Oñate, J. Arteaga and R. Flores. A weigthed recovered residual method for error estimation and mesh adaptivity in convection-diffusion problems. Research Report CIMNE, Barcelona, September 2004.
 GiD. The personal pre/postprocessing system. CIMNE, Barcelona, www.gidhome.com, 2004.
 M.A. Cruchaga and E. Oñate. A finite element formulation for incompressible flow problems using a generalized streamline operator. Comput. Methods Appl. Mech. Engng., 143, 49–67, 1997.
 T.J.R. Hughes and M. Mallet. A new finite element formulations for computational fluid dynamics: III. The generalized streamline operator for multidimensional advective-diffusive systems. Comput Methods Appl. Mech. Engrg., 58, 305–328, (1986a).
 T.J.R. Hughes, L.P. Franca and G.M. Hulbert. A new finite element formulation for computational fluid dynamics: VIII. The Galerkin/least-squares method for advective-diffusive equations. Comput. Methods Appl. Mech. Engrg., 73, 173–189, (1989).
 T.J.R. Hughes. Multiscale phenomena: Green functions, subgrid scale models, bubbles and the origins of stabilized methods. Comput. Methods Appl. Mech. Engrg, 127, 387–401, (1995).
 R. Codina. Comparison of some finite element methods for solving the diffusion-convection-reaction equation. Comput. Methods Appl. Mech. Engrg., 156, 185–210, (1998).
 R. Codina. Stabilization of incompressibility and convection through orthogonal sub-scales in finite element method. Comput. Methods Appl. Mech. Engrg., 190, 1579–1599, (2000).
 O.C. Zienkiewicz and R. Codina. A general algorithm for compressible and incompressible flow. Part I: The split characteristic based scheme. Int. J. Num. Meth. in Fluids, 20, 869-85, (1995).
 Tdyn. An unstructured finite element code for fluid dynamic analysis, COMPASS Ingeniería y Sistemas SA, www.compassis.com, 2003.