FIC-FEM formulation for the multidimensional transient advection-diffusion-absorption equation (preprint)

<!-- metadata commented in wiki content
==FIC-FEM formulation for the multidimensional transient advection-diffusion-absorption equation==

'''Albert Puigferrat<math>^{1}</math>,  Ignasi de-Pouplana<math>^{1,2}</math> and Eugenio Oñate<math>^{1,2}</math>'''

{|
|-
|<math>^{1}</math> Centre Internacional de Metodes Numerics a l'Enginyeria (CIMNE), Barcelona, Spain
|-
| <math>^{2}</math> Universitat Politecnica de Catalunya (UPC), Barcelona, Spain
|}
-->
==Abstract==

In this paper we present  a  stabilized FIC-FEM formulation for the multidimensional transient advection-diffusion-absorption equation. The starting point is the non-local form of the governing equations for the multidimensional transient advection-diffusion-absorption problems obtained via the Finite Increment Calculus (FIC) procedure. The FIC governing equations have a residual form that introduces a characteristic length vector that depends on streamline, absorption and shock capturing stabilization parameters, as well as on a characteristic element size  that ensures a stabilized numerical solution using a standard Galerkin FEM. The value of the stabilization parameters is obtained as an extension of the steady-state form. The accuracy of the FIC-FEM formulation is verified in the solution of several transient advection-diffusion-absorption problems using regular meshes of 3-noded triangles and 4-noded quadrilaterals.

'''Keywords''': Advection-diffusion-absorption, transient solution, multidimensional, Finite element method, Finite increment calculus, FIC

==1 Introduction==

The numerical solution of the advection-diffusion-absorption equation by the Galerkin FEM can exhibit oscillations due to a number of reasons. For the stationary convection-dominated case global oscillations in the analysis can be found, as well as local Gibbs oscillations along the characteristic layers for 2D and 3D problems. Gibbs oscillations near the Dirichlet boundaries and in regions where the source term has an irregular distribution can also appear for stationary absorption-dominated problems. For transient problems dispersive oscillations can appear when the initial solution and/or the distributed source term are non regular <span id='citeF-12'></span><span id='citeF-68'></span>[[#cite-12|[12,68]]].

A number of numerical procedures has been derived in recent decades aiming to overcome the above limitations of Galerkin FEM for solving  advection-diffusion-absorption problems. The goal in all cases has been to obtain “stabilized” numerical solution that are free of spurious oscillations and have a physical meaning. Among the many procedures we can mention the family of streamline-upward Petrov-Galerkin (SUPG) methods. The original SUPG methods for advection-diffusion problems <span id='citeF-5'></span><span id='citeF-25'></span><span id='citeF-26'></span><span id='citeF-33'></span><span id='citeF-35'></span>[[#cite-5|[5,25,26,33,35]]] where enhanced with the addition of non linear shock-capturing terms aiming to control the Gibbs oscillations across internal/boundary layers for advection-diffusion problems <span id='citeF-7'></span><span id='citeF-10'></span><span id='citeF-13'></span><span id='citeF-14'></span><span id='citeF-20'></span><span id='citeF-32'></span><span id='citeF-34'></span><span id='citeF-36'></span><span id='citeF-51'></span>[[#cite-7|[7,10,13,14,20,32,34,36,51]]]. Other families of stabilized methods for this type of problems were based on Galerkin Least Squares (GLS) <span id='citeF-27'></span>[[#cite-27|[27]]], Variational Multiscale  (VMS) <span id='citeF-28'></span>[[#cite-28|[28]]] and Characteristic split procedure <span id='citeF-66'></span><span id='citeF-68'></span>[[#cite-66|[66,68]]]. Many of these methods are reviewed in <span id='citeF-31'></span>[[#cite-31|[31]]].

Control over the dispersive oscillations for  transient advection-diffusion problems via SUPG methods and space-time FEM were reported in <span id='citeF-30'></span>[[#cite-30|[30]]] and <span id='citeF-65'></span>[[#cite-65|[65]]], respectively.

Extensions of the SUPG method for advection-diffusion-reaction problems were reported in <span id='citeF-6'></span><span id='citeF-29'></span><span id='citeF-64'></span>[[#cite-6|[6,29,64]]]. Other stabilization techniques for these problems were based in the GLS method <span id='citeF-17'></span><span id='citeF-21'></span>[[#cite-17|[17,21]]], the addition of internal bubbles <span id='citeF-1'></span><span id='citeF-2'></span><span id='citeF-3'></span><span id='citeF-4'></span><span id='citeF-18'></span><span id='citeF-19'></span>[[#cite-1|[1,2,3,4,18,19]]] and on variations of VMS procedure <span id='citeF-9'></span><span id='citeF-22'></span><span id='citeF-23'></span><span id='citeF-24'></span><span id='citeF-63'></span>[[#cite-9|[9,22,23,24,63]]]. A review of some of these methods is presented in <span id='citeF-8'></span>[[#cite-8|[8]]].

Oñate ''et al.'' <span id='citeF-52'></span><span id='citeF-55'></span>[[#cite-52|[52,55]]] derived a non linear stabilization procedure for the steady state advection-diffusion-absorption problem using a single stabilization parameter via a Finite Increment Calculus (FIC) approach. Oñate and Felippa <span id='citeF-15'></span>[[#cite-15|[15]]] used the variational FIC method for obtaining exact nodal solutions for 1D diffusion-reaction problems (including Helmholtz problems) using also a single linear stabilization parameter.

The homogeneous steady advection-diffusion-reaction equation has typically two fundamental solutions. This has lead the way to the derivation of linear stabilized methods that provide nodally exact solutions in 1D using ''two stabilization parameters'' via internal bubbles <span id='citeF-4'></span>[[#cite-4|[4]]], SUPG <span id='citeF-21'></span>[[#cite-21|[21]]] and VMS <span id='citeF-24'></span>[[#cite-24|[24]]] procedures.

A particular class of two-parameter stabilized models for 1D and multidimensional steady and transient advection-diffusion production problems using the FIC approach and the FEM vas presented by Nadukandi, Oñate and García-Espinosa <span id='citeF-37'></span><span id='citeF-38'></span>[[#cite-37|[37,38]]].

Starting from a different perspective, but still within the FIC framework, Oñate, Miquel and Nadukandi <span id='citeF-61'></span>[[#cite-61|[61]]] presented  a FIC-FEM formulation for the 1D steady-state and transient advection-diffusion-absorption equations  using two linear stabilization parameters. The FIC-FEM formulation yielded exact nodal solution for 1D steady-state  problems using regular meshes of 2-noded linear elements. Very accurate results were also obtained in <span id='citeF-61'></span>[[#cite-61|[61]]] for 1D stationary problems solved with irregular meshes, as well as for 1D transient problems. The FIC-FEM formulation presented in <span id='citeF-61'></span>[[#cite-61|[61]]] was extended by the same authors in  <span id='citeF-62'></span>[[#cite-62|[62]]] to multidimensional steady-state advection-diffusion-absorption problems. Good results were obtained for the problems analyzed using linear 3-noded triangles and bilinear 4-noded quadrilaterals.

In this work we extend to the multidimensional transient case the 1D FIC-FEM formulation presented in <span id='citeF-62'></span>[[#cite-62|[62]]]. The lay-out of the paper is the following. In the next section we formulate the FIC form of the equations governing multidimensional transient convection-diffusion-absorption problems. The finite element discretization is then presented. The stabilization parameters are obtained as an extension of the expression for the stationary case <span id='citeF-61'></span>[[#cite-61|[61]]]. The  accuracy of the multidimensional transient FIC-FEM formulation is verified in the solution of a number of 2D transient  advection-diffusion-absorption problems using uniform meshes of 3-noded triangles and 4-noded quadrilateral elements. Accurate solutions are obtained in all cases.

==2 The multidimensional transient advection-diffusion-reaction problem==

===2.1 Governing equations===

==='''Transport balance'''===

The transport balance equation in a domain of area/volume <math display="inline">\Omega </math> can be expressed as

<span id="eq-1.a"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>r_t=0\quad \hbox{in } \Omega  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (1.a)
|}

with

<span id="eq-1.b"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>r_t:=\rho c \left(\frac{\partial \phi }{\partial t}+\mathbf{v}^T {\boldsymbol \nabla } {\phi }\right)- {\boldsymbol \nabla }^T {\boldsymbol D}{\boldsymbol \nabla }\phi + s\phi -Q  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (1.b)
|}

For 3D problems,

<span id="eq-2"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{v} = [ v_1 ~,~ v_2~,~ v_3]T \quad , \quad  {\boldsymbol D} = \left[\begin{array}{ccc}k_1 & 0&0 \\                       0 & k_2&0\\                       0 &0 & k_3                     \end{array}  \right]\quad , \quad {\boldsymbol \nabla } = \left[\frac{\partial }{\partial x_1} ~,~\frac{\partial }{\partial x_2}~,~\frac{\partial }{\partial x_3}\right]^T                      </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
|}

In Eqs.([[#eq-1.a|1.a]])&#8211;([[#eq-1.b|1.b]]) <math display="inline">\phi </math> is the transported variable (i.e., the temperature in a heat transfer problem or the concentration in a mass transfer problem), <math display="inline">v_i</math> is the <math display="inline">i</math>th component of the velocity vector <math display="inline">{\mathbf{}v}</math>; <math display="inline">\rho ,c</math> and <math display="inline">k_i</math> are the  density, the specific flux parameter and the conductivity of the material along the <math display="inline">i</math>th global direction, respectively and <math display="inline">s</math> is the reaction parameter. In the following, and unless otherwise specified, we will assume that the velocity field is solenoidal and that the problem parameters (<math display="inline">\rho ,  c ,k,s</math>) are constant over the analysis domain <math display="inline">\Omega </math>.

==='''Boundary conditions'''===

<span id="eq-3"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\phi - \phi ^p =0 \quad \hbox{on }\Gamma _\phi  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
|}

<span id="eq-4"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>r_{\scriptscriptstyle \Gamma } =0\quad \hbox{on } \Gamma _q </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
|}

with

<span id="eq-5"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>r_\Gamma := - q_n + q^p_n   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
|}

where

<span id="eq-6"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>q_n = {\boldsymbol q}^T {\boldsymbol n} \quad ,\quad {\boldsymbol q} =  - {\boldsymbol D}  {\boldsymbol \nabla }\phi  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
|}

In Eqs.([[#eq-3|3]])&#8211;[[#eq-6|6]] <math display="inline">\phi ^p</math> and <math display="inline">q^p</math> are the prescribed values of the transported variable and the outgoing diffusive flux at the Dirichlet and Neumann boundaries <math display="inline">\Gamma _\phi </math> and <math display="inline">\Gamma _q</math>, respectively, with <math display="inline">\Gamma _\phi \cup \Gamma _q = \Gamma </math>, <math display="inline">\Gamma </math> being the total boundary of the domain and '''n''' is the unit vector normal to the boundary.

The definition of the problem is completed with the ''initial conditions''

<span id="eq-7"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\phi ({\boldsymbol x},t_0) = \phi _0 ({\boldsymbol x}) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
|}

where <math display="inline">\phi _0 ({\boldsymbol x})</math> is the value of the transported variable at time <math display="inline">t=t_0</math>.

In this work we will consider cases for which <math display="inline">s\ge 0</math> only. This includes the following particular problems:

* Advection-diffusion-absorption (<math display="inline">|{\boldsymbol v}|\not =0</math>, <math display="inline">K \not = 0</math>, <math display="inline">s>0</math>).
* Advection-diffusion (<math display="inline">|{\boldsymbol v}|\not =0</math>, <math display="inline">K \not = 0</math>, <math display="inline">s=0</math>).
* Diffusion-absorption (<math display="inline">|{\boldsymbol v}| =0</math>, <math display="inline">K \not = 0</math>, <math display="inline">s>0</math>).
* Advection-absorption (<math display="inline">|{\boldsymbol v}|\not =0</math>, <math display="inline">K  = 0</math>, <math display="inline">s>0</math>).

In the above <math display="inline">K</math> is the  average difussion given by <math display="inline">K = \left[\frac{1}{n_d}\sum \limits _{i=1}^{n_d}(k_i)^2 \right]^{1/2}</math>, where <math display="inline">n_d</math> is the number of space dimensions (i.e. <math display="inline">n_d =2</math> for 2D problems).

===2.2 Finite increment calculus (FIC) expressions===

The governing equations [[#eq-1.a|1.a]] and [[#eq-1.b|1.b]] and the boundary conditions [[#eq-3|3]]-[[#eq-6|6]] are expressed using the FIC theory as <span id='citeF-55'></span>[[#cite-55|[55]]].

===Transport balance===

<span id="eq-8"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>r_t-\underline{\frac{1}{2} {\boldsymbol h}^T {\boldsymbol \nabla } r_t}\frac{1}{2} {\boldsymbol h}^T {\boldsymbol \nabla } r_t = 0\qquad \hbox{in } \Omega   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
|}

with <math display="inline">{\boldsymbol h}= [h_1,h_2,h_3]^T</math> in 3D.

===Boundary conditions===

<span id="eq-9.a"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\phi - \phi ^p =0 \quad \hbox{on }\Gamma _\phi  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (9.a)
|}

<span id="eq-9.b"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>r_\Gamma + \underline{\frac{1}{2}h_n r_t}\frac{1}{2}h_n r_t =0 \quad \hbox{on }\Gamma _q \quad ,  \quad  \hbox{with } h_n= {\boldsymbol h}^T {\boldsymbol n}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (9.b)
|}

Eqs.([[#eq-8|8]]) and ([[#eq-9.b|9.b]]) are obtained by expressing the balance of fluxes in an arbitrary prismatic space domain of size <math display="inline">h_1 \times h_2 \times h_3</math> within the global problem domain and at the Neumann boundary, respectively. The distances <math display="inline">h_i</math> are termed ''characteristic lengths'' of the FIC method. The variations of the transported variable within the balance domain are approximated by Taylor series expansions retaining one order higher terms than in the infinitesimal theory <span id='citeF-40'></span>[[#cite-40|[40]]]. This higher order expansions lead to the underlined terms in Eqs.([[#eq-8|8]]) and ([[#Boundary conditions|2.2]]). These terms lead naturally to stabilized numerical schemes.

Note that  as the characteristic length vector <math display="inline">\boldsymbol h</math> tends to zero the FIC governing equations gradually recover the standard infinitesimal form, giving in the limit (for <math display="inline">\mathbf{h}=\mathbf{0}</math>) <math display="inline">r_t=0</math> in <math display="inline">\Omega </math> and <math display="inline">r_{\scriptscriptstyle \Gamma } =0\quad \hbox{on }\Gamma _q</math>.

The stability and accuracy of the numerical solution depends on the components of the characteristic length  vector <math display="inline">\boldsymbol h</math>. At the discretization level  the length of <math display="inline">\boldsymbol h</math> is usually expressed as a proportion of a typical grid dimension <span id='citeF-40'></span>[[#cite-40|[40]]].

The FIC equations are the starting point for deriving different stabilized numerical methods. Combining the FIC equations with the Galerkin FEM leads to the so-called FIC-FEM procedure that has been successfully applied to the solution of many problems in convective transport, fluid and solid mechanics such as advection-diffusion <span id='citeF-40'></span><span id='citeF-42'></span><span id='citeF-43'></span><span id='citeF-51'></span>[[#cite-40|[40,42,43,51]]], diffusion-absorption and Helmholtz <span id='citeF-15'></span>[[#cite-15|[15]]], advection-diffusion-absorption <span id='citeF-52'></span><span id='citeF-55'></span>[[#cite-52|[52,55]]], advection-diffusion-reaction <span id='citeF-61'></span>[[#cite-61|[61]]], incompressible fluid flow <span id='citeF-53'></span><span id='citeF-54'></span><span id='citeF-56'></span><span id='citeF-58'></span><span id='citeF-60'></span>[[#cite-53|[53,54,56,58,60]]], fluid-structure-interaction <span id='citeF-45'></span><span id='citeF-49'></span><span id='citeF-57'></span>[[#cite-45|[45,49,57]]], particle-laden flows and standard and incompressible solid mechanics <span id='citeF-47'></span><span id='citeF-50'></span><span id='citeF-59'></span>[[#cite-47|[47,50,59]]]. The FIC approach has also been applied to the solution of a variety of problems in mechanics using the meshless finite point method  <span id='citeF-39'></span><span id='citeF-41'></span><span id='citeF-44'></span><span id='citeF-46'></span>[[#cite-39|[39,41,44,46]]].

==3 Definition of the characteristic length vector==

The characteristic length vector <math display="inline">\boldsymbol h</math> is designed so that the expression for all the stabilization matrices and vectors reduce to those given for the 1D case as reported in <span id='citeF-61'></span>[[#cite-61|[61]]].

With this objective in mind the following expression for the characteristic length vector has been chosen in this work

<span id="eq-10"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol h}=  {\boldsymbol h}_v + {\boldsymbol h}_r +  {\boldsymbol h}_{sc} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
|}

The characteristic length vectors in the r.h.s. of Eq.[[#eq-10|10]] are defined as follows.

==='''Streamline characteristic length vector''', h<sub>v</sub>===

In Eq.[[#eq-10|10]] <math display="inline">{\boldsymbol h}_v</math> is a length vector along the velocity direction defined as

<span id="eq-11"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol h}_v = \alpha _v l_v \frac{\mathbf{v}}{|\mathbf{v}|} = \alpha _v l_v \hat {\boldsymbol v}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
|}

where <math display="inline">\alpha _v</math> is a streamline stabilization parameter, <math display="inline">l_v</math> is a characteristic element dimension and <math display="inline">\hat {\boldsymbol v}</math> is a unit velocity vector.

==='''Absorption characteristic length vector''', h<sub>r</sub>===

In Eq.[[#eq-10|10]] <math display="inline">{\boldsymbol h}_r</math> is a characteristic length vector induced by absorption effects and defined as

<span id="eq-12.a"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol h}_r = {\boldsymbol H}_r{\boldsymbol \nabla }\phi \qquad \hbox{with} \quad {\boldsymbol H}_r = \frac{2sgn (r_t)}{r_s}\left[{\boldsymbol D}_s + \alpha _r D \hat {\boldsymbol v}\hat {\boldsymbol v}^T\right] </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (12.a)
|}

where <math display="inline">\alpha _r</math> is an absorption stabilization parameter, <math display="inline">{D} = \hat {\boldsymbol v}^T {\boldsymbol D} \hat{\boldsymbol v}</math> and <math display="inline">r_s</math> is the space residual defined as

<span id="eq-12.b"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>r_s :=r_t - \rho c {\partial \phi  \over \partial t} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (12.b)
|}

Matrix <math display="inline">{\boldsymbol D}_s</math> in Eq.[[#eq-12.a|12.a]] is defined for different element types as follows.

===<u>3-noded triangles and 4-noded tetrahedra</u>:===

<span id="eq-13"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol D}_s=\frac{s}{(n+1)}\sum \limits _{i=1}^{n} {\boldsymbol l}_i {\boldsymbol l}_i^T </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (13)
|}

where <math display="inline">{\boldsymbol l}_i</math> is the vector joining the baricenter of the element and the <math display="inline">i</math>th node, and <math display="inline">n</math> is the number of nodes of the element.

The diffusion introduced by matrix <math display="inline">{\boldsymbol D}_s</math> takes care of the instabilities induced by the irregularity of the triangular mesh near boundaries that develop parabolic layes <span id='citeF-62'></span>[[#cite-62|[62]]].

===<u>Any other element:</u> Dₛ=[0]===

==='''Shock capturing characteristic length vector''' h<sub>s<sub>c</sub></sub>===

In Eq.[[#eq-10|10]] <math display="inline">{\boldsymbol h}_{sc}</math> is a shock-capturing characteristic length vector in the direction of the gradient of the solution. This vector accounts for the Gibbs oscillations across characteristic internal/boundary layers for advection-diffusion problems. It is defined as

<span id="eq-14.a"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol h}_{sc} = {h}_{sc}  \widehat{{\boldsymbol \nabla }\phi }  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (14.a)
|}

where <math display="inline">\widehat{{\boldsymbol \nabla }\phi }=\frac{{\boldsymbol \nabla }\phi }{|{\boldsymbol \nabla }\phi |}</math> is the unit gradient vector, and

<span id="eq-14.b"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{h}_{sc} = (1-\beta ^2)\left[l_{sc} sgn (r_s) - \frac{2|{\boldsymbol \nabla }\phi |}{r_s} ({\boldsymbol D} + {\boldsymbol D}_s): ({\boldsymbol I}- \hat{\boldsymbol v}\hat {\boldsymbol v}^T )\right] </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (14.b)
|}

In Eq.[[#eq-14.b|14.b]] <math display="inline">\beta </math> is a non linear parameter that depends on  the angle <math display="inline">\theta </math> between the velocity vector and the gradient vector. It is defined as

<span id="eq-14.c"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\beta = \left\{\begin{array}{lll}1 & \hbox{ if } & \theta < \theta _c\\ \displaystyle \hat {\boldsymbol v}^T  \widehat{{\boldsymbol \nabla }\phi } & \hbox{ if } & \theta \ge \theta _c \end{array}\right. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (14.c)
|}

where <math display="inline">\theta _c</math> is a critical angle. In this work we have taken <math display="inline">\theta _c =20^\circ </math>.

The parameter <math display="inline">\beta </math> controls the amount of shock-capturing nonlinear diffusion active at any point of the domain <span id='citeF-62'></span>[[#cite-62|[62]]]. When the  gradient vector is parallel to the velocity vector then <math display="inline">\beta =1</math>, <math display="inline">h_{sc}</math> vanishes and the linear stabilization terms suffice to diminish spurious numerical oscillations about the layers. The term <math display="inline">1-\beta ^{2}</math> in Eq.[[#eq-14.b|14.b]] gradually increases the magnitude of the shock-capturing term from zero, when the  gradient vector is aligned to velocity vector, to a maximum value when the layer gradient is orthogonal to the velocity.

==4 FIC-FEM formulation for the multidimensional transient advection-diffusion-absorption problem==

===4.1 Weighted residual form of the FIC equations===

The weighted residual form of the FIC governing equations [[#eq-8|8]] and [[#Boundary conditions|2.2]] is written as

<span id="eq-15"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int _\Omega W (r_t - \frac{1}{2} {\boldsymbol h}^T {\boldsymbol \nabla } r_t)d\Omega + \int _{\Gamma _q}W (-q_n +  q_n^p+\frac{1}{2}h_n r_t) d\Gamma =0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
|}

where <math display="inline">W</math> are arbitrary weighting functions.

Integrating by parts the FIC term in the first integral of [[#eq-15|15]] gives

<span id="eq-16"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int _\Omega \left(W + \frac{1}{2} ({\boldsymbol \nabla }^T W) {\boldsymbol h}   \right)r_t  d\Omega + \oint _{\Gamma _q} W(-q_n + q_n^p) d\Gamma =0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (16)
|}

Note that the FIC term has vanished from the boundary integral, as it usual in the FIC-FEM approach <span id='citeF-40'></span><span id='citeF-48'></span>[[#cite-40|[40,48]]].

Let us substitute the expression for the characteristic vector <math display="inline">{\boldsymbol h}</math> of Eq.[[#eq-10|10]] into Eq.[[#eq-16|16]]. This gives (using Eqs.[[#eq-11|11]], [[#eq-12.a|12.a]] and [[#eq-14.a|14.a]])

<span id="eq-17"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int _\Omega \left[W r_t + \frac{1}{2} ({\boldsymbol \nabla }^T W) \left(\alpha _v l_v \hat{\boldsymbol v} + {\boldsymbol H}_r {\boldsymbol \nabla }\phi + {h}_{sc} \frac{{\boldsymbol \nabla }\phi }{|{\boldsymbol \nabla }\phi |}\right)r_t \right]d\Omega + \oint _{\Gamma _q} W (-q_n + q_n^p) d\Gamma =0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (17)
|}

The final step is the integration by parts of the diffusive term in the expression of <math display="inline">r_t</math> in the first term of the first integral of Eq.[[#eq-17|17]]. This gives, after replacing the definition of <math display="inline">{\boldsymbol H}_r </math> in Eq.[[#eq-17|17]] and grouping some terms, the following expression for the weak variational form of the FIC governing equations

<span id="eq-18"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\displaystyle \int _\Omega \left[\rho c \bar W \frac{\partial \phi }{\partial t} + \rho c W  {\boldsymbol v}^T {\boldsymbol \nabla } \phi + ({\boldsymbol \nabla }^T W) {\boldsymbol D}_T {\boldsymbol \nabla }\phi + W s\phi  + \frac{1}{2} ({\boldsymbol \nabla }^T W){\boldsymbol h}_v s\phi \right]d\Omega + </math>
|-
| style="text-align: center;" | <math>  \displaystyle \int _\Omega \frac{1}{2} ({\boldsymbol \nabla }^T W) {\boldsymbol h} \left[-  {\boldsymbol \nabla }^T ( {\boldsymbol D} {\boldsymbol \nabla }\phi ) - Q\right]d\Omega  - \int _\Omega WQ d\Omega + \oint _{\Gamma _q} W q_n^p d\Gamma =0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (18)
|}

The space weighting function <math display="inline">\bar W</math> is given by

<span id="eq-19"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\bar W = W +\frac{1}{2} {\boldsymbol h}_v^T  {\boldsymbol \nabla } W   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (19)
|}

The expression of the total diffusivity matrix <math display="inline">{\boldsymbol D}_T </math> in Eq.[[#eq-18|18]] is

<span id="eq-20"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\displaystyle {\boldsymbol D}_T  = {\boldsymbol D} + \alpha _v {\boldsymbol D}_v + {\boldsymbol D}_s + \alpha _r {D} \hat {\boldsymbol v} \hat {\boldsymbol v}^T + {\boldsymbol D}_{sc}{\boldsymbol I} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (20)
|}

where <math display="inline">{\boldsymbol D}</math> and <math display="inline">{\boldsymbol D}_s</math> are defined in Eqs.[[#eq-2|2]] and [[#eq-13|13]], respectively, '''I''' is the unit matrix and

<span id="eq-21.a"></span>
<span id="eq-21.b"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol D}_v \!\!\!\!=\!\! \frac{l_v}{2} \hat {\boldsymbol v} {\boldsymbol v}^T</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (21.a)
|-
| style="text-align: center;" | <math> {\boldsymbol D}_{sc} \!\!\!\!=\!\! \left(\frac{1}{2}  l_{s_c} \frac{|r_t|}{|{\boldsymbol \nabla }\phi |} - ({\boldsymbol D} + {\boldsymbol D}_s): ({\boldsymbol I}- \hat{\boldsymbol v}\hat {\boldsymbol v}^T )\right)(1-\beta ^2 ) </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (21.b)
|}
|}

===4.2 FIC-FEM equations===

We interpolate the transported variable <math display="inline">\phi </math> in the standard FEM fashion over a mesh of elements with <math display="inline">n</math> nodes <span id='citeF-67'></span><span id='citeF-68'></span>[[#cite-67|[67,68]]] as

<span id="eq-22"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\phi \simeq \hat \phi = \sum \limits _{i=1}^n N_i \phi _i </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (22)
|}

where <math display="inline">N_i</math> are the space shape functions and <math display="inline">\phi _i</math> the nodal variables.

Introducing Eq.[[#eq-22|22]] into the FIC variational form [[#eq-18|18]] and using a Galerkin approach (<math display="inline">W_i =N_i</math>) gives the final system of discretized equations as

<span id="eq-23"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{M} \dot{\boldsymbol \phi } + [\mathbf{K} + \mathbf{C} + \mathbf{S}] {\boldsymbol \phi } = {\boldsymbol f} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (23)
|}

where <math display="inline">{\boldsymbol \phi } = [\phi _1,\phi _2,\cdots , \phi _N]^T</math> is the vector of nodal unknowns, with <math display="inline">N</math> being total number of nodes in the mesh.

The rest of matrices and vector <math display="inline">{\boldsymbol f}</math> in Eq.[[#eq-23|23]] are obtained in the standard FEM fashion by assembling the element contribution given by

<span id="eq-24"></span>
<span id="eq-25"></span>
<span id="eq-26"></span>
<span id="eq-27"></span>
<span id="eq-28"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{M}_{ij}^e = \int _{\Omega ^e} \rho c N_i N_j d\Omega + \int _{\Omega ^e} \rho c  \frac{1}{2} ({\boldsymbol \nabla }^T N_i ){\boldsymbol h}_v N_j d\Omega </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (24)
|-
| style="text-align: center;" | <math> \mathbf{K}_{ij}^e = \int _{\Omega ^e} ({\boldsymbol \nabla }^T N_i ) {\boldsymbol D}_T {\boldsymbol \nabla } N_j d\Omega - \int _{\Omega ^e} \frac{1}{2} ({\boldsymbol \nabla }^T N_i) {\boldsymbol h} {\boldsymbol \nabla }^T ({\boldsymbol D} {\boldsymbol \nabla } N_j)d\Omega </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (25)
|-
| style="text-align: center;" | <math> \mathbf{C}^e_{ij} =  \int _{\Omega ^e} \rho c  N_i {\boldsymbol v}^T {\boldsymbol \nabla }    N_j d\Omega </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (26)
|-
| style="text-align: center;" | <math> \mathbf{S}^e_{ij} =  \int _{\Omega ^e} s N_i N_j d\Omega + \int _{\Omega ^e}  \frac{1}{2} ({\boldsymbol \nabla }^T N_i ){\boldsymbol h}_v s N_j d\Omega </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (27)
|-
| style="text-align: center;" | <math> {\boldsymbol f}_i^e = \int _{\Omega ^e} \left(N_i + \frac{1}{2}( {\boldsymbol \nabla }^T N_i) {\boldsymbol h} \right)Q d\Omega - \oint _{\Gamma _q^e} N_i q_n^p d\Gamma  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (28)
|}
|}

We note that the second integral in Eq.[[#eq-25|25]] vanishes for linear finite element approximations, such as those used in this work.

Good results have been obtained for the solution of the problems presented in this paper using a lumped form of matrix <math display="inline">{\boldsymbol M}</math>. An exception is the problem of the evolution of a concentration field presented in Section [[#7.6 Advection-diffusion of a concentration field |7.6]] for which using the consistent form of matrix <math display="inline">{\boldsymbol M}</math> is essential to reduce excessive diffusion.

==5 Transient solution scheme==

===5.1 Implicit time integration scheme===

We discretize in time the system of equations ([[#eq-23|23]]) using a Generalized Trapezoidal rule <span id='citeF-11'></span><span id='citeF-67'></span>[[#cite-11|[11,67]]]. The solution for the nodal values at a time instant is found using an incremental iterative strategy as

<span id="eq-29"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{}^{i}{\boldsymbol H}^{n} \Delta {\boldsymbol \phi }= - {}^i{\boldsymbol r}_t^n </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (29)
|}

where <math display="inline">\Delta {\boldsymbol \phi }</math> is the increment of the nodal variables, <math display="inline">\theta </math> is a non dimensional time parameter (<math display="inline">0.5 <\theta \le 1</math> is required for the integration scheme to be unconditionally stable <span id='citeF-11'></span><span id='citeF-67'></span><span id='citeF-68'></span>[[#cite-11|[11,67,68]]]), <math display="inline">(\cdot )^n</math> denotes values at time <math display="inline">t=t_n</math>, <math display="inline">{}^i(\cdot )</math> denotes values at the <math display="inline">i</math>th iteration and

<span id="eq-30"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{}^{i}{\boldsymbol H}^{n} = \frac{1}{\theta \Delta t}{\boldsymbol M} + {}^i{\boldsymbol K}^{n} + {\boldsymbol C}+{\boldsymbol S} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (30)
|}

<span id="eq-31"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{}^i{\boldsymbol r}_t^n := \mathbf{M} \dot{\boldsymbol \phi } + [{}^i\mathbf{K}^n + \mathbf{C} + \mathbf{S}] {}^{i}{\boldsymbol \phi }^{n+\theta } - {}^{i}{\boldsymbol f}^{n} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (31)
|}

In Eq.[[#eq-31|31]] we define <math display="inline">\dot {\boldsymbol \phi } = \frac{{\boldsymbol \phi }^{n+\theta }-{\boldsymbol \phi }^n}{\theta \Delta t}</math>.

From the value of <math display="inline">\Delta {\boldsymbol \phi }</math> obtained from solving Eq.[[#eq-29|29]] we compute the value of <math display="inline">{\boldsymbol \phi }^{n+\theta }</math> at the <math display="inline">i+1</math> iteration as

<span id="eq-32"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{}^{i+1}{\boldsymbol \phi }^{n+\theta } = {}^{i}{\boldsymbol \phi }^{n+\theta }+ \Delta {\boldsymbol \phi } </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (32)
|}

The iterative solution at <math display="inline">t_{n+1}</math>, <math display="inline">{}^{i+1}{\boldsymbol \phi }^{n+1}</math>, is  obtained as

<span id="eq-33"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{}^{i+1}{\boldsymbol \phi }^{n+1} = \frac{1}{\theta } {}^{i+1}{\boldsymbol \phi }^{n+\theta } + \left(1-\frac{1}{\theta }\right){\boldsymbol \phi }^n  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (33)
|}

The non-linearity in the expression of matrix <math display="inline">\boldsymbol H</math> in Eq.[[#eq-30|30]] is due to the dependence of the shock-capturing term in matrix <math display="inline">\boldsymbol K</math> via matrix <math display="inline">{\boldsymbol D}_{sc}</math> (see Eq.[[#eq-21.b|21.b]]).

The iterations proceed until convergence is achieved for both the unknown <math display="inline">\phi </math> and the residual measured in a <math display="inline">L_2</math> norm. In the transient problems solved in this work, convergence within each time step was typically achieved in 2&#8211;3 iterations.

===5.2 Explicit time integration scheme===

If we rewrite the residual in Eq.[[#eq-31|31]] choosing <math display="inline">\theta =0</math> and defining now <math display="inline">\dot{\boldsymbol \phi }</math> as <math display="inline">\frac{{\boldsymbol \phi }^{n+1}-{\boldsymbol \phi }^n}{\Delta t}</math> leads to the following explicit solution scheme

<span id="eq-34"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{\boldsymbol \phi }^{n+1} ={\boldsymbol \phi }^n + \Delta t {\boldsymbol M}^{-1}_d\left[{\boldsymbol f}^n - ({\boldsymbol K}^n + {\boldsymbol C} - {\boldsymbol S}) {\boldsymbol \phi }^n\right] </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (34)
|}

where <math display="inline">{\boldsymbol M}_d</math> denotes the lumped diagonal form of matrix <math display="inline">{\boldsymbol M}</math>.

In the expression of <math display="inline">{\boldsymbol K}</math> of Eq.[[#eq-34|34]] the shock-capturing terms are taken as constant within a time step.

==6 Computation of the stabilization parameters==

===6.1 General expression of the stabilization parameters===

The optimal value of the stabilization parameters <math display="inline">\alpha _v</math> and <math display="inline">\alpha _r</math> giving nodally exact solutions is quite difficult in the transient case due to the multiple forms that the solution can take as it evolves in time. We  present a procedure for computing a quasi-optimal value of <math display="inline">\alpha _v</math> and <math display="inline">\alpha _r</math> that has proved to yield  accurate results for transient advection-diffusion-absorption problems evolving towards a steady state solution, as well as for solutions involving the pure advection of a discontinuous function. The method is an extension of the approach proposed and successfully tested for 1D transient problems in <span id='citeF-61'></span>[[#cite-61|[61]]].

The transient equation ([[#eq-1.b|1.b]]) can be written as follows

<span id="eq-35"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\rho c {\boldsymbol v}^T {\boldsymbol \nabla }\phi - {\boldsymbol \nabla }^T ( {\boldsymbol D} {\boldsymbol \nabla }\phi ) + \bar s \phi - Q =0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (35)
|}

where

<span id="eq-36"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\bar s =s + s_t \quad \hbox{with } s_t = \rho c \frac{\dot \phi }{\phi }  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (36)
|}

Eq.([[#eq-35|35]]) defines a pseudo-stationary problem in which a non linear reaction term  has been introduced.

The non linear reaction term <math display="inline">s_t</math> can be approximated as

<span id="eq-37"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>s_t = \frac{\rho c}{\theta \Delta t} f(\kappa )  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (37)
|}

with

<span id="eq-38"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>f(\kappa ) \ge 0 \quad \hbox{and}\quad \kappa  = \frac{\phi ^{n+\theta }  - \phi ^n }{\phi ^{n+\theta }}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (38)
|}

The value of <math display="inline">s_t</math> in Eq.[[#eq-37|37]] has been limited in this work to a ten percent of the value of the actual absorption parameter <math display="inline">s</math>.

From Eq.([[#eq-36|36]]) we can define an equivalent Damköhler number <math display="inline">\bar \sigma _v</math> as

<span id="eq-39"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\bar \sigma  =\sigma + \sigma _t \quad \hbox{with } \quad \sigma _t = \frac{1}{\theta C} f(\kappa )  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (39)
|}

where <math display="inline">\sigma </math> is the  Damköhler number defined in Eq.[[#eq-44|44]] and <math display="inline">C= \frac{v \Delta t}{l^e}</math> is the element Courant number where <math display="inline">l^e</math> is a characteristic element dimension. For 1D problems <math display="inline">l^e</math> is taken as the element length. For 2D problems we have chosen <math display="inline">l^e =\sqrt{2\Omega ^e}</math> where <math display="inline">\Omega ^e</math> is the area of the element.

The following definition for <math display="inline">f(\kappa )</math> has been chosen in the examples solved in this work

<span id="eq-40"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>f(\kappa )= \frac{2}{\tanh{1}}\tanh \left(\delta \frac{\big|\phi ^{n}  - \phi ^{n-1} \big|_\infty ^e}{\big|\phi ^{n}+\phi ^{n-1} \big|_\infty ^e}\right)  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (40)
|}

where <math display="inline">|a|_\infty ^e</math> denotes the maximum value of <math display="inline">a</math> within an element and <math display="inline">\delta </math> is a positive number that controls the slope of the function <math display="inline">\tanh (\cdot )</math> that ranges from zero to one.

Function <math display="inline">f(\kappa )</math> has been designed so that <math display="inline">f(\kappa )=0</math> (and <math display="inline">\sigma _t =0</math>) for a steady-state problem (or in zones where <math display="inline">\dot \phi =0</math>), and <math display="inline">f(\kappa ) = \frac{2}{\tanh{1}}</math> (and <math display="inline">\sigma _t = \frac{2}{\theta C}</math>) for cases when <math display="inline">\phi (x,t)</math> suddenly changes from a zero value to a finite value at a node.

Good results were obtained in the transient problems solved using Eq.([[#eq-40|40]]) with <math display="inline">\delta = 1</math>. The optimal choice of <math display="inline">\delta </math> in terms of the nature of the transient solution is a matter that deserves further research.

We highlight that for the explicit time integration scheme we have used <math display="inline">\theta = 1</math> in the definition of <math display="inline">s_t</math> and <math display="inline">\sigma _t</math>.

The stabilization parameters <math display="inline">\alpha _v</math> and <math display="inline">\alpha _r</math> in Eqs.[[#eq-11|11]] and [[#eq-12.a|12.a]] for the transient convection-diffusion-radiation problem are defined as follows

<span id="eq-41"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\alpha _v = \begin{cases}\displaystyle \frac{2}{\bar \sigma } \left(1 - \frac{\bar \sigma \tanh \gamma }{\xi -1} \right)& , \quad  \bar \sigma \ge 2^{-12}\\ \displaystyle \frac{\bar \sigma }{3} + \bar \alpha _v \left(1 - \frac{\sigma }{\gamma } \right)& , \quad  \bar \sigma < 2^{-12} \end{cases} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (41)
|}

with

<span id="eq-42"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\bar \alpha _v = \coth \gamma - \frac{1}{\gamma } </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (42)
|}

and

<span id="eq-43"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\alpha _r =\gamma \left[\frac{\sigma }{\varphi } \left(\frac{\xi -1 +\varphi }{\xi - 1} \right)- \alpha _v\right] - 1-\frac{1}{D}\hat {\boldsymbol v}^T  {\boldsymbol D}_s \hat{\boldsymbol v} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (43)
|}

In the above expressions

<span id="eq-44"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{l}\displaystyle \xi = \frac{\cosh \lambda }{\cosh \gamma }\quad \hbox{with } \quad \lambda = \left(\gamma ^2 + w^2   \right)^{1/2}\\ \displaystyle \gamma =  \frac{\rho c |{\boldsymbol v}| l_v}{2D}\quad , \quad w = \frac{\rho c s(l_v)^2}{D} \quad , \quad \sigma = \frac{sl_v}{|{\boldsymbol v}|}= \frac{w}{2\gamma } \end{array}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (44)
|}

The expressions  of <math display="inline">\alpha _v</math> and <math display="inline">\alpha _r</math> in Eqs.[[#eq-41|41]] and [[#eq-43|43]] are an extension of the values for the steady-state problem presented in <span id='citeF-62'></span>[[#cite-62|[62]]].

In Eq.[[#eq-43|43]] <math display="inline">\varphi </math> is a constant such that <math display="inline">2\le \varphi \le 3</math>. The “exact” expression of <math display="inline">\alpha _r</math> for 1D problems requires choosing <math display="inline">\varphi =3</math> <span id='citeF-62'></span>[[#cite-62|[62]]]. In our computations for 2D problems we have obtained good results using <math display="inline">\varphi =2</math>.

The characteristic length <math display="inline">l_v</math> is defined as the size of the element in the direction of the velocity <math display="inline">{\boldsymbol v}</math>. Representing the element's edge joining nodes <math display="inline">a</math> and <math display="inline">b</math> with the vector <math display="inline">{\boldsymbol l}_{ab}</math>, <math display="inline">l_v</math> is computed as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>l_v = \underset{\hbox{edges}}{\max } \left\{{\boldsymbol v}^T {\boldsymbol l}_{ab}\right\} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (45)
|}

The characteristic length <math display="inline">l_{sc}</math> is generally defined as

<span id="eq-46"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>l_{sc} =\sqrt{2\Omega ^e} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (46)
|}

where <math display="inline">\Omega ^e</math> is the element area. For elements belonging to the boundary in meshes of 3-noded triangles we have chosen

<span id="eq-47"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>l_{sc} =2\sqrt{\Omega ^e} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (47)
|}

The expression of <math display="inline">l_{sc}</math> in Eq.[[#eq-47|47]] introduces a higher value of the shock-capturing term at the boundaries where sharp gradients exist in directions different from the velocity directions.

==7 Examples==

In the examples shown next the density <math display="inline">\rho </math> and the specific flux c are chosen such that <math display="inline">\rho c = 1</math>. We have also assumed an isotropic diffusion of value <math display="inline">k</math>.

===7.1 Transient advection-diffusion-absorption problem===

The analysis domain <math display="inline">(x,y) = [0,8] \times [0,8]</math> is discretised into a regular mesh of <math display="inline">(2\times 8) \times (2\times 8)</math> 3-noded triangles of unit rectangular side (<math display="inline">l_v =1</math>) (fig:dac_scheme). The advection, diffusion and reaction coefficients are chosen as <math display="inline">v_1 = 8</math>, <math display="inline">k = 2</math> and <math display="inline">s = 2</math>. The transported variable is the temperature. The schematics of the problem can be seen in fig:dac_scheme. The problem data yields the dimensionless numbers <math display="inline">\gamma _v = 2</math> and <math display="inline">\omega _v = 1</math>.

The Dirichlet boundary conditions <math display="inline">\phi (x=0)=3</math> and <math display="inline">\phi (x=8)=8</math> are employed. The initial solution is chosen to have a linear profile. The transient solution was obtained iteratively using the implicit scheme of Eqs.[[#eq-29|29]]&#8211;[[#eq-33|33]] with <math display="inline">\theta = 1.0</math> and  a time step of <math display="inline">\Delta t = 0.0625s</math>  (fig:dac_results). This corresponds to an element Courant number <math display="inline">C = 0.5</math>. An exponential layer gradually develops at the right boundary which triggers a global instability in the Galerkin FEM. As there is not significant dispersive phenomena this global instability is successfully controlled by the SUPG method (well-known result) and the FIC-FEM method developed in this work. Figure [[#img-3|3]] show a perspective view of the quasi-steady-state solution for the temperature for <math display="inline">t=2s</math>. The steady state solution matches the 1D result reported in <span id='citeF-61'></span>[[#cite-61|[61]]].

<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_437326896-dac_scheme.png|420px|Transient advection-absorption problem. Square domain with linear velocity and zero source.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 1:''' Transient advection-absorption problem. Square domain with linear velocity and zero source.
|}

<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_437326896-dac_results.png|420px|Transient advection-absorption problem. Transient FIC-FEM solution obtained with a structured mesh of 2 ×8 ×8 three-noded triangles. The transient solutions are plotted along line  A-A' at times 0.0 s, 0.125 s, 0.25 s, 0.5 s, 1 s and 2 s.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 2:''' Transient advection-absorption problem. Transient FIC-FEM solution obtained with a structured mesh of <math>2 \times 8 \times 8</math> three-noded triangles. The transient solutions are plotted along line  A-A' at times 0.0 s, 0.125 s, 0.25 s, 0.5 s, 1 s and 2 s.
|}

<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_437326896-dac_image.png|420px|Transient advection-diffusion-absorption problem. Solution obtained at t = 2 s.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 3:''' Transient advection-diffusion-absorption problem. Solution obtained at t = 2 s.
|}

===7.2 Skewed flow is a square domain with non-uniform boundary conditions===

The domain is <math display="inline">(x,y) = [0,1]\times [0,1]</math>. The problem data is: <math display="inline">{\boldsymbol v}=[-5\cdot{10}^6,-9\cdot{10}^6]^T</math>, <math display="inline">k=1</math>, <math display="inline">\rho c =1</math>, <math display="inline">s=0</math> and <math display="inline">Q=0</math>. The transported variable is again the temperature. The boundary conditions are: <math display="inline">\phi=1</math> on (<math display="inline">x=1,y>0.7)\cup (x<1,y=1)</math>, <math display="inline">\phi=0.5</math> at (<math display="inline">x=1,y=0.7</math>) and <math display="inline">\phi=0</math> on the rest of the boundary (fig:SC_scheme). The transient solution was obtained using the implicit scheme here presented with <math display="inline">\theta = 0.8</math> and a time step of <math display="inline">\Delta t = 1</math>e-09 s. Convergence within each time step was found in 2-3 iterations. The solution develops an exponential boundary layer at the outflow boundary and an internal characteristic layer (parabolic) which is skewed to the mesh and the boundary. Both layers are subgrid phenomena for the considered mesh resolution. Both the exponential and characteristic layers are reproduced in the FIC-FEM solution without spurious oscillations near the layers. The FIC-FEM result of the evolution of the solution in time towards steady state are shown in fig:SC_results. Figure [[#img-6|6]] shows a perspective view of the quasi steady-state solution for <math display="inline">t=2</math>e-07s. The steady-state solution agrees with that reported in <span id='citeF-51'></span><span id='citeF-55'></span><span id='citeF-62'></span>[[#cite-51|[51,55,62]]].

<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_437326896-SC_scheme.png|420px|Square domain with non-uniform Dirichlet conditions, downwards diagonal velocity and zero source. Mesh of 2×20×20 three-noded triangles]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 4:''' Square domain with non-uniform Dirichlet conditions, downwards diagonal velocity and zero source. Mesh of <math>2\times 20\times 20</math> three-noded triangles
|}

<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_437326896-SC_results.png|420px|Square domain with non-uniform Dirichlet conditions. FIC-FEM solution along line A-A'. The transient solutions are plotted at times 5e-08 s, 6e-08 s, 7e-08 s, 1e-07 s and 2e-07 s.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 5:''' Square domain with non-uniform Dirichlet conditions. FIC-FEM solution along line A-A'. The transient solutions are plotted at times 5e-08 s, 6e-08 s, 7e-08 s, 1e-07 s and 2e-07 s.
|}

<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_437326896-SC_image.png|420px|Square domain with non-uniform Dirichlet conditions. Solution at time 2e-07 s.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 6:''' Square domain with non-uniform Dirichlet conditions. Solution at time 2e-07 s.
|}

===7.3 Uniform advection problem with constant source===

We study the uniform advection of a temperature field with a constant source term. The problem data is: <math display="inline">{\boldsymbol v}=[1,0]^T</math>, <math display="inline">k=10^{-8}</math>, <math display="inline">\rho c=1</math>, <math display="inline">s=0</math> and <math display="inline">Q=1</math>. The homogeneous boundary condition <math display="inline">\phi =0</math> is imposed everywhere as shown in  fig:Q_source_scheme. The solution develops an exponential layer at the outflow boundary <math display="inline">x=1</math> and parabolic layers at the boundaries <math display="inline">y=0</math> and <math display="inline">y=1</math>. The FIC-FEM results of the transient solution towards steady-state are shown in fig:Q_source_results. A perspective of the quasi steady-state solution at time <math display="inline">t=1.5 s</math> is shown in Figure [[#img-9|9]]. Both layers are subgrid phenomena for the considered mesh resolution. The transient solution was obtained iteratively using the implicit scheme presented in the paper with <math display="inline">\theta = 0.8</math> and a time step of <math display="inline">\Delta t = 0.005 s</math>. Convergence within each time step was found in 2-3 iterations. The steady-state solution agrees with that reported in <span id='citeF-62'></span>[[#cite-62|[62]]].

<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_437326896-Q_source.png|420px|Square domain with linear velocity and constant source. Structured mesh of 2×20×20 three node triangles]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 7:''' Square domain with linear velocity and constant source. Structured mesh of <math>2\times 20\times 20</math> three node triangles
|}

<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_437326896-Q_source_results.png|420px|Square domain with linear velocity and constant source. FIC solution obtained along line A-A'. The transient solutions are plotted at times 0.25 s, 0.5 s, 0.75 s, 1 s and 1.5 s.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 8:''' Square domain with linear velocity and constant source. FIC solution obtained along line A-A'. The transient solutions are plotted at times 0.25 s, 0.5 s, 0.75 s, 1 s and 1.5 s.
|}

<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_437326896-Q_source_image.png|420px|Square domain with linear velocity and constant source. Solution at time 1.5 s.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 9:''' Square domain with linear velocity and constant source. Solution at time 1.5 s.
|}

===7.4 Advective-diffusive transport of the temperature in a rectangular domain with Neumann and non-uniform Dirichlet conditions, rotational velocity field and zero source===

The rectangular domain of sides 1 <math display="inline">\times </math> 2 units and the boundary conditions are shown in fig:rotational_scheme. The rotational velocity field is defined as <math display="inline">\mathbf{u}=10^4 \times [y(1-x^2),-x(1-y^2)]^T</math>. A unit isotropic diffusion is assumed.  The value of the temperature (<math display="inline">\phi </math>) is prescribed at part of the domain sides as shown in fig:rotational_scheme. In the rest of the sides, the Neumann condition (<math display="inline">q^p_n=0</math>) is imposed. The problem was solved with the uniform mesh of 3-noded triangles shown in fig:rotational_scheme. The correct solution is a uniform plateau of <math display="inline">\phi =100</math> with a boundary layer at the right hand side where the solution drops towards the prescribed value of <math display="inline">\phi =0</math> and a circular sharp gradient region around the lower circular zone where the solution takes a zero value. The FIC-FEM results of the transient solution towards steady-state are shown in fig:rotational_results. A perspective of the quasi steady-state solution at time <math display="inline">t=1</math>e-03s is shown in Figure [[#img-12|12]]. The transient solution was obtained using the implicit scheme presented in this work with <math display="inline">\theta = 0.8</math> and a time step of <math display="inline">\Delta t = 1</math>e-05s. Convergence within each time step was found in 2-3 iterations. The steady-state solution agrees with that reported in <span id='citeF-51'></span>[[#cite-51|[51]]].

<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_437326896-rotational_scheme.png|420px|Advective-diffusive transport of the temperature in a rectangular domain with Neumann and non-uniform Dirichlet conditions, rotational velocity and zero source. Structured mesh of 2 x 40 x 20 three-noded triangles]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 10:''' Advective-diffusive transport of the temperature in a rectangular domain with Neumann and non-uniform Dirichlet conditions, rotational velocity and zero source. Structured mesh of 2 x 40 x 20 three-noded triangles
|}

<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_437326896-rotational_results.png|420px|FIC solution of the problem in Figure [[#img-10|10]]. The transient solutions are plotted along line A-A' at times 5e-05 s, 1e-04 s, 3e-04 s, 4e-04, 5e-04 and 1e-03 s.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 11:''' FIC solution of the problem in Figure [[#img-10|10]]. The transient solutions are plotted along line A-A' at times 5e-05 s, 1e-04 s, 3e-04 s, 4e-04, 5e-04 and 1e-03 s.
|}

<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_437326896-rotational_image.png|420px|Solution at 1e-03 s.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 12:''' Solution at 1e-03 s.
|}

===7.5 Pure convection of a temperature plateau===

We study the pure convection of an initial temperature plateau. The domain dimensions are <math display="inline">(x,y) = [0,1]\times [0,1]</math>. The problem data is: <math display="inline">{\boldsymbol v}=[1,0]^T</math>, <math display="inline">k=0</math>, <math display="inline">s=0</math>, <math display="inline">\rho c =1</math> and <math display="inline">Q=0</math>. The homogeneous boundary condition <math display="inline">\phi =0</math> is imposed at the left wall as shown in  fig:explicit_wave_scheme. The temperature plateau is generated by imposing a value of the temperature <math display="inline">\phi =1.0</math> from <math display="inline">x=0.1</math> to <math display="inline">x=0.2</math>. The problem has been solved using a regular mesh of <math display="inline">50\times 50</math> four-noded quadrilateral elements. Figures [[#img-14|14]] and [[#img-15|15]] show that the initial solution maintains its shape through the domain. The transient solution was obtained using an explicit scheme and  a time step of <math display="inline">\Delta t = 0.2 s</math>, giving an element Courant number of <math display="inline">C= 1</math>.

<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_437326896-explicit_wave_scheme.png|420px|Pure convection of a temperature plateau in a square domain with horizontal velocity. Structured mesh of 50 x 50 4-noded quadrilaterals.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 13:''' Pure convection of a temperature plateau in a square domain with horizontal velocity. Structured mesh of 50 x 50 4-noded quadrilaterals.
|}

<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_437326896-explicit_wave_results.png|420px|Pure convection of a temperature plateau. FIC-FEM solution plotted along line A-A' at times 0 s, 0.2 s and 0.7 s.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 14:''' Pure convection of a temperature plateau. FIC-FEM solution plotted along line A-A' at times 0 s, 0.2 s and 0.7 s.
|}

<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_437326896-explicit_wave_image.png|420px|Pure convective case. FIC-FEM solution at time 0.7 s.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 15:''' Pure convective case. FIC-FEM solution at time 0.7 s.
|}

===7.6 Advection-diffusion of a concentration field ===

In this problem we study the advection and diffusion of a concentration that is instantaneously  released as a Dirac source in a point of a stream. This problem is generally known as the plume transport problem <span id='citeF-16'></span>[[#cite-16|[16]]].

The domain dimensions are <math display="inline">(x,y) = [0,35]\times [0,10]</math>. The problem data is: <math display="inline">{\boldsymbol v} =[1,0]^T</math>, <math display="inline">\rho c =1</math>, <math display="inline">k=0.1</math>, <math display="inline">s=0</math> and <math display="inline">Q=0</math>. The homogeneous boundary condition <math display="inline">\phi =0</math> is imposed on the left wall as shown in  fig:plume_scheme. The value <math display="inline">\phi _0 =1000</math> is imposed as a initial value of the concentration at the node <math display="inline">(x,y)=(2,5)</math>.

The analytical solution for this 2D problem is <span id='citeF-16'></span>[[#cite-16|[16]]]

<span id="eq-48.a"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\phi (x,y,t)= \frac{\phi _0}{L_3 4\pi Dt} e^{-A} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (48.a)
|}

with

<span id="eq-48.b"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>A = - \frac{1}{4Dt} \left\{[x_1 - (x_1^0 +v_1 t)]^2 + [x_2 - (x_2^0 + v_2 t)]^2  \right\} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (48.b)
|}

In Eqs.[[#7.6 Advection-diffusion of a concentration field |7.6]] <math display="inline">\phi </math> is the initial value of the concentration that is dropped at the point with coordinates <math display="inline">{\boldsymbol x}^0 = [x_1^0,x_2^0]^T</math>, <math display="inline">L_3</math> is the vertical dimension of the analysis domain (in this problem we have taken <math display="inline">L_3 =1</math>), <math display="inline">v_1</math>  and <math display="inline">v_2</math>  are the horizontal and vertical components of the velocity vector, respectively and <math display="inline">D=\frac{K}{\rho c}</math> where <math display="inline">K</math> is the average diffusion defined in Section 2.1.

The transient solution was obtained using an implicit scheme with <math display="inline">\theta = 0.5</math> and a time step of <math display="inline">\Delta t = 0.25 s</math>, giving a Courant number of <math display="inline">C= 0.5</math>. We highlight that the shock-capturing terms were not taken into account for solving this problem. The use of the consistent form of the mass matrix was essential for obtaining good results for this problem. fig:plume_results_1 shows how the solution is advected and diffuses as time increases. The numerical results for the distribution of the contraction along line A-A' are in agreement with the analytical solution shown in fig:plume_results_2. The graphs are taken along the A-A' line, which goes from (0,5) to (10,5). A comparison of the contours of the numerical and analytical distributions of the concentration at <math display="inline">t=15s</math> can be seen in Figures [[#img-19|19]] and [[#img-20|20]], respectively.

<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_437326896-plume_scheme.png|540px|Advection-diffusion of a concentration field. Rectangular domain with horizontal velocity. Structured mesh of 2 x 70 x 20 three-noded triangles]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 16:''' Advection-diffusion of a concentration field. Rectangular domain with horizontal velocity. Structured mesh of 2 x 70 x 20 three-noded triangles
|}

<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_437326896-plume_results_1.png|420px|Advection-diffusion of a concentration field. FIC-FEM solution along line A-A' obtained with the mesh of Figure [[#img-16|16]]. Structured mesh of 2 x 70 x 20 three-noded triangles. Solution at times 0 s, 0.25 s, 0.5 s, 0.75 s, 1.0 s, 1.25 s, 1.5 s, 1.75 s, 2.0 s, 4.5 s and 5.5 s.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 17:''' Advection-diffusion of a concentration field. FIC-FEM solution along line A-A' obtained with the mesh of Figure [[#img-16|16]]. Structured mesh of 2 x 70 x 20 three-noded triangles. Solution at times 0 s, 0.25 s, 0.5 s, 0.75 s, 1.0 s, 1.25 s, 1.5 s, 1.75 s, 2.0 s, 4.5 s and 5.5 s.
|}

<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_437326896-plume_results_2.png|420px|Advection-diffusion of a concentration field. Analytical solution along line A-A' plotted on the structured mesh of Figure [[#img-16|16]]. Solution at times 0 s, 0.25 s, 0.5 s, 0.75 s, 1.0 s, 1.25 s, 1.5 s, 1.75 s, 2.0 s, 4.5 s and 5.5 s.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 18:''' Advection-diffusion of a concentration field. Analytical solution along line A-A' plotted on the structured mesh of Figure [[#img-16|16]]. Solution at times 0 s, 0.25 s, 0.5 s, 0.75 s, 1.0 s, 1.25 s, 1.5 s, 1.75 s, 2.0 s, 4.5 s and 5.5 s.
|}

<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_437326896-plume_num.png|510px|Advection-diffusion of a concentration field. Contours of the FIC-FEM results at t = 15s.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 19:''' Advection-diffusion of a concentration field. Contours of the FIC-FEM results at t = 15s.
|}

<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_437326896-plume_ana.png|510px|Advection-diffusion of a concentration field. Contours of the analytical results at t = 15s.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 20:''' Advection-diffusion of a concentration field. Contours of the analytical results at t = 15s.
|}

==8 Concluding remarks==

In this work we have extended the steady-state multidimensional FIC-FEM formulation reported in <span id='citeF-62'></span>[[#cite-62|[62]]] to the transient case. The stabilized variational expression has the standard residual form typical of the FIC-FEM procedure. The stabilized terms introduced by the FIC approach depend on a characteristic element length and two stabilization parameters which expression is obtained as an extension of the steady-state forms presented in <span id='citeF-62'></span>[[#cite-62|[62]]]. A shock capturing term is introduced to account for the Gibbs oscillations across internal/boundary layers.

The good behavior of the FIC-FEM formulation has been verified for transient advection-diffusion-absorption problems with numerical solutions evolving to a steady-state, as well as to two fully transient problems.

The numerical results were obtained in most cases using an implicit scheme. An exception is the pure convection problem (Section 7.5) that was solved with an explicit scheme and a Courant number of <math display="inline">C =1</math>. The shock capturing terms were not taken into account in the solution of the fully transient problems (Sections 7.5 and 7.6). We also remark the importance of using a consistent stabilized mass matrix in the concentration field example.

==Acknowledgements==

This research was partially funded by the PRECISE project (BIA2017-83805-R) of the Natural Research Plan of the Spanish Government. Support for this work was also provided from the Office for Naval Research Global (ONRG) of the US Navy through the NICE-SHIP project. We also acknowledge the financial support of the CERCA programme of the Generalitat de Catalunya.

==References==

===BIBLIOGRAPHY===

<div id="cite-1"></div>
'''[[#citeF-1|[1]]]'''  Baiocchi C, Brezzi F, Franca LP. Virtual bubbles and the GaLS. ''Computer Methods in Applied Mechanics and Engineering'' 1993; '''105''':125&#8211;141.

<div id="cite-2"></div>
'''[[#citeF-2|[2]]]'''  Brezzi F, Bristeau MO, Franca LP, Mallet M, Rogé G. A relationship between stabilized finite element methods and the Galerkin method with bubble functions. ''Computer Methods in Applied Mechanics and Engineering'' 1992;  '''96''':117&#8211;129.

<div id="cite-3"></div>
'''[[#citeF-3|[3]]]'''  Brezzi F, Russo A. Choosing bubbles for advection-diffusion problems. ''Mathematical Models and Methods in Applied Sciences'' 1994;  '''4''':571&#8211;587.

<div id="cite-4"></div>
'''[[#citeF-4|[4]]]'''  Brezzi F, Hauke G, Marini LD, Sangalli G. Link-cutting bubbles for the stabilization of convection&#8211;diffusion&#8211;reaction problems. ''Mathematical Models and Methods in Applied Sciences'' 2003;  '''13'''(3):445&#8211;461.

<div id="cite-5"></div>
'''[[#citeF-5|[5]]]'''  Brooks AN, Hughes TJR. Streamline Upwind/Petrov&#8211;Galerkin formulations for the convective dominated flows with particular emphasis on the incompressible Navier-Stokes equations. ''Computer Methods in Applied Mechanics and Engineering'' 1982;  '''32''' (1&#8211;3):199&#8211;259.

<div id="cite-6"></div>
'''[[#citeF-6|[6]]]'''  Burman E, Ern A. Nonlinear diffusion and discrete maximum principle for stabilized Galerkin approximations of the convection&#8211;diffusion&#8211;reaction equation,'' ''Computer Methods in Applied Mechanics and Engineering'' 2002,  191:3833-3855.

<div id="cite-7"></div>
'''[[#citeF-7|[7]]]'''  Codina R. A discontinuity-capturing crosswind-dissipation for the finite element solution of the convection-diffusion equation,'' ''Computer Methods in Applied Mechanics and Engineering'' 1993,  110:325-342.

<div id="cite-8"></div>
'''[[#citeF-8|[8]]]'''  Codina R. Comparison of some finite element methods for solving the diffusion-convection-reaction equation, ''Computer Methods in Applied Mechanics and Engineering'' 1998,  156:185-210.

<div id="cite-9"></div>
'''[[#citeF-9|[9]]]'''  Codina R. On stabilized finite element methods for linear systems of convection&#8211;diffusion&#8211;reaction equations, ''Computer Methods in Applied Mechanics and Engineering'' 2000a,  190:2681-2706.

<div id="cite-10"></div>
'''[[#citeF-10|[10]]]'''  de Sampaio PAB, Coutinho ALGA. A natural derivation of discontinuity capturing operator for convection-diffusion problems, ''Computer Methods in Applied Mechanics and Engineering'' 2001,  190:6291-6308.

<div id="cite-11"></div>
'''[[#citeF-11|[11]]]'''  Donea J. A Taylor-Galerkin method for convective transport problems, ''International Journal for Numerical Methods in Engineering'' 1984,  20:101-119.

<div id="cite-12"></div>
'''[[#citeF-12|[12]]]'''   J. Donea,  A. Huerta,   ''Finite element method for flow problems''. J. Wiley & Sons, 2003.

<div id="cite-13"></div>
'''[[#citeF-13|[13]]]'''  Douglas J, Russell TF. Numerical methods for convection-dominated diffusion problems based on combining the method of characteristics with the finite element or finite difference procedures, ''SIAM Journal of Numerical Analysis'' 1982,  19:871-885.

<div id="cite-14"></div>
'''[[#citeF-14|[14]]]'''  Dutra do Carmo EG, Alvarez GB. A new upwind function in stabilized finite element formulations, using linear and quadratic elements for the scalar convection-diffusion problems, ''Computer Methods in Applied Mechanics and Engineering'' 2004,  193:2383-2402.

<div id="cite-15"></div>
'''[[#citeF-15|[15]]]'''  Felippa C, Oñate E. Nodally exact Ritz discretizations of the 1D diffusion-absorption and Helmholtz equations by variational FIC and modified equation methods, ''Computational Mechanics'' 2007,  39(2):91-111.

<div id="cite-16"></div>
'''[[#citeF-16|[16]]]'''  Fisher H, List J, Koh C, Imberger J, Brooks N. Mixing in Inland and Coastal Waters, Elsevier, 2013.

<div id="cite-17"></div>
'''[[#citeF-17|[17]]]'''  Franca LP, Dutra do Carmo EG. The Galerkin gradient least-squares method, ''Computer Methods in Applied Mechanics and Engineering'' 1989,  74:41-54.

<div id="cite-18"></div>
'''[[#citeF-18|[18]]]'''  Franca LP, Farhat C. Bubble functions prompt unusual stabilized finite element methods, ''Computer Methods in Applied Mechanics and Engineering'' 1995,  123:299-308.

<div id="cite-19"></div>
'''[[#citeF-19|[19]]]'''  Franca LP, Valentin F. On an improved unusual stabilized finite element method for the advective-reactive-diffusive equation, ''Computer Methods in Applied Mechanics and Engineering'' 2001,  190:1785-1800.

<div id="cite-20"></div>
'''[[#citeF-20|[20]]]'''  Galeo AC, Dutra do Carmo EG. A consistent approximate upwind Petrov&#8211;Galerkin method for the convection-dominated problems. ''Computer Methods in Applied Mechanics and Engineering'' 1988,  68(1):83-95.

<div id="cite-21"></div>
'''[[#citeF-21|[21]]]'''  Harari I, Hughes TJR. Stabilized finite element methods for steady advection-diffusion with production, ''Computer Methods in Applied Mechanics and Engineering'' 1994,  115:165-191.

<div id="cite-22"></div>
'''[[#citeF-22|[22]]]'''  Hauke G, García-Olivares A. Variational subgrid scale formulations for the advection-diffusion&#8211;reaction equation, ''Computer Methods in Applied Mechanics and Engineering'' 2001,  190:6847-6865.

<div id="cite-23"></div>
'''[[#citeF-23|[23]]]'''  Hauke G. A simple subgrid scale stabilized method for the advection-diffusion&#8211;reaction equation, ''Computer Methods in Applied Mechanics and Engineering'' 2002,  191:2925-2947.

<div id="cite-24"></div>
'''[[#citeF-24|[24]]]'''  Hauke G, Sangalli G, Doweidar MH. Combining adjoint stabilized methods for the advection-diffusion&#8211;reaction problem, ''Mathematical Models and Methods in Applied Sciences'' 2007,  17(2):305-326.

<div id="cite-25"></div>
'''[[#citeF-25|[25]]]'''  Hughes TJR, Brooks AN. A theoretical framework for Petrov&#8211;Galerkin methods, with discontinuous weighting functions: application to the streamline upwind procedure, in Gallagher RH, Norrie DM, Oden JT, Zienkiewicz OC (eds.), ''Finite Elements in Fluids'' 1982,  IV,Wiley, Chichester.

<div id="cite-26"></div>
'''[[#citeF-26|[26]]]'''  Hughes TJR,  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.'' 1986a,  58:305&#8211;328.

<div id="cite-27"></div>
'''[[#citeF-27|[27]]]'''  Hughes TJR, Franca LP, Hulbert GM. A new finite element formulation for computational fluid dynamics: VIII. The Galerkin/least-squares method for advective-diffusive equations, ''Computer Methods in Applied Mechanics and Engineering'' 1989,  73:173-189.

<div id="cite-28"></div>
'''[[#citeF-28|[28]]]'''  Hughes TJR, Feijoo GR, Mazzei L, Quincy JB. The variational multiscale method: a paradigm for computational mechanics, ''Computer Methods in Applied Mechanics and Engineering'' 1998,  166:3-24.

<div id="cite-29"></div>
'''[[#citeF-29|[29]]]'''  Idelsohn SR, Nigro N, Storti M, Buscaglia G. A Petrov&#8211;Galerkin formulation for advective-reaction-diffusion problems, ''Computer Methods in Applied Mechanics and Engineering'' 1996,  136:27-46.

<div id="cite-30"></div>
'''[[#citeF-30|[30]]]'''  Idelsohn SR, Heinrich JC, Oñate E. Petrov&#8211;Galerkin methods for the transient advective-diffusive equation with sharp gradients, ''International Journal for Numerical Methods in Engineering'' 1996,  39:1455-1473.

<div id="cite-31"></div>
'''[[#citeF-31|[31]]]'''  John V, Knobloch P. On spurious oscillations at layers diminishing (SOLD) methods for convection-diffusion equations: Part I - A review, ''Computer Methods in Applied Mechanics and Engineering'' 2007,  196:2197&#8211;2215.

<div id="cite-32"></div>
'''[[#citeF-32|[32]]]'''  Johnson C, Szepessy A, Hansbo P. On the convergence of shock-capturing streamline diffusion finite element methods for hyperbolic conservation laws, ''Mathematics of Computation'' 1990,  54:107-129.

<div id="cite-33"></div>
'''[[#citeF-33|[33]]]'''  Kikuchi F, Ushijima T. Theoretical analysis of some finite element schemes   for convective diffusion equations, in: R.&nbsp;Gallagher, D.&nbsp;Norrie, J.&nbsp;Oden,   O.&nbsp;C. Zienkiewicz (Eds.), Finite Elements in Fluids, Vol.&nbsp;IV, John Wiley and   Sons Ltd, Chichester, 1982.

<div id="cite-34"></div>
'''[[#citeF-34|[34]]]'''  Knobloch P. Improvements of the Mizukami-Hughes method for convection-diffusion equations, ''Computer Methods in Applied Mechanics and Engineering'' 2006,  196:579-594.

<div id="cite-35"></div>
'''[[#citeF-35|[35]]]'''  Löhner R, Morgan K, Zienkiewicz OC. The solution of non-linear hyperbolic equation systems by the finite element method, ''International Journal for Numerical Methods in Fluids'' 1984,  4:1043-1063.

<div id="cite-36"></div>
'''[[#citeF-36|[36]]]'''  Mizukami A, Hughes TJR. A Petrov&#8211;Galerkin finite element method for convection-dominated flows: An accurate upwinding technique for satisfying the maximum principle, ''Computer Methods in Applied Mechanics and Engineering'' 1985,  50(2):181-193.

<div id="cite-37"></div>
'''[[#citeF-37|[37]]]'''  Nadukandi P, Oñate E, García Espinosa J. A high-resolution Petrov-Galerkin method for the 1D convection-diffusion-reaction problem, ''Computer Methods in Applied Mechanics and Engineering'' 2010,  199(9&#8211;12):525-546.

<div id="cite-38"></div>
'''[[#citeF-38|[38]]]'''  Nadukandi P, Oñate E, García Espinosa J. A high-resolution Petrov-Galerkin method for the convection-diffusion-reaction problem. Part II. A multidimensional extension, ''Computer Methods in Applied Mechanics and Engineering'' 2012,  213-216:327-352.

<div id="cite-39"></div>
'''[[#citeF-39|[39]]]'''  Oñate E, Idelsohn SR, Zienkiewicz OC, TaylorRL, Sacco C. A stabilized finite point method for analysis of fluid mechanics problems. ''Computer Methods in Applied Mechanics and Engineering'' 1996,  139:315&#8211;346.

<div id="cite-40"></div>
'''[[#citeF-40|[40]]]'''  Oñate E. Derivation of stabilized equations for numerical solution of advective-diffusive transport and fluid flow problems, ''Computer Methods in Applied Mechanics and Engineering'' 1998,  151:233-265.

<div id="cite-41"></div>
'''[[#citeF-41|[41]]]'''  Oñate E, Idelsohn SR. A mesh-free finite point method for advective-diffusive transport and fluid flow problems. ''Comput. Mech.'' 1998,  23:283&#8211;292.

<div id="cite-42"></div>
'''[[#citeF-42|[42]]]''' Oñate E, Manzan M. A general procedure for deriving stabilized space-time finite element methods for advective-diffusive problems. ''Int. J. Num. Meth. Fluids'' 1999,  31:203&#8211;221.

<div id="cite-43"></div>
'''[[#citeF-43|[43]]]'''  Oñate E, Manzan M. Stabilization techniques for finite element analysis of convection-diffusion problems, in Sundén B, Comini G. (eds.), ''Computational Analysis of Convection Heat Transfer'' 2000, 71-117, WIT Press, Southampton (UK), 2000.

<div id="cite-44"></div>
'''[[#citeF-44|[44]]]'''  Oñate E, Sacco C, Idelsohn SR. A finite point method for incompressible flow problems. ''Computing and Visualization in Science'' 2000,  2:67&#8211;75.

<div id="cite-45"></div>
'''[[#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. ''Computer Methods in Applied Mechanics and Engineering'' 2001,  191:635&#8211;660.

<div id="cite-46"></div>
'''[[#citeF-46|[46]]]'''  Oñate E, Perazzo F, Miquel J. A finite point method for elasticity problems. ''Computers and Structures'' 2001,  79:2151&#8211;2163.

<div id="cite-47"></div>
'''[[#citeF-47|[47]]]'''   Oñate E,  Taylor RL, Zienkiewicz OC, Rojek J. A residual correction method based on finite calculus. ''Engineering Computations'' 2003,   20(5/6):629&#8211;658.

<div id="cite-48"></div>
'''[[#citeF-48|[48]]]'''  Oñate E. Possibilities of finite calculus in computational mechanics.  ''International Journal for Numerical Methods in Engineering'' 2004,  60(1), 255&#8211;281.

<div id="cite-49"></div>
'''[[#citeF-49|[49]]]'''  Oñate E, Idelsohn SR, Del Pin F, Aubry R. The particle finite element method. An overview. ''International Journal of Computational Methods'' 2004,  1(2):267&#8211;307.

<div id="cite-50"></div>
'''[[#citeF-50|[50]]]'''   Oñate E,  Rojek J,  Taylor RL  Zienkiewicz OC. Finite calculus formulation for incompressible solids using linear triangles and tetrahedra. ''International Journal for Numerical Methods in  Engineering'' 2004,   59:1473&#8211;1500.

<div id="cite-51"></div>
'''[[#citeF-51|[51]]]'''  Oñate E, Zarate F, Idelsohn SR. Finite element formulation for the convective-diffusive problems with sharp gradients using finite calculus, ''Computer Methods in Applied Mechanics and Engineering'' 2006,  195:1793-1825.

<div id="cite-52"></div>
'''[[#citeF-52|[52]]]'''  Oñate E, Miquel J, Hauke G. Stabilized formulation for the advection-diffusion-absorption equation using finite calculus and linear finite elements, ''Computer Methods in Applied Mechanics and Engineering'' 2006,  195(33&#8211;36):3926&#8211;3946.

<div id="cite-53"></div>
'''[[#citeF-53|[53]]]'''  Oñate E, Valls A, García J. FIC/FEM formulation with matrix   stabilizing terms for incompressible flows at low and high Reynolds   numbers. ''Computational Mechanics'' 2006,  38(4-5):440&#8211;455.

<div id="cite-54"></div>
'''[[#citeF-54|[54]]]'''  Oñate E,  García J,  Idelsohn SR, Del Pin F. FIC   formulations for finite element analysis of incompressible flows. Eulerian,   ALE and Lagrangian approaches.  ''Computer Methods in Applied Mechanics and   Engineering'' 2006,  195(23-24):3001&#8211;3037.

<div id="cite-55"></div>
'''[[#citeF-55|[55]]]'''  Oñate E, Miquel J, Zarate F. Stabilized solution of the multidimensional advection-diffusion-absorption equation using linear finite elements, ''Computers and Fluids'' 2007,  36:92&#8211;112.

<div id="cite-56"></div>
'''[[#citeF-56|[56]]]'''  Oñate E, Valls A, García J. Modeling incompressible flows at low and high Reynolds numbers via a finite calculus-finite element approach. ''Journal of Computational Physics'' 2007,  224:332&#8211;351.

<div id="cite-57"></div>
'''[[#citeF-57|[57]]]'''  Oñate E, Idelsohn SR, Celigueta MA, Rossi R. Advances in the particle finite element method for the analysis of fluid-multibody interaction and bed erosion in free surface flows. ''Comput. Methods Appl. Mech. Engrg.'' 2008,  197(19-20):1777–-1800.

<div id="cite-58"></div>
'''[[#citeF-58|[58]]]'''  Oñate E, Nadukandi P, Idelsohn S, García J, Felippa C. A family of residual-based stabilized finite element methods for Stokes flows. ''Int. J. Num. Meth. in Fluids'' 2011,  65 (1-3):106&#8211;134.

<div id="cite-59"></div>
'''[[#citeF-59|[59]]]'''   Oñate E, Idelsohn SR, Felippa C. Consistent pressure Laplacian stabilization for incompressible continua via higher order finite calculus. ''International Journal for Numerical Methods in  Engineering'' 2011,  87 (1-5):171–-195.

<div id="cite-60"></div>
'''[[#citeF-60|[60]]]'''  Oñate E, Nadukandi P, Idelsohn S. P1/P0+ elements for incompressible  flows  with discontinuous material properties. ''Comput. Methods Appl. Mech. Engrg.'' 2014,  271:185&#8211;209

<div id="cite-61"></div>
'''[[#citeF-61|[61]]]'''  Oñate E, Miquel J, Nadukandi P.  An accurate FIC-FEM formulation for the 1D advection–diffusion–reaction equation. ''Comput. Methods Appl. Mech. Engrg.'' 2016, 298:373–406.

<div id="cite-62"></div>
'''[[#citeF-62|[62]]]'''  Oñate E, Nadukandi P, Miquel J. Accurate FIC-FEM formulation for the multidimensional steady-state advection–diffusion–absorption equation. ''Comput. Methods Appl. Mech. Engrg.'' 2017, 327:352-368.

<div id="cite-63"></div>
'''[[#citeF-63|[63]]]'''  J.&nbsp;Principe, R.&nbsp;Codina,   [ On the   stabilization parameter in the subgrid scale approximation of scalar   convection–diffusion–reaction equations on distorted meshes], Computer   Methods in Applied Mechanics and Engineering 2010, 199(21-22):1386&#8211;1402.

<div id="cite-64"></div>
'''[[#citeF-64|[64]]]'''  Tezduyar TE, Park YJ. Discontinuity-capturing finite element formulations for nonlinear convection&#8211;diffusion&#8211;reaction equations. ''Computer Methods in Applied Mechanics and Engineering'' 1986; '''59''':307&#8211;325.

<div id="cite-65"></div>
'''[[#citeF-65|[65]]]'''  Yu CC, Heinrich JC. Petrov&#8211;Galerkin methods for the time-dependent convective transport equations, ''International Journal for Numerical Methods in Engineering'' 1986,  23:883-901.

<div id="cite-66"></div>
'''[[#citeF-66|[66]]]'''  Zienkiewicz OC, Codina R. A general algorithm for compressible and incompressible flows. Part I: the split, characteristic based scheme, ''International Journal for Numerical Methods in Fluids'' 1995,  20:869-885.

<div id="cite-67"></div>
'''[[#citeF-67|[67]]]'''  Zienkiewicz OC,  Taylor RL, Zhu JZ. ''The finite element method. The basis''.  6th Ed., Elsevier, 2005.

<div id="cite-68"></div>
'''[[#citeF-68|[68]]]'''  Zienkiewicz OC,  Taylor RL, Nithiarasu P.   ''The finite element   method for fluid dynamics''.  6th Ed.,  Elsevier, 2005.