Finite element formulation for convective-diffusive problems with sharp gradients using finite calculus

Latest revision as of 11:15, 12 November 2018

Published in Comput. Meth. Appl. Mech. Engng., Vol. 195 (13-16), pp. 1793–1825, 2006
doi: 10.1016/j.cma.2005.05.036

Abstract

A finite element method (FEM) for steady-state convective-diffusive problems presenting sharp gradients of the solution both in the interior of the domain and in boundary layers is presented. The necessary stabilization of the numerical solution is provided by the Finite Calculus (FIC) approach. The FIC method is based in the solution by the Galerkin FEM of a modified set of governing equations which include characteristic length parameters. It is shown that the FIC balance equation for the multidimensional convection-diffusion problem written in the principal curvature axes of the solution, introduces an orthotropic diffusion which stabilizes the numerical solution both in smooth regions as well in the vicinity of sharp gradients. The dependence of the stabilization terms with the principal curvature directions of the solution makes the method non linear. Details of the iterative scheme to obtain stabilized results are presented together with examples of application which show the efficiency and accuracy of the approach.

1 INTRODUCTION

It is well known that the standard Galerkin FEM solution of the steady-state convective-diffusive equation is unstable for values of the Peclet number greater than one (see Volume 3 in [1]). A number of numerical schemes have been proposed in order to guarantee that the numerical solution is stable, that is, that the solution has a physical meaning. In the first attempts to solve this problem the underdiffusive character of the Galerkin FEM (and the analogous central finite difference scheme) for convective-diffusive problems was corrected by adding “artificial diffusion terms” to the governing equations [1,2]. The relationship of this approach with the upwind finite difference method [2] lead to the derivation of a variety of Petrov-Galerkin FEM. All these methods can be interpreted as extensions of the standard Galerkin variational form of the FEM by adding residual-based integral terms computed over the element domains. Among the many stabilization methods of this or similar kind we name the Upwind FEM [3,4], the Streamline Upwind Petrov-Galerkin (SUPG) method [5-9], the Taylor-Galerkin method [10,11, the generalized Galerkin method [12,13], the Galerkin Least Square method and related approaches [14-16], the Characteristic Galerkin method [17,18], the Characteristic Based Split method [19], the Subgrid Scale method [20-23], the Residual Free Bubbles method [24], the Discontinuous Enrichment Method [25] and the Streamline Upwind with Boundary Terms method [26]. Basically all the methods make use of a single stabilization parameter which suffices to stabilize the numerical solution along the velocity (streamline) direction. The computation of the streamline stabilization parameter for multidimensional problems is usually based on extensions of the optimal value of the parameter for the simpler 1D case. Specific attempts to design the stability parameter for multidimensional problems in the context of the Petrov-Galerkin formulation have been recently reported [27-31]. In general all the stabilized methods above mentioned yield good results for 1D-type problems where the velocity vector is aligned with the direction of the gradient of the solution. However, the use of a single stabilization parameter is insufficient to provide stabilized solutions in the vicinity of sharp gradients not aligned with the velocity direction, which may appear at the interior of the domain or at boundary layers. The usual remedy for these situations is to use the so called “shock capturing” or “discontinuity-capturing” schemes [32-36]. These methods basically add additional transverse diffusion terms in selected elements in order to correct the undershoots and overshoots yielded by the Petrov-Galerkin solution in the sharp gradient zones. Usually the new shock capturing diffusion terms depend on the gradient of the solution and the scheme becomes non linear.

The aim of this work is to develop a general finite element formulation which can provide stabilized numerical results for all range of convective-diffusive problems. The new formulation therefore has the necessary intrinsic features to deal with streamline-type instabilities, as well as with the typical undershoots and overshoots in the vicinity of sharp gradients at different angles with the velocity direction. The new formulation thus provides a unified theoretical and computational framework incorporating all the ingredients of the traditional Petrov-Galerkin and shock-capturing methods.

The formulation proposed is based in the so called Finite Calculus (FIC) approach [37,38]. The FIC method is based in expressing the equation of balance of fluxes in a domain of finite size. This introduces additional stabilizing terms in the differential equations of the infinitessimal theory which are a function of the balance domain dimensions. These dimensions are termed characteristic length parameters and play a key role in the stabilization process.

The modified governing equations lead to stabilized numerical schemes using whatever numerical method. It is interesting that many of the stabilized FEM can be recovered using the FIC formulation. The FIC method has been successfully applied to problems of convection-diffusion [37-41], convection-diffusion-absorption [42,43], incompressible fluid flow [44-47] and incompressible solid mechanics [48,49].

The key to the stabilization in the FIC method is to choose adequately the characteristic length parameters for each element. For 1D problems the standard optimal and critical values of the single stabilization parameter can be found, as shown in Section 3. For 2D/3D problems the characteristic length parameters along each space direction are grouped in a characteristic length vector ${\textstyle h}$ . Here the key issue is to choose correctly the direction of this vector. It is interesting that if the direction of ${\textstyle h}$ is taken parallel to that of the velocity, then the resulting FIC stabilized equations coincide with that of the standard SUPG method [37,38].

In previous work of the authors, transverse sharp gradients to the velocity direction were accurately captured by computing iteratively the direction and modulus of vector ${\textstyle h}$ which minimizes a residual norm [37-41]. In [38] accurate stabilized solutions for convection-diffusion problems with sharp gradients were obtaining by spliting the characteristic lenght vector ${\textstyle {h}}$ as the sum of two vectors. The first vector is chosen aligned with the velocity direction, hence yielding the standard SUPG terms, while the second vector is taken parallel to the direction of the gradient of the solution.

In this work a slightly different and more consistent method is proposed. We have found that the key to the general stabilization algorithm for convection-diffusion problems is to express the FIC balance equation taking as coordinate axes the principal curvature directions of the solution. These equations contain the necessary additional diffusion to stabilize the numerical solution in all situations. It is interesting that when one of the principal curvature directions coincides with the velocity direction, then the classical SUPG method is recovered.

The content of the paper is organized as follows. In the next section, the basic FIC equations for the multidimensional steady-state convective-diffusive problem are given. The general equation is particularized for the simple 1D problem and the optimal and critical values of the 1D stabilization parameter are obtained. Then the form of the FIC balance equation written in the principal curvature axes is presented. The computation of the stabilization parameters is detailed and the general iterative solution scheme is explained. The formulation is particularized for linear elements (four node quadrilaterals and three node triangles) where the computation of the curvature directions requires a derivative recovery algorithm. Very good results have been obtained by approximating the main principal curvature direction within the element by the gradient direction at the element center, which is a much more economical approach. A number of examples showing the accuracy and efficiency of the method are presented. Convergence of the iterative scheme in problems with arbitrary sharp gradients was found in most cases in just two iterations.

2 FIC GOVERNING EQUATIONS FOR STEADY-STATE CONVECTION-DIFFUSIVE PROBLEMS

The basic FIC equations are obtained by expressing the balance of fluxes in interior and boundary domains of finite size and retaining one order higher terms in the Taylor expansions than those retained in the infinitessimal theory. The resulting governing equations are

(1)

with the boundary conditions

(2a)

(2b)

where ${\textstyle \Gamma _{\phi }}$ and ${\textstyle \Gamma _{q}}$ are the Dirichlet and Neumann boundaries where the variable ${\textstyle \phi }$ and the outgoing normal diffusive flux are prescribed to values ${\textstyle \phi ^{p}}$ and ${\textstyle q_{n}^{p}}$ , respectively. The modified equation (2b) is obtained by invoking higher order balance of fluxes in a finite domain next to the Neumann boundary [37,38]. In above equations

(3)

where ${\textstyle u}$ is the divergence-free velocity vector, ${\textstyle {\textbf {D}}}$ is the diffusion matrix, ${\textstyle {\boldsymbol {\nabla }}}$ is the gradient operator, ${\textstyle Q}$ is the external source term and ${\textstyle {\textbf {n}}}$ is the normal vector. Vector ${\textstyle {\textbf {h}}}$ is the characteristic length vector. For 2D problems ${\textstyle {\textbf {h}}=[h_{x},h_{y}]^{T}}$ , where ${\textstyle h_{x}}$ and ${\textstyle h_{y}}$ are characteristic distances along the sides of the rectangular domain where higher order balance of fluxes is enforced [37,38].

The underlined terms in Equations (1) and (2b) introduce the necessary stabilization in the numerical solution of the convective-diffusive problem [37-41]. The finite element formulation will be presented next.

2.1 Finite element discretization

A finite element interpolation of the unknown ${\textstyle \phi }$ can be written as

(4)

where ${\textstyle N_{i}}$ are the shape functions and ${\textstyle \phi _{i}}$ are the nodal values of the approximate function ${\textstyle {\hat {\phi }}}$ [1].

Application of the Galerkin FE method to Equations (1) and (2b) gives, after integrating by parts the term ${\textstyle {h}^{T}{\boldsymbol {\nabla }}{\hat {r}}}$

(5)

The last integral in Equation (5) has been expressed as a sum of the element contributions to allow for interelement discontinuities in the term ${\textstyle {\boldsymbol {\nabla }}{\hat {r}}}$ of Eq.(1), where ${\textstyle {\hat {r}}=r({\hat {\phi }})}$ is the residual of the FE approximation of the infinitesimal governing equation.

Note that the residual terms have disappeared from the Neumann boundary ${\textstyle \Gamma _{q}}$ . This is due to the consistency between the FIC terms in Eqs.(1) and (2b).

Integrating by parts the diffusive terms in the first integral of Eq.(5) leads to

\int _{\Omega }N_{i}[{\textbf {u}}^{T}\nabla {\hat {\phi }}+{\boldsymbol {\nabla }}^{T}N_{i}{\textbf {D}}{\boldsymbol {\nabla }}{\hat {\phi }}]d\Omega -\sum \limits _{e}{1 \over 2}\int _{\Omega ^{e}}({\textbf {h}}^{T}{\boldsymbol {\nabla }}N_{i}+N_{i}{\boldsymbol {\nabla }}^{T}{\textbf {h}}){\hat {r}}d\Omega -\int _{\Omega }N_{i}Qd\Omega +\int _{\Gamma _{q}}N_{i}q_{n}^{p}d\Gamma =0

(6)

In matrix form

(7)

Matrix ${\textbf {K}}$ and vector ${\textbf {f}}$ are assembled from the element contributions given by

$K_{ij}^{e}=\int _{\Omega ^{e}}[N_{i}{\textbf {u}}^{T}{\boldsymbol {\nabla }}N_{j}+{\boldsymbol {\nabla }}^{T}N_{i}({\textbf {D}}+{1 \over 2}{\textbf {h}}{\textbf {u}}^{T}){\boldsymbol {\nabla }}N_{j}+N_{i}({\boldsymbol {\nabla }}^{T}{\textbf {h}}){\textbf {u}}^{T}{\boldsymbol {\nabla }}N_{j}]d\Omega -$
$-{1 \over 2}\int _{\Omega ^{e}}({\textbf {h}}^{T}{\boldsymbol {\nabla }}N_{i}+N_{i}{\boldsymbol {\nabla }}^{T}{\textbf {h}}){\boldsymbol {\nabla }}^{T}({D}{\boldsymbol {\nabla }}N_{j})d\Omega$	(8.a)
$f_{i}^{e}=\int _{\Omega ^{e}}[N_{i}+{1 \over 2}({\textbf {h}}^{T}{\boldsymbol {\nabla }}N_{i}+N_{i}{\boldsymbol {\nabla }}^{T}{\textbf {h}})]Qd\Omega -\int _{\Gamma _{q}^{e}}N_{i}q_{n}^{p}d\Gamma$	(8.b)

Note that the method is equivalent to modifying the original diffusion matrix ${\textstyle D}$ by ${\textstyle {D}_{G}}$ with

(9)

In Eq.(9) ${\textstyle {\bar {\textbf {D}}}}$ is the balancing diffusion matrix introduced by the FIC method.

We note that the Galerkin/FIC formulation is residual based, i.e. the stabilization term is a function of the FEM residual which progressively vanishes as the numerical results approaches the “exact” solution (see Eq.(5)). This preserves the consistency of the Galerkin method in order to obtain improved convergence rates with finite elements of any interpolation order. In the following we will restrict the application of the method to linear elements only.

Note that when linear elements are used the second integral of Equation (8a) vanishes. The same happens with the term ${\textstyle {\textbf {h}}^{T}{\boldsymbol {\nabla }}N_{i}Q}$ of the first integral of Equation (8b) when ${\textstyle N_{i}}$ is linear and ${\textstyle Q}$ is constant. Finally, all terms involving the derivatives of ${\textstyle {\textbf {h}}}$ vanish if ${\textstyle {\textbf {h}}}$ is constant over the element. The evaluation of these integrals is mandatory in any other case.

2.2 Equivalence with the SUPG method

Let us now assume that the direction of vector ${\textstyle h}$ is parallel to that of the velocity ${\textstyle {\textbf {u}}}$ , i.e. ${\textstyle {\textbf {h}}=h{{\textbf {u}} \over \vert {\textbf {u}}\vert }}$ where ${\textstyle h}$ is a characteristic length. Under these conditions, Eq.(5) reads (assuming ${\textstyle {\boldsymbol {\nabla }}^{T}{\textbf {h}}=0}$ )

(10)

Equation (10) coincides precisely with the SUPG method. The ratio ${\textstyle \displaystyle {h \over 2\vert {\textbf {u}}\vert }}$ has dimensions of time and it is usually termed element intrinsic time parameter ${\textstyle \tau }$ .

The balancing diffusion matrix ${\textstyle {\bar {\textbf {D}}}}$ of Eq.(9) is now given by

(11)

It can be shown that the definition of ${\textstyle {\bar {\textbf {D}}}}$ of Eq.(11) is equivalent to introducing an artificial diffusion of value ${\textstyle {h\vert {\textbf {u}}\vert \over 2}}$ along the streamlines [1,7,9] and this explains the name of the SUPG method. The element characteristic length ${\textstyle h}$ (or the value of ${\textstyle \tau }$ ) is computed in practice by heuristic linear and non-linear extensions of the optimal expression for the 1D problem [5-9,27-31].

In the following sections we will denote by the standard SUPG method to the stabilization procedure leading to a balancing diffusion matrix ${\textstyle {\bar {\textbf {D}}}}$ given by Eq.(11).

It is important to note that the SUPG expression is a particular case of the more general FIC formulation. This explains the limitations of the SUPG method to provide stabilized numerical results in the vicinity of sharp gradients of the solution transverse to the flow direction and the need to introduce in these cases additional shock capturing terms (see Section 5). However, in the FIC formulation the direction of ${\textstyle {\textbf {h}}}$ is arbitrary and not necessarily coincident with that of ${\textstyle {\textbf {u}}}$ . The components of ${\textstyle {\textbf {h}}}$ introduce the necessary stabilization along the streamlines and the transverse directions to the flow. In this manner, the FIC method reproduces the best-features of the so-called stabilized discontinuity-capturing schemes [32-36].

3 COMPUTATION OF THE CHARACTERISTIC LENGTH FOR THE 1D CONVECTION-DIFFUSION EQUATION

The computation of the characteristic length is a crucial step as its value affect to the stability and accuracy of the numerical solution.

Before we face the general multidimensional problem, some based concepts of the computation of the characteristic length for the 1D convection-diffusion equation are briefly given.

The FIC equation for the 1D steady-state convection-diffusion problem is written as

(12)

with

(13)

where primes denote differentiation with respect to the independent space variable ${\textstyle x}$ .

The element matrices of Eqs.(8) are now written for two node linear elements as

$K_{ij}^{e}=\int _{l^{e}}\left(N_{i}uN_{j}^{\prime }+N_{i}^{\prime }\left(k+{hu \over 2}\right)N_{j}^{\prime }\right)dx$	(14.a)
$f_{i}^{e}=\int _{l^{e}}\left(N_{i}+{h \over 2}N_{j}^{\prime }\right)Qdx-[N_{i}q^{p}]_{x_{L}}$	(14.b)

where ${\textstyle l^{e}}$ is the element length and ${\textstyle x_{L}}$ is the coordinate of the end point of the 1D domain where the outgoing flux is prescribed to a value ${\textstyle q^{p}}$ . In Eq.(14a) the space derivatives of ${\textstyle h}$ have been neglected as ${\textstyle h}$ is assumed to be constant within each element.

Note that the FIC method introduces a diffusion of value ${\textstyle {hu \over 2}}$ . This term is analogous to that provided by artificial diffusion and upwinding methods [1,2].

The characteristic length parameter ${\textstyle h}$ for each element can be made proportional to the element length as ${\textstyle h=\alpha l^{e}}$ where ${\textstyle \alpha }$ is an element stabilization parameter. A typical stencil for a mesh of uniform size 1D elements of length ${\textstyle l^{e}=l}$ is written as

(15)

Injecting into Eq.(15) the exact analytical solution ${\textstyle \phi _{i}=A+Be^{{u \over k}x_{i}}}$ where ${\textstyle A}$ and ${\textstyle B}$ are contants, gives the optimal value of ${\textstyle \alpha =\alpha _{opt}}$ yielding the exact solution at the nodes with

(16)

where ${\textstyle \gamma ={ul \over 2k}}$ is the element Peclet number.

A critical value of ${\textstyle \alpha }$ ensuring physically correct results can be found by writing the stencil (15) for a two element mesh with nodes 1, 2 and 3 and Dirichlet boundary conditions ${\textstyle \phi _{1}=0}$ and ${\textstyle \phi _{3}={\bar {\phi }}}$ . The resulting equation reads

(17)

An unphysical (unstable) solution will imply ${\textstyle \phi _{2}<0}$ . This is avoided if the following value of ${\textstyle \alpha \geq \alpha _{c}}$ is chosen with

(18)

where ${\textstyle \alpha _{c}}$ is termed the critical stabilization parameter. It can be verified that that ${\textstyle \alpha _{c}\simeq \alpha _{opt}}$ for ${\textstyle \alpha >2}$ . Note that both ${\textstyle \alpha }$ and ${\textstyle \alpha _{c}}$ are constant within each element, as initially assumed.

Above well know results from the finite element and finite difference literature will be useful for the general multidimensional case described next.

Remark 1. The stabilizing diffusion of Eq.(15) can be expressed in terms of the equivalent intrinsic time parameter  ${\textstyle \tau }$  by replacing the term  ${\textstyle \alpha {ul \over 2}}$  by  ${\textstyle \tau u^{2}}$  where  ${\textstyle \tau =\alpha {l \over 2u}}$ . The optimal and critical values of  ${\textstyle \tau }$  are therefore obtained by dividing by  ${\textstyle {l \over 2u}}$  the expressions of  ${\textstyle \alpha _{opt}}$  and  ${\textstyle \alpha _{c}}$  of Eqs.(16) and (18), respectively.

4 COMPUTATION OF THE CHARACTERISTIC LENGTH VECTOR

We present here a general procedure to compute the characteristic length vector ${\textstyle h}$ for convection-diffusion problems. For the sake of preciseness the method is explained for 2D problems although it is equally applicable to 3D problems.

Figure 1: Global axes (

) and principal curvature axes (

)

Let us write down the FIC balance equation in the principal curvature axes of the solution ${\textstyle \xi ,\eta }$ (Figure 1). For simplicity we consider the 2D sourceless case ( ${\textstyle Q=0}$ ) with an isotropic diffusion defined by a constant diffusion parameter ${\textstyle k}$ . The FIC balance equation is

-u_{\xi }{\partial \phi  \over \partial \xi }-u_{\eta }{\partial \phi  \over \partial \eta }+k\left({\partial ^{2}\phi  \over \partial \xi ^{2}}+{\partial ^{2}\phi  \over \partial \eta ^{2}}\right)-{h_{\xi } \over 2}{\partial  \over \partial \xi }\left[-u_{\xi }{\partial \phi  \over \partial \xi }-u_{\eta }{\partial \phi  \over \partial \eta }+k\left({\partial ^{2}\phi  \over \partial \xi ^{2}}+{\partial ^{2}\phi  \over \partial \eta ^{2}}\right)\right]

(19)

where ${\textstyle u_{\xi },u_{\eta }}$ are the velocities along the principal axes of curvature ${\textstyle \xi }$ and ${\textstyle \eta }$ , respectively.

As ${\textstyle \xi }$ and ${\textstyle \eta }$ are the principal curvature axes of the solution then

(20)

Introducing this simplification into Eq.(19) we can rewrite this equation as

-u_{\xi }{\partial \phi  \over \partial \xi }-u_{\eta }{\partial \phi  \over \partial \eta }+\left(k+{u_{\xi }h_{\xi } \over 2}\right){\partial ^{2}\phi  \over \partial \xi ^{2}}+\left(k+{u_{\eta }h_{\eta } \over 2}\right){\partial ^{2}\phi  \over \partial \eta ^{2}}-{k \over 2}\left(h_{\xi }{\partial ^{3}\phi  \over \partial \xi ^{3}}+h_{\eta }{\partial ^{3}\phi  \over \partial \eta ^{3}}\right)=0

(21)

We can see clearly from Eq.(21) that the FIC governing equations introduce orthotopic diffusion parameters of values ${\textstyle {u_{\xi }h_{\xi } \over 2}}$ and ${\textstyle {u_{\eta }h_{\eta } \over 2}}$ along the ${\textstyle \xi }$ and ${\textstyle \eta }$ axes, respectively. Also note that the last term of Eq.(21) will vanish after discretization for linear elements.

Eq.(21) can be rewritten in matrix form (neglecting the last term) as

(22)

where ${\textstyle {\textbf {u}}^{\prime }=[u_{\xi },u_{\eta }]^{T}}$ , ${\textstyle {\boldsymbol {\nabla }}^{\prime }=\left[{\partial \over \partial \xi },{\partial \over \partial \eta }\right]^{T}}$ , ${\textstyle {\textbf {D}}}$ is the “physical” isotropic diffusion matrix and ${\textstyle {\bar {\textbf {D}}}'}$ is the balancing diffusion matrix in the local axes ${\textstyle \xi }$ and ${\textstyle \eta }$ . The form of these matrices is

(23)

The velocities along the principal curvature axes ${\textstyle u_{\xi }}$ and ${\textstyle u_{\eta }}$ can be obtained by projecting the cartesian velocities into the principal curvature axes ${\textstyle \xi }$ and ${\textstyle \eta }$ as

(24)

where ${\textstyle c_{\alpha }=\cos \alpha }$ , ${\textstyle s_{\alpha }=\sin \alpha }$ and ${\textstyle \alpha }$ is the angle which the ${\textstyle \xi }$ axis forms with the ${\textstyle x}$ axis (Figure 1). Note that as the solution is continuous the principal curvature directions ${\textstyle \xi }$ and ${\textstyle \eta }$ are orthogonal.

The characteristic length distances ${\textstyle h_{\xi }}$ and ${\textstyle h_{\eta }}$ are defined as

(25)

where ${\textstyle l_{\xi }}$ and ${\textstyle l_{\eta }}$ are typical element dimensions along the ${\textstyle \xi }$ and ${\textstyle \eta }$ axes, respectively and ${\textstyle \alpha _{\xi }}$ and ${\textstyle \alpha _{\eta }}$ are the corresponding stabilization parameters.

The values of ${\textstyle \alpha _{\xi }}$ and ${\textstyle \alpha _{\eta }}$ are computed by considering the solution of two uncoupled 1D problems along the ${\textstyle \xi }$ and ${\textstyle \eta }$ directions. This gives from Eq.(16)

$\alpha _{\xi }=\coth \gamma _{\xi }-{1 \over \gamma _{\xi }}\quad ,\quad \gamma _{\xi }={u_{\xi }l_{\xi } \over 2k}$	(26.a)
$\alpha _{\eta }=\coth \gamma _{\eta }-{1 \over \gamma _{\eta }}\quad ,\quad \gamma _{\eta }={u_{\eta }l_{\eta } \over 2k}$	(26.b)

The lengths ${\textstyle l_{\xi }}$ and ${\textstyle l_{\eta }}$ are taken as the maximum projection of the velocities ${\textstyle u_{\xi }}$ and ${\textstyle u_{\eta }}$ along the element sides (for triangles) and the element diagonals (for quadrilaterals), i.e.

(27.a)

with

(27.b)

In Eq.(27a) ${\textstyle {\textbf {u}}_{\xi }}$ and ${\textstyle {\textbf {u}}_{\eta }}$ contain the global components of the velocity vectors ${\textstyle {\vec {u}}_{\xi }}$ and ${\textstyle {\vec {u}}_{\eta }}$ , respectively. For triangles ${\textstyle {\textbf {d}}_{j}}$ are the element sides vectors, whereas for quadrilaterals ${\textstyle {\textbf {d}}_{j}}$ are the element diagonals vectors.

The next step is to transform Eq.(22) to global axes ${\textstyle x,y}$ . The resulting equation is written as

(28)

where the global diffusion matrix ${\textstyle {\textbf {D}}_{G}}$ is

(29a)

The isotropic diffusion matrix ${\textstyle {\textbf {D}}}$ is given by Eq.(23) and the global balancing diffusion matrix ${\textstyle {\bar {\textbf {D}}}}$ is

(29b)

where the transformation matrix T is given in Eq.(24).

Remark 2. We can write the local balancing diffusion matrix ${\textstyle {\bar {\textbf {D}}}'}$ of Eq.(23) for ${\textstyle k_{\xi }{-}k_{\eta }>0}$ as

(30a)

where

(30b)

If ${\textstyle k_{\eta }}$ is greater than ${\textstyle k_{\xi }}$ a similar split is performed with ${\textstyle (k_{\eta }-k_{\xi })}$ in the position 22 of matrix ${\textstyle {\bar {\textbf {D}}}_{\xi }}$ and ${\textstyle k_{\eta }}$ is replaced by ${\textstyle k_{\xi }}$ in ${\textstyle {\bar {\textbf {D}}}_{iso}}$ .

Eq.(30a) shows that the diffusivity matrix ${\textstyle {\bar {\textbf {D}}}'}$ is equivalent to the sum of a balancing diffusion along the principal curvature direction ${\textstyle {\vec {\xi }}}$ (or ${\textstyle {\vec {\eta }}}$ if ${\textstyle k_{\eta }>k_{\xi }}$ ) and an isotropic diffusion matrix ${\textstyle {\bar {\textbf {D}}}_{iso}}$

The global balancing diffusion matrix of Eq.(29a) can therefore be written as

(31)

Clearly, if the principal curvature direction ${\textstyle {\vec {\xi }}}$ is parallel to the velocity direction, then ${\textstyle u_{\eta }=0}$ , ${\textstyle k_{\eta }=0}$ , ${\textstyle {\bar {\textbf {D}}}_{iso}={0}}$ and

(32)

where ${\textstyle u_{\xi }=\vert {u}\vert }$ and ${\textstyle h_{\xi }}$ is computed by Eq.(26a). Note that the method coincides with the standard SUPG approach in this case.

Remark 3. The global balance diffusion matrix ${\textstyle {\bar {\textbf {D}}}}$ can be also computed from the expression of vector ${\textstyle {\textbf {h}}}$ in global axes as

(33)

where ${\textstyle {\textbf {h}}'=[h_{\xi },h_{\eta }]^{T}}$ and ${\textstyle {\textbf {T}}}$ is the transformation matrix of Eq.(24). The proof of this equivalence is given in the Appendix. In our computation however we have chosen to compute the balancing diffusion matrix ${\textstyle {\bar {\textbf {D}}}}$ via Eqs.(23)–(29).

4.1 Computation of the principal curvature axes for linear elements

The principal curvature axes ${\textstyle \xi ,\eta }$ can be estimated as the eigenvectors of the curvature matrix

(34)

In general matrix C varies within each element.

Computation of the second derivatives field for linear elements is not straightforward and it requires the nodal recovery of the first derivative field. This can be performed by using a simple nodal averaging procedure or more sophisticated superconvergence patch recovery techniques for the first derivative field [50,51].

Excellent results have been obtained by the authors in this work by approximating the principal curvature direction ${\vec {\xi }}$ by the direction of the gradient vector ${\textstyle {\boldsymbol {\nabla }}\phi }$ .

This simplification allows us to estimate the direction ${\textstyle {\vec {\xi }}}$ in a very economical manner as the gradient vector ${\textstyle {\boldsymbol {\nabla }}\phi }$ can be directly computed at any point of a linear element. Direction ${\textstyle {\vec {\eta }}}$ is taken orthogonal to that of ${\textstyle {\vec {\xi }}}$ in an anti-clockwise sense.

For linear triangles ${\textstyle {\boldsymbol {\nabla }}\phi }$ is constant within the element. For four node quadrilaterals ${\textstyle {\boldsymbol {\nabla }}\phi }$ varies linearly. We have assumed in this case that the direction of ${\textstyle {\vec {\xi }}}$ is constant within the element and equal to that of vector ${\textstyle {\boldsymbol {\nabla }}\phi }$ computed at the element center.

The dependence of the balancing diffusion matrix ${\textstyle {\bar {\textbf {D}}}}$ with the principal curvature directions ${\textstyle {\vec {\xi }}}$ and ${\textstyle {\vec {\eta }}}$ introduces a non linearity in the solution process. A simple and effective iterative algorithm is described next.

4.2 General iterative scheme

A stabilized numerical solution can be found by the following algorithm.

Step 0 (SUPG step). At each integration point choose ${\textstyle {}^{0}{\boldsymbol {\xi }}={\textbf {u}}}$ , i.e. the gradient direction is taken coincident with the velocity direction. Compute ${\textstyle {}^{0}{\boldsymbol {\eta }},{}^{1}{\bar {\textbf {D}}}'}$ , ${\textstyle {}^{0}{\bar {\textbf {D}}}}$ and ${\textstyle {}^{0}{\textbf {D}}_{G}}$ from Eqs.(23–29). The expression of the balancing diffusion matrix coincides now precisely with the standard (linear) SUPG form.

Solve for ${\textstyle {}^{0}{\bar {\boldsymbol {\phi }}}}$ .

Verify that the solution is stable. This can be performed by verifying that there are not undershoots or overshoots in the numerical results with respect to the expected physical values. If the SUPG solution is unstable, then implement the following iterative scheme.

For each iteration:

Step 1 Compute at the element center. ${\textstyle {}^{1}{\boldsymbol {\xi }}={\boldsymbol {\nabla }}^{0}{\bar {\phi }}}$ . Then compute ${\textstyle {}^{1}{\boldsymbol {\eta }},{}^{1}{\bar {\textbf {D}}}'}$ , ${\textstyle {}^{1}{\bar {\textbf {D}}}}$ and ${\textstyle {}^{1}{\textbf {D}}_{G}}$ .

Solve for ${\textstyle {}^{1}{\bar {\boldsymbol {\phi }}}}$ .

Step 2 Estimate the convergence of the process. We have chosen the following convergence norm

(35)

where ${\textstyle N}$ is the total number of nodes in the mesh and ${\textstyle \phi _{max}}$ is the maximum prescribed value at the Dirichlet boundary (if ${\textstyle {\bar {\phi }}_{max}=0}$ then ${\textstyle {\bar {\phi }}_{max}=1}$ ). In above steps the left upper indices denote the iteration number.

In the examples shown in the next section ${\textstyle \varepsilon =10^{-3}}$ has been taken in Eq.(35).

If condition (35) is not satisfied, start a new iteration and repeat steps 1 and 2 until convergence. Indexes 0 and 1 are replaced now by ${\textstyle i-1}$ and ${\textstyle i}$ , respectively.

Remark 4. The convergence of the iterative scheme of the previous section obviously depends on the “quality” of the standard SUPG solution obtained in Step 0. Clearly when boundary layers or sharp transverse internal layers dominate the solution, the non linearity of the processes increases. The authors have found that the convergence of the scheme can be substantially improved in some problems if the following expression for the stabilizing dissipation matrix is used in Step 1

(36)

where ${\textstyle \beta }$ is a relaxation factor such that ${\textstyle 0\leq \beta \leq 1}$ .

Convergence of the iterative scheme was achieved in just two iteration in Examples 6.1–6.3 for values of ${\textstyle \beta }$ between 0.8 and 1. The best solution in Example 6.4 was found in five iterations for ${\textstyle \beta =0.3}$ . This is due to the higher non linearity of the process induced by the transverse boundary layers in this case.

Remark 5. The direct iterative scheme presented in Section 4.2 is obviously not the only option to solve the non linear equations resulting from the dependence of the stabilizing diffusion matrix with the principal curvatures of the solution.

A simple and economical alternative is to use a time relaxation technique based in the solution of a pseudo transient problem with a forward Euler scheme and a diagonal mass matrix. This technique was used by Codina in [9,33] for solving a similar class of convection-diffusion problems.

5 EQUIVALENCE WITH SUPG AND SHOCK-CAPTURING TECHNIQUES

As mentioned earlier, the SUPG formulation based on the use of a stabilization term of the form shown in Eq.(10) has severe limitations to provide stabilized solutions (i.e. without overshoots and undershoots) in the vicinity of sharp gradients of the solution transverse to the flow direction. The usual remedy in these cases is to introduce an additional stabilization term (the shock capturing term).

In esence all shock capturing formulations are equivalent to an split of the stabilization term in two parts. One part is taken equal to the standard linear SUPG term. This term suffices to stabilize the effect of the gradient of the solution in the direction of the velocity. The second part is introduced in order to account for gradients of the solution along a direction not aligned with the velocity. This second term (the shock capturing or transverse diffusion term) usually involves the velocity projection along the gradient direction. This term introduces therefore a non linearity in the solution process. The usual approach is to start by computing the linear SUPG solution, check if there are sharp transverse gradients which deserve further stabilization and, if so, introduce the shock capturing term and iterate until a stabilized solution is found in the whole domain.

Despite their effectiveness, most shock capturing methods lack of a general theoretical framework and can be viewed as “ad hoc” extensions of the SUPG method. A serious drawback of these methods is that, in some cases, the amount of transverse diffusion can not be properly controlled and hence the final solution is either underdiffusive or overdiffusive.

The general stabilized formulation we propose here includes the best features of the SUPG and shock capturing techniques. The iterative scheme in fact follows the lines of the standard shock capturing methods [32-36]. The innovation in the proposed formulation lays in the definition of the shock capturing terms, which are deduced form the FIC governing equations written in the principal curvature directions. Numerical results obtained in all examples show that the method preserves the quality of the numerical solution in all situations by introducing the necessary streamline and transverse dissipation.

A similar shock capturing method of this kind was proposed by Oñate [38] in the context of the FIC approach. There the characteristic length vector was split as the sum of a vector along the velocity direction (the linear SUPG term) and a vector in the direction of the solution gradient (the non linear schock capturing term). Good results were obtained with this method for a variety of advection-diffusion problems [38]. The formulation we present here is more general as it encompasses the SUPG method and the shock capturing strategy in a unified approach.

We have shown in Remark 2 that the proposed approach coincides with the standard SUPG method (with no discontinuity capturing) when the principal direction ${\textstyle {\vec {\xi }}}$ is parallel to the velocity direction. This analogy is discussed further next.

Let us split the characteristic vector ${\textstyle h}$ as the sum of two vectors along the velocity and the gradient directions respectively as

(37)

where ${\textstyle h_{s}}$ and ${\textstyle h_{t}}$ are the characteristic lengths along the velocity and gradient directions, respectively. When ${\textstyle h_{s}}$ is linear (i.e., independent of the solution) the first term in the r.h.s of Eq.(37) is the classical SUPG term. The second term depending on the solution gradient is the shock-capturing term.

Eq.(37) allows to solve for ${\textstyle h_{s}}$ and ${\textstyle h_{t}}$ in terms of the components of ${\textstyle {u},{\boldsymbol {\nabla }}{\hat {\phi }}}$ and ${\textstyle {h}}$ as

(38.a)

(38.b)

where ${\textstyle h_{x}}$ and ${\textstyle h_{y}}$ are computed from ${\textstyle h_{\xi }}$ and ${\textstyle h_{\eta }}$ by Eq.(33).

Eqs.(38) clearly show that in the general case the expressions of ${\textstyle h_{s}}$ and ${\textstyle h_{t}}$ are a function of the numerical solution ${\textstyle {\hat {\phi }}}$ . Hence, for arbitrary orientations of the principal curvature direction ${\textstyle {\vec {\xi }}}$ , the stabilized method proposed can be interpreted as the combination of a non linear SUPG term and a non linear shock capturing term. The introduction of this gradient-dependent non linearity in the SUPG term is one of the distinct features of the new formulation.

When the principal curvature direction ${\textstyle {\vec {\xi }}}$ is parallel to the velocity direction then ${\textstyle h_{\eta }=0}$ and

(39a)

Hence

(39b)

and the method coincides in this case with the standard linear SUPG approach, as expected.

A case of practical interest is when the directions of ${\textstyle {\textbf {u}}}$ and ${\textstyle {\boldsymbol {\nabla }}\phi }$ are orthogonal. This occurs for some specific distributions of the velocity and the solution fields or in the vicinity of sharp layers when the mesh is fine enough, or it has some specific orientation. In these cases we have

(40)

Hence from Eq.(33)

(41.a)

and from Eqs.(38)

(41.b)

Once again, the linear SUPG method is recovered and the stabilized formulation does not introduce any additional transverse dissipation, as should be expected.

Remark 6. From above arguments we conclude that the simplified expression ${\textstyle {\textbf {h}}=h_{s}{{\textbf {u}} \over \vert {\textbf {u}}\vert }}$ with ${\textstyle h_{s}}$ being a non linear function of the solution leads to a non linear SUPG term. This has a similarity with the non linear SUPG formulation reported in [28,31]. We have found that in general this assumption does not suffice to stabilize the solution in presence of high gradients transverse to the solution. In these cases the introduction of the term involving ${\textstyle h_{t}}$ in Eq.(37) is essential in order to obtain an improved solution.

6 EXAMPLES

The efficiency of the algorithm is verified in a number of examples of application including sharp gradients at the interior of the domain as well as in boundary layers. The examples show the ability of the algorithm to provide a stabilized and accurate solution in all cases in very few iterations. All units in the examples are in the International System.

In the examples solved a solution is consider unstable if ${\textstyle (1+{\bar {\varepsilon }})\phi _{max}^{p}\leq {\bar {\phi }}\leq (1+{\bar {\varepsilon }})\phi _{min}^{p}}$ where ${\textstyle \phi _{max}^{p}}$ and ${\textstyle \phi _{min}^{p}}$ are the maximum and minimum prescribed values of ${\textstyle {\bar {\phi }}}$ at the Dirichlet boundaries. The parameter ${\textstyle {\bar {\epsilon }}}$ has been taken equal to ${\textstyle 10^{-2}}$ . This condition basically ensures that the solution has no overshoots or undershoots with respect to the prescribed boundary values.

In all examples presented the classical SUPG results (with no discontinuity capturing) are compared with those obtained with the FIC iterative scheme. Examples 6.1–6.3 were solved for a value of ${\textstyle \beta =1}$ in Eq.(36). A converged and accurate solution was obtained in all cases in just two iterations. The convergence process was quite insensitive to the value of ${\textstyle \beta }$ for values of ${\textstyle \beta }$ ranging between 0.8 and 1.0. We note that good stabilized results were also obtained in these examples using just one iteration. However, a slight increase in the accuracy of the results was obtained using two iterations, while no significant increase in accuracy was reached by iterating any further.

The FIC solution reported for Example 6.4 was obtained in five iterations for a value of ${\textstyle \beta =0.3}$ in Eq.(36).

We finally note that most of the examples presented have been successfully solved in the past using different shock capturing techniques [32-36]. In this paper the FIC results are compared with values obtained using the standard (linear) SUPG formulation. The basic aim of the examples presented here is to show that the FIC formulation integrates naturally the best features of both the SUPG and the shock capturing formulations.

6.1 Square domain with uniform Dirichlet conditions, upward diagonal velocity and zero source

Figure 2 shows the geometry of the analysis domain and the boundary conditions. The domain is a square of side equal to 10. Zero values of the variable ${\textstyle \phi }$ are prescribed at the inlet boundaries ${\textstyle x=0}$ and ${\textstyle y=0}$ , whereas ${\textstyle \phi =10}$ is prescribed at the outlet boundaries ${\textstyle x=10}$ and ${\textstyle y=10}$ . The isotropic diffusion ${\textstyle k=0.01}$ and the velocity vectors are taken parallel to the diagonal of the domain with ${\textstyle {u}={3{\sqrt {2}} \over 2}[1,1]^{T}}$ . The problem has been solved with the structured mesh of ${\textstyle 10\times 10}$ four node rectangles shown in the figure. The Peclet number along the velocity direction is ${\textstyle \gamma =300}$ .

Figure 2: Square domain with uniform Dirichlet conditions, upward diagonal velocity and zero source. SUPG and FIC solutions obtained with a structured mesh of

linear four node square elements

Results in Figure 2 show contours of the numerical solution for the SUPG and the FIC iterative scheme proposed. FIC ${\textstyle ^{2}}$ in the figures denotes results obtained after two iterations of the scheme described in Section 4.2. A value of ${\textstyle \beta =1.0}$ in Eq.(36) has been chosen. Note that the SUPG method yields incorrect negative values of ${\textstyle \phi }$ in the vicinity of the upper right corner, as expected, whereas the FIC ${\textstyle ^{2}}$ results are physically correct over the domain. A plot of the distribution of the solution along the horizontal center line and the diagonal line is shown. The difference between the FIC and SUPG results along the diagonal line is clearly seen.

Figure 3 shows results for the same problem solved with a structurerd mesh of ${\textstyle 2\times 10\times 10}$ three node triangular elements. The effect of the stabilization process in more visible in this case as the SUPG solution has more undershoots than for four node square elements.

Figure 3: Square domain with uniform Dirichlet conditions, upward diagonal velocity and zero source. SUPG and FIC solutions obtained with a structured mesh of

linear three node triangles

The convergence of the iterative process using square and triangular elements is plotted in Figure 4 where the evolution of the values of ${\textstyle {\bar {\phi }}}$ at three different points within the analysis domain is plotted. The evolution of the norm of Eq.(35) with the number of iterations is also shown. Note that a converged solution ( ${\textstyle \varepsilon <10^{-3}}$ ) is obtained in just two iterations.

Figure 4: Square domain with uniform Dirichlet conditions, upward diagonal velocity and zero source. Convergence of the FIC solution for square and triangular elements

Figures 5 and 6 show the solution of a similar problem using the same number of four node rectangular elements and linear triangles with a geometrical aspect ratio of 2:1. The SUPG solution of both cases is clearly non physical and the FIC iterative process provides again a stabilized and accurate solution in just two iterations.

Figure 5: Square domain with uniform Dirichlet conditions, upward diagonal velocity and zero source. SUPG and FIC solutions obtained with a structured mesh of

linear four node rectangular elements. Geometrical aspect ratio 2:1

Figure 6: Square domain with uniform Dirichlet conditions, upward diagonal velocity and zero source. SUPG and FIC solutions obtained with a structured mesh of

linear triangles. Geometrical aspect ratio 2:1

6.2 Square domain with non uniform Dirichlet conditions, downwards diagonal velocity and zero source

The geometry of the problem and the boundary conditions are shown in Figure 7. The side of the square domain has a length equal to one. The velocity vector is ${\textstyle {u}=[5,-9]^{T}\times 10^{6}}$ . An isotropic diffusion of unit value is assumed. The correct solution is a uniform plateau of value ${\textstyle \phi =100}$ with an internal sharp layer in the vicinity of the left lower region where the solution drops to ${\textstyle {\bar {\phi }}=0}$ , a boundary layer at the right hand side and at partial boundary layer at the bottom side where the value of ${\textstyle \phi }$ is prescribed to zero.

Figure 7: Square domain with non uniform Dirichlet conditions, downwards diagonal velocity and zero source. SUPG and FIC solutions obtained with a structured mesh of

linear four node square elements

Figures 7 and 8 show results obtained with structured meshes of ${\textstyle 20\times 20}$ four node square elements and ${\textstyle 2\times 20\times 20}$ linear triangles, respectively. The SUPG solution in both cases has overshoots in the vicinity of the sharp gradient zones, whereas the FIC iterative scheme yields a stabilized and accurate solution in two iterations. Again ${\textstyle \beta =1}$ in Eq.(36) was chosen. Note that the boundary layers and the internal sharp gradient zone are well captured without overshoots or undershoots with the FIC method using both types of elements.

Figure 8: Square domain with non uniform Dirichlet conditions, downwards diagonal velocity and zero source. SUPG and FIC solutions obtained with a structured mesh of

three node triangles

The convergence of the process for both elements is shown in Figure 9.

Figure 9: Square domain with uniform Dirichlet conditions, upward diagonal velocity and zero source. Convergence of the FIC solution for square and triangular elements

Figures 10 and 11 show results for the same problem using coarse and refined unstructured meshes of linear quadrilaterals and triangles, respectively. The performance of the FIC algorithm is very good in both cases despite the difficulty to capture the sharp layers with coarse unstructured meshes.

Figure 10: Square domain with non uniform Dirichlet conditions, downwards diagonal velocity and zero source. SUPG and FIC solutions obtained with an unstructured mesh of 432 linear four node quadrilateral elements

Figure 11: Square domain with non uniform Dirichlet conditions, downwards diagonal velocity and zero source. SUPG and FIC solutions obtained with an unstructured mesh of 780 three node triangles

The solution is improved in both cases by using a refined unstructured mesh as shown in Figures 12 and 13. Note that the overall solutions obtained with the linear triangle and unstructured meshes are slightly better than those obtained with linear quadrilaterals for a similar number of elements. The improvements are clearer in the vicinity of the boundary layer at the right hand side of the domain.

Figure 12: Square domain with non uniform Dirichlet conditions, downwards diagonal velocity and zero source. SUPG and FIC solutions obtained with a refined unstructured mesh of 1826 linear four node quadrilaterals

Figure 13: Square domain with non uniform Dirichlet conditions, downwards diagonal velocity and zero source. SUPG and FIC solutions obtained with a refined unstructured mesh of 1898 linear triangles

Figures 14 and 15 show results for a similar problem solved now with structured meshes of linear rectangles and triangles with a geometrical aspect ratio of 2:1. The correct solution in this case is a central plateau of ${\textstyle \phi =10}$ with a boundary layer at the right hand side where the solution drops to ${\textstyle \phi =0}$ , and two sharp gradient zones at both sides of the plateau within the domain. The figures show the contours of ${\textstyle {\bar {\phi }}}$ over the analysis domain and the distribution of ${\textstyle {\bar {\phi }}}$ along the vertical and horizontal central lines. As expected the SUPG solution yields over and undershoots. Two iterations of the FIC iterative scheme again suffice to provide a correct solution over the whole domain.

Figure 14: Square domain with non uniform Dirichlet conditions, upwards diagonal velocity and zero source. SUPG and FIC solutions obtained with a structured mesh of

linear four node rectangular elements. Geometrical aspect ratio 2:1

Figure 15: Square domain with non uniform Dirichlet conditions, upwards diagonal velocity and zero source. SUPG and FIC solutions obtained with a structured mesh of

three node triangles. Geometrical aspect ratio 2:1

This problem has been chosen by many authors as a test for stabilized formulations. A comparison of the solution obtained with different shock capturing techniques was reported by Codina [9,33] using a time relaxation scheme. The FIC results presented here reproduce the best results obtained for this problem with the shock capturing methods analyzed in [9,33].

6.3 Rectangular domain with Neumann and non uniform Dirichlet conditions, rotational velocity field and zero source

The rectangular domain of sides ${\textstyle 1\times 2}$ units and the boundary conditions are shown in Figure 16. The rotational velocity field is defined by ${\textstyle {u}=10^{4}\times [y(1-x^{2}),-x(1-y^{2})]^{T}}$ . A unit isotropic diffusion is assumed. The velocity is prescribed at the domain sides shown. In the rest of the sides the Neumann boundary condition with ${\textstyle q_{n}^{p}=0}$ is assumed. The correct solution is a uniform plateau of ${\textstyle \phi =100}$ with a boundary layer at the right hand side where the solution drops towards the prescribed value of ${\textstyle \phi =0}$ and a circular sharp gradient region around the lower circular zone where the solution takes a zero value.

Figures 16 and 17 present the results obtained with the structured meshes of linear square and triangular elements respectively shown. The contours of ${\textstyle {\bar {\phi }}}$ over the domain and the distribution of the solution along an horizontal line are plotted. Results clearly show the incapacity of the SUPG method to correctly capture the boundary layer and the internal gradient region. The solution obtained with two iterations of the FIC scheme presented is physically sound and accurate over the whole analysis domain as shown in the figures. FIC results were once more obtained for a value of ${\textstyle \beta =1}$ in Eq.(36).

Figure 16: Rectangular domain with Neumann and non uniform Dirichlet conditions, rotational velocity field and zero source. SUPG and FIC solutions obtained with a structured mesh of

four node square elements

Figure 17: Rectangular domain with Neumann and non uniform Dirichlet conditions, rotational velocity field and zero source. SUPG and FIC solutions obtained with a structured mesh of

three node triangles

Figure 18 shows results for the solution of the same problem using coarse and refined unstructured meshes of linear quadrilaterals. Note the accuracy of the FIC stabilized solution, even for the coarse mesh.


Figure 18: Rectangular domain with Neumann and non uniform Dirichlet conditions, rotational velocity field and zero source. SUPG and FIC solutions obtained with two unstructured meshes of four node quadrilaterals. (a) Coarse mesh of 892 elements. (b) Refined mesh of 2054 elements

The same problem was finally solved with a unstructured mesh of 1554 linear triangles. The number of elements chosen is quite similar to that used for the structured solution (1600 elements). The solution obtained shown in Figure 19 is quite accurate over the whole domain, despite the relative coarseness of the mesh. This is another prove of the best performance of the linear triangle in unstructured meshes.

Figure 19: Rectangular domain with Neumann and non uniform Dirichlet conditions, rotational velocity field and zero source. SUPG and FIC solutions obtained with an unstructured mesh of 1554 three node triangles

6.4 Square domain with homogeneous Dirichlet condition and constant source

The computational domain for this last example is again the unit square, discretized now with uniform meshes of 20 ${\textstyle \times }$ 20 linear quadrilateral and ${\textstyle 2\times 20\times 20}$ three node triangles. The value of ${\textstyle \phi }$ is prescribed to zero on the whole boundary. The velocity vector is ${\textstyle {u}=[1,0]^{T}}$ and the isotropic diffusion is ${\textstyle k=10^{-8}}$ . A uniformly distributed source ${\textstyle Q=1}$ over the whole domain has been taken.

The correct solution of this problem is a linear distribution of ${\textstyle \phi }$ increasing uniformly from a value of zero at the inlet boundary until it reaches a maximum value in the vecinity of the outlet boundary. Three boundary layers exist for this problem, one at the exit boundary and two at the two lateral sides, where the value of ${\textstyle \phi }$ drops drastically to zero. Figures 20 and 21 show results obtained with the two structured meshes of quadrilateral and triangles, respectively. The linear SUPG solution yields oscillations in the direction normal to the lateral layers, although a good resolution of the layer normal to the velocity field is obtained, as expected. The FIC solution for ${\textstyle \beta =1}$ in Eq.(36) yielded a smoooth solution near the boundary layers, free of the lateral oscillations which are present in the SUPG solution. The best FIC solution in this case was however obtained for a value of ${\textstyle \beta =0.3}$ in Eq.(36). Five iterations of the scheme were needed to obtain convergence in this case. The numerical results obtained for both meshes are shown in Figures 20 and 21.

Figure 20: Square domain with homogeneous Dirichlet conditions and constant source. Solution for a structural mesh of

four node quadrilaterals

Figure 21: Square domain with homogeneous Dirichlet conditions and constant source. Solution for an unstructured mesh of

three node triangles

This problem was also studied by Codina [9,33] with four node square elements as another benchmark to study the performance of different shock capturing methods. The FIC method proposed here has yielded accurate solutions for both quadrilateral and triangular elements, which are difficult to obtain for many of the existing stabilization procedures, as reported in [9,33].

7 CONCLUDING REMARKS

The FIC governing equations for the convection-diffusion problem written in the principal curvature directions introduce naturally the necessary balancing diffusion terms to accurately solve convection dominated problems with sharp internal gradients and boundary layers. In our computations over four node quadrilaterals and three node triangles very good results have been obtained by assuming that the principal curvature direction ${\textstyle \xi }$ is constant over each element and equal to the direction of the gradient vector at the element center.

The dependence of the balancing diffusion with the solution gradient makes the stabilization process non linear. An efficient iterative scheme has been proposed where the standard SUPG formulation is used to obtain an initial solution. Convergence of the iterative scheme has proved to be very fast in all problems solved.

The balancing diffusion terms can be interpreted as the sum of a non linear SUPG diffusion term and a non linear shock capturing term. If the velocity is parallel to the gradient vector, then the standard (linear) SUPG term is recovered.

The efficiency and accuracy of the method was tested in a number of problems with arbitrary sharp gradients solved with structured and unstructured meshes of linear quadrilaterals and triangles.

Very accurate results were obtained in just two iterations for most of the problems solved (Examples 6.1–6.3). One iteration sufficed to provide good quality stabilized results in many cases.

As a general rule, both linear triangles and bilinear quadrilaterals yielded similar good results for structured meshes. The linear triangle shows a slightly superior performance for unstructured meshes.

The examples presented show that the FIC formulation incorporates in a unified theoretical and computational framework the best features of the classical streamline stabilization methods and shock capturing techniques.

ACKNOWLEDGEMENTS

The authors thank S. Badia, C. Felippa, R. Löhner, R.L. Taylor and O.C. Zienkiewicz for many useful discussions.

APPENDIX A

Computation of the balancing diffusion matrix ${\overline {\textbf {D}}}$ in local and global axesComputation of the balancing diffusion matrix D ¯ {\displaystyle {\overline {\textbf {D}}}} in local and global axes

The balancing diffusion matrix ${\textstyle {\bar {\textbf {D}}}}$ can be computed by transforming the local matrix ${\textstyle {\bar {\textbf {D}}}'}$ of Eq.(23) to global axes, or by using directly Eq.(33). We will prove next that both expressions yield the same system of governing differential equations, i.e. we will prove that

(A.1)

where

(A.2)

and

${\textbf {h}}=\left\{{\begin{array}{c}h_{x}\\h_{y}\end{array}}\right\}={\textbf {T}}^{T}\left\{{\begin{array}{c}h_{\xi }\\h_{\eta }\end{array}}\right\}={\textbf {T}}^{T}{\textbf {h}}'$	(A.3)
${\textbf {u}}=\left\{{\begin{array}{c}u\\v\end{array}}\right\}={\textbf {T}}^{T}\left\{{\begin{array}{c}u_{\xi }\\u_{\eta }\end{array}}\right\}={\textbf {T}}^{T}{\textbf {u}}'$	(A.4)

In above ${\textstyle c_{\alpha }=\cos \alpha ,s_{\alpha }=\sin \alpha }$ where ${\textstyle \alpha }$ is the angle that the principal curvature axis ${\textstyle {\vec {\xi }}}$ forms with the global direction ${\textstyle x}$ . The characteristic length distances ${\textstyle h_{\xi }}$ and ${\textstyle h_{\eta }}$ are computed as described in Eqs.(25–27).

In order to prove the identity in Eq.(A.1) we proceed as follows

{\bar {\textbf {D}}}={\textbf {T}}^{T}{\bar {\textbf {D}}}'{\textbf {T}}={1 \over 2}\left[{\begin{array}{ccc}\displaystyle {u_{\xi }h_{\xi }c_{\alpha }^{2}+u_{\eta }h_{\eta }s_{\alpha }^{2} \over {\hbox{Sym.}}}&\!\!\!\!{\Bigg |}&\!\!\!\!\displaystyle {(u_{\xi }h_{\xi }-u_{\eta }h_{\eta })s_{\alpha }c_{\alpha } \over u_{\xi }h_{\xi }s_{\alpha }^{2}+u_{\eta }h_{\eta }c_{\alpha }^{2}}\end{array}}\right]

(A.5)

{\begin{array}{ll}{\boldsymbol {\nabla }}^{T}({\textbf {T}}^{T}{\bar {\textbf {D}}}'{\textbf {T}}){\boldsymbol {\nabla }}\phi =&{1 \over 2}\left[(u_{\xi }h_{\xi }c_{\alpha }^{2}+u_{\eta }h_{\eta }s_{\alpha }^{2})\phi _{x}^{''}+(u_{\xi }h_{\xi }s_{\alpha }^{2}+u_{\eta }h_{\eta }c_{\alpha }^{2})\phi _{y}^{''}+\right.\\&\left.+2(u_{\xi }h_{\xi }-u_{\eta }h_{\eta })s_{\alpha }c_{\alpha }\phi _{xy}^{''}\right]\end{array}}

(A.6)

On the other hand

(A.7)

Finally

{\begin{array}{ll}{\boldsymbol {\nabla }}^{T}\left({1 \over 2}{\textbf {h}}{\textbf {u}}^{T}\right){\boldsymbol {\nabla }}\phi =&{1 \over 2}{\boldsymbol {\nabla }}^{T}[{\textbf {T}}^{T}({\textbf {h}}'{\textbf {u}}^{'T}){\textbf {T}}]{\boldsymbol {\nabla }}\phi ]={1 \over 2}\left[(u_{\xi }h_{\xi }c_{\alpha }^{2}+\right.\\&\left.u_{\eta }h_{\eta }s_{\alpha }^{2})\phi _{x}^{''}+(u_{\xi }h_{\xi }s_{\alpha }^{2}+u_{\eta }h_{\eta }c_{\alpha }^{2})\phi _{y}^{''}+2(u_{\xi }h_{\xi }-u_{\eta }h_{\eta })s_{\alpha }c_{\alpha }\phi _{xy}^{''}\right]+\\&\left[s_{\alpha }c_{\alpha }(\phi _{y}^{''}-\phi _{x}^{''})+(c_{\alpha }^{2}-s_{\alpha }^{2})\phi _{xy}^{''}\right](u_{\xi }h_{\xi }+u_{\eta }h_{\eta })\end{array}}

(A.8)

Clearly Eqs.(A.6) and (A.8) coincide if the last term in Eq.(A.8) is zero. This happens when ${\textstyle {\vec {\xi }}}$ and ${\textstyle {\vec {\eta }}}$ are the principal curvature directions. In order to prove this we will consider the usual case for which ${\textstyle \phi _{xy}^{''}=\phi _{yx}^{''}}$ and both principal curvature directions are orthogonal. Then it exists a transformation matrix ${\textstyle {\textbf {T}}}$ such that

{\textbf {T}}\left[{\begin{array}{cc}\phi _{x}^{''}&\phi _{xy}^{''}\\\phi _{xy}^{''}&\phi _{y}^{''}\end{array}}\right]{\textbf {T}}^{T}=\left[{\begin{array}{cc}\phi _{\xi }^{''}&0\\0&\phi _{\eta }^{''}\end{array}}\right]

(A.9)

with ${\textstyle {\textbf {T}}}$ defined as in Eq.(A.2).

The condition which ensures a diagonal matrix in the r.h.s. of Eq.(A.9) is

(A.10)

which is precisely the condition which cancels out the last term in Eq.(A.8).

Hence the identity of Eq.(A.1) is satisfied.

REFERENCES

[1] Zienkiewicz O.C. and Taylor R.L. The finite element method. 5th Edition, 3 Volumes, Butterworth–Heinemann, 2000.

[2] Hirsch C. Numerical computation of internal and external flow, J. Wiley, Vol. 1 1988, Vol. 2, 1990.

[3] Heinrich J.C., Huyakorn P.S. and Zienkiewicz O.C. An upwind finite element scheme for two dimensional convective transport equations, Int. J. Num. Meth. Engng., 11, 131–143, 1977.

[4] Kelly D.W., Nakazawa S., Zienkiewicz O.C. and Heinrich J.C. A note on upwind and anisotropic balancing dissipation in finite element approximation to convective diffusion problems. Int. J. Num. Meth. Engng., 15, pp. 1705–1711, 1980.

[5] Brooks A.N. and Hughes T.J.R. Streamline upwind/Petrov-Galerkin formulation for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations. Comput. Methods Appl. Mech. Engrg., 32, 199–259, 1982.

[6] Hughes T.J.R and Tezduyar, T.E. Finite element method for first order hyperbolic systems with particular emphasis on the compressible Euler equations. Comput. Meth. Appl. Mech. Engng., 45, 217–284, 1984.

[7] Hughes T.J.R and Mallet M. 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.

[8] Codina R., Oñate E. and Cervera M. The intrinsic time for the streamline-upwind Petrov-Galerkin formulation using quadratic elements. Int. J. Num. Meth. Engrg., 94, 239–262, 1992.

[9] Codina R. A finite element formulation for the numerical solution of the convection-diffusion equations. Monograph no. 14, CIMNE, Barcelona 1993.

[10] Donea J. A Taylor-Galerkin method for convective transport problems. Int. J. Num. Meth. Engng., 20, pp. 101–119, 1984.

[11] Donea J. and Huerta A. Finite element methods for flow problems. J. Wiley, 2003.

[12] Donea J., Belytschko T. and Smolinski P. A generalized Galerkin method for steady state convection-diffusion problems with applications to quadratic shape functions elements. Comp. Meth. Appl. Mech. Engng., 48, 25–43, 1985.

[13] Goldschmit M.B. and Dvorkin E.N. On the solution of the steady-state convection-diffusion equation using quadratic elements. A generalized Galerkin technique also reliable with distorted meshes. Engng. Comput., 11, 6, 565–579, 1994.

[14] Hughes T.J.R, Franca L.P. and Hulbert G.M. 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.

[15] Franca L.P., Frey S.L. and Hughes T.J.R. Stabilized finite element methods: I. Application to the advective-diffusive model. Comput. Meth. Appl. Mech. Engn, Vol. 95, pp. 253–276, 1992.

[16] Franca L.P. and Dutra do Carmo E.G. The Galerkin Gradient Least Square method. Comput. Methods Appl. Mech. Engrg., 74, 41–54, 1989.

[17] Douglas J. and Russell T.F. Numerical methods for convection dominated diffusion problems based on combining the method of characteristics with finite element or finite difference procedures. SIAM J. Numer. Anal., 19, 871, 1982.

[18] Löhner R., Morgan K. and Zienkiewicz O.C. The solution of non-linear hyperbolic equation systems by the finite element method. Int. J. Num. Meth. in Fluids, 4, 1043, 1984.

[19] Zienkiewicz O.C. and Codina R. 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.

[20] Hughes T.J.R. Multiscale phenomena: Green functions, subgrid scale models, bubbles and the origins of stabilized methods. Comput. Methods Appl. Mech. Engrg 1995, 127, 387–401, 1995.

[21] Brezzi F., Franca L.P., Hughes T.J.R. and Russo A. $b=\int g$ . Comput. Methods Appl. Mech. Engrg., 145, 329–339, 1997.

[22] Codina R. Stabilization of incompressibility and convection through orthogonal sub-scales in finite element method. Comput. Methods Appl. Mech. Engrg., 190, 1579–1599, 2000.

[23] Hauke G. A simple subgrid scale stabilized method for the advection-diffusion-reaction equation. Comput. Methods Appl. Mech. Engrg., 191, 2925–2948, 2002.

[24] Brezzi, F. and Russo A. Choosing bubbles for advection-diffusion problems. Math Models Methods Appl. Sci., 4, 571–587, 1994.

[25] Farhat C., Harari I. and Franca L.P. The discontinuous enrichment method. Comput. Methods Appl. Mech. Engrg., 190, 6455–6479, 2001.

[26] Ilinca F., Hétu J.F. and Pelletier D. On stabilized finite element formulation for incompressible advective-diffusive transport and fluid flow problems. Comput. Methods Appl. Mech. Engrg., 188, 235–257, 2000.

[27] Tezduyar T.E. and Osawa Y. Finite element stabilization parameters computed from element matrices and vectors. Comput. Methods Appl. Mech. Engrg., 190, 411–430, 2000.

[28] Tezduyar T.E. Adaptive determination of the finite element stabilization parameters. In Proceedings of the ECCOMAS Computational Fluid Dynamics Conference 2001 (CD-ROM), Swansea, Wales, United Kingdom 2001.

[29] Harari I., Franca L.P. and Oliveira S.P. Streamline design of stability parameters for advection-diffusion problems. J. Comput. Physics, 171, 115–131, 2001.

[30] Nesliturk A. and Harari I. The nearly-optimal Petrov-Galerkin method for convection-diffusion problems. Comput Methods Appl. Mech. Engrg., 192, 2501–2519, 2003.

[31] Tezduyar T.E. Computation of moving boundaries and stabilization parameters. Int. J. Num. Meth. Fluids, 43, 555–575, 2003.

[32] Hughes T.J.R. and Mallet M. A new finite element formulations for computational fluid dynamics: IV. A discontinuity capturing operator for multidimensional advective-diffusive system. Comput. Methods Appl. Mech. Engrg., 58, 329–336, 1986b.

[33] Codina R. A discontinuity-capturing crosswind dissipation for the finite element solution of the convection-diffusion equation. Comput. Methods Appl. Mech. Engrg., 110, 325–342, 1993.

[34] Hughes T.J.R., Mallet M. and Mizukami A. A new finite element formulation for comptutational fluid dynamics: II Beyond SUPG. Comput. Methods Appl. Mech. Engrg., 54, 341–355, 1986.

[35] Galeo A.C. and Dutra do Carmo E.G. A consistent approximate upwind Petrov-Galerkin method for convection-dominated problems. Comput. Methods Appl. Mech. Engrg., 68, 83–95, 1988.

[36] Tezduyar T.E. and Park Y.J. Discontinuity-capturing finite element formulations for nonlinear convection-diffusion-reaction equations. Comput. Methods Appl. Mech. Engrg., 59, 307–325, 1986

[37] Oñate E. Derivation of stabilized equations for advective-diffusive transport and fluid flow problems. Comput. Methods Appl. Mech. Engrg., 151:1-2, 233–267, 1998.

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

[39] Oñate E. and Manzan M. Stabilization techniques for finite element analysis of convection diffusion problems. In Comput. Anal. of Heat Transfer, WIT Press, G. Comini and B. Sunden (Eds.), 2000.

[40] Oñate E., García J. and Idelsohn S.R. Computation of the stabilization parameter for the finite element solution of advective-diffusive problems. Int. J. Num. Meth. Fluids, 25, 1385–1407, 1997.

[41] Oñate E. and Manzan M. A general procedure for deriving stabilized space-time finite element methods for advective-diffusive problems. Int. J. Num. Meth. Fluids, 31, 203–221, 1999.

[42] Felippa C.A. and Oñate E. Nodally exact Ritz discretization of 1D diffusion-absorption and Helmholtz equations by variational FIC and modified equation methods. Research Report No. PI 237, CIMNE, Barcelona 2004. Submitted to Computational Mechanics.

[43] Oñate E., Miquel J. and Hauke G. Stabilized formulation for the 1D advection-diffusion-absorption equation using finite calculus and the finite element method. Submitted to Comput. Meth. Appl. Mech. Engng., Dec. 2004.

[44] Oñate E. A stabilized finite element method for incompressible viscous flows using a finite increment calculus formulation. Comput. Methods Appl. Mech. Engrg. 2000; 182:1–2:355–370.

[45] Oñate E., García J. A finite element method for fluid-structure interaction with surface waves using a finite calculus formulation. in Comput. Methods Appl. Mech. Engrg., 191:6-7, 635-660, 2001.

[46] García J. and Oñate E. An unstructured finite element solver for ship hydrodynamic problems. Journal of Appl. Mechanics, 70, 18–26, 2003.

[47] Oñate E., García J., Idelsohn S.R. Ship Hydrodynamics. Encyclopedia of Computational Mechanics, T. Hughes, R. de Borst and E. Stein (Eds.), J. Wiley, 2004.

[48] Oñate E., Taylor R.L., Zienkiewicz O.C. and Rojek J. A residual correction method based on finite calculus. Engineering Computations, Vol. 20, No. 5/6, 629–658, 2003.

[49] Oñate E., Rojek J., Taylor R.L. and Zienkiewicz O.C. Finite calculus formulation for analysis of incompressible solids using linear triangles and tetrahedra. Int. J. Num. Meth. Engng., 59, (11), 1473–1500, 2004.

[50] Zienkiewicz OC, Zhu JZ. The Superconvergent patch recovery (SPR) and adaptive finite element refinement. Comput. Methods Appl. Mech. Engrg. 1992; 101:207–224.

[51] Wiberg NW, Abdulwahab F, Li XD. (1997). Error estimation and adaptive procedures based on superconvergent patch recovery. Archives Comput. Meth. Engng. 1997; 4:3:203–242.

Latest revision as of 11:15, 12 November 2018

Abstract

1 INTRODUCTION

2 FIC GOVERNING EQUATIONS FOR STEADY-STATE CONVECTION-DIFFUSIVE PROBLEMS

2.1 Finite element discretization

2.2 Equivalence with the SUPG method

3 COMPUTATION OF THE CHARACTERISTIC LENGTH FOR THE 1D CONVECTION-DIFFUSION EQUATION

4 COMPUTATION OF THE CHARACTERISTIC LENGTH VECTOR

4.1 Computation of the principal curvature axes for linear elements

4.2 General iterative scheme

5 EQUIVALENCE WITH SUPG AND SHOCK-CAPTURING TECHNIQUES

6 EXAMPLES

6.1 Square domain with uniform Dirichlet conditions, upward diagonal velocity and zero source

6.2 Square domain with non uniform Dirichlet conditions, downwards diagonal velocity and zero source

6.3 Rectangular domain with Neumann and non uniform Dirichlet conditions, rotational velocity field and zero source

6.4 Square domain with homogeneous Dirichlet condition and constant source

7 CONCLUDING REMARKS

ACKNOWLEDGEMENTS

APPENDIX A

Computation of the balancing diffusion matrix ${\overline {\textbf {D}}}$ in local and global axesComputation of the balancing diffusion matrix D ¯ {\displaystyle {\overline {\textbf {D}}}} in local and global axes

REFERENCES

Document information

Document Score

Share this document

@@ Line 1: / Line 1: @@
-<!-- metadata commented in wiki content
+Published in ''Comput. Meth. Appl. Mech. Engng.'', Vol. 195 (13-16), pp. 1793–1825, 2006<br />
-==FINITE ELEMENT FORMULATION FOR CONVECTIVE-DIFFUSIVE PROBLEMS WITH SHARP GRADIENTS USING FINITE CALCULUS==
+doi: 10.1016/j.cma.2005.05.036
-'''Eugenio Oñate, Francisco Zárate and Sergio R. Idelsohn'''
-{|
-|-
-|International Center for Numerical Methods in Engineering (CIMNE)
-|-
-| Universidad Politécnica de Cataluña
-|-
-| Edificio C1, Gran Capitán s/n, 08034, Barcelona, Spain
-|-
-| e-mail: [mailto:onate@cimne.upc.edu onate@cimne.upc.edu]
-|}
--->
 ==Abstract==
 A finite element method (FEM) for steady-state convective-diffusive problems presenting sharp gradients of the solution both in the interior of the domain and in boundary layers is presented. The necessary stabilization of the numerical solution is provided by the Finite Calculus (FIC) approach. The FIC method is based in the solution by the Galerkin FEM of a modified set of governing equations which include characteristic length parameters. It is shown that the FIC balance equation for the multidimensional convection-diffusion problem written in the principal curvature axes of the solution, introduces an orthotropic diffusion which stabilizes the numerical solution both in smooth regions as well in the vicinity of sharp gradients. The dependence of the stabilization terms with the principal curvature directions of the solution makes the method non linear. Details of the iterative scheme to obtain stabilized results are presented together with examples of application which show the efficiency and accuracy of the approach.
-<br/><br/>
 ==1 INTRODUCTION==
-It is well known that the standard Galerkin FEM solution of the steady-state convective-diffusive equation is unstable for values of the Peclet number greater than one (see Volume 3 in <span id='citeF-1'></span>[[#cite-1|1]]). A number of numerical schemes have been proposed in order to guarantee that the numerical solution is stable, that is, that the solution has a physical meaning. In the first attempts to solve this problem the underdiffusive character of the Galerkin FEM (and the analogous central finite difference scheme) for convective-diffusive problems was corrected by adding “artificial diffusion terms” to the governing equations <span id='citeF-1'></span>[[#cite-1|1]]. The relationship of this approach with the upwind finite difference method <span id='citeF-2'></span>[[#cite-2|2]] lead to the derivation of a variety of Petrov-Galerkin FEM. All these methods can be interpreted as extensions of the standard Galerkin variational form of the FEM by adding residual-based integral terms computed over the element domains. Among the many stabilization methods of this or similar kind we name the Upwind FEM <span id='citeF-3'></span>[[#cite-3|3]], the Streamline Upwind Petrov-Galerkin (SUPG) method <span id='citeF-5'></span>[[#cite-5|5]], the Taylor-Galerkin method <span id='citeF-12'></span>[[#cite-12|12]], the generalized Galerkin method <span id='citeF-17'></span>[[#cite-17|17]], the Galerkin Least Square method and related approaches <span id='citeF-19'></span>[[#cite-19|19]], the Characteristic Galerkin method <span id='citeF-22'></span>[[#cite-22|22]], the Characteristic Based Split method <span id='citeF-24'></span>[[#cite-24|24]], the  Subgrid Scale method <span id='citeF-25'></span>[[#cite-25|25]], the Residual Free Bubbles method <span id='citeF-37'></span>[[#cite-37|37]], the Discontinuous Enrichment Method <span id='citeF-39'></span>[[#cite-39|39]] and the Streamline Upwind with Boundary Terms method <span id='citeF-40'></span>[[#cite-40|40]].  Basically all the methods make use of a single stabilization parameter which suffices to stabilize the numerical solution along the velocity (streamline) direction. The computation of the streamline stabilization parameter for multidimensional problems is usually based on extensions of the optimal value of the parameter for the simpler 1D case. Specific attempts to design the stability parameter for multidimensional problems in the context of the Petrov-Galerkin formulation have been recently reported <span id='citeF-41'></span>[[#cite-41|41]]. In general all the stabilized methods above mentioned yield good results for 1D-type problems where the velocity vector is aligned with the direction of the gradient of the solution. However, the use of a single stabilization parameter is  insufficient to provide stabilized solutions in the vicinity of sharp gradients not aligned with the velocity direction, which may appear at the interior of the domain or at boundary layers. The usual remedy for these situations is to use the so called “shock capturing” or “discontinuity-capturing” schemes <span id='citeF-49'></span>[[#cite-49|49]]. These methods basically add additional transverse diffusion terms in selected elements in order to correct the undershoots and overshoots yielded by the Petrov-Galerkin solution in the sharp gradient zones. Usually the new shock capturing diffusion terms depend on the gradient of the solution and the scheme becomes non linear.
+It is well known that the standard Galerkin FEM solution of the steady-state convective-diffusive equation is unstable for values of the Peclet number greater than one (see Volume 3 in <span id='citeF-1'></span>[[#cite-1|[1]]]). A number of numerical schemes have been proposed in order to guarantee that the numerical solution is stable, that is, that the solution has a physical meaning. In the first attempts to solve this problem the underdiffusive character of the Galerkin FEM (and the analogous central finite difference scheme) for convective-diffusive problems was corrected by adding “artificial diffusion terms” to the governing equations <span id='citeF-1'></span>[[#cite-1|[1,2]]]. The relationship of this approach with the upwind finite difference method <span id='citeF-2'></span>[[#cite-2|[2]]] lead to the derivation of a variety of Petrov-Galerkin FEM. All these methods can be interpreted as extensions of the standard Galerkin variational form of the FEM by adding residual-based integral terms computed over the element domains. Among the many stabilization methods of this or similar kind we name the Upwind FEM <span id='citeF-3'></span>[[#cite-3|[3,4]]], the Streamline Upwind Petrov-Galerkin (SUPG) method <span id='citeF-5'></span>[[#cite-5|[5-9]]], the Taylor-Galerkin method <span id='citeF-12'></span>[[#cite-12|[10,11]], the generalized Galerkin method <span id='citeF-17'></span>[[#cite-17|[12,13]]], the Galerkin Least Square method and related approaches <span id='citeF-19'></span>[[#cite-19|[14-16]]], the Characteristic Galerkin method <span id='citeF-22'></span>[[#cite-22|[17,18]]], the Characteristic Based Split method <span id='citeF-24'></span>[[#cite-24|[19]]], the  Subgrid Scale method <span id='citeF-25'></span>[[#cite-25|[20-23]]], the Residual Free Bubbles method <span id='citeF-37'></span>[[#cite-37|[24]]], the Discontinuous Enrichment Method <span id='citeF-39'></span>[[#cite-39|[25]]] and the Streamline Upwind with Boundary Terms method <span id='citeF-40'></span>[[#cite-40|[26]]].  Basically all the methods make use of a single stabilization parameter which suffices to stabilize the numerical solution along the velocity (streamline) direction. The computation of the streamline stabilization parameter for multidimensional problems is usually based on extensions of the optimal value of the parameter for the simpler 1D case. Specific attempts to design the stability parameter for multidimensional problems in the context of the Petrov-Galerkin formulation have been recently reported <span id='citeF-41'></span>[[#cite-41|[27-31]]]. In general all the stabilized methods above mentioned yield good results for 1D-type problems where the velocity vector is aligned with the direction of the gradient of the solution. However, the use of a single stabilization parameter is  insufficient to provide stabilized solutions in the vicinity of sharp gradients not aligned with the velocity direction, which may appear at the interior of the domain or at boundary layers. The usual remedy for these situations is to use the so called “shock capturing” or “discontinuity-capturing” schemes <span id='citeF-49'></span>[[#cite-49|[32-36]]]. These methods basically add additional transverse diffusion terms in selected elements in order to correct the undershoots and overshoots yielded by the Petrov-Galerkin solution in the sharp gradient zones. Usually the new shock capturing diffusion terms depend on the gradient of the solution and the scheme becomes non linear.
 The aim of this work is to develop a general finite element formulation which can provide stabilized numerical results for all range of convective-diffusive problems. The new formulation therefore has the necessary intrinsic features to deal with streamline-type instabilities, as well as with the typical undershoots and overshoots in the vicinity of sharp gradients at different angles with the velocity direction. The new formulation thus provides a unified theoretical and computational framework incorporating all the ingredients of the traditional Petrov-Galerkin and shock-capturing methods.
-The formulation proposed is based in the so called Finite Calculus (FIC) approach <span id='citeF-50'></span>[[#cite-50|50]]. The FIC method is based in expressing the equation of balance of fluxes in a domain of finite size. This introduces additional stabilizing terms in the differential equations of the infinitessimal theory which are a function of the balance domain dimensions. These dimensions are termed ''characteristic length parameters'' and play a key role in the stabilization process.
+The formulation proposed is based in the so called Finite Calculus (FIC) approach <span id='citeF-50'></span>[[#cite-50|[37,38]]]. The FIC method is based in expressing the equation of balance of fluxes in a domain of finite size. This introduces additional stabilizing terms in the differential equations of the infinitessimal theory which are a function of the balance domain dimensions. These dimensions are termed ''characteristic length parameters'' and play a key role in the stabilization process.
-The modified governing equations lead to stabilized numerical schemes using whatever numerical method. It is interesting that many of the stabilized FEM can be recovered using the FIC formulation. The FIC method has been successfully applied to problems of convection-diffusion <span id='citeF-50'></span>[[#cite-50|50]], convection-diffusion-absorption [42,43], incompressible fluid flow <span id='citeF-47'></span>[[#cite-47|47]] and incompressible solid mechanics [48,49].
+The modified governing equations lead to stabilized numerical schemes using whatever numerical method. It is interesting that many of the stabilized FEM can be recovered using the FIC formulation. The FIC method has been successfully applied to problems of convection-diffusion <span id='citeF-50'></span>[[#cite-50|[37-41]]], convection-diffusion-absorption [42,43], incompressible fluid flow <span id='citeF-47'></span>[[#cite-47|[44-47]]] and incompressible solid mechanics [48,49].
-The key to the stabilization in the FIC method is to choose adequately the characteristic length parameters  for each element. For 1D problems the standard optimal  and critical values of the single stabilization parameter can be found, as shown in Section 3. For 2D/3D problems the characteristic length parameters along each space direction are grouped in a characteristic length vector <math display="inline">h</math>. Here the key issue is to choose correctly the direction of this vector. It is interesting that if the direction of <math display="inline">h</math> is taken parallel to that of the velocity, then the resulting FIC stabilized equations ''coincide'' with that of the standard SUPG method <span id='citeF-50'></span>[[#cite-50|50]].
+The key to the stabilization in the FIC method is to choose adequately the characteristic length parameters  for each element. For 1D problems the standard optimal  and critical values of the single stabilization parameter can be found, as shown in Section 3. For 2D/3D problems the characteristic length parameters along each space direction are grouped in a characteristic length vector <math display="inline">h</math>. Here the key issue is to choose correctly the direction of this vector. It is interesting that if the direction of <math display="inline">h</math> is taken parallel to that of the velocity, then the resulting FIC stabilized equations ''coincide'' with that of the standard SUPG method <span id='citeF-50'></span>[[#cite-50|[37,38]]].
-In previous work of the authors, transverse sharp gradients to the velocity direction were accurately captured  by computing iteratively the direction and modulus of vector <math display="inline">h</math> which minimizes a residual norm <span id='citeF-50'></span>[[#cite-50|50]]. In <span id='citeF-51'></span>[[#cite-51|51]] accurate stabilized solutions for convection-diffusion problems with sharp gradients were obtaining by spliting the characteristic lenght vector <math display="inline">{h}</math> as the sum of two vectors. The first vector is chosen aligned with the velocity direction, hence yielding the standard SUPG terms, while the second vector is taken parallel to the direction of the gradient of the solution.
+In previous work of the authors, transverse sharp gradients to the velocity direction were accurately captured  by computing iteratively the direction and modulus of vector <math display="inline">h</math> which minimizes a residual norm <span id='citeF-50'></span>[[#cite-50|[37-41]]]. In <span id='citeF-51'></span>[[#cite-51|[38]]] accurate stabilized solutions for convection-diffusion problems with sharp gradients were obtaining by spliting the characteristic lenght vector <math display="inline">{h}</math> as the sum of two vectors. The first vector is chosen aligned with the velocity direction, hence yielding the standard SUPG terms, while the second vector is taken parallel to the direction of the gradient of the solution.
 In this work a slightly different and more consistent method is proposed. We have found that the key to the general stabilization algorithm for convection-diffusion problems is to express the  FIC balance equation taking as  coordinate axes the principal curvature directions of the solution. These equations contain the necessary additional diffusion to stabilize the numerical solution ''in all situations''. It is interesting that when one of the principal curvature directions coincides with the velocity direction, then the  classical SUPG method is recovered.
@@ Line 76: / Line 61: @@
 |}
-where <math display="inline">\Gamma _\phi </math> and <math display="inline">\Gamma _q</math> are the Dirichlet and Neumann boundaries where the variable <math display="inline">\phi </math> and the outgoing normal diffusive flux are prescribed to values <math display="inline">\phi ^p</math> and <math display="inline">q_n^p</math>, respectively. The modified equation (2b) is obtained by invoking higher order balance of fluxes in a finite domain next to the Neumann boundary <span id='citeF-50'></span>[[#cite-50|50]]. In above equations
+where <math display="inline">\Gamma _\phi </math> and <math display="inline">\Gamma _q</math> are the Dirichlet and Neumann boundaries where the variable <math display="inline">\phi </math> and the outgoing normal diffusive flux are prescribed to values <math display="inline">\phi ^p</math> and <math display="inline">q_n^p</math>, respectively. The modified equation (2b) is obtained by invoking higher order balance of fluxes in a finite domain next to the Neumann boundary <span id='citeF-50'></span>[[#cite-37|[37,38]]]. In above equations
 {| class="formulaSCP" style="width: 100%; text-align: left;"
@@ Line 88: / Line 73: @@
 |}
-where <math display="inline">u</math> is the divergence-free velocity vector, <math display="inline">{\textbf D}</math> is the diffusion matrix, <math display="inline">{\boldsymbol \nabla }</math> is the gradient operator, <math display="inline">Q</math> is the external source term and <math display="inline">{\textbf n}</math> is the normal vector. Vector <math display="inline">{\textbf h}</math>  is the ''characteristic length'' vector. For 2D problems <math display="inline">{\textbf h}=[h_x, h_y]^T</math>, where <math display="inline">h_x</math> and <math display="inline">h_y</math> are characteristic distances along the sides of the rectangular domain where higher order balance of fluxes is enforced <span id='citeF-50'></span>[[#cite-50|50]].
+where <math display="inline">u</math> is the divergence-free velocity vector, <math display="inline">{\textbf D}</math> is the diffusion matrix, <math display="inline">{\boldsymbol \nabla }</math> is the gradient operator, <math display="inline">Q</math> is the external source term and <math display="inline">{\textbf n}</math> is the normal vector. Vector <math display="inline">{\textbf h}</math>  is the ''characteristic length'' vector. For 2D problems <math display="inline">{\textbf h}=[h_x, h_y]^T</math>, where <math display="inline">h_x</math> and <math display="inline">h_y</math> are characteristic distances along the sides of the rectangular domain where higher order balance of fluxes is enforced <span id='citeF-50'></span>[[#cite-37|[37,38]]].
-The underlined terms in Equations (1) and (2b) introduce the necessary stabilization in the numerical solution of the convective-diffusive problem  <span id='citeF-50'></span>[[#cite-50|50]]. The finite element formulation will be presented next.
+The underlined terms in Equations (1) and (2b) introduce the necessary stabilization in the numerical solution of the convective-diffusive problem  <span id='citeF-50'></span>[[#cite-37|[37-41]]]. The finite element formulation will be presented next.
 ===2.1 Finite element discretization===
@@ Line 106: / Line 91: @@
 |}
-where <math display="inline">N_i</math> are the shape functions and <math display="inline"> \phi _i</math> are the nodal values of the approximate function <math display="inline">\hat \phi </math> <span id='citeF-1'></span>[[#cite-1|1]].
+where <math display="inline">N_i</math> are the shape functions and <math display="inline"> \phi _i</math> are the nodal values of the approximate function <math display="inline">\hat \phi </math> <span id='citeF-1'></span>[[#cite-1|[1]]].
 Application of the Galerkin FE method to Equations (1) and (2b) gives, after integrating by parts the term <math display="inline"> {h}^T{\boldsymbol \nabla }\hat r</math>
@@ Line 211: / Line 196: @@
 |}
-It can be shown that the definition of <math display="inline">\bar {\textbf D}</math> of Eq.(11) is equivalent to introducing an artificial diffusion of value <math display="inline">{h\vert {\textbf u}\vert \over 2}</math> along the streamlines <span id='citeF-1'></span>[[#cite-1|1]] and this explains the name of the SUPG method. The element characteristic length <math display="inline">h</math> (or the value of <math display="inline">\tau </math>) is computed in practice  by  heuristic linear and non-linear extensions of the optimal expression for the 1D problem <span id='citeF-5'></span>[[#cite-5|5]].
+It can be shown that the definition of <math display="inline">\bar {\textbf D}</math> of Eq.(11) is equivalent to introducing an artificial diffusion of value <math display="inline">{h\vert {\textbf u}\vert \over 2}</math> along the streamlines <span id='citeF-1'></span>[[#cite-1|[1,7,9]]] and this explains the name of the SUPG method. The element characteristic length <math display="inline">h</math> (or the value of <math display="inline">\tau </math>) is computed in practice  by  heuristic linear and non-linear extensions of the optimal expression for the 1D problem <span id='citeF-5'></span>[[#cite-5|[5-9,27-31]]].
 In the following sections we will denote by the standard SUPG method to the stabilization procedure leading to a balancing diffusion matrix <math display="inline">\bar {\textbf D}</math> given by Eq.(11).
-It is important to note that the SUPG expression is a ''particular case'' of the more general FIC formulation. This explains the limitations of the SUPG method to provide stabilized numerical results in the vicinity of sharp gradients of the solution transverse to the flow direction and the need to introduce in these cases additional shock capturing terms (see Section 5). However, in the FIC formulation the direction of <math display="inline">{\textbf h}</math> is arbitrary and not necessarily coincident with that of <math display="inline">{\textbf u}</math>. The components of <math display="inline">{\textbf h}</math>  introduce the necessary stabilization along the streamlines and the transverse directions to the flow. In this manner, the FIC method reproduces the best-features of the so-called stabilized discontinuity-capturing schemes <span id='citeF-49'></span>[[#cite-49|49]].
+It is important to note that the SUPG expression is a ''particular case'' of the more general FIC formulation. This explains the limitations of the SUPG method to provide stabilized numerical results in the vicinity of sharp gradients of the solution transverse to the flow direction and the need to introduce in these cases additional shock capturing terms (see Section 5). However, in the FIC formulation the direction of <math display="inline">{\textbf h}</math> is arbitrary and not necessarily coincident with that of <math display="inline">{\textbf u}</math>. The components of <math display="inline">{\textbf h}</math>  introduce the necessary stabilization along the streamlines and the transverse directions to the flow. In this manner, the FIC method reproduces the best-features of the so-called stabilized discontinuity-capturing schemes <span id='citeF-49'></span>[[#cite-32|[32-36]]].
 ==3 COMPUTATION OF THE CHARACTERISTIC LENGTH FOR THE 1D CONVECTION-DIFFUSION EQUATION==
@@ Line 266: / Line 251: @@
 where <math display="inline">l^e</math> is the element length and <math display="inline">x_L</math> is the coordinate of the end point of the 1D domain where the outgoing flux is prescribed to a value <math display="inline">q^p</math>. In Eq.(14a) the space derivatives of <math display="inline">h</math> have been neglected as <math display="inline">h</math> is assumed to be constant within each element.
-Note that the  FIC method introduces a diffusion of value <math display="inline">{hu\over 2}</math>. This term is analogous to that provided by artificial diffusion and upwinding methods <span id='citeF-1'></span>[[#cite-1|1]].
+Note that the  FIC method introduces a diffusion of value <math display="inline">{hu\over 2}</math>. This term is analogous to that provided by artificial diffusion and upwinding methods <span id='citeF-1'></span>[[#cite-1|[1,2]]].
 The characteristic length parameter <math display="inline">h</math> for each element can  be made proportional to the element length as <math display="inline">h= \alpha l^e</math> where <math display="inline">\alpha </math> is an element stabilization parameter. A typical stencil for a mesh of uniform size 1D elements of length <math display="inline">l^e =l</math> is written as
@@ Line 322: / Line 307: @@
 Above well know results from the finite element and finite difference literature will be useful for the general multidimensional case described next.
-===Remark 1===
+ '''Remark 1'''. The stabilizing diffusion of Eq.(15) can be expressed in terms of the equivalent intrinsic time parameter <math display="inline">\tau </math> by replacing the term <math display="inline">\alpha {ul\over 2}</math> by <math display="inline">\tau u^2</math> where <math display="inline">\tau = \alpha {l\over 2u}</math>. The optimal and critical values of <math display="inline">\tau </math> are therefore obtained by dividing by <math display="inline">{l\over 2u}</math> the expressions of <math display="inline">\alpha _{opt}</math> and <math display="inline">\alpha _c</math> of Eqs.(16) and (18), respectively.
-The stabilizing diffusion of Eq.(15) can be expressed in terms of the equivalent intrinsic time parameter <math display="inline">\tau </math> by replacing the term <math display="inline">\alpha {ul\over 2}</math> by <math display="inline">\tau u^2</math> where <math display="inline">\tau = \alpha {l\over 2u}</math>. The optimal and critical values of <math display="inline">\tau </math> are therefore obtained by dividing by <math display="inline">{l\over 2u}</math> the expressions of <math display="inline">\alpha _{opt}</math> and <math display="inline">\alpha _c</math> of Eqs.(16) and (18), respectively.
 ==4 COMPUTATION OF THE CHARACTERISTIC LENGTH VECTOR==
 We present here a general procedure to compute the characteristic length vector <math display="inline">h</math> for convection-diffusion problems. For the sake of preciseness the method is explained for 2D problems although it is equally applicable to 3D problems.
+<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_426828278-Figure1.png|351px|Global axes (x,y) and principal curvature axes (ξ,η)]]
+|- style="text-align: center; font-size: 75%;"
+| colspan="1" | '''Figure 1:''' Global axes (<math>x,y</math>) and principal curvature axes (<math>\xi ,\eta </math>)
+|}
 Let us write down the FIC balance equation  in the principal curvature axes of the solution <math display="inline">\xi ,\eta </math> (Figure 1). For simplicity we consider the 2D sourceless case (<math display="inline">Q=0</math>) with an isotropic diffusion defined by a constant diffusion parameter <math display="inline">k</math>.  The FIC balance equation is
@@ Line 365: / Line 357: @@
 {| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>-u_\xi {\partial \phi  \over \partial \xi }-u_\eta {\partial \phi  \over \partial \eta }+\left(k + {u_\xi h_\xi \over 2}\right){\partial ^2\phi \over \partial \xi ^2} + \left(k + {u_\eta h_\eta \over 2}\right){\partial ^2\phi \over \partial \eta ^2}  - {k\over 2} \left({h_\xi {h_\xi \over 2}{\partial ^3 \phi \over \partial \xi ^3}+  h_\eta {h_\eta \over 2} {\partial ^3 \phi \over \partial \eta ^3}\right)=0 }</math>
+| style="text-align: center;" | <math> -u_\xi {\partial \phi  \over \partial \xi }-u_\eta {\partial \phi  \over \partial \eta }+\left(k + {u_\xi h_\xi \over 2}\right){\partial ^2\phi \over \partial \xi ^2} + \left(k + {u_\eta h_\eta \over 2}\right){\partial ^2\phi \over \partial \eta ^2}  - {k\over 2} \left(h_\xi {\partial ^3 \phi \over \partial \xi ^3}+  h_\eta  {\partial ^3 \phi \over \partial \eta ^3}\right)=0</math>
 |}
 | style="width: 5px;text-align: right;white-space: nowrap;" | (21)
 |}
-}
 We can see clearly from Eq.(21) that the FIC governing equations introduce orthotopic diffusion parameters  of values <math display="inline">{u_\xi h_\xi \over 2}</math> and <math display="inline"> {u_\eta h_\eta \over 2}</math> along the <math display="inline">\xi </math> and <math display="inline">\eta </math> axes, respectively. Also note that the last term of Eq.(21) will vanish after discretization for linear elements.
@@ Line 381: / Line 372: @@
 {| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>- {u}^{\prime T} {\boldsymbol \nabla }^\prime \phi + {\boldsymbol \nabla }^{\prime T} ({D}+\bar {D}^\prime ) {\boldsymbol \nabla }^{\prime }\phi =0 </math>
+| style="text-align: center;" | <math>- {\textbf u}^{\prime T} {\boldsymbol \nabla }^\prime \phi + {\boldsymbol \nabla }^{\prime T} ({\textbf D}+\bar {\textbf D}^\prime ) {\boldsymbol \nabla }^{\prime }\phi =0 </math>
 |}
 | style="width: 5px;text-align: right;white-space: nowrap;" | (22)
 |}
-where <math display="inline">{u}^\prime =[u_\xi , u_\eta ]^T</math>, <math display="inline">{\boldsymbol \nabla }^\prime = \left[{\partial \over \partial \xi }, {\partial \over \partial \eta }\right]^T</math>, <math display="inline">D</math> is the “physical” isotropic diffusion matrix and <math display="inline">\bar {D}'</math> is the balancing diffusion matrix in the local axes <math display="inline">\xi </math> and <math display="inline">\eta </math>. The form of these matrices is
+where <math display="inline">{\textbf u}^\prime =[u_\xi , u_\eta ]^T</math>, <math display="inline">{\boldsymbol \nabla }^\prime = \left[{\partial \over \partial \xi }, {\partial \over \partial \eta }\right]^T</math>, <math display="inline">\textbf D</math> is the “physical” isotropic diffusion matrix and <math display="inline">\bar {\textbf D}'</math> is the balancing diffusion matrix in the local axes <math display="inline">\xi </math> and <math display="inline">\eta </math>. The form of these matrices is
 {| class="formulaSCP" style="width: 100%; text-align: left;"
@@ Line 393: / Line 384: @@
 {| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>{D}= \left[\begin{array}{cc}k &0\\0&k\end{array}\right]\qquad \hbox{and} \qquad \bar {D}^\prime = \left[\begin{array}{cc}\displaystyle{u_\xi h_\xi \over 2} & 0\\ 0 & \displaystyle{u_\eta h_\eta \over 2}\end{array}\right] </math>
+| style="text-align: center;" | <math>{\textbf D}= \left[\begin{array}{cc}k &0\\0&k\end{array}\right]\qquad \hbox{and} \qquad \bar {\textbf D}^\prime = \left[\begin{array}{cc}\displaystyle{u_\xi h_\xi \over 2} & 0\\ 0 & \displaystyle{u_\eta h_\eta \over 2}\end{array}\right] </math>
 |}
 | style="width: 5px;text-align: right;white-space: nowrap;" | (23)
@@ Line 405: / Line 396: @@
 {| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>{u}'=\left\{\begin{array}{c}u_\xi \\ u_\eta \end{array}\right\}= {T} {u}  \quad \hbox{with}\quad {T}= \left[\begin{array}{cc}c_\alpha & s_\alpha \\ - s_\alpha & c_\alpha \end{array}\right]\quad ,\quad {u}=\left\{\begin{array}{c}u\\v\end{array}\right\} </math>
+| style="text-align: center;" | <math>{\textbf u}'=\left\{\begin{array}{c}u_\xi \\ u_\eta \end{array}\right\}= {\textbf T} {\textbf u}  \quad \hbox{with}\quad {\textbf T}= \left[\begin{array}{cc}c_\alpha & s_\alpha \\ - s_\alpha & c_\alpha \end{array}\right]\quad ,\quad {\textbf u}=\left\{\begin{array}{c}u\\v\end{array}\right\} </math>
 |}
 | style="width: 5px;text-align: right;white-space: nowrap;" | (24)
@@ Line 465: / Line 456: @@
 |}
-In Eq.(27a) <math display="inline">{u}_\xi </math> and <math display="inline">{u}_\eta </math> contain the global components of the  velocity vectors <math display="inline">\vec u_\xi </math> and <math display="inline">\vec u_\eta </math>, respectively. For triangles <math display="inline">{d}_j</math> are the element sides vectors, whereas for quadrilaterals <math display="inline">{d}_j</math> are the element diagonals vectors.
+In Eq.(27a) <math display="inline">{\textbf u}_\xi </math> and <math display="inline">{\textbf u}_\eta </math> contain the global components of the  velocity vectors <math display="inline">\vec u_\xi </math> and <math display="inline">\vec u_\eta </math>, respectively. For triangles <math display="inline">{\textbf d}_j</math> are the element sides vectors, whereas for quadrilaterals <math display="inline">{\textbf d}_j</math> are the element diagonals vectors.
 The next step is to transform Eq.(22) to global axes <math display="inline">x,y</math>. The resulting equation is written as
@@ Line 474: / Line 465: @@
 {| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>-{u}^T {\boldsymbol \nabla }\phi +{\boldsymbol \nabla }^T {D}_G {\boldsymbol \nabla }\phi=0 </math>
+| style="text-align: center;" | <math>-{\textbf u}^T {\boldsymbol \nabla }\phi +{\boldsymbol \nabla }^T {\textbf D}_G {\boldsymbol \nabla }\phi=0 </math>
 |}
 | style="width: 5px;text-align: right;white-space: nowrap;" | (28)
 |}
-where the global diffusion matrix <math display="inline">{D}_G</math> is
+where the global diffusion matrix <math display="inline">{\textbf D}_G</math> is
 {| class="formulaSCP" style="width: 100%; text-align: left;"
@@ Line 486: / Line 477: @@
 {| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>{D}_G={D}+\bar{D}</math>
+| style="text-align: center;" | <math>{\textbf D}_G={\textbf D}+\bar{\textbf D}</math>
 |}
 | style="width: 5px;text-align: right;white-space: nowrap;" | (29a)
 |}
-The isotropic diffusion matrix <math display="inline">D</math> is given by Eq.(23) and  the global balancing diffusion matrix <math display="inline">\bar{D}</math> is
+The isotropic diffusion matrix <math display="inline">\textbf D</math> is given by Eq.(23) and  the global balancing diffusion matrix <math display="inline">\bar{\textbf D}</math> is
 {| class="formulaSCP" style="width: 100%; text-align: left;"
@@ Line 498: / Line 489: @@
 {| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>\bar {D} = {T}^T \bar {D}^\prime {T}</math>
+| style="text-align: center;" | <math>\bar {\textbf D} = {\textbf T}^T \bar {\textbf D}^\prime {\textbf T}</math>
 |}
 | style="width: 5px;text-align: right;white-space: nowrap;" | (29b)
@@ Line 505: / Line 496: @@
 where the transformation matrix '''T''' is given in Eq.(24).
-===Remark 2===
+'''Remark 2'''. We can write the local balancing diffusion matrix <math display="inline">\bar{\textbf D}'</math> of Eq.(23) for <math display="inline">k_\xi{-}k_\eta >0</math>  as<br />
-We can write the local balancing diffusion matrix <math display="inline">\bar{D}'</math> of Eq.(23) for <math display="inline">k_\xi{-}k_\eta >0</math>  as
 {| class="formulaSCP" style="width: 100%; text-align: left;"
 |-
@@ Line 514: / Line 502: @@
 {| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>\bar{D}' ={1\over 2}\left[\begin{array}{cc} (k_\xi -k_\eta ) &0\\ 0&0\end{array}\right]+ {1\over 2}\left[\begin{array}{cc} k_\eta  &0\\ 0&k_\eta \end{array}\right]=\bar{D}_\xi + \bar{D}_{iso}</math>
+| style="text-align: center;" | <math>\bar{\textbf D}' ={1\over 2}\left[\begin{array}{cc} (k_\xi -k_\eta ) &0\\ 0&0\end{array}\right]+ {1\over 2}\left[\begin{array}{cc} k_\eta  &0\\ 0&k_\eta \end{array}\right]=\bar{\textbf D}_\xi + \bar{\textbf D}_{iso}</math>
 |}
 | style="width: 5px;text-align: right;white-space: nowrap;" | (30a)
@@ Line 531: / Line 519: @@
 |}
-If <math display="inline">k_\eta </math> is greater than <math display="inline"> k_\xi </math> a similar split is performed with <math display="inline">(k_\eta - k_\xi )</math> in the position 22 of matrix <math display="inline">\bar{D}_\xi </math> and <math display="inline">k_\eta </math> is replaced by <math display="inline">k_\xi </math>  in <math display="inline">\bar{D}_{iso}</math>.
+If <math display="inline">k_\eta </math> is greater than <math display="inline"> k_\xi </math> a similar split is performed with <math display="inline">(k_\eta - k_\xi )</math> in the position 22 of matrix <math display="inline">\bar{\textbf D}_\xi </math> and <math display="inline">k_\eta </math> is replaced by <math display="inline">k_\xi </math>  in <math display="inline">\bar{\textbf D}_{iso}</math>.
-Eq.(30a) shows that the diffusivity matrix <math display="inline">\bar{D}'</math> is equivalent to the sum of a balancing diffusion along the principal curvature direction <math display="inline">\vec {\xi }</math> (or <math display="inline">\vec {\eta }</math> if <math display="inline">k_\eta > k_\xi </math>) and an isotropic diffusion matrix <math display="inline">\bar{D}_{iso}</math>
+Eq.(30a) shows that the diffusivity matrix <math display="inline">\bar{\textbf D}'</math> is equivalent to the sum of a balancing diffusion along the principal curvature direction <math display="inline">\vec {\xi }</math> (or <math display="inline">\vec {\eta }</math> if <math display="inline">k_\eta > k_\xi </math>) and an isotropic diffusion matrix <math display="inline">\bar{\textbf D}_{iso}</math>
 The global balancing diffusion matrix of Eq.(29a) can therefore be written as
@@ Line 542: / Line 530: @@
 {| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>\bar {D} = \bar{D}_{iso} + {T}^T\bar{D}_\xi {T} </math>
+| style="text-align: center;" | <math>\bar {\textbf D} = \bar{\textbf D}_{iso} + {\textbf T}^T\bar{\textbf D}_\xi {\textbf T} </math>
 |}
 | style="width: 5px;text-align: right;white-space: nowrap;" | (31)
 |}
-Clearly, if the principal curvature direction <math display="inline">\vec {\xi }</math> is parallel to the velocity direction, then <math display="inline">u_\eta =0</math>, <math display="inline">k_\eta =0</math>, <math display="inline"> \bar{D}_{iso}={0}</math> and
+Clearly, if the principal curvature direction <math display="inline">\vec {\xi }</math> is parallel to the velocity direction, then <math display="inline">u_\eta =0</math>, <math display="inline">k_\eta =0</math>, <math display="inline"> \bar{\textbf D}_{iso}={0}</math> and
 {| class="formulaSCP" style="width: 100%; text-align: left;"
@@ Line 554: / Line 542: @@
 {| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>\bar{D}= \bar{D}_\xi =\left[\begin{array}{cc}{u_\xi h_\xi \over 2} &0\\0&0\end{array}\right] </math>
+| style="text-align: center;" | <math>\bar{\textbf D}= \bar{\textbf D}_\xi =\left[\begin{array}{cc}{u_\xi h_\xi \over 2} &0\\0&0\end{array}\right] </math>
 |}
 | style="width: 5px;text-align: right;white-space: nowrap;" | (32)
@@ Line 561: / Line 549: @@
 where <math display="inline">u_\xi = \vert {u}\vert </math> and <math display="inline">h_\xi </math> is computed by Eq.(26a). Note that the method coincides with the standard SUPG approach in this case.
-===Remark 3===
+'''Remark 3'''. The global balance diffusion matrix <math display="inline">\bar  {\textbf D}</math> can be also computed from the expression of vector <math display="inline">\textbf h</math> in global axes as
-The global balance diffusion matrix <math display="inline">\bar  {D}</math> can be also computed from the expression of vector <math display="inline">h</math> in global axes as
 {| class="formulaSCP" style="width: 100%; text-align: left;"
@@ Line 570: / Line 556: @@
 {| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>\bar  {D} ={1\over 2} {h}^T {u}\quad \hbox{with }{h}=\left\{\begin{array}{c}h_x\\h_y \end{array}\right\}  = {T}^T {h}' </math>
+| style="text-align: center;" | <math>\bar  {\textbf D} ={1\over 2} {\textbf h}^T {\textbf u}\quad \hbox{with }{\textbf h}=\left\{\begin{array}{c}h_x\\h_y \end{array}\right\}  = {\textbf T}^T {\textbf h}' </math>
 |}
 | style="width: 5px;text-align: right;white-space: nowrap;" | (33)
 |}
-where <math display="inline">{h}'=[h_\xi , h_\eta ]^T</math> and <math display="inline">{T}</math> is the transformation matrix of Eq.(24). The proof of this equivalence is given in the Appendix. In our computation however we have chosen to compute the balancing diffusion matrix <math display="inline">\bar  {D}</math> via Eqs.(23)&#8211;(29).
+where <math display="inline">{\textbf h}'=[h_\xi , h_\eta ]^T</math> and <math display="inline">{\textbf T}</math> is the transformation matrix of Eq.(24). The proof of this equivalence is given in the Appendix. In our computation however we have chosen to compute the balancing diffusion matrix <math display="inline">\bar  {\textbf D}</math> via Eqs.(23)&#8211;(29).
 ===4.1 Computation of the principal curvature axes for linear elements===
@@ Line 586: / Line 572: @@
 {| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>{C}= \left[\begin{array}{cc}\phi _x^{''} & \phi ^{''}_{xy}\\ \phi ^{''}_{xy} & \phi _y^{''}  \end{array}\right]\quad ,\quad \hbox{where } \phi _x^{''}={\partial ^2\phi \over \partial x^2}, \hbox{etc.} </math>
+| style="text-align: center;" | <math>{\textbf C}= \left[\begin{array}{cc}\phi _x^{''} & \phi ^{''}_{xy}\\ \phi ^{''}_{xy} & \phi _y^{''}  \end{array}\right]\quad ,\quad \hbox{where } \phi _x^{''}={\partial ^2\phi \over \partial x^2}, \hbox{etc.} </math>
 |}
 | style="width: 5px;text-align: right;white-space: nowrap;" | (34)
@@ Line 601: / Line 587: @@
 For linear triangles <math display="inline">{\boldsymbol \nabla } \phi </math> is constant within the element. For four node quadrilaterals <math display="inline">{\boldsymbol \nabla } \phi </math> varies linearly. We have assumed in this case that the direction of <math display="inline">\vec \xi </math> is constant within the element and equal to that of vector <math display="inline">{\boldsymbol \nabla } \phi </math> computed at the element center.
-The dependence of the balancing diffusion matrix <math display="inline">\bar{D}</math> with the principal  curvature directions <math display="inline">\vec \xi </math> and <math display="inline">\vec {\eta }</math> introduces a non linearity in the solution process. A simple and effective iterative algorithm is described next.
+The dependence of the balancing diffusion matrix <math display="inline">\bar{\textbf D}</math> with the principal  curvature directions <math display="inline">\vec \xi </math> and <math display="inline">\vec {\eta }</math> introduces a non linearity in the solution process. A simple and effective iterative algorithm is described next.
 ===4.2 General iterative scheme===
@@ Line 607: / Line 593: @@
 A stabilized numerical solution can be found by the following algorithm.
-'''Step 0'''  (SUPG step). At each integration point choose <math display="inline">{}^0{\boldsymbol \xi } ={u}</math>, i.e. the gradient direction is taken coincident with the velocity direction. Compute <math display="inline">{}^0{\boldsymbol \eta }, {}^1\bar  {D}'</math>, <math display="inline">{}^0\bar  {D}</math> and <math display="inline">{}^0{D}_G</math> from Eqs.(23&#8211;29). ''The  expression of the balancing diffusion matrix coincides  now precisely with the standard (linear) SUPG form''.
+'''Step 0'''  (SUPG step). At each integration point choose <math display="inline">{}^0{\boldsymbol \xi } ={\textbf u}</math>, i.e. the gradient direction is taken coincident with the velocity direction. Compute <math display="inline">{}^0{\boldsymbol \eta }, {}^1\bar  {\textbf D}'</math>, <math display="inline">{}^0\bar  {\textbf D}</math> and <math display="inline">{}^0{\textbf D}_G</math> from Eqs.(23&#8211;29). ''The  expression of the balancing diffusion matrix coincides  now precisely with the standard (linear) SUPG form''.
 Solve for <math display="inline">{}^0\bar{\boldsymbol \phi }</math>.
@@ Line 613: / Line 599: @@
 Verify that the solution is stable. This can be performed by verifying that there are not undershoots or overshoots in the numerical results with respect to the expected physical values. If the SUPG solution is unstable, then implement the following iterative scheme.
-For each iteration: '''Step 1'''  Compute at the element center. <math display="inline">{}^1{\boldsymbol \xi } = {\boldsymbol \nabla }^0\bar \phi </math>. Then compute <math display="inline">{}^1{\boldsymbol \eta } , {}^1\bar  {D}'</math>, <math display="inline">{}^1\bar  {D}</math> and <math display="inline">{}^1{D}_G</math>.
+For each iteration:
+'''Step 1'''  Compute at the element center. <math display="inline">{}^1{\boldsymbol \xi } = {\boldsymbol \nabla }^0\bar \phi </math>. Then compute <math display="inline">{}^1{\boldsymbol \eta } , {}^1\bar  {\textbf D}'</math>, <math display="inline">{}^1\bar  {\textbf D}</math> and <math display="inline">{}^1{\textbf D}_G</math>.
 Solve for <math display="inline">{}^1\bar{\boldsymbol \phi }</math>.
@@ Line 635: / Line 623: @@
 If condition (35) is not satisfied, start a new iteration and repeat steps 1 and 2 until convergence. Indexes 0 and 1 are replaced now by <math display="inline">i-1</math> and <math display="inline">i</math>, respectively.
-===Remark 4===
+'''Remark 4'''. The convergence of the iterative scheme of the previous section obviously depends on the “quality” of the standard SUPG solution obtained in Step 0. Clearly when boundary layers or  sharp transverse internal layers dominate the solution, the non linearity of the processes increases. The authors have found that the convergence of the scheme can be substantially improved in some problems if the following expression for the stabilizing dissipation matrix is used in Step 1
-The convergence of the iterative scheme of the previous section obviously depends on the “quality” of the standard SUPG solution obtained in Step 0. Clearly when boundary layers or  sharp transverse internal layers dominate the solution, the non linearity of the processes increases. The authors have found that the convergence of the scheme can be substantially improved in some problems if the following expression for the stabilizing dissipation matrix is used in Step 1
 {| class="formulaSCP" style="width: 100%; text-align: left;"
@@ Line 644: / Line 630: @@
 {| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>{D}_G = \beta ^i \bar{D}_G + (1-\beta ) {}^{i-1}\bar{D}_G  </math>
+| style="text-align: center;" | <math>{\textbf D}_G = \beta ^i \bar{\textbf D}_G + (1-\beta ) {}^{i-1}\bar{\textbf D}_G  </math>
 |}
 | style="width: 5px;text-align: right;white-space: nowrap;" | (36)
@@ Line 653: / Line 639: @@
 Convergence of the iterative scheme was achieved in just two iteration in Examples 6.1&#8211;6.3 for values of <math display="inline">\beta </math> between 0.8 and 1. The best solution in Example 6.4 was found in five iterations for <math display="inline">\beta =0.3</math>. This is due to the higher non linearity of the process induced by the transverse boundary layers in this case.
-===Remark 5===
+'''Remark 5'''. The direct iterative scheme presented in Section 4.2 is obviously not the only option to solve the non linear equations resulting from the dependence of the stabilizing diffusion matrix with the principal curvatures of the solution.
-The direct iterative scheme presented in Section 4.2 is obviously not the only option to solve the non linear equations resulting from the dependence of the stabilizing diffusion matrix with the principal curvatures of the solution.
-A simple and economical alternative is to use a time relaxation technique based in the solution of a pseudo transient problem with a forward Euler scheme and a diagonal mass matrix. This technique was used by Codina in <span id='citeF-11'></span>[[#cite-11|11]] for solving a similar class of convection-diffusion problems.
+A simple and economical alternative is to use a time relaxation technique based in the solution of a pseudo transient problem with a forward Euler scheme and a diagonal mass matrix. This technique was used by Codina in <span id='citeF-11'></span>[[#cite-9|[9,33]]] for solving a similar class of convection-diffusion problems.
 ==5 EQUIVALENCE WITH SUPG AND SHOCK-CAPTURING TECHNIQUES==
@@ Line 667: / Line 651: @@
 Despite their effectiveness, most shock capturing methods lack of a general theoretical framework and can be viewed as “ad hoc” extensions of the SUPG method. A serious drawback of these methods is that, in some cases, the amount of transverse diffusion can not be properly controlled and    hence the final solution is either underdiffusive or overdiffusive.
-The general stabilized formulation we propose here includes the best features of the SUPG and shock capturing techniques. The iterative scheme in fact follows the lines of the standard shock capturing methods <span id='citeF-49'></span>[[#cite-49|49]]. The innovation in the proposed formulation  lays in the definition of the shock capturing terms, which are deduced form the FIC governing equations written in the principal curvature directions. Numerical results obtained in all examples show that the method preserves the quality of the numerical solution in all situations by introducing the necessary  streamline and transverse dissipation.
+The general stabilized formulation we propose here includes the best features of the SUPG and shock capturing techniques. The iterative scheme in fact follows the lines of the standard shock capturing methods <span id='citeF-49'></span>[[#cite-32|[32-36]]]. The innovation in the proposed formulation  lays in the definition of the shock capturing terms, which are deduced form the FIC governing equations written in the principal curvature directions. Numerical results obtained in all examples show that the method preserves the quality of the numerical solution in all situations by introducing the necessary  streamline and transverse dissipation.
-A similar shock capturing method of this kind was proposed by Oñate <span id='citeF-51'></span>[[#cite-51|51]] in the context of the FIC approach. There the characteristic length vector was split as the sum of a  vector along the velocity direction (the linear SUPG term) and a vector in the direction of the solution gradient (the non linear schock capturing term). Good results were obtained with this method for a variety of advection-diffusion problems <span id='citeF-51'></span>[[#cite-51|51]]. The formulation we present here is more general as it encompasses the SUPG method and the shock capturing strategy in a unified approach.
+A similar shock capturing method of this kind was proposed by Oñate <span id='citeF-51'></span>[[#cite-38|[38]]] in the context of the FIC approach. There the characteristic length vector was split as the sum of a  vector along the velocity direction (the linear SUPG term) and a vector in the direction of the solution gradient (the non linear schock capturing term). Good results were obtained with this method for a variety of advection-diffusion problems <span id='citeF-51'></span>[[#cite-38|[38]]]. The formulation we present here is more general as it encompasses the SUPG method and the shock capturing strategy in a unified approach.
 We have shown in Remark 2 that the proposed approach coincides with the standard SUPG method (with no discontinuity capturing) when the principal direction <math display="inline">\vec {\xi }</math> is parallel to the velocity direction. This analogy is discussed further next.
@@ Line 680: / Line 664: @@
 {| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>{h}= h_s {{u}\over \vert {u}\vert } +h_t {{\boldsymbol \nabla }\hat \phi \over \vert {\boldsymbol \nabla } \hat \phi \vert } </math>
+| style="text-align: center;" | <math>{\textbf h}= h_s {{\textbf u}\over \vert {\textbf u}\vert } +h_t {{\boldsymbol \nabla }\hat \phi \over \vert {\boldsymbol \nabla } \hat \phi \vert } </math>
 |}
 | style="width: 5px;text-align: right;white-space: nowrap;" | (37)
@@ Line 694: / Line 678: @@
 {| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>h_s = {h_x \hat \phi _y^\prime - h_y \hat \phi _x^\prime \over u \hat \phi _y^\prime - v \hat \phi _x^\prime }\vert {u}\vert </math>
+| style="text-align: center;" | <math>h_s = {h_x \hat \phi _y^\prime - h_y \hat \phi _x^\prime \over u \hat \phi _y^\prime - v \hat \phi _x^\prime }\vert {\textbf u}\vert </math>
 |}
 | style="width: 5px;text-align: right;white-space: nowrap;" | (38.a)
@@ Line 720: / Line 704: @@
 {| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>c_\alpha = {u\over \vert {u}\vert }, \quad s_\alpha ={v\over \vert {v}\vert }, \quad h_x = {uh_\xi  \over \vert {u}\vert }, \quad h_y ={vh_\xi  \over \vert {u}\vert }</math>
+| style="text-align: center;" | <math>c_\alpha = {u\over \vert {\textbf u}\vert }, \quad s_\alpha ={v\over \vert {\textbf v}\vert }, \quad h_x = {uh_\xi  \over \vert {\textbf u}\vert }, \quad h_y ={vh_\xi  \over \vert {\textbf u}\vert }</math>
 |}
 | style="width: 5px;text-align: right;white-space: nowrap;" |  (39a)
@@ Line 739: / Line 723: @@
 and the method coincides in this case with the standard linear SUPG approach, as expected.
-A case of practical interest is when the directions of <math display="inline">u</math> and <math display="inline">{\boldsymbol \nabla }\phi </math> are ''orthogonal''. This occurs for some specific distributions of the velocity and the solution fields or in the vicinity of sharp layers when the mesh is fine enough, or it has some specific orientation. In these cases we have
+A case of practical interest is when the directions of <math display="inline">\textbf u</math> and <math display="inline">{\boldsymbol \nabla }\phi </math> are ''orthogonal''. This occurs for some specific distributions of the velocity and the solution fields or in the vicinity of sharp layers when the mesh is fine enough, or it has some specific orientation. In these cases we have
 {| class="formulaSCP" style="width: 100%; text-align: left;"
@@ Line 746: / Line 730: @@
 {| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>h_\xi = 0\quad ,\quad h_\eta \not =0 \,\, \hbox{(it is computed by Eqs.(25)--(27))}\,\, ,\,\, c_\alpha = - {v\over \vert {u}\vert }\quad \hbox{and} \quad s_\alpha ={u\over \vert {u}\vert } </math>
+| style="text-align: center;" | <math>h_\xi = 0\quad ,\quad h_\eta \not =0 \,\, \hbox{(it is computed by Eqs.(25)--(27))}\,\, ,\,\, c_\alpha = - {v\over \vert {\textbf u}\vert }\quad \hbox{and} \quad s_\alpha ={u\over \vert {\textbf u}\vert } </math>
 |}
 | style="width: 5px;text-align: right;white-space: nowrap;" | (40)
@@ Line 758: / Line 742: @@
 {| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>h_x = - h_\eta   {u\over \vert {u}\vert } \quad ,\quad  h_y = - h_\eta  {v\over \vert {u}\vert } </math>
+| style="text-align: center;" | <math>h_x = - h_\eta   {u\over \vert {\textbf u}\vert } \quad ,\quad  h_y = - h_\eta  {v\over \vert {\textbf u}\vert } </math>
 |}
 | style="width: 5px;text-align: right;white-space: nowrap;" | (41.a)
@@ Line 777: / Line 761: @@
 Once again, the linear SUPG method is  recovered and the stabilized formulation does not introduce any additional transverse dissipation, as should be expected.
-===Remark 6===
+'''Remark 6'''. From above arguments we conclude that the simplified expression <math display="inline">{\textbf h}=h_s {{\textbf u} \over \vert {\textbf u}\vert }</math> with <math display="inline">h_s</math> being a non linear function of the solution leads to a non linear SUPG term. This has a similarity with the non linear SUPG formulation reported in <span id='citeF-28'></span>[[#cite-28|[28,31]]]. We have found that in general this assumption does not suffice to stabilize the solution in presence of high gradients transverse to the solution. In these cases the introduction of the term involving <math display="inline">h_t</math> in Eq.(37) is essential in order to obtain an improved solution.
-From above arguments we conclude that the simplified expression <math display="inline">{h}=h_s {{u} \over \vert {u}\vert }</math> with <math display="inline">h_s</math> being a non linear function of the solution leads to a non linear SUPG term. This has a similarity with the non linear SUPG formulation reported in <span id='citeF-28'></span>[[#cite-28|28]]. We have found that in general this assumption does not suffice to stabilize the solution in presence of high gradients transverse to the solution. In these cases the introduction of the term involving <math display="inline">h_t</math> in Eq.(37) is essential in order to obtain an improved solution.
 ==6 EXAMPLES==
@@ Line 791: / Line 773: @@
 The FIC solution reported for Example 6.4 was obtained in five iterations for a value of <math display="inline">\beta =0.3</math> in Eq.(36).
-We finally note that most of the examples presented have been successfully solved in the past using different shock capturing techniques <span id='citeF-49'></span>[[#cite-49|49]]. In this paper the FIC results are compared with values obtained using the standard (linear) SUPG formulation. The basic aim of the examples presented here is to show that the FIC formulation integrates naturally the best features of both the SUPG and the shock capturing formulations.
+We finally note that most of the examples presented have been successfully solved in the past using different shock capturing techniques <span id='citeF-49'></span>[[#cite-32|[32-36]]]. In this paper the FIC results are compared with values obtained using the standard (linear) SUPG formulation. The basic aim of the examples presented here is to show that the FIC formulation integrates naturally the best features of both the SUPG and the shock capturing formulations.
 ===6.1 Square domain with uniform Dirichlet conditions, upward diagonal velocity and zero source===
 Figure 2 shows the geometry of the analysis domain and the boundary conditions. The domain is a square of side equal to 10. Zero values of the variable <math display="inline">\phi </math> are prescribed at the inlet boundaries <math display="inline">x=0</math> and <math display="inline">y=0</math>, whereas <math display="inline">\phi =10</math> is prescribed at the outlet boundaries <math display="inline">x=10</math> and <math display="inline">y=10</math>. The isotropic diffusion <math display="inline">k=0.01</math> and the velocity vectors are taken parallel to the diagonal of the domain with <math display="inline">{u}={3\sqrt 2\over 2}[1,1]^T</math>. The problem has been solved with the structured mesh of <math display="inline">10\times 10</math> four node rectangles shown in the figure. The Peclet number along the velocity direction is <math display="inline">\gamma =300</math>.
+<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_426828278-Figure2.png|551px|Square domain with uniform Dirichlet conditions, upward diagonal velocity and zero source. SUPG and FIC solutions obtained with a structured mesh of 10×10 linear four node square elements]]
+|- style="text-align: center; font-size: 75%;"
+| colspan="1" | '''Figure 2:''' Square domain with uniform Dirichlet conditions, upward diagonal velocity and zero source. SUPG and FIC solutions obtained with a structured mesh of <math>10\times 10</math> linear four node square elements
+|}
 Results in Figure 2 show contours of the numerical solution for the SUPG and the FIC iterative scheme proposed. FIC<math display="inline">^2</math> in the figures denotes results obtained after two iterations of the scheme described in Section 4.2. A value of <math display="inline">\beta =1.0</math> in Eq.(36) has been chosen. Note that the SUPG method yields incorrect negative values of <math display="inline">\phi </math> in the vicinity of the upper right corner, as expected, whereas the FIC<math display="inline">^2</math> results are physically correct over the domain. A plot of the distribution of the solution along the horizontal center line and the diagonal line is shown. The difference between the FIC and SUPG results along the diagonal line is clearly seen.
 Figure 3 shows results  for the same problem solved with a structurerd mesh of <math display="inline">2\times 10 \times 10</math> three  node triangular elements. The effect of the stabilization process in more visible in this case as the SUPG solution has more undershoots than for four node square elements.
+<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_426828278-Figure3.png|551px|Square domain with uniform Dirichlet conditions, upward diagonal velocity and zero source. SUPG and FIC solutions obtained with a structured mesh of 2×10×10 linear three node triangles]]
+|- style="text-align: center; font-size: 75%;"
+| colspan="1" | '''Figure 3:''' Square domain with uniform Dirichlet conditions, upward diagonal velocity and zero source. SUPG and FIC solutions obtained with a structured mesh of <math>2\times 10\times 10</math> linear three node triangles
+|}
 The convergence of the iterative process using square and triangular elements is plotted in Figure 4 where the evolution of the values of <math display="inline">\bar \phi </math> at three different points within the analysis domain is plotted. The evolution of the norm of Eq.(35) with the number of iterations is also shown. Note that a converged solution (<math display="inline">\varepsilon <10^{-3}</math>) is obtained in just two iterations.
+<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_426828278-Figure4.png|551px|Square domain with uniform Dirichlet conditions, upward diagonal velocity and zero source. Convergence of the FIC solution for square and triangular elements]]
+|- style="text-align: center; font-size: 75%;"
+| colspan="1" | '''Figure 4:''' Square domain with uniform Dirichlet conditions, upward diagonal velocity and zero source. Convergence of the FIC solution for square and triangular elements
+|}
 Figures 5 and 6 show the solution of a similar problem using the same number of four node rectangular elements and linear triangles with a geometrical aspect ratio of 2:1. The SUPG solution of both cases is clearly non physical and the FIC iterative process provides again a stabilized and accurate solution in just two iterations.
+<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_426828278-Figure5.png|551px|Square domain with uniform Dirichlet conditions, upward diagonal velocity and zero source. SUPG and FIC solutions obtained with a structured mesh of 10×20 linear four node rectangular elements. Geometrical aspect ratio 2:1]]
+|- style="text-align: center; font-size: 75%;"
+| colspan="1" | '''Figure 5:''' Square domain with uniform Dirichlet conditions, upward diagonal velocity and zero source. SUPG and FIC solutions obtained with a structured mesh of <math>10\times 20</math> linear four node rectangular elements. Geometrical aspect ratio 2:1
+|}
+<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_426828278-Figure6.png|551px|Square domain with uniform Dirichlet conditions, upward diagonal velocity and zero source. SUPG and FIC solutions obtained with a structured mesh of 2×10×20 linear triangles. Geometrical aspect ratio 2:1]]
+|- style="text-align: center; font-size: 75%;"
+| colspan="1" | '''Figure 6:''' Square domain with uniform Dirichlet conditions, upward diagonal velocity and zero source. SUPG and FIC solutions obtained with a structured mesh of <math>2\times 10\times 20</math> linear triangles. Geometrical aspect ratio 2:1
+|}
 ===6.2 Square domain with non uniform Dirichlet conditions, downwards diagonal velocity and zero source===
 The geometry of the problem and the boundary conditions are shown in Figure 7. The side of the square domain has a length equal to one. The velocity vector is <math display="inline">{u} =[5,-9]^T \times 10^{6}</math>. An isotropic diffusion of unit value is assumed. The correct solution is a uniform plateau of value <math display="inline">\phi =100</math> with an internal sharp layer in the vicinity of the left lower region where the solution drops to <math display="inline">\bar \phi =0</math>, a boundary layer at the right hand side and at partial boundary layer at the bottom side where the value of <math display="inline">\phi </math> is prescribed to zero.
+<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_426828278-Figure7.png|551px|Square domain with non uniform Dirichlet conditions, downwards diagonal velocity and zero source. SUPG and FIC solutions obtained with a structured mesh of 20×20 linear four node square elements]]
+|- style="text-align: center; font-size: 75%;"
+| colspan="1" | '''Figure 7:''' Square domain with non uniform Dirichlet conditions, downwards diagonal velocity and zero source. SUPG and FIC solutions obtained with a structured mesh of <math>20\times 20</math> linear four node square elements
+|}
 Figures 7 and 8 show results obtained with structured meshes of <math display="inline">20\times 20</math> four node square elements and <math display="inline">2\times 20\times 20</math> linear triangles, respectively. The SUPG solution in both cases has overshoots in the vicinity of the sharp gradient zones, whereas the FIC iterative scheme yields a stabilized and accurate solution in two iterations. Again <math display="inline">\beta =1</math> in Eq.(36) was chosen. Note that the boundary layers and the internal sharp gradient zone are well captured without overshoots or undershoots with the FIC method using both types of elements.
+<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_426828278-Figure8.png|551px|Square domain with non uniform Dirichlet conditions, downwards diagonal velocity  and zero source. SUPG and FIC solutions obtained with a structured mesh of 2×20×20 three node triangles]]
+|- style="text-align: center; font-size: 75%;"
+| colspan="1" | '''Figure 8:''' Square domain with non uniform Dirichlet conditions, downwards diagonal velocity  and zero source. SUPG and FIC solutions obtained with a structured mesh of <math>2\times 20\times 20</math> three node triangles
+|}
 The convergence of the process for both elements is shown in Figure 9.
+<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_426828278-Figure9.png|551px|Square domain with uniform Dirichlet conditions, upward diagonal velocity and zero source. Convergence of the FIC solution for square and triangular elements]]
+|- style="text-align: center; font-size: 75%;"
+| colspan="1" | '''Figure 9:''' Square domain with uniform Dirichlet conditions, upward diagonal velocity and zero source. Convergence of the FIC solution for square and triangular elements
+|}
 Figures 10 and 11 show results for the same problem using coarse and refined unstructured meshes of  linear quadrilaterals and triangles, respectively. The performance of the FIC algorithm is very good  in both cases despite the difficulty to capture the sharp layers with coarse unstructured meshes.
+<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_426828278-Figure10.png|551px|Square domain with non uniform Dirichlet conditions, downwards diagonal velocity and zero source. SUPG and FIC solutions obtained with an unstructured mesh of 432 linear four node quadrilateral elements]]
+|- style="text-align: center; font-size: 75%;"
+| colspan="1" | '''Figure 10:''' Square domain with non uniform Dirichlet conditions, downwards diagonal velocity and zero source. SUPG and FIC solutions obtained with an unstructured mesh of 432 linear four node quadrilateral elements
+|}
+<div id='img-11'></div>
+{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
+|-
+|[[Image:Draft_Samper_426828278-Figure11.png|551px|Square domain with non uniform Dirichlet conditions, downwards diagonal velocity  and zero source. SUPG and FIC solutions obtained with an unstructured mesh of 780 three node triangles]]
+|- style="text-align: center; font-size: 75%;"
+| colspan="1" | '''Figure 11:''' Square domain with non uniform Dirichlet conditions, downwards diagonal velocity  and zero source. SUPG and FIC solutions obtained with an unstructured mesh of 780 three node triangles
+|}
 The solution is improved in both cases by using a refined unstructured mesh as shown in Figures 12 and 13. Note that the overall solutions obtained with the linear triangle and unstructured meshes are slightly better than those obtained with linear quadrilaterals for a similar number of elements. The improvements are clearer in the vicinity of the boundary layer at the right hand side of the domain.
+<div id='img-12'></div>
+{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
+|-
+|[[Image:Draft_Samper_426828278-Figure12.png|551px|Square domain with non uniform Dirichlet conditions, downwards diagonal velocity and zero source. SUPG and FIC solutions obtained with a refined unstructured mesh of 1826 linear four node quadrilaterals]]
+|- style="text-align: center; font-size: 75%;"
+| colspan="1" | '''Figure 12:''' Square domain with non uniform Dirichlet conditions, downwards diagonal velocity and zero source. SUPG and FIC solutions obtained with a refined unstructured mesh of 1826 linear four node quadrilaterals
+|}
+<div id='img-13'></div>
+{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
+|-
+|[[Image:Draft_Samper_426828278-Figure13.png|551px|Square domain with non uniform Dirichlet conditions, downwards diagonal velocity and zero source. SUPG and FIC solutions obtained with a refined unstructured mesh of 1898 linear triangles]]
+|- style="text-align: center; font-size: 75%;"
+| colspan="1" | '''Figure 13:''' Square domain with non uniform Dirichlet conditions, downwards diagonal velocity and zero source. SUPG and FIC solutions obtained with a refined unstructured mesh of 1898 linear triangles
+|}
 Figures 14 and 15 show results for a similar problem solved now with structured meshes of linear rectangles and triangles with a geometrical aspect ratio of 2:1. The correct solution in this case is a central plateau of <math display="inline"> \phi =10</math> with a boundary layer at the right hand side where the solution drops to <math display="inline">\phi =0</math>, and two sharp gradient zones at both sides of the plateau within the domain. The figures show the contours of <math display="inline">\bar \phi </math> over the analysis domain and the distribution of <math display="inline">\bar \phi </math> along the vertical and horizontal central lines. As expected the SUPG solution yields over and undershoots. Two iterations of the FIC iterative scheme again suffice to provide a correct solution  over the whole domain.
-This problem has been chosen by many authors as a test for stabilized formulations. A comparison of the solution obtained with different shock capturing techniques was reported by Codina <span id='citeF-11'></span>[[#cite-11|11]] using a time relaxation scheme. The FIC results presented here reproduce the best results obtained for this problem with the shock capturing methods analyzed in <span id='citeF-11'></span>[[#cite-11|11]].
+<div id='img-14'></div>
+{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
+|-
+|[[Image:Draft_Samper_426828278-Figure14.png|551px|Square domain with non uniform Dirichlet conditions, upwards diagonal velocity  and zero source. SUPG and FIC solutions obtained with a structured mesh of 10×20 linear four node rectangular elements. Geometrical aspect ratio 2:1]]
+|- style="text-align: center; font-size: 75%;"
+| colspan="1" | '''Figure 14:''' Square domain with non uniform Dirichlet conditions, upwards diagonal velocity  and zero source. SUPG and FIC solutions obtained with a structured mesh of <math>10\times 20</math> linear four node rectangular elements. Geometrical aspect ratio 2:1
+|}
+<div id='img-15'></div>
+{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
+|-
+|[[Image:Draft_Samper_426828278-Figure15.png|551px|Square domain with non uniform Dirichlet conditions, upwards diagonal velocity  and zero source. SUPG and FIC solutions obtained with a structured mesh of 2×10×20 three node triangles. Geometrical aspect ratio 2:1]]
+|- style="text-align: center; font-size: 75%;"
+| colspan="1" | '''Figure 15:''' Square domain with non uniform Dirichlet conditions, upwards diagonal velocity  and zero source. SUPG and FIC solutions obtained with a structured mesh of <math>2\times 10\times 20</math> three node triangles. Geometrical aspect ratio 2:1
+|}
+This problem has been chosen by many authors as a test for stabilized formulations. A comparison of the solution obtained with different shock capturing techniques was reported by Codina <span id='citeF-11'></span>[[#cite-9|[9,33]]] using a time relaxation scheme. The FIC results presented here reproduce the best results obtained for this problem with the shock capturing methods analyzed in <span id='citeF-11'></span>[[#cite-9|[9,33]]].
 ===6.3 Rectangular domain with Neumann and non uniform Dirichlet conditions, rotational  velocity field and zero source===
@@ Line 826: / Line 930: @@
 Figures 16 and 17 present the results obtained with the structured meshes of linear square and triangular elements respectively shown. The contours of <math display="inline">\bar \phi </math> over the domain and the distribution of the solution along an horizontal line are plotted. Results clearly show the incapacity of the SUPG method to correctly capture the boundary layer and the internal gradient region. The solution obtained with two iterations of the FIC scheme presented is physically sound and accurate over the whole analysis domain as shown in the figures. FIC results were once more obtained for a value of <math display="inline">\beta =1</math> in Eq.(36).
+<div id='img-16'></div>
+{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
+|-
+|[[Image:Draft_Samper_426828278-Figure16.png|551px|Rectangular domain with Neumann and  non uniform Dirichlet conditions, rotational velocity field and zero source. SUPG and FIC solutions obtained with a structured mesh of 40×20 four node square elements]]
+|- style="text-align: center; font-size: 75%;"
+| colspan="1" | '''Figure 16:''' Rectangular domain with Neumann and  non uniform Dirichlet conditions, rotational velocity field and zero source. SUPG and FIC solutions obtained with a structured mesh of <math>40\times 20</math> four node square elements
+|}
+<div id='img-17'></div>
+{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
+|-
+|[[Image:Draft_Samper_426828278-Figure17.png|551px|Rectangular domain with Neumann and  non uniform Dirichlet conditions, rotational velocity field and zero source. SUPG and FIC solutions obtained with  a structured mesh of 2×40×20   three node triangles]]
+|- style="text-align: center; font-size: 75%;"
+| colspan="1" | '''Figure 17:''' Rectangular domain with Neumann and  non uniform Dirichlet conditions, rotational velocity field and zero source. SUPG and FIC solutions obtained with  a structured mesh of <math>2\times 40\times 20</math>   three node triangles
+|}
 Figure 18 shows results for the solution of the same problem using coarse and refined unstructured meshes of linear quadrilaterals. Note the accuracy of the FIC stabilized solution, even for the coarse mesh.
+<div id='img-18'></div>
+{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
+|-
+|[[Image:Draft_Samper_426828278-Figure18.png|551px|]]
+|[[Image:Draft_Samper_426828278-Figure18b.png|551px|Rectangular domain with Neumann and  non uniform Dirichlet conditions, rotational velocity field and zero source. SUPG and FIC solutions obtained with two unstructured meshes of four node quadrilaterals. (a) Coarse mesh of 892 elements. (b) Refined mesh of 2054 elements]]
+|- style="text-align: center; font-size: 75%;"
+| colspan="2" | '''Figure 18:''' Rectangular domain with Neumann and  non uniform Dirichlet conditions, rotational velocity field and zero source. SUPG and FIC solutions obtained with two unstructured meshes of four node quadrilaterals. (a) Coarse mesh of 892 elements. (b) Refined mesh of 2054 elements
+|}
 The same problem was finally solved with a unstructured mesh of 1554 linear triangles. The number of elements chosen is quite similar to that used for the structured solution (1600 elements). The solution obtained shown in Figure 19 is quite accurate over the whole domain, despite the relative coarseness of the mesh. This is another prove of the best performance of the linear triangle in unstructured meshes.
+<div id='img-19'></div>
+{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
+|-
+|[[Image:Draft_Samper_426828278-Figure19.png|551px|Rectangular domain with Neumann and  non uniform Dirichlet conditions, rotational velocity field and zero source. SUPG and FIC solutions obtained with  an unstructured mesh of 1554 three node triangles]]
+|- style="text-align: center; font-size: 75%;"
+| colspan="1" | '''Figure 19:''' Rectangular domain with Neumann and  non uniform Dirichlet conditions, rotational velocity field and zero source. SUPG and FIC solutions obtained with  an unstructured mesh of 1554 three node triangles
+|}
 ===6.4 Square domain with homogeneous Dirichlet condition and constant source===
@@ Line 837: / Line 977: @@
 The correct solution of this problem is a linear distribution of <math display="inline">\phi </math> increasing uniformly from a value of zero at the inlet boundary until it reaches a maximum value in the vecinity of the outlet boundary. Three boundary layers exist for this problem, one at the exit boundary and two at the two lateral sides, where the value of <math display="inline">\phi </math> drops drastically to zero. Figures 20 and 21 show results obtained with the two structured meshes of quadrilateral and triangles, respectively. The linear SUPG solution yields oscillations in the direction normal to the lateral layers, although a good resolution of the layer normal to the velocity field is obtained, as expected. The FIC solution for <math display="inline">\beta =1</math> in Eq.(36) yielded a smoooth solution near the boundary layers, free of the lateral oscillations which are present in the SUPG solution. The best FIC solution in this case was however obtained for a value of <math display="inline">\beta =0.3</math> in Eq.(36). Five iterations of the scheme were needed to obtain convergence in this case. The numerical results obtained for both meshes are shown in Figures 20 and 21.
-This problem was also studied by Codina <span id='citeF-11'></span>[[#cite-11|11]] with four node square elements as another benchmark to study the performance of different shock capturing methods. The FIC method proposed here has yielded accurate solutions for both quadrilateral and triangular elements, which are difficult to obtain for many of the existing stabilization procedures, as reported in <span id='citeF-11'></span>[[#cite-11|11]].
+<div id='img-20'></div>
+{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
+|-
+|[[Image:Draft_Samper_426828278-Figura20.png|551px|Square domain with homogeneous Dirichlet conditions and constant source. Solution for a structural mesh of 20×20 four node quadrilaterals]]
+|- style="text-align: center; font-size: 75%;"
+| colspan="1" | '''Figure 20:''' Square domain with homogeneous Dirichlet conditions and constant source. Solution for a structural mesh of <math>20\times 20</math> four node quadrilaterals
+|}
+<div id='img-21'></div>
+{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
+|-
+|[[Image:Draft_Samper_426828278-Figura21.png|551px|Square domain with homogeneous Dirichlet conditions and constant source. Solution for an unstructured mesh of 2 ×20×20 three node triangles]]
+|- style="text-align: center; font-size: 75%;"
+| colspan="1" | '''Figure 21:''' Square domain with homogeneous Dirichlet conditions and constant source. Solution for an unstructured mesh of <math>2 \times 20\times 20</math> three node triangles
+|}
+This problem was also studied by Codina <span id='citeF-11'></span>[[#cite-9|[9,33]]] with four node square elements as another benchmark to study the performance of different shock capturing methods. The FIC method proposed here has yielded accurate solutions for both quadrilateral and triangular elements, which are difficult to obtain for many of the existing stabilization procedures, as reported in <span id='citeF-11'></span>[[#cite-9|[9,33]]].
 ==7 CONCLUDING REMARKS==
@@ Line 859: / Line 1,016: @@
 The authors thank S. Badia, C. Felippa, R. Löhner, R.L. Taylor and O.C. Zienkiewicz for many useful discussions.
-==APPENDIX==
+==APPENDIX A==
-===Computation of the balancing diffusion matrix ̄D in local and global axes===
+===Computation of the balancing diffusion matrix <math>\overline{\textbf D}</math> in local and global axes===
-The balancing diffusion matrix <math display="inline">\bar {D}</math> can be computed by transforming the local matrix <math display="inline">\bar{D}'</math> of Eq.(23) to global axes, or by using directly Eq.(33). We will prove next that both expressions yield the same system of governing differential equations, i.e. we will prove that
+The balancing diffusion matrix <math display="inline">\bar {\textbf D}</math> can be computed by transforming the local matrix <math display="inline">\bar{\textbf D}'</math> of Eq.(23) to global axes, or by using directly Eq.(33). We will prove next that both expressions yield the same system of governing differential equations, i.e. we will prove that
 {| class="formulaSCP" style="width: 100%; text-align: left;"
@@ Line 870: / Line 1,027: @@
 {| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>{\boldsymbol \nabla }^T ( {T}^T \bar{D}' {T}) {\boldsymbol \nabla }\phi = {\boldsymbol \nabla }^T \left({1\over 2} {h} {u}^T\right){\boldsymbol \nabla }\phi  </math>
+| style="text-align: center;" | <math>{\boldsymbol \nabla }^T ( {\textbf T}^T \bar{\textbf D}' {\textbf T}) {\boldsymbol \nabla }\phi = {\boldsymbol \nabla }^T \left({1\over 2} {\textbf h} {\textbf u}^T\right){\boldsymbol \nabla }\phi  </math>
 |}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (42)
+| style="width: 5px;text-align: right;white-space: nowrap;" | (A.1)
 |}
@@ Line 882: / Line 1,039: @@
 {| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>\bar{D}' = {1\over 2} \left[\begin{array}{cc}u_\xi h_\xi & 0\\ 0 & u_\eta h_\eta \end{array}\right]\quad ,\quad {T} =\left[\begin{array}{cc}c_\alpha & s_\alpha \\ - s_\alpha & c_\alpha \end{array}\right] </math>
+| style="text-align: center;" | <math>\bar{\textbf D}' = {1\over 2} \left[\begin{array}{cc}u_\xi h_\xi & 0\\ 0 & u_\eta h_\eta \end{array}\right]\quad ,\quad {\textbf T} =\left[\begin{array}{cc}c_\alpha & s_\alpha \\ - s_\alpha & c_\alpha \end{array}\right] </math>
 |}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (43)
+| style="width: 5px;text-align: right;white-space: nowrap;" | (A.2)
 |}
@@ Line 894: / Line 1,051: @@
 {| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>{h}= \left\{\begin{array}{c}h_x \\h_y\end{array} \right\}= {T}^T \left\{\begin{array}{c}h_\xi \\h_\eta \end{array} \right\}= {T}^T{h}'</math>
+| style="text-align: center;" | <math>{\textbf h}= \left\{\begin{array}{c}h_x \\h_y\end{array} \right\}= {\textbf T}^T \left\{\begin{array}{c}h_\xi \\h_\eta \end{array} \right\}= {\textbf T}^T{\textbf h}'</math>
-| style="width: 5px;text-align: right;white-space: nowrap;" | (46)
+| style="width: 5px;text-align: right;white-space: nowrap;" | (A.3)
 |-
-| style="text-align: center;" | <math> {u}= \left\{\begin{array}{c}u\\ v \end{array} \right\}= {T}^T \left\{\begin{array}{c}u_\xi \\u_\eta \end{array} \right\}= {T}^T{u}' </math>
+| style="text-align: center;" | <math> {\textbf u}= \left\{\begin{array}{c}u\\ v \end{array} \right\}= {\textbf T}^T \left\{\begin{array}{c}u_\xi \\u_\eta \end{array} \right\}= {\textbf T}^T{\textbf u}' </math>
-| style="width: 5px;text-align: right;white-space: nowrap;" | (49)
+| style="width: 5px;text-align: right;white-space: nowrap;" | (A.4)
 |}
 |}
@@ Line 911: / Line 1,068: @@
 {| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>\bar{D} = {T}^T \bar{D}'{T} = {1\over 2}  \left[\begin{array}{ccc}\displaystyle{u_\xi h_\xi c_\alpha ^2 + u_\eta h_\eta s_\alpha ^2  \over \hbox{Sym.}} & \!\!\!\!\Bigg|& \!\!\!\!\displaystyle{(u_\xi h_\xi - u_\eta h_\eta ) s_\alpha c_\alpha \over  u_\xi h_\xi s_\alpha ^2 + u_\eta h_\eta c_\alpha ^2}\end{array}\right] </math>
+| style="text-align: center;" | <math>\bar{\textbf D} = {\textbf T}^T \bar{\textbf D}'{\textbf T} = {1\over 2}  \left[\begin{array}{ccc}\displaystyle{u_\xi h_\xi c_\alpha ^2 + u_\eta h_\eta s_\alpha ^2  \over \hbox{Sym.}} & \!\!\!\!\Bigg|& \!\!\!\!\displaystyle{(u_\xi h_\xi - u_\eta h_\eta ) s_\alpha c_\alpha \over  u_\xi h_\xi s_\alpha ^2 + u_\eta h_\eta c_\alpha ^2}\end{array}\right] </math>
 |}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (50)
+| style="width: 5px;text-align: right;white-space: nowrap;" | (A.5)
 |}
@@ Line 921: / Line 1,078: @@
 {| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>\begin{array}{ll}{\boldsymbol \nabla }^T ({T}^T \bar{D}'{T}) {\boldsymbol \nabla }\phi = & {1\over 2} \left[(u_\xi h_\xi c_\alpha ^2  + u_\eta h_\eta s_\alpha ^2) \phi _x^{''} + (u_\xi h_\xi s_\alpha ^2 + u_\eta h_\eta c_\alpha ^2) \phi _y^{''} +\right.\\ &\left.+  2 (u_\xi h_\xi - u_\eta h_\eta ) s_\alpha c_\alpha \phi _{xy}^{''}\right] \end{array} </math>
+| style="text-align: center;" | <math>\begin{array}{ll}{\boldsymbol \nabla }^T ({\textbf T}^T \bar{\textbf D}'{\textbf T}) {\boldsymbol \nabla }\phi = & {1\over 2} \left[(u_\xi h_\xi c_\alpha ^2  + u_\eta h_\eta s_\alpha ^2) \phi _x^{''} + (u_\xi h_\xi s_\alpha ^2 + u_\eta h_\eta c_\alpha ^2) \phi _y^{''} +\right.\\ &\left.+  2 (u_\xi h_\xi - u_\eta h_\eta ) s_\alpha c_\alpha \phi _{xy}^{''}\right] \end{array} </math>
 |}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (51)
+| style="width: 5px;text-align: right;white-space: nowrap;" | (A.6)
 |}
@@ Line 933: / Line 1,090: @@
 {| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>{1\over 2} {h }{u}^T = {1\over 2} {T}^T ({h}' {u}^{'T}) {T} </math>
+| style="text-align: center;" | <math>{1\over 2} {\textbf h }{\textbf u}^T = {1\over 2} {\textbf T}^T ({\textbf h}' {\textbf u}^{'T}) {\textbf T} </math>
 |}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (52)
+| style="width: 5px;text-align: right;white-space: nowrap;" | (A.7)
 |}
@@ Line 945: / Line 1,102: @@
 {| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>\begin{array}{ll}{\boldsymbol \nabla }^T \left({1\over 2}{h }{u}^T\right){\boldsymbol \nabla } \phi = & {1\over 2} {\boldsymbol \nabla }^T [{T}^T ({h}' {u}^{'T}){T}]  {\boldsymbol \nabla } \phi ]= {1\over 2} \left[(u_\xi h_\xi c_\alpha ^2  + \right.\\ &\left.u_\eta h_\eta s_\alpha ^2) \phi _x^{''} + (u_\xi h_\xi s_\alpha ^2 + u_\eta h_\eta c_\alpha ^2) \phi _y^{''} + 2 (u_\xi h_\xi - u_\eta h_\eta ) s_\alpha c_\alpha \phi _{xy}^{''}\right]+ \\ & \left[s_\alpha c_\alpha (\phi _y^{''}-\phi _x^{''}) + (c_\alpha ^2 - s_\alpha ^2 ) \phi _{xy}^{''}\right](u_\xi h_\xi + u_\eta h_\eta ) \end{array} </math>
+| style="text-align: center;" | <math>\begin{array}{ll}{\boldsymbol \nabla }^T \left({1\over 2}{\textbf h }{\textbf u}^T\right){\boldsymbol \nabla } \phi = & {1\over 2} {\boldsymbol \nabla }^T [{\textbf T}^T ({\textbf h}' {\textbf u}^{'T}){\textbf T}]  {\boldsymbol \nabla } \phi ]= {1\over 2} \left[(u_\xi h_\xi c_\alpha ^2  + \right.\\ &\left.u_\eta h_\eta s_\alpha ^2) \phi _x^{''} + (u_\xi h_\xi s_\alpha ^2 + u_\eta h_\eta c_\alpha ^2) \phi _y^{''} + 2 (u_\xi h_\xi - u_\eta h_\eta ) s_\alpha c_\alpha \phi _{xy}^{''}\right]+ \\ & \left[s_\alpha c_\alpha (\phi _y^{''}-\phi _x^{''}) + (c_\alpha ^2 - s_\alpha ^2 ) \phi _{xy}^{''}\right](u_\xi h_\xi + u_\eta h_\eta ) \end{array} </math>
 |}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (53)
+| style="width: 5px;text-align: right;white-space: nowrap;" | (A.8)
 |}
-Clearly Eqs.(A.6) and (A.8) coincide if the last term in Eq.(A.8) is zero. This happens when <math display="inline">\vec{\xi }</math> and <math display="inline">\vec{\eta }</math> are the principal curvature directions. In order to prove this we will consider the usual case for which <math display="inline">\phi _{xy}^{''} = \phi _{yx}^{''}</math> and both principal curvature directions are orthogonal. Then it exists a transformation matrix <math display="inline">T</math> such that
+Clearly Eqs.(A.6) and (A.8) coincide if the last term in Eq.(A.8) is zero. This happens when <math display="inline">\vec{\xi }</math> and <math display="inline">\vec{\eta }</math> are the principal curvature directions. In order to prove this we will consider the usual case for which <math display="inline">\phi _{xy}^{''} = \phi _{yx}^{''}</math> and both principal curvature directions are orthogonal. Then it exists a transformation matrix <math display="inline">\textbf T</math> such that
 {| class="formulaSCP" style="width: 100%; text-align: left;"
@@ Line 957: / Line 1,114: @@
 {| style="text-align: left; margin:auto;width: 100%;"
 |-
-| style="text-align: center;" | <math>{T}\left[\begin{array}{cc}\phi _x^{''} & \phi _{xy}^{''}\\ \phi _{xy}^{''}&  \phi _y^{''}  \end{array}\right]{T} ^T= \left[\begin{array}{cc}\phi _{\xi }^{''} &0\\ 0 & \phi _{\eta }^{''} \end{array}\right] </math>
+| style="text-align: center;" | <math>{\textbf T}\left[\begin{array}{cc}\phi _x^{''} & \phi _{xy}^{''}\\ \phi _{xy}^{''}&  \phi _y^{''}  \end{array}\right]{\textbf T} ^T= \left[\begin{array}{cc}\phi _{\xi }^{''} &0\\ 0 & \phi _{\eta }^{''} \end{array}\right] </math>
 |}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (54)
+| style="width: 5px;text-align: right;white-space: nowrap;" | (A.9)
 |}
-with <math display="inline">{T}</math> defined as in Eq.(A.2).
+with <math display="inline">{\textbf T}</math> defined as in Eq.(A.2).
 The condition which ensures a diagonal matrix in the r.h.s. of Eq.(A.9) is
@@ Line 973: / Line 1,130: @@
 | style="text-align: center;" | <math>s_\alpha c_\alpha (\phi _y^{''}-\phi _x^{''}) +  ( c_\alpha ^2 - s_\alpha ^2 )\phi _{xy}^{''} =0 </math>
 |}
-| style="width: 5px;text-align: right;white-space: nowrap;" | (55)
+| style="width: 5px;text-align: right;white-space: nowrap;" | (A.10)
 |}
@@ Line 980: / Line 1,137: @@
 Hence the identity of Eq.(A.1) is satisfied.
-===BIBLIOGRAPHY===
+==REFERENCES==
 <div id="cite-1"></div>
@@ Line 992: / Line 1,149: @@
 <div id="cite-4"></div>
-'''[4]'''  Kelly D.W.,  Nakazawa S.,  Zienkiewicz O.C. and  Heinrich  J.C. A  note on upwind and anisotropic balancing dissipation in finite element  approximation to convective diffusion problems. ''Int. J. Num. Meth.  Engng.'', '''15''', pp. 1705&#8211;1711, 1980.
+'''[[#citeF-4|[4]]]'''  Kelly D.W.,  Nakazawa S.,  Zienkiewicz O.C. and  Heinrich  J.C. A  note on upwind and anisotropic balancing dissipation in finite element  approximation to convective diffusion problems. ''Int. J. Num. Meth.  Engng.'', '''15''', pp. 1705&#8211;1711, 1980.
 <div id="cite-5"></div>
@@ Line 998: / Line 1,155: @@
 <div id="cite-6"></div>
-'''[6]''' Hughes T.J.R and Tezduyar, T.E. Finite element method for first order hyperbolic systems with particular emphasis on the compressible Euler equations. ''Comput. Meth. Appl. Mech. Engng.'', '''45''', 217&#8211;284, 1984.
+'''[[#citeF-6|[6]]]''' Hughes T.J.R and Tezduyar, T.E. Finite element method for first order hyperbolic systems with particular emphasis on the compressible Euler equations. ''Comput. Meth. Appl. Mech. Engng.'', '''45''', 217&#8211;284, 1984.
 <div id="cite-7"></div>
-'''[7]''' Hughes T.J.R and Mallet M. 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&#8211;328, 1986a.
+'''[[#citeF-7|[7]]]''' Hughes T.J.R and Mallet M. 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&#8211;328, 1986a.
 <div id="cite-8"></div>
-'''[8]'''  Codina R., Oñate E. and Cervera M. The intrinsic time for the streamline-upwind Petrov-Galerkin formulation using quadratic elements. ''Int. J. Num. Meth. Engrg.'', '''94''', 239&#8211;262, 1992.
+'''[[#citeF-8|[8]]]'''  Codina R., Oñate E. and Cervera M. The intrinsic time for the streamline-upwind Petrov-Galerkin formulation using quadratic elements. ''Int. J. Num. Meth. Engrg.'', '''94''', 239&#8211;262, 1992.
 <div id="cite-9"></div>
-'''[9]'''  Codina R. ''A finite element formulation for the numerical solution of the convection-diffusion equations''. Monograph no. 14, CIMNE, Barcelona 1993.
+'''[[#citeF-9|[9]]]'''  Codina R. ''A finite element formulation for the numerical solution of the convection-diffusion equations''. Monograph no. 14, CIMNE, Barcelona 1993.
 <div id="cite-10"></div>
-'''[10]'''  Donea J. A Taylor-Galerkin method for convective transport  problems. ''Int. J. Num. Meth. Engng.'', '''20''', pp. 101&#8211;119,  1984.
+'''[[#citeF-10|[10]]]'''  Donea J. A Taylor-Galerkin method for convective transport  problems. ''Int. J. Num. Meth. Engng.'', '''20''', pp. 101&#8211;119,  1984.
 <div id="cite-11"></div>
@@ Line 1,019: / Line 1,176: @@
 <div id="cite-13"></div>
-'''[13]'''   Goldschmit M.B. and Dvorkin E.N. On the solution of the  steady-state convection-diffusion equation using quadratic elements. A  generalized Galerkin technique also reliable with distorted meshes. '' Engng. Comput.'', '''11''', 6, 565&#8211;579, 1994.
+'''[[#citeF-13|[13]]]'''   Goldschmit M.B. and Dvorkin E.N. On the solution of the  steady-state convection-diffusion equation using quadratic elements. A  generalized Galerkin technique also reliable with distorted meshes. '' Engng. Comput.'', '''11''', 6, 565&#8211;579, 1994.
 <div id="cite-14"></div>
-'''[14]''' Hughes T.J.R, Franca L.P. and Hulbert G.M. 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&#8211;189, 1989.
+'''[[#citeF-14|[14]]]''' Hughes T.J.R, Franca L.P. and Hulbert G.M. 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&#8211;189, 1989.
 <div id="cite-15"></div>
-'''[15]'''   Franca L.P.,  Frey S.L. and  Hughes T.J.R. Stabilized finite element methods: I. Application to the  advective-diffusive model. ''Comput. Meth. Appl. Mech. Engn'', Vol. '''95''', pp. 253&#8211;276, 1992.
+'''[[#citeF-15|[15]]]'''   Franca L.P.,  Frey S.L. and  Hughes T.J.R. Stabilized finite element methods: I. Application to the  advective-diffusive model. ''Comput. Meth. Appl. Mech. Engn'', Vol. '''95''', pp. 253&#8211;276, 1992.
 <div id="cite-16"></div>
-'''[16]''' Franca  L.P. and Dutra do Carmo E.G. The Galerkin Gradient Least Square method. ''Comput. Methods Appl. Mech. Engrg.'', '''74''', 41&#8211;54, 1989.
+'''[[#citeF-16|[16]]]''' Franca  L.P. and Dutra do Carmo E.G. The Galerkin Gradient Least Square method. ''Comput. Methods Appl. Mech. Engrg.'', '''74''', 41&#8211;54, 1989.
 <div id="cite-17"></div>
@@ Line 1,034: / Line 1,191: @@
 <div id="cite-18"></div>
-'''[18]'''  Löhner R., Morgan K. and  Zienkiewicz O.C. The solution of non-linear hyperbolic equation systems by the finite element method. ''Int. J. Num. Meth. in Fluids'', '''4''', 1043, 1984.
+'''[[#citeF-18|[18]]]'''  Löhner R., Morgan K. and  Zienkiewicz O.C. The solution of non-linear hyperbolic equation systems by the finite element method. ''Int. J. Num. Meth. in Fluids'', '''4''', 1043, 1984.
 <div id="cite-19"></div>
@@ Line 1,040: / Line 1,197: @@
 <div id="cite-20"></div>
-'''[20]''' Hughes T.J.R. Multiscale phenomena: Green functions,  subgrid scale models, bubbles and the origins of stabilized methods.  ''Comput. Methods Appl. Mech. Engrg'' 1995, '''127''', 387&#8211;401, 1995.
+'''[[#citeF-20|[20]]]''' Hughes T.J.R. Multiscale phenomena: Green functions,  subgrid scale models, bubbles and the origins of stabilized methods.  ''Comput. Methods Appl. Mech. Engrg'' 1995, '''127''', 387&#8211;401, 1995.
 <div id="cite-21"></div>
-'''[21]''' Brezzi F., Franca L.P., Hughes T.J.R. and  Russo A. <math>b=\int  g</math>. ''Comput. Methods Appl. Mech. Engrg.'', '''145''', 329&#8211;339, 1997.
+'''[[#citeF-21|[21]]]''' Brezzi F., Franca L.P., Hughes T.J.R. and  Russo A. <math>b=\int  g</math>. ''Comput. Methods Appl. Mech. Engrg.'', '''145''', 329&#8211;339, 1997.
 <div id="cite-22"></div>
@@ Line 1,049: / Line 1,206: @@
 <div id="cite-23"></div>
-'''[23]'''  Hauke G. A simple subgrid scale stabilized method for the advection-diffusion-reaction equation. ''Comput. Methods Appl. Mech. Engrg.'', '''191''', 2925&#8211;2948, 2002.
+'''[[#citeF-23|[23]]]'''  Hauke G. A simple subgrid scale stabilized method for the advection-diffusion-reaction equation. ''Comput. Methods Appl. Mech. Engrg.'', '''191''', 2925&#8211;2948, 2002.
 <div id="cite-24"></div>
@@ Line 1,058: / Line 1,215: @@
 <div id="cite-26"></div>
-'''[26]'''  Ilinca F., Hétu J.F. and Pelletier D. On stabilized finite element formulation for incompressible advective-diffusive transport and fluid flow problems. ''Comput. Methods  Appl. Mech. Engrg.'', '''188''', 235&#8211;257, 2000.
+'''[[#citeF-26|[26]]]'''  Ilinca F., Hétu J.F. and Pelletier D. On stabilized finite element formulation for incompressible advective-diffusive transport and fluid flow problems. ''Comput. Methods  Appl. Mech. Engrg.'', '''188''', 235&#8211;257, 2000.
 <div id="cite-27"></div>
-'''[27]'''  Tezduyar T.E. and Osawa Y. Finite element stabilization parameters computed from element matrices and vectors. ''Comput. Methods Appl. Mech. Engrg.'', '''190''', 411&#8211;430, 2000.
+'''[[#citeF-27|[27]]]'''  Tezduyar T.E. and Osawa Y. Finite element stabilization parameters computed from element matrices and vectors. ''Comput. Methods Appl. Mech. Engrg.'', '''190''', 411&#8211;430, 2000.
 <div id="cite-28"></div>
@@ Line 1,067: / Line 1,224: @@
 <div id="cite-29"></div>
-'''[29]'''  Harari I.,  Franca L.P. and  Oliveira S.P. Streamline design of stability parameters for advection-diffusion problems. ''J. Comput. Physics'', '''171''', 115&#8211;131, 2001.
+'''[[#citeF-29|[29]]]'''  Harari I.,  Franca L.P. and  Oliveira S.P. Streamline design of stability parameters for advection-diffusion problems. ''J. Comput. Physics'', '''171''', 115&#8211;131, 2001.
 <div id="cite-30"></div>
-'''[30]'''  Nesliturk A. and Harari I. The nearly-optimal Petrov-Galerkin method for convection-diffusion problems. ''Comput Methods  Appl. Mech. Engrg.'', '''192''', 2501&#8211;2519, 2003.
+'''[[#citeF-30|[30]]]'''  Nesliturk A. and Harari I. The nearly-optimal Petrov-Galerkin method for convection-diffusion problems. ''Comput Methods  Appl. Mech. Engrg.'', '''192''', 2501&#8211;2519, 2003.
 <div id="cite-31"></div>
-'''[31]'''  Tezduyar T.E. Computation of moving boundaries and stabilization parameters. ''Int. J. Num. Meth. Fluids'', '''43''', 555&#8211;575, 2003.
+'''[[#citeF-31|[31]]]'''  Tezduyar T.E. Computation of moving boundaries and stabilization parameters. ''Int. J. Num. Meth. Fluids'', '''43''', 555&#8211;575, 2003.
 <div id="cite-32"></div>
-'''[32]'''  Hughes T.J.R. and Mallet M. A new finite element formulations for computational fluid dynamics: IV. A discontinuity capturing operator for multidimensional advective-diffusive system. ''Comput. Methods Appl. Mech. Engrg.'', '''58''', 329&#8211;336, 1986b.
+'''[[#citeF-32|[32]]]'''  Hughes T.J.R. and Mallet M. A new finite element formulations for computational fluid dynamics: IV. A discontinuity capturing operator for multidimensional advective-diffusive system. ''Comput. Methods Appl. Mech. Engrg.'', '''58''', 329&#8211;336, 1986b.
 <div id="cite-33"></div>
-'''[33]'''  Codina R. A discontinuity-capturing crosswind dissipation for the finite element solution of the convection-diffusion equation. ''Comput. Methods Appl. Mech. Engrg.'', '''110''', 325&#8211;342, 1993.
+'''[[#citeF-33|[33]]]'''  Codina R. A discontinuity-capturing crosswind dissipation for the finite element solution of the convection-diffusion equation. ''Comput. Methods Appl. Mech. Engrg.'', '''110''', 325&#8211;342, 1993.
 <div id="cite-34"></div>
-'''[34]'''  Hughes T.J.R.,  Mallet M. and Mizukami A. A new finite element formulation for comptutational fluid dynamics: II Beyond SUPG. ''Comput. Methods Appl. Mech. Engrg.'', '''54''', 341&#8211;355, 1986.
+'''[[#citeF-34|[34]]]'''  Hughes T.J.R.,  Mallet M. and Mizukami A. A new finite element formulation for comptutational fluid dynamics: II Beyond SUPG. ''Comput. Methods Appl. Mech. Engrg.'', '''54''', 341&#8211;355, 1986.
 <div id="cite-35"></div>
-'''[35]'''  Galeo A.C. and  Dutra do Carmo E.G. A consistent approximate upwind Petrov-Galerkin method for convection-dominated problems. ''Comput. Methods Appl. Mech. Engrg.'', '''68''', 83&#8211;95, 1988.
+'''[[#citeF-35|[35]]]'''  Galeo A.C. and  Dutra do Carmo E.G. A consistent approximate upwind Petrov-Galerkin method for convection-dominated problems. ''Comput. Methods Appl. Mech. Engrg.'', '''68''', 83&#8211;95, 1988.
 <div id="cite-36"></div>
-'''[36]'''  Tezduyar T.E. and  Park Y.J. Discontinuity-capturing finite element formulations for nonlinear convection-diffusion-reaction equations. ''Comput. Methods Appl. Mech. Engrg.'', '''59''', 307&#8211;325, 1986
+'''[[#citeF-36|[36]]]'''  Tezduyar T.E. and  Park Y.J. Discontinuity-capturing finite element formulations for nonlinear convection-diffusion-reaction equations. ''Comput. Methods Appl. Mech. Engrg.'', '''59''', 307&#8211;325, 1986
 <div id="cite-37"></div>
@@ Line 1,094: / Line 1,251: @@
 <div id="cite-38"></div>
-'''[38]'''  Oñate E. Possibilities of finite calculus in computational mechanics.  ''Int. J. Num. Meth. Engng.'' 2004; '''60''':255&#8211;281.
+'''[[#citeF-38|[38]]]'''  Oñate E. Possibilities of finite calculus in computational mechanics.  ''Int. J. Num. Meth. Engng.'' 2004; '''60''':255&#8211;281.
 <div id="cite-39"></div>
@@ Line 1,106: / Line 1,263: @@
 <div id="cite-42"></div>
-'''[42]'''  Felippa C.A. and  Oñate E. Nodally exact Ritz discretization of 1D diffusion-absorption and Helmholtz equations by variational FIC and modified equation methods. Research Report No. PI 237, CIMNE, Barcelona 2004. Submitted to ''Computational Mechanics''.
+'''[[#citeF-42|[42]]]'''  Felippa C.A. and  Oñate E. Nodally exact Ritz discretization of 1D diffusion-absorption and Helmholtz equations by variational FIC and modified equation methods. Research Report No. PI 237, CIMNE, Barcelona 2004. Submitted to ''Computational Mechanics''.
 <div id="cite-43"></div>
-'''[43]'''  Oñate E., Miquel J. and Hauke G. Stabilized formulation for the 1D advection-diffusion-absorption equation using finite calculus and the finite element method. Submitted to ''Comput. Meth. Appl. Mech. Engng.'', Dec. 2004.
+'''[[#citeF-43|[43]]]'''  Oñate E., Miquel J. and Hauke G. Stabilized formulation for the 1D advection-diffusion-absorption equation using finite calculus and the finite element method. Submitted to ''Comput. Meth. Appl. Mech. Engng.'', Dec. 2004.
 <div id="cite-44"></div>
-'''[44]'''  Oñate E. A stabilized finite element method for incompressible viscous flows using a finite increment calculus formulation. ''Comput. Methods Appl. Mech. Engrg.'' 2000; '''182:1&#8211;2''':355&#8211;370.
+'''[[#citeF-44|[44]]]'''  Oñate E. A stabilized finite element method for incompressible viscous flows using a finite increment calculus formulation. ''Comput. Methods Appl. Mech. Engrg.'' 2000; '''182:1&#8211;2''':355&#8211;370.
 <div id="cite-45"></div>
-'''[45]'''  Oñate E.,  García J. A finite element method for fluid-structure interaction with surface waves using a finite calculus formulation. in ''Comput. Methods Appl. Mech. Engrg.'', '''191:6-7''', 635-660, 2001.
+'''[[#citeF-45|[45]]]'''  Oñate E.,  García J. A finite element method for fluid-structure interaction with surface waves using a finite calculus formulation. in ''Comput. Methods Appl. Mech. Engrg.'', '''191:6-7''', 635-660, 2001.
 <div id="cite-46"></div>
-'''[46]'''  García J. and Oñate E. An unstructured finite element solver for ship hydrodynamic problems. ''Journal of Appl. Mechanics'', '''70''', 18&#8211;26, 2003.
+'''[[#citeF-46|[46]]]'''  García J. and Oñate E. An unstructured finite element solver for ship hydrodynamic problems. ''Journal of Appl. Mechanics'', '''70''', 18&#8211;26, 2003.
 <div id="cite-47"></div>
@@ Line 1,124: / Line 1,281: @@
 <div id="cite-48"></div>
-'''[48]'''  Oñate E., Taylor R.L., Zienkiewicz O.C. and  Rojek J. A residual correction method based on finite calculus. ''Engineering Computations'', Vol. '''20''', No. '''5/6''', 629&#8211;658, 2003.
+'''[[#citeF-48|[48]]]'''  Oñate E., Taylor R.L., Zienkiewicz O.C. and  Rojek J. A residual correction method based on finite calculus. ''Engineering Computations'', Vol. '''20''', No. '''5/6''', 629&#8211;658, 2003.
 <div id="cite-49"></div>
@@ Line 1,134: / Line 1,291: @@
 <div id="cite-51"></div>
 '''[[#citeF-51|[51]]]'''  Wiberg NW, Abdulwahab F,  Li XD. (1997). Error estimation  and adaptive procedures based on superconvergent patch recovery. ''Archives Comput. Meth. Engng.'' 1997; '''4:3''':203&#8211;242.
-<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_426828278-Figure1.png|351px|Global axes (x,y) and principal curvature axes (ξ,η)]]
-|- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 1:''' Global axes (<math>x,y</math>) and principal curvature axes (<math>\xi ,\eta </math>)
-|}
-<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_426828278-Figure2.png|351px|Square domain with uniform Dirichlet conditions, upward diagonal velocity and zero source. SUPG and FIC solutions obtained with a structured mesh of 10×10 linear four node square elements]]
-|- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 2:''' Square domain with uniform Dirichlet conditions, upward diagonal velocity and zero source. SUPG and FIC solutions obtained with a structured mesh of <math>10\times 10</math> linear four node square elements
-|}
-<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_426828278-Figure3.png|351px|Square domain with uniform Dirichlet conditions, upward diagonal velocity and zero source. SUPG and FIC solutions obtained with a structured mesh of 2×10×10 linear three node triangles]]
-|- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 3:''' Square domain with uniform Dirichlet conditions, upward diagonal velocity and zero source. SUPG and FIC solutions obtained with a structured mesh of <math>2\times 10\times 10</math> linear three node triangles
-|}
-<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_426828278-Figure4.png|351px|Square domain with uniform Dirichlet conditions, upward diagonal velocity and zero source. Convergence of the FIC solution for square and triangular elements]]
-|- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 4:''' Square domain with uniform Dirichlet conditions, upward diagonal velocity and zero source. Convergence of the FIC solution for square and triangular elements
-|}
-<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_426828278-Figure5.png|351px|Square domain with uniform Dirichlet conditions, upward diagonal velocity and zero source. SUPG and FIC solutions obtained with a structured mesh of 10×20 linear four node rectangular elements. Geometrical aspect ratio 2:1]]
-|- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 5:''' Square domain with uniform Dirichlet conditions, upward diagonal velocity and zero source. SUPG and FIC solutions obtained with a structured mesh of <math>10\times 20</math> linear four node rectangular elements. Geometrical aspect ratio 2:1
-|}
-<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_426828278-Figure6.png|351px|Square domain with uniform Dirichlet conditions, upward diagonal velocity and zero source. SUPG and FIC solutions obtained with a structured mesh of 2×10×20 linear triangles. Geometrical aspect ratio 2:1]]
-|- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 6:''' Square domain with uniform Dirichlet conditions, upward diagonal velocity and zero source. SUPG and FIC solutions obtained with a structured mesh of <math>2\times 10\times 20</math> linear triangles. Geometrical aspect ratio 2:1
-|}
-<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_426828278-Figure7.png|351px|Square domain with non uniform Dirichlet conditions, downwards diagonal velocity and zero source. SUPG and FIC solutions obtained with a structured mesh of 20×20 linear four node square elements]]
-|- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 7:''' Square domain with non uniform Dirichlet conditions, downwards diagonal velocity and zero source. SUPG and FIC solutions obtained with a structured mesh of <math>20\times 20</math> linear four node square elements
-|}
-<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_426828278-Figure8.png|351px|Square domain with non uniform Dirichlet conditions, downwards diagonal velocity  and zero source. SUPG and FIC solutions obtained with a structured mesh of 2×20×20 three node triangles]]
-|- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 8:''' Square domain with non uniform Dirichlet conditions, downwards diagonal velocity  and zero source. SUPG and FIC solutions obtained with a structured mesh of <math>2\times 20\times 20</math> three node triangles
-|}
-<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_426828278-Figure9.png|351px|Square domain with uniform Dirichlet conditions, upward diagonal velocity and zero source. Convergence of the FIC solution for square and triangular elements]]
-|- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 9:''' Square domain with uniform Dirichlet conditions, upward diagonal velocity and zero source. Convergence of the FIC solution for square and triangular elements
-|}
-<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_426828278-Figure10.png|351px|Square domain with non uniform Dirichlet conditions, downwards diagonal velocity and zero source. SUPG and FIC solutions obtained with an unstructured mesh of 432 linear four node quadrilateral elements]]
-|- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 10:''' Square domain with non uniform Dirichlet conditions, downwards diagonal velocity and zero source. SUPG and FIC solutions obtained with an unstructured mesh of 432 linear four node quadrilateral elements
-|}
-<div id='img-11'></div>
-{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
-|-
-|[[Image:Draft_Samper_426828278-Figure11.png|351px|Square domain with non uniform Dirichlet conditions, downwards diagonal velocity  and zero source. SUPG and FIC solutions obtained with an unstructured mesh of 780 three node triangles]]
-|- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 11:''' Square domain with non uniform Dirichlet conditions, downwards diagonal velocity  and zero source. SUPG and FIC solutions obtained with an unstructured mesh of 780 three node triangles
-|}
-<div id='img-12'></div>
-{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
-|-
-|[[Image:Draft_Samper_426828278-Figure12.png|351px|Square domain with non uniform Dirichlet conditions, downwards diagonal velocity and zero source. SUPG and FIC solutions obtained with a refined unstructured mesh of 1826 linear four node quadrilaterals]]
-|- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 12:''' Square domain with non uniform Dirichlet conditions, downwards diagonal velocity and zero source. SUPG and FIC solutions obtained with a refined unstructured mesh of 1826 linear four node quadrilaterals
-|}
-<div id='img-13'></div>
-{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
-|-
-|[[Image:Draft_Samper_426828278-Figure13.png|351px|Square domain with non uniform Dirichlet conditions, downwards diagonal velocity and zero source. SUPG and FIC solutions obtained with a refined unstructured mesh of 1898 linear triangles]]
-|- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 13:''' Square domain with non uniform Dirichlet conditions, downwards diagonal velocity and zero source. SUPG and FIC solutions obtained with a refined unstructured mesh of 1898 linear triangles
-|}
-<div id='img-14'></div>
-{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
-|-
-|[[Image:Draft_Samper_426828278-Figure14.png|351px|Square domain with non uniform Dirichlet conditions, upwards diagonal velocity  and zero source. SUPG and FIC solutions obtained with a structured mesh of 10×20 linear four node rectangular elements. Geometrical aspect ratio 2:1]]
-|- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 14:''' Square domain with non uniform Dirichlet conditions, upwards diagonal velocity  and zero source. SUPG and FIC solutions obtained with a structured mesh of <math>10\times 20</math> linear four node rectangular elements. Geometrical aspect ratio 2:1
-|}
-<div id='img-15'></div>
-{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
-|-
-|[[Image:Draft_Samper_426828278-Figure15.png|351px|Square domain with non uniform Dirichlet conditions, upwards diagonal velocity  and zero source. SUPG and FIC solutions obtained with a structured mesh of 2×10×20 three node triangles. Geometrical aspect ratio 2:1]]
-|- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 15:''' Square domain with non uniform Dirichlet conditions, upwards diagonal velocity  and zero source. SUPG and FIC solutions obtained with a structured mesh of <math>2\times 10\times 20</math> three node triangles. Geometrical aspect ratio 2:1
-|}
-<div id='img-16'></div>
-{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
-|-
-|[[Image:Draft_Samper_426828278-Figure16.png|351px|Rectangular domain with Neumann and  non uniform Dirichlet conditions, rotational velocity field and zero source. SUPG and FIC solutions obtained with a structured mesh of 40×20 four node square elements]]
-|- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 16:''' Rectangular domain with Neumann and  non uniform Dirichlet conditions, rotational velocity field and zero source. SUPG and FIC solutions obtained with a structured mesh of <math>40\times 20</math> four node square elements
-|}
-<div id='img-17'></div>
-{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
-|-
-|[[Image:Draft_Samper_426828278-Figure17.png|351px|Rectangular domain with Neumann and  non uniform Dirichlet conditions, rotational velocity field and zero source. SUPG and FIC solutions obtained with  a structured mesh of 2×40×20   three node triangles]]
-|- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 17:''' Rectangular domain with Neumann and  non uniform Dirichlet conditions, rotational velocity field and zero source. SUPG and FIC solutions obtained with  a structured mesh of <math>2\times 40\times 20</math>   three node triangles
-|}
-<div id='img-18'></div>
-{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
-|-
-|[[Image:Draft_Samper_426828278-Figure18.png|351px|]]
-|[[Image:Draft_Samper_426828278-Figure18b.png|351px|Rectangular domain with Neumann and  non uniform Dirichlet conditions, rotational velocity field and zero source. SUPG and FIC solutions obtained with two unstructured meshes of four node quadrilaterals. (a) Coarse mesh of 892 elements. (b) Refined mesh of 2054 elements]]
-|- style="text-align: center; font-size: 75%;"
-| colspan="2" | '''Figure 18:''' Rectangular domain with Neumann and  non uniform Dirichlet conditions, rotational velocity field and zero source. SUPG and FIC solutions obtained with two unstructured meshes of four node quadrilaterals. (a) Coarse mesh of 892 elements. (b) Refined mesh of 2054 elements
-|}
-<div id='img-19'></div>
-{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
-|-
-|[[Image:Draft_Samper_426828278-Figure19.png|351px|Rectangular domain with Neumann and  non uniform Dirichlet conditions, rotational velocity field and zero source. SUPG and FIC solutions obtained with  an unstructured mesh of 1554 three node triangles]]
-|- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 19:''' Rectangular domain with Neumann and  non uniform Dirichlet conditions, rotational velocity field and zero source. SUPG and FIC solutions obtained with  an unstructured mesh of 1554 three node triangles
-|}
-<div id='img-20'></div>
-{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
-|-
-|[[Image:Draft_Samper_426828278-Figura20.png|351px|Square domain with homogeneous Dirichlet conditions and constant source. Solution for a structural mesh of 20×20 four node quadrilaterals]]
-|- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 20:''' Square domain with homogeneous Dirichlet conditions and constant source. Solution for a structural mesh of <math>20\times 20</math> four node quadrilaterals
-|}
-<div id='img-21'></div>
-{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
-|-
-|[[Image:Draft_Samper_426828278-Figura21.png|351px|Square domain with homogeneous Dirichlet conditions and constant source. Solution for an unstructured mesh of 2 ×20×20 three node triangles]]
-|- style="text-align: center; font-size: 75%;"
-| colspan="1" | '''Figure 21:''' Square domain with homogeneous Dirichlet conditions and constant source. Solution for an unstructured mesh of <math>2 \times 20\times 20</math> three node triangles
-|}

Latest revision as of 11:15, 12 November 2018

Abstract

1 INTRODUCTION

2 FIC GOVERNING EQUATIONS FOR STEADY-STATE CONVECTION-DIFFUSIVE PROBLEMS

2.1 Finite element discretization

2.2 Equivalence with the SUPG method

3 COMPUTATION OF THE CHARACTERISTIC LENGTH FOR THE 1D CONVECTION-DIFFUSION EQUATION

4 COMPUTATION OF THE CHARACTERISTIC LENGTH VECTOR

4.1 Computation of the principal curvature axes for linear elements

4.2 General iterative scheme

5 EQUIVALENCE WITH SUPG AND SHOCK-CAPTURING TECHNIQUES

6 EXAMPLES

6.1 Square domain with uniform Dirichlet conditions, upward diagonal velocity and zero source

6.2 Square domain with non uniform Dirichlet conditions, downwards diagonal velocity and zero source

6.3 Rectangular domain with Neumann and non uniform Dirichlet conditions, rotational velocity field and zero source

6.4 Square domain with homogeneous Dirichlet condition and constant source

7 CONCLUDING REMARKS

ACKNOWLEDGEMENTS

APPENDIX A

Computation of the balancing diffusion matrix D ¯ {\displaystyle {\overline {\textbf {D}}}} in local and global axesComputation of the balancing diffusion matrix D ¯ {\displaystyle {\overline {\textbf {D}}}} in local and global axes

REFERENCES

Document information

Document Score

Share this document

Computation of the balancing diffusion matrix ${\overline {\textbf {D}}}$ in local and global axesComputation of the balancing diffusion matrix D ¯ {\displaystyle {\overline {\textbf {D}}}} in local and global axes