doi: 10.1016/j.cma.2004.07.054

​

​

==Abstract==

​

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.

​

==1 INTRODUCTION==

​

The estimation of the error in the finite element solution of the Navier-Stokes equations for an incompressible fluid is a difficult problem <span id='citeF-1'></span>[[#cite-1|[1]]]. 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 <span id='citeF-2'></span>[[#cite-2|[2-4]]]. Extension of these methods to fluid mechanics problems have been reported in <span id='citeF-5'></span>[[#cite-5|[5]]]. 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 <span id='citeF-6'></span>[[#cite-6|[6-10]]]. Alternative error estimation approaches for fluid flow problems  can be found in <span id='citeF-11'></span>[[#cite-11|[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 <span id='citeF-17'></span>[[#cite-17|[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 <span id='citeF-22'></span>[[#cite-22|[22]]].

​

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.

​

==2 STABILIZED FINITE ELEMENT SOLUTION OF THE INCOMPRESSIBLE NAVIER-STOKES EQUATIONS==

​

We consider here the solution of the Navier-Stokes equations for an incompressible fluid written as

​

'''Momentum'''

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>r_{m_i}=0 \qquad \hbox{in }\Omega  </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (1)

|}

​

'''Mass balance'''

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>r_v=0 \qquad \hbox{in }\Omega  </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (2)

|}

​

For the transient case

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>r_{m_i}=\rho \left[{\partial u_i \over \partial t} + u_j {\partial u_i \over \partial x_j}\right]+{\partial p\over \partial x_i}-{\partial \tau _{ij}\over \partial x_j}-b_i </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (3)

|}

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>r_v = {\partial u_i\over \partial x_i} </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (4)

|}

​

with <math display="inline">i,j=1,n_d</math> when <math display="inline">n_d</math> is the number of space dimensions (<math display="inline">n_d =2</math> for 2D flow problems).

​

In eq.(3) <math display="inline">\rho </math> is the fluid density (here assumed to be constant), <math display="inline">u_i</math> is the velocity component in the <math display="inline">i</math>-th direction, <math display="inline">p</math> the pressure, <math display="inline">b_i</math> the body forces and <math display="inline">\tau _{ij}</math> the viscous stress tensor components related to the velocity gradients through the fluid viscosity <math display="inline">\mu </math> by

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\tau _{ij}=2\mu \left(\varepsilon _{ij}-{1\over 3} {\partial u_k\over \partial x_k}\delta _{ij}\right)</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" |  (5a)

|}

​

with

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\varepsilon _{ij}={1\over 2}\left({\partial u_i\over \partial x_j}+{\partial u_j\over \partial x_i}\right)</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" |  (5b)

|}

​

Summation convention for repeated indexes in products and derivatives is used unless otherwise specified.

​

The field equations are completed by the standard boundary conditions

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\sigma _{ij}n_j - t_i^p = 0 \quad \hbox{ on }  \Gamma _t</math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (6)

|-

| style="text-align: center;" | <math> u_i-u_i^p  = 0 \quad \hbox{ on }  \Gamma _u </math>

| style="width: 5px;text-align: right;white-space: nowrap;" | (7)

|}

|}

​

where <math display="inline">t_i^p</math> and <math display="inline">u_i^p</math> are prescribed tractions and the prescribed displacement at the Neuman and Dirichlet boundaries <math display="inline">\Gamma _t</math> and <math display="inline">\Gamma _u</math>, respectively.

​

We will choose an equal order finite element (FE) approximation for the velocity and pressure

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>u_i\simeq \hat u_i =  \sum \limits _{j=1}^n N_j (\hat u_i)_j = {\bf N}\hat{\bf u}_i </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" |  (8a)

|}

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>p\simeq \hat p =   \sum \limits _{j=1}^n N_j \hat p_j = {\bf N}\hat{\bf p}</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" |  (8b)

|}

​

where <math display="inline">\hat {(\cdot )}</math> denotes approximate values, <math display="inline">N_j</math> are the shape functions, <math display="inline">(\cdot )_j</math> are nodal values and <math display="inline">n</math> is the number of nodes in the mesh <span id='citeF-1'></span>[[#cite-1|[1]]]. 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

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\int _\Omega N_k \hat r_{m_i} d\Omega +\int _{\Gamma _t} N_k (\hat \sigma _{ij} n_j - t_i^p)d\Gamma +\sum _e \int _{\Omega ^e}\tau _{m_i}^e P_{m_i} (N_k)\hat r_{m_i} d\Omega =0 </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (9)

|}

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\int _\Omega N_k \hat r_v d\Omega + \sum _e \int _{\Omega ^e}\tau _v ^e P_v (N_k)\hat r_v d\Omega =0 </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (10)

|}

​

In above <math display="inline">\hat r_{m_i}\equiv r_{m_i}(\hat u_j, \hat p)</math> and <math display="inline">\hat r_v =r_v(\hat u_j, \hat p)</math> denote the approximate finite element residuals for the momentum and mass balance equations, respectively, <math display="inline"> P_{m_i} (N_k)</math> and <math display="inline">P_v (N_k)</math> are appropriate functions of the shape functions <math display="inline">N_k</math> and <math display="inline">\tau _{m_i}^e</math> and <math display="inline">\tau _v ^e</math> 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 <math display="inline">P_{m_i}</math> and <math display="inline">P_v</math> and in the values of the stabilization parameters <span id='citeF-1'></span>[[#cite-1|[1,24-30]]].

​

A somehow different class of stabilized method developed by the authors is based on a finite calculus (FIC) formulation <span id='citeF-17'></span>[[#cite-17|[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.

​

==3 RESIDUAL POWER FORMS BASED ON THE RECOVERED BALANCE EQUATIONS==

​

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:

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\begin{array}{c}\displaystyle R_{m_i}= \bar{\hat r}_{m_i} \\  ~~\\ \displaystyle R_v=\bar{\hat r}_v \\  ~~\\ \displaystyle R_{\Gamma _i}=n_j \bar{\hat \sigma }_{ij}-t_i^p \end{array} </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (11)

|}

​

where <math display="inline">\bar{\hat{(\cdot )}}</math> denote values computed using the recovered derivative fields. In Eqs.(11) <math display="inline">R_{m_i}</math> can be interpreted as the ith component of an out-of-balance body force, whereas <math display="inline">R_v</math> represents an spurious non-zero volumetric strain rate. Similarly, the recovered residuals <math display="inline">R_{\Gamma _i}</math> are out-of-balance tractions at the Neuman boundary.

​

The ''power'' due to the recovered residuals 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>\displaystyle P_{m_i}= \int _\Omega W_i^u  R_{m_i} d\Omega + \int _{\Gamma _t} \bar W_i^u  R_{\Gamma _i} d\Gamma \qquad \hbox{no sum in }i</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" |  (12a)

|}

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\displaystyle P_v= \int _\Omega \bar W^p  R_v d\Omega </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" |  (12b)

|}

​

where the weights <math display="inline">(W_i^u,\bar W_i^u)</math> and <math display="inline">W^p</math> are functions of the approximate velocity vector field and the pressure, respectively. Indeed the simplest choice is <math display="inline">W_i^u=\bar W_i^u =\hat u_i</math> and <math display="inline">W^p =\hat p</math>. However, many other options are possible. In this work we have chosen

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>W_i^u=\hat u_i + {h_j\over 2} {\partial \hat u_i\over \partial x_j}\quad ,\quad \bar W_i^u  = 0\quad \hbox{and} \quad  W^p =\hat p + {h_j\over 2} {\partial \hat p\over \partial x_j} </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (13)

|}

​

where <math display="inline">h_j</math> 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 <math display="inline">\hat {u}=0</math>  and  at the exterior boundary <math display="inline">\hat \sigma _{ij}n_j\simeq t_i</math> on <math display="inline">\Gamma _t</math>.  This justifies neglecting the  boundary integral in Eq.(12a) by choosing <math display="inline">\bar W_i^u  = 0</math>.

​

The weights <math display="inline">W_i^u</math> and <math display="inline">W^p</math> 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 <math display="inline">\hat u_i</math> and <math display="inline">\hat p</math> 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

​

{| 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 P_{m_i} = \sum \limits _{k=1}^{\bar N} [P_{m_i}]_k= \sum \limits _{k=1}^{\bar N} \left|\left[W_i^u R_{m_i}\right]_k \right|\Omega _k \quad \hbox{no sum in } i\\ \\ \displaystyle P_v= \sum \limits _{k=1}^{\bar N} [P_v]_k= \sum \limits _{k=1}^{\bar N} \left|\left[W_i^p  R_v\right]_k\right|\Omega _k  \end{array} </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (14)

|}

​

where <math display="inline">[\cdot ]_k</math> denotes values computed at node <math display="inline">k</math>, <math display="inline">\bar N</math> is the number of ''internal'' nodes in the mesh and <math display="inline">\Omega _k </math> is the tributary area of node <math display="inline">k</math> (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

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\Vert P\Vert _\Omega = \sum \limits _{k=1}^{\bar N} \Vert P\Vert _{\Omega ^k}  </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (15)

|}

​

where <math display="inline">\Vert P\Vert _{\Omega ^k} </math> is the estimated error value for an individual node defined as

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\Vert P\Vert _{\Omega ^k} = \sum \limits _{i=1}^{n_d} [P_{m_i}]_k +[P_v]_k </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (16)

|}

​

===Remark===

​

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 <math display="inline">P_{m_i}</math> and <math display="inline">P_v</math>, 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 <span id='citeF-15'></span>[[#cite-15|[15]]] and references herein.

​

​

==4 RECOVERY OF THE NODAL DERIVATIVES OF THE VELOCITIES AND THE PRESSURE==

​

The recovered derivative field <math display="inline">\left(\displaystyle \overline{\partial \hat u_i\over \partial x_j}\right)</math> is defined so as to satisfy the following variational form

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\int _\Omega N_k \left[\left(\overline{\partial \hat u_i\over \partial x_j}\right)- {\partial \hat u_i\over \partial x_j}\right]d\Omega =0\quad \hbox{no sum in }i </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (17)

|}

​

where <math display="inline">N_k</math> are the same shape functions used to approximate the velocity field and <math display="inline">\displaystyle{\partial \hat u_i\over \partial x_j}</math> 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.

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\left(\overline{\partial \hat u_i\over \partial x_j}\right)=\sum \limits _k N_k \left(\overline{\partial \hat u_i\over \partial x_j}\right)_k </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (18)

|}

​

Substituting Eq.(18) into (17) gives the following system of equations from which the recovered nodal derivative values can be obtained

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>{\bf M} \bar{\hat {\bf d}}_{ij}= {\bf f}_{ij} </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (19)

|}

​

with

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\bar{\hat {\bf d}}_{ij}= \left[\left(\overline{\partial \hat u_i\over \partial x_j}\right)_1,\left(\overline{\partial \hat u_i\over \partial x_j}\right)_2,\cdots , \left(\overline{\partial \hat u_i\over \partial x_j}\right)_N\right]^T </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (20)

|}

​

where <math display="inline">N</math> is the number of nodes in the mesh.

​

The element contributions to <math display="inline">{\bf M}</math> and <math display="inline">{\bf f}</math> are

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>M_{kl}^e =\int _{\Omega ^e} N_k N_l d\Omega \quad , \quad (f_{ij}^e)_k = \int _{\Omega ^e} N_k {\partial \hat u_i\over \partial x_j} d\Omega  </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (21)

|}

​

We recall that in above <math display="inline">\displaystyle{\partial \hat u_i\over \partial x_j}</math> is the (discontinuous) derivative field obtained from the finite element solution and <math display="inline">\left(\displaystyle \overline{\partial \hat u_i\over \partial x_j}\right)_k</math> are the recovered derivative values at node <math display="inline">k</math>. The process is repeated for <math display="inline">i,j=1,n_d</math>.

​

The full form of matrix <math display="inline">{\bf M}</math> 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 <math display="inline">\left(\displaystyle{\partial ^2 \hat u_k\over \partial x_i\partial x_j}\right)</math> needed for computation of all the terms in the recovered momentum equations <math display="inline">R_{m_i}</math>.

​

==5 ERROR INDICATORS AND MESH ADAPTATION PROCEDURES==

​

A ''global error parameter'' is defined as

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\xi _g = {\Vert P\Vert _\Omega \over \mu U} </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (22)

|}

​

where <math display="inline">\mu </math> is the prescribed global error value and <math display="inline">U</math> is the total power in the mesh defined here as

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>U= \int _\Omega \left[\tau _{ij} \varepsilon _{ij} + {\partial p\over \partial x_i}u_i\right]d\Omega  </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (23)

|}

​

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 <math display="inline">\xi _g =1</math>) 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.

​

===5.1 Mesh adaptation based on the equidistribution of the global error===

​

A ''local error parameter'' <math display="inline"> \xi _g</math> is defined for each element as

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\xi ^e= {\Vert P\Vert _{\Omega ^e} N\over \Vert P\Vert _\Omega } </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (24)

|}

​

where <math display="inline">\Vert P\Vert _{\Omega ^e}</math> is the residual power error for element <math display="inline">e</math> and <math display="inline">N</math> 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 <math display="inline">\xi ^e= 1</math> for all elements then the total residual power error is uniformly distributed in the mesh <span id='citeF-4'></span>[[#cite-4|[4]]].

​

The error indicator for the element is defined as

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\gamma ^e = (\xi ^e)^\alpha (\xi _g)^\beta  </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (25)

|}

​

Parameters <math display="inline">\alpha </math> and <math display="inline">\beta </math> are found by analyzing the convergence rate of the global and local error parameters as described in <span id='citeF-4'></span>[[#cite-4|[4]]]. In the examples presented in the paper the following values of <math display="inline">\alpha </math> and <math display="inline">\beta </math> have been used

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\alpha = {1\over p+n_d}\quad ,\quad \beta = {1\over p+1} </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (26)

|}

​

where <math display="inline">p</math> is the order of the finite element interpolation (i.e. <math display="inline">p=1</math> for linear elements).

​

The new element size <math display="inline">h_{new}</math> is computed in terms of the error indicator <math display="inline">\gamma ^e</math> and the current element size as

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>h_{new}= h_{old} [\gamma ^e]^{-1} </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (27)

|}

​

===5.2 Mesh adaptation based on the equidistribution of the error density===

​

An alternative local error parameter based on the equidistribution of the density of the residual power can be defined as

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\xi ^e= {\Omega \Vert P\Vert _{\Omega ^e} \over \Omega ^e \Vert P\Vert _\Omega } </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (28)

|}

​

where <math display="inline">\Omega </math> and <math display="inline">\Omega ^e</math> are the area of the analysis domain and the area of an individual element, respectively. Now <math display="inline">\xi ^e = 1</math> 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 <math display="inline">\alpha </math> and <math display="inline">\beta </math> are now chosen as deduced following the procedure described in <span id='citeF-4'></span>[[#cite-4|[4]]]

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\alpha = {1\over p} \quad ,\quad \beta = {1\over p} </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (29)

|}

​

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 <span id='citeF-4'></span>[[#cite-4|[4]]]. 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 <math display="inline">\xi ^e=1</math>. 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 <span id='citeF-4'></span>[[#cite-4|[4]]].

​

===5.3 Mesh refinement strategy===

​

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 <math display="inline">\xi _g</math> and <math display="inline">\xi ^e</math>. 

​

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

​

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.

​

==6 STABILIZED FEM FOR THE NAVIER-STOKES EQUATIONS USING FINITE CALCULUS==

​

===6.1 Finite calculus form of the Navier-Stokes equations===

​

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 <span id='citeF-17'></span>[[#cite-17|[17-21]]]:<br/>

​

''Momentum''

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>r_{m_i} - \underline{{1\over 2} h_{j} {\partial r_{m_i}\over \partial x_ j}}{1\over 2} h_{j} {\partial r_{m_i}\over \partial x_ j} =0 \qquad \hbox{in } \Omega  </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (30)

|}

​

''Mass balance''

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>r_v - \underline{{1\over 2} h_{j}{\partial r_v \over \partial x_j}}{1\over 2} h_{j}{\partial r_v \over \partial x_j} =0 \qquad \hbox{in } \Omega  </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (31)

|}

​

where <math display="inline">r_{m_i}</math> and <math display="inline">r_v</math> 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 <math display="inline">{h}</math> is defined as (for 2D problems)

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>{h}=\left\{\begin{matrix} h_{1}\\ h_{2}\\\end{matrix}\right\} </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (32)

|}

​

where <math display="inline">h_{1}</math> and <math display="inline">h_{2}</math> are characteristic dimensions of the finite control domain where balance of momentum and mass conservation is enforced. The components of vector <math display="inline">{h}</math> introduce the necessary stabilization along the streamline and transverse directions to the flow in the discrete problem <span id='citeF-17'></span>[[#cite-17|[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 <math display="inline">{h}\to 0</math> the standard infinitesimal form of the momentum and mass balance equations is recovered <span id='citeF-1'></span>[[#cite-1|[1]]].

​

Eqs.(30) and (31) are completed by the following FIC boundary conditions <span id='citeF-17'></span>[[#cite-17|[17]]].

​

<br/>

​

''Equilibrium of surface tractions  at the boundary'' <math display="inline">\Gamma _t</math>

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>n_j \sigma _{ij}-t_i^p + \underline{{1\over 2} h_{j}n_j r_{m_i}}{1\over 2} h_{j}n_j r_{m_i} = 0\qquad \hbox{on}~~\Gamma _t  </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (33)

|}

​

''Prescribed velocity at  the boundaries ''

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>u_t =u_t^p \qquad \hbox{on}~~\Gamma _{u_t} </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (34)

|}

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>u_n - \underline{{1\over 2}h_{i}n_ir_v}{1\over 2}h_{i}n_ir_v=u_n^p \qquad \hbox{on} ~~\Gamma _{u_n} </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (35)

|}

​

In eq.(34) <math display="inline">u_t</math> and <math display="inline">u_t^p</math> 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 <math display="inline">u_n</math> and <math display="inline">u_n^p</math> denote the velocity normal to the boundary and its prescribed value, respectively. The value of <math display="inline">u_n^p</math> is  zero on solid walls and stationary free surfaces.

​

Also in eqs.(34) and (35) <math display="inline">\Gamma _{u_t}</math> and <math display="inline">\Gamma _{u_n}</math> are the parts of the boundary <math display="inline">\Gamma </math> of <math display="inline">\Omega </math> where the tangential and normal velocities are prescribed, respectively. The Dirichlet boundary is defined as <math display="inline">\Gamma _u=\Gamma _{u_t}\cup \Gamma _{u_n}</math>.

​

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&#8211;35) can be found in <span id='citeF-17'></span>[[#cite-17|[17]]] whereas the derivation of eq.(9) is shown in <span id='citeF-18'></span>[[#cite-18|[18]]].

​

===6.2 Residual power forms of FIC equations===

​

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'' <span id='citeF-19'></span>[[#cite-19|[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):

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\begin{array}{c}\displaystyle R_{m_i}= \bar{\hat r}_{m_i} - {h_{j} \over 2} {\partial \bar{\hat r}_{m_i}\over \partial x_j}\\  \\ \displaystyle R_v=\bar{\hat r}_v - {h_{j}\over 2} {\partial \bar{\hat r}_v \over \partial x_j}  \end{array} </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (36)

|}

​

where once more <math display="inline">\bar{\hat{(\cdot )}}</math> denote values computed using the recovered derivative fields.

​

The total  power due to the FIC residuals 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>\begin{array}{c}\displaystyle P_{m_i}= \int _\Omega \hat u_i R_{m_i} d\Omega + \int _{\Gamma _t} \hat u_i  \left(n_j \bar{\hat \sigma }_{ij} - t_i + {h_{j}\over 2} n_j \bar{\hat r}_{m_i}\right)d\Gamma \\ \\ \displaystyle P_v= \int _\Omega \hat p R_v d\Omega - \int _{\Gamma _{u_n}} \hat p \left(\hat u_n -u_n^p - {h_{j}\over 2}n_j  \bar{\hat r}_v\right)d\Gamma  \end{array} </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (37)

|}

​

Substituting the expressions of <math display="inline">R_{m_i}</math> and <math display="inline">R_v</math> of Eqs.(36) into (37) and integrating by parts some derivative terms gives

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>P_{m_i}= \int _\Omega \left[\hat u_i + {h_j\over 2}\left({\partial \hat u_i \over  \partial x_j} \right)\right]\bar{\hat r}_{m_i} d\Omega +  \int _{\Gamma _t} \hat u_i \left(\bar{\hat \sigma }_{ij}n_j  -t_i\right)d\Gamma  - \int _{\Gamma _u}  {h_j\over 2} n_j \hat u_i \bar{\hat r}_{m_i} d\Gamma \quad \hbox{no sum in } i</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" |  (38a)

|}

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>P_v= \int _\Omega \left[\hat p  +{h_j\over 2} \left({\partial \hat p\over \partial x_j}\right)\right] \bar{\hat r}_v d\Omega - \int _{\Gamma _{u_n}} \hat p (\hat u_n - u_n^p)d\Gamma - \int _{\Gamma -\Gamma _{u_n}}\hat p {h_j\over 2}n_j \bar{\hat r}_v d\Gamma </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" |  (38b)

|}

​

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 <math display="inline">\hat u_j=0</math> and hence the line integrals in Eq.(38a) vanish there. Note also that at the exterior boundary <math display="inline">\hat \sigma _{ij}n_j\simeq t_i</math> on <math display="inline">\Gamma _t</math> and <math display="inline">r_{m_i}\simeq 0</math> on <math display="inline">\Gamma _u</math> and this justifies neglecting the line integral at those boundaries. In Eq.(38b) the first line integral vanishes if <math display="inline">\hat u_n</math> is made equal to the prescribed normal velocity <math display="inline">u_n^p</math> when solving the system of discretized equations.  The second integral in Eq.(38b) can be neglected by accepting that either <math display="inline">p=0</math> (external boundary) or <math display="inline">\bar{\hat r}_v \simeq 0</math> at the <math display="inline">\Gamma - \Gamma _{u_n}</math> boundary.

​

In summary, the relevant residual power terms for the FIC momentum and mass balance equations are

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\displaystyle P_{m_i}= \int _\Omega \left[\hat u_i + {h_j\over 2} \left({\partial \hat u_i \over  \partial x_j}\right)\right]\bar{\hat r}_{m_i} d\Omega \qquad \hbox{no sum in }i</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" |  (39a)

|}

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;width: 100%;" 

|-

| style="text-align: center;" | <math>\displaystyle P_v= \int _\Omega \left[\hat p + {h_j\over 2} \left({\partial \hat p\over \partial x_j}\right)\right]\bar{\hat r}_v  d\Omega </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" |  (39b)

|}

​

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

​

{| 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 P_{m_i}= \sum \limits _{k=1}^{\bar N}\left|\left[\left(\hat u_i + {h_j\over 2}\left(\overline{\partial \hat u_i \over  \partial x_j}\right)\right)\bar{\hat r}_{m_i}\right]_k\right|\Omega _k\\ \\ \displaystyle P_v= \sum \limits _{k=1}^{\bar N} \left|\left[\left(\hat p  +{h_j\over 2} \left(\overline{\partial \hat p\over \partial x_j}\right)\right)\bar{\hat r}_v \right]_k \right|\Omega _k  \end{array} </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (40)

|}

​

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

​

<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_249588187-figura0.png|300px|Tributary domain for a node k and definition of characteristic lengths.]]

|- style="text-align: center; font-size: 75%;"

| colspan="1" | '''Figure 1:''' Tributary domain for a node <math>k</math> 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.

​

===Remark===

​

Eqs.(40) define the weighting functions in the general residual power forms of Eqs.(12) and (13) as

​

{| 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_i^u = \hat u_i + {h_j\over 2} \left(\overline{\partial \hat u_i \over  \partial x_j}\right)} \quad , \quad \bar W_i^u  = 0 \\ \\ \displaystyle{W^p  = \hat p +  {h_j\over 2}  \left(\overline{\partial \hat p\over \partial x_j}\right)} \end{array} </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | (41)

|}

​

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

​

We note again however that Eqs.(41) are not the only option for <math display="inline">W_i^u</math> and <math display="inline">W^p</math> and some other choices are explored in <span id='citeF-22'></span>[[#cite-22|[22]]] to derive a variety of error estimation and mesh refinement procedures for convection-diffusion problems. Numerical results reported in <span id='citeF-22'></span>[[#cite-22|[22]]] show that the expressions given by Eqs.(41) have a better performance.

​

==7 EXAMPLES==

​

Typical parameters used for the  examples presented are: prescribed global  error <math display="inline">\mu = 1%</math>, maximum and minimum sizes of the new elements generated <math display="inline">h_{max}=0.3</math>m and <math display="inline">h_{min}=0.001</math>m.

​

===7.1 Flow past a NACA airfoil===

​

This a transient problem. The mesh refinement process has been carried out at a time <math display="inline">t=2s</math>. 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.

​

<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_249588187-fig1.png|600px|Flow past a NACA airfoil. Characteristic length =1m, Re=100. Original mesh]]

|- style="text-align: center; font-size: 75%;"

| colspan="1" | '''Figure 2:''' Flow past a NACA airfoil. Characteristic length =1m, <math>Re=100</math>. 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.

​

<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_249588187-fig2.png|600px|Mesh obtained after five refinement steps using the equidistribution of the global error.]]

|- style="text-align: center; font-size: 75%;"

| colspan="1" | '''Figure 3:''' Mesh obtained after five refinement steps using the equidistribution of the global error.

|}

​

<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_249588187-fig3.png|468px|]]

|style="padding-left:10px;"|[[Image:Draft_Samper_249588187-fig3_b.png|534px|Details of the refined mesh of Figure 3 in the vecinity of the airfoil]]

|- style="text-align: center; font-size: 75%;"

| colspan="2" | '''Figure 4:''' Details of the refined mesh of Figure 3 in the vecinity of the airfoil

|}

​

<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_249588187-naca_line_1.png|600px|Velocity field after five refinement steps]]

|- style="text-align: center; font-size: 75%;"

| colspan="1" | '''Figure 5:''' Velocity field after five refinement steps

|}

​

<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_249588187-naca_line_2.png|510px|Detail of the velocity field]]

|- style="text-align: center; font-size: 75%;"

| colspan="1" | '''Figure 6:''' Detail of the velocity field

|}

​

​

{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"

|+ style="font-size: 75%;" |Table. 1 Flow past a NACA airfoil (equidistribution of global error). Convergence of the global error with the mesh refinement. Minimum element size = 0.001m

|- style="border-top: 2px solid;"

| style="border-left: 2px solid;border-right: 2px solid;" |  Mesh 

| style="border-left: 2px solid;border-right: 2px solid;" | Global error <math display="inline">\left({\Vert P\Vert _\Omega \over U}\right)</math>(%)  

| style="border-left: 2px solid;border-right: 2px solid;" | Nodes 

| style="border-left: 2px solid;border-right: 2px solid;" | Elements

| style="border-left: 2px solid;border-right: 2px solid;" | Drag error %

| style="border-left: 2px solid;border-right: 2px solid;" | Lift error% 

|- style="border-top: 2px solid;"

| style="border-left: 2px solid;border-right: 2px solid;" |  1

| style="border-left: 2px solid;border-right: 2px solid;" | 5.55

| style="border-left: 2px solid;border-right: 2px solid;" | 2574

| style="border-left: 2px solid;border-right: 2px solid;" | 4784

| style="border-left: 2px solid;border-right: 2px solid;" | 29.1810942 

| style="border-left: 2px solid;border-right: 2px solid;" | 2.99902143

|- style="border-top: 2px solid;"

| style="border-left: 2px solid;border-right: 2px solid;" |  2

| style="border-left: 2px solid;border-right: 2px solid;" | 2.93

| style="border-left: 2px solid;border-right: 2px solid;" | 5340

| style="border-left: 2px solid;border-right: 2px solid;" | 10680

| style="border-left: 2px solid;border-right: 2px solid;" | 28.7443059

| style="border-left: 2px solid;border-right: 2px solid;" | 2.5500113

|- style="border-top: 2px solid;"

| style="border-left: 2px solid;border-right: 2px solid;" |  3

| style="border-left: 2px solid;border-right: 2px solid;" | 2.87

| style="border-left: 2px solid;border-right: 2px solid;" | 7942

| style="border-left: 2px solid;border-right: 2px solid;" | 15884

| style="border-left: 2px solid;border-right: 2px solid;" | 22.1078024

| style="border-left: 2px solid;border-right: 2px solid;" | 2.00851784

|- style="border-top: 2px solid;"

| style="border-left: 2px solid;border-right: 2px solid;" |  4

| style="border-left: 2px solid;border-right: 2px solid;" | 2.77

| style="border-left: 2px solid;border-right: 2px solid;" | 11948

| style="border-left: 2px solid;border-right: 2px solid;" | 23896

| style="border-left: 2px solid;border-right: 2px solid;" | 15.2893198

| style="border-left: 2px solid;border-right: 2px solid;" | 1.54898957

|- style="border-top: 2px solid;"

| style="border-left: 2px solid;border-right: 2px solid;" |  5

| style="border-left: 2px solid;border-right: 2px solid;" | 1.81

| style="border-left: 2px solid;border-right: 2px solid;" | 17142

| style="border-left: 2px solid;border-right: 2px solid;" | 34284

| style="border-left: 2px solid;border-right: 2px solid;" | 10.6478917

| style="border-left: 2px solid;border-right: 2px solid;" | 1.17033653

|- style="border-top: 2px solid;"

| style="border-left: 2px solid;border-right: 2px solid;" |  6

| style="border-left: 2px solid;border-right: 2px solid;" | 1.81

| style="border-left: 2px solid;border-right: 2px solid;" | 17255

| style="border-left: 2px solid;border-right: 2px solid;" | 34510

| style="border-left: 2px solid;border-right: 2px solid;" | 8.54878424

| style="border-left: 2px solid;border-right: 2px solid;" | 1.02567233

|- style="border-top: 2px solid;"

| style="border-left: 2px solid;border-right: 2px solid;" |  7

| style="border-left: 2px solid;border-right: 2px solid;" | 1.71

| style="border-left: 2px solid;border-right: 2px solid;" | 17372

| style="border-left: 2px solid;border-right: 2px solid;" | 34744

| style="border-left: 2px solid;border-right: 2px solid;" | 6.93290243

| style="border-left: 2px solid;border-right: 2px solid;" | 0.86514698

|- style="border-top: 2px solid;border-bottom: 2px solid;"

| style="border-left: 2px solid;border-right: 2px solid;" |  8

| style="border-left: 2px solid;border-right: 2px solid;" | 1.72

| style="border-left: 2px solid;border-right: 2px solid;" | 16729

| style="border-left: 2px solid;border-right: 2px solid;" | 33458

| style="border-left: 2px solid;border-right: 2px solid;" | 6.40377809

| style="border-left: 2px solid;border-right: 2px solid;" | 0.97782406

​

|}

​

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., <math display="inline">h_{max}=0.3</math>m and <math display="inline">h_{min}=0.001</math>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.

​

<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_249588187-fig6.png|600px|Mesh obtained after two refinement steps using the equidistribution of the error density (min size = 0.001m)]]

|- style="text-align: center; font-size: 75%;"

| colspan="1" | '''Figure 7:''' Mesh obtained after two refinement steps using the equidistribution of the error density (min size = 0.001m)

|}

​

<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_249588187-fig7.png|474px|]]

|style="padding-left:10px;"|[[Image:Draft_Samper_249588187-fig7_b.png|534px|Detail of the refined mesh shown in the  previous figure]]

|- style="text-align: center; font-size: 75%;"

| colspan="2" | '''Figure 8:''' Detail of the refined mesh shown in the  previous figure

|}

​

​

{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"

|+ style="font-size: 75%;" |Table. 2 Flow past a NACA profile (equidistribution of error density). Convergence of the global error with the mesh refinement (min. element size = 0.001m)

|- style="border-top: 2px solid;"

| style="border-left: 2px solid;border-right: 2px solid;" |  Mesh 

| style="border-left: 2px solid;border-right: 2px solid;" | Global error <math display="inline">\left({\Vert P\Vert _\Omega \over U}\right)</math>(%)  

| style="border-left: 2px solid;border-right: 2px solid;" | Nodes 

| style="border-left: 2px solid;border-right: 2px solid;" | Elements

| style="border-left: 2px solid;border-right: 2px solid;" | Drag error %

| style="border-left: 2px solid;border-right: 2px solid;" | Lift error% 

|- style="border-top: 2px solid;"

| style="border-left: 2px solid;border-right: 2px solid;" |  1

| style="border-left: 2px solid;border-right: 2px solid;" | 5.55

| style="border-left: 2px solid;border-right: 2px solid;" | 2574

| style="border-left: 2px solid;border-right: 2px solid;" | 4784

| style="border-left: 2px solid;border-right: 2px solid;" | 29.1810942 

| style="border-left: 2px solid;border-right: 2px solid;" | 201.915023

|- style="border-top: 2px solid;"

| style="border-left: 2px solid;border-right: 2px solid;" |  2

| style="border-left: 2px solid;border-right: 2px solid;" | 2.18

| style="border-left: 2px solid;border-right: 2px solid;" | 10594

| style="border-left: 2px solid;border-right: 2px solid;" | 21188

| style="border-left: 2px solid;border-right: 2px solid;" | 10.4796861

| style="border-left: 2px solid;border-right: 2px solid;" | 1.16034921

|- style="border-top: 2px solid;border-bottom: 2px solid;"

| style="border-left: 2px solid;border-right: 2px solid;" |  3

| style="border-left: 2px solid;border-right: 2px solid;" | 2.14

| style="border-left: 2px solid;border-right: 2px solid;" | 39895

| style="border-left: 2px solid;border-right: 2px solid;" | 79790

| style="border-left: 2px solid;border-right: 2px solid;" | 3,17339951

| style="border-left: 2px solid;border-right: 2px solid;" | 0.5988364

​

|}

​

<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_249588187-Fig8.png|600px|Mesh obtained after four refinement steps using the equidistribution of the error density (min size = 0.002m).]]

|- style="text-align: center; font-size: 75%;"

| colspan="1" | '''Figure 9:''' Mesh obtained after four refinement steps using the equidistribution of the error density (min size = 0.002m).

|}

​

<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_249588187-Fig9.png|480px|]]

|style="padding-left:10px;"|[[Image:Draft_Samper_249588187-Fig9_b.png|540px|Detail of the refined mesh of Figure 9.]]

|- style="text-align: center; font-size: 75%;"

| colspan="2" | '''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.

​

​

{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"

|+ style="font-size: 75%;" |Table. 3 Flow past a NACA airfoil (equidistribution of error density). Convergence of the global error with the mesh refinement. Minimum  element size = 0.002m

|- style="border-top: 2px solid;"

| style="border-left: 2px solid;border-right: 2px solid;" |  Mesh 

| style="border-left: 2px solid;border-right: 2px solid;" | Global error <math display="inline">\left({\Vert P\Vert _\Omega \over U}\right)</math>(%)  

| style="text-align: right;border-left: 2px solid;border-right: 2px solid;" | Nodes 

| style="text-align: right;border-left: 2px solid;border-right: 2px solid;" | Elements

| style="border-left: 2px solid;border-right: 2px solid;" | Drag error %

| style="border-left: 2px solid;border-right: 2px solid;" | Lift error% 

|- style="border-top: 2px solid;"

| style="border-left: 2px solid;border-right: 2px solid;" |  1

| style="border-left: 2px solid;border-right: 2px solid;" | 6.38

| style="text-align: right;border-left: 2px solid;border-right: 2px solid;" | 2574 

| style="text-align: right;border-left: 2px solid;border-right: 2px solid;" | 4784 

| style="border-left: 2px solid;border-right: 2px solid;" | 29.1810942

| style="border-left: 2px solid;border-right: 2px solid;" | 2.99902143

|- style="border-top: 2px solid;"

| style="border-left: 2px solid;border-right: 2px solid;" |  2 

| style="border-left: 2px solid;border-right: 2px solid;" | 3.42  

| style="text-align: right;border-left: 2px solid;border-right: 2px solid;" | 9871 

| style="text-align: right;border-left: 2px solid;border-right: 2px solid;" | 19742