Revision as of 10:55, 19 March 2020 (view source) Cinmemj (talk | contribs) ← Older edit		Latest revision as of 13:47, 2 June 2020 (view source) Cinmemj (talk | contribs)
Line 1:		Line 1:
−	In ''Computer Methods in Applied Mechanics and Engineering'', Article in Press, 2020<br />	+	In ''Computer Methods in Applied Mechanics and Engineering'', Vol. 365, pp. 112984, 2020<br />
	DOI: 10.1016/j.cma.2020.112984		DOI: 10.1016/j.cma.2020.112984
	<!-- metadata commented in wiki content		<!-- metadata commented in wiki content

---

Latest revision as of 13:47, 2 June 2020

In Computer Methods in Applied Mechanics and Engineering, Vol. 365, pp. 112984, 2020
DOI: 10.1016/j.cma.2020.112984

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 [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 [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 [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) [27], Variational Multiscale (VMS) [28] and Characteristic split procedure [66,68]. Many of these methods are reviewed in [31].

Control over the dispersive oscillations for transient advection-diffusion problems via SUPG methods and space-time FEM were reported in [30] and [65], respectively.

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

Oñate et al. [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 [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 [4], SUPG [21] and VMS [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 [37,38].

Starting from a different perspective, but still within the FIC framework, Oñate, Miquel and Nadukandi [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 [61] for 1D stationary problems solved with irregular meshes, as well as for 1D transient problems. The FIC-FEM formulation presented in [61] was extended by the same authors in [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 [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 [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 ${\textstyle \Omega }$ can be expressed as

${\displaystyle r_{t}=0\quad {\hbox{in }}\Omega }$

(1.a)

with

${\displaystyle 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}$

(1.b)

For 3D problems,

${\displaystyle \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}}$

(2)

In Eqs.(1.a)–(1.b) ${\textstyle \phi }$ is the transported variable (i.e., the temperature in a heat transfer problem or the concentration in a mass transfer problem), ${\textstyle v_{i}}$ is the ${\textstyle i}$th component of the velocity vector ${\textstyle {\mathbf {} v}}$; ${\textstyle \rho ,c}$ and ${\textstyle k_{i}}$ are the density, the specific flux parameter and the conductivity of the material along the ${\textstyle i}$th global direction, respectively and ${\textstyle s}$ 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 (${\textstyle \rho ,c,k,s}$) are constant over the analysis domain ${\textstyle \Omega }$.

Boundary conditions

${\displaystyle \phi -\phi ^{p}=0\quad {\hbox{on }}\Gamma _{\phi }}$

(3)

${\displaystyle r_{\scriptscriptstyle \Gamma }=0\quad {\hbox{on }}\Gamma _{q}}$

(4)

with

${\displaystyle r_{\Gamma }:=-q_{n}+q_{n}^{p}}$

(5)

where

${\displaystyle q_{n}={\boldsymbol {q}}^{T}{\boldsymbol {n}}\quad ,\quad {\boldsymbol {q}}=-{\boldsymbol {D}}{\boldsymbol {\nabla }}\phi }$

(6)

In Eqs.(3)–6 ${\textstyle \phi ^{p}}$ and ${\textstyle q^{p}}$ are the prescribed values of the transported variable and the outgoing diffusive flux at the Dirichlet and Neumann boundaries ${\textstyle \Gamma _{\phi }}$ and ${\textstyle \Gamma _{q}}$, respectively, with ${\textstyle \Gamma _{\phi }\cup \Gamma _{q}=\Gamma }$, ${\textstyle \Gamma }$ 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

${\displaystyle \phi ({\boldsymbol {x}},t_{0})=\phi _{0}({\boldsymbol {x}})}$

(7)

where ${\textstyle \phi _{0}({\boldsymbol {x}})}$ is the value of the transported variable at time ${\textstyle t=t_{0}}$.

In this work we will consider cases for which ${\textstyle s\geq 0}$ only. This includes the following particular problems:

• Advection-diffusion-absorption (${\textstyle |{\boldsymbol {v}}|\not =0}$, ${\textstyle K\not =0}$, ${\textstyle s>0}$).
• Advection-diffusion (${\textstyle |{\boldsymbol {v}}|\not =0}$, ${\textstyle K\not =0}$, ${\textstyle s=0}$).
• Diffusion-absorption (${\textstyle |{\boldsymbol {v}}|=0}$, ${\textstyle K\not =0}$, ${\textstyle s>0}$).
• Advection-absorption (${\textstyle |{\boldsymbol {v}}|\not =0}$, ${\textstyle K=0}$, ${\textstyle s>0}$).

In the above ${\textstyle K}$ is the average difussion given by ${\textstyle K=\left[{\frac {1}{n_{d}}}\sum \limits _{i=1}^{n_{d}}(k_{i})^{2}\right]^{1/2}}$, where ${\textstyle n_{d}}$ is the number of space dimensions (i.e. ${\textstyle n_{d}=2}$ for 2D problems).

2.2 Finite increment calculus (FIC) expressions

The governing equations 1.a and 1.b and the boundary conditions 3-6 are expressed using the FIC theory as [55].

Transport balance

${\displaystyle 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 }$

(8)

with ${\textstyle {\boldsymbol {h}}=[h_{1},h_{2},h_{3}]^{T}}$ in 3D.

Boundary conditions

${\displaystyle \phi -\phi ^{p}=0\quad {\hbox{on }}\Gamma _{\phi }}$

(9.a)

${\displaystyle 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}}}$

(9.b)

Eqs.(8) and (9.b) are obtained by expressing the balance of fluxes in an arbitrary prismatic space domain of size ${\textstyle h_{1}\times h_{2}\times h_{3}}$ within the global problem domain and at the Neumann boundary, respectively. The distances ${\textstyle h_{i}}$ 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 [40]. This higher order expansions lead to the underlined terms in Eqs.(8) and (2.2). These terms lead naturally to stabilized numerical schemes.

Note that as the characteristic length vector ${\textstyle {\boldsymbol {h}}}$ tends to zero the FIC governing equations gradually recover the standard infinitesimal form, giving in the limit (for ${\textstyle \mathbf {h} =\mathbf {0} }$) ${\textstyle r_{t}=0}$ in ${\textstyle \Omega }$ and ${\textstyle r_{\scriptscriptstyle \Gamma }=0\quad {\hbox{on }}\Gamma _{q}}$.

The stability and accuracy of the numerical solution depends on the components of the characteristic length vector ${\textstyle {\boldsymbol {h}}}$. At the discretization level the length of ${\textstyle {\boldsymbol {h}}}$ is usually expressed as a proportion of a typical grid dimension [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 [40,42,43,51], diffusion-absorption and Helmholtz [15], advection-diffusion-absorption [52,55], advection-diffusion-reaction [61], incompressible fluid flow [53,54,56,58,60], fluid-structure-interaction [45,49,57], particle-laden flows and standard and incompressible solid mechanics [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 [39,41,44,46].

3. Definition of the characteristic length vector

The characteristic length vector ${\textstyle {\boldsymbol {h}}}$ is designed so that the expression for all the stabilization matrices and vectors reduce to those given for the 1D case as reported in [61].

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

${\displaystyle {\boldsymbol {h}}={\boldsymbol {h}}_{v}+{\boldsymbol {h}}_{r}+{\boldsymbol {h}}_{sc}}$

(10)

The characteristic length vectors in the r.h.s. of Eq.10 are defined as follows.

Streamline characteristic length vector, h_v

In Eq.10 ${\textstyle {\boldsymbol {h}}_{v}}$ is a length vector along the velocity direction defined as

${\displaystyle {\boldsymbol {h}}_{v}=\alpha _{v}l_{v}{\frac {\mathbf {v} }{|\mathbf {v} |}}=\alpha _{v}l_{v}{\hat {\boldsymbol {v}}}}$

(11)

where ${\textstyle \alpha _{v}}$ is a streamline stabilization parameter, ${\textstyle l_{v}}$ is a characteristic element dimension and ${\textstyle {\hat {\boldsymbol {v}}}}$ is a unit velocity vector.

Absorption characteristic length vector, h_r

In Eq.10 ${\textstyle {\boldsymbol {h}}_{r}}$ is a characteristic length vector induced by absorption effects and defined as

${\displaystyle {\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]}$

(12.a)

where ${\textstyle \alpha _{r}}$ is an absorption stabilization parameter, ${\textstyle {D}={\hat {\boldsymbol {v}}}^{T}{\boldsymbol {D}}{\hat {\boldsymbol {v}}}}$ and ${\textstyle r_{s}}$ is the space residual defined as

${\displaystyle r_{s}:=r_{t}-\rho c{\partial \phi \over \partial t}}$

(12.b)

Matrix ${\textstyle {\boldsymbol {D}}_{s}}$ in Eq.12.a is defined for different element types as follows.

3-noded triangles and 4-noded tetrahedra:

${\displaystyle {\boldsymbol {D}}_{s}={\frac {s}{(n+1)}}\sum \limits _{i=1}^{n}{\boldsymbol {l}}_{i}{\boldsymbol {l}}_{i}^{T}}$

(13)

where ${\textstyle {\boldsymbol {l}}_{i}}$ is the vector joining the baricenter of the element and the ${\textstyle i}$th node, and ${\textstyle n}$ is the number of nodes of the element.

The diffusion introduced by matrix ${\textstyle {\boldsymbol {D}}_{s}}$ takes care of the instabilities induced by the irregularity of the triangular mesh near boundaries that develop parabolic layes [62].

Any other element: Dₛ=[0]

Shock capturing characteristic length vector h_sc

In Eq.10 ${\textstyle {\boldsymbol {h}}_{sc}}$ 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

${\displaystyle {\boldsymbol {h}}_{sc}={h}_{sc}{\widehat {{\boldsymbol {\nabla }}\phi }}}$

(14.a)

where ${\textstyle {\widehat {{\boldsymbol {\nabla }}\phi }}={\frac {{\boldsymbol {\nabla }}\phi }{|{\boldsymbol {\nabla }}\phi |}}}$ is the unit gradient vector, and

${\displaystyle {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]}$

(14.b)

In Eq.14.b ${\textstyle \beta }$ is a non linear parameter that depends on the angle ${\textstyle \theta }$ between the velocity vector and the gradient vector. It is defined as

${\displaystyle \beta =\left\{{\begin{array}{lll}1&{\hbox{ if }}&\theta <\theta _{c}\\\displaystyle {\hat {\boldsymbol {v}}}^{T}{\widehat {{\boldsymbol {\nabla }}\phi }}&{\hbox{ if }}&\theta \geq \theta _{c}\end{array}}\right.}$

(14.c)

where ${\textstyle \theta _{c}}$ is a critical angle. In this work we have taken ${\textstyle \theta _{c}=20^{\circ }}$.

The parameter ${\textstyle \beta }$ controls the amount of shock-capturing nonlinear diffusion active at any point of the domain [62]. When the gradient vector is parallel to the velocity vector then ${\textstyle \beta =1}$, ${\textstyle h_{sc}}$ vanishes and the linear stabilization terms suffice to diminish spurious numerical oscillations about the layers. The term ${\textstyle 1-\beta ^{2}}$ in Eq.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 8 and 2.2 is written as

${\displaystyle \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}$

(15)

where ${\textstyle W}$ are arbitrary weighting functions.

Integrating by parts the FIC term in the first integral of 15 gives

${\displaystyle \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}$

(16)

Note that the FIC term has vanished from the boundary integral, as it usual in the FIC-FEM approach [40,48].

Let us substitute the expression for the characteristic vector ${\textstyle {\boldsymbol {h}}}$ of Eq.10 into Eq.16. This gives (using Eqs.11, 12.a and 14.a)

${\displaystyle \int _{\Omega }\left[Wr_{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}$

(17)

The final step is the integration by parts of the diffusive term in the expression of ${\textstyle r_{t}}$ in the first term of the first integral of Eq.17. This gives, after replacing the definition of ${\textstyle {\boldsymbol {H}}_{r}}$ in Eq.17 and grouping some terms, the following expression for the weak variational form of the FIC governing equations

${\displaystyle \displaystyle \int _{\Omega }\left[\rho c{\bar {W}}{\frac {\partial \phi }{\partial t}}+\rho cW{\boldsymbol {v}}^{T}{\boldsymbol {\nabla }}\phi +({\boldsymbol {\nabla }}^{T}W){\boldsymbol {D}}_{T}{\boldsymbol {\nabla }}\phi +Ws\phi +{\frac {1}{2}}({\boldsymbol {\nabla }}^{T}W){\boldsymbol {h}}_{v}s\phi \right]d\Omega +}$ ${\displaystyle \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 }WQd\Omega +\oint _{\Gamma _{q}}Wq_{n}^{p}d\Gamma =0}$

(18)

The space weighting function ${\textstyle {\bar {W}}}$ is given by

${\displaystyle {\bar {W}}=W+{\frac {1}{2}}{\boldsymbol {h}}_{v}^{T}{\boldsymbol {\nabla }}W}$

(19)

The expression of the total diffusivity matrix ${\textstyle {\boldsymbol {D}}_{T}}$ in Eq.18 is

${\displaystyle \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}}}$

(20)

where ${\textstyle {\boldsymbol {D}}}$ and ${\textstyle {\boldsymbol {D}}_{s}}$ are defined in Eqs.2 and 13, respectively, I is the unit matrix and

${\displaystyle {\boldsymbol {D}}_{v}\!\!\!\!=\!\!{\frac {l_{v}}{2}}{\hat {\boldsymbol {v}}}{\boldsymbol {v}}^{T}}$ (21.a) ${\displaystyle {\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})}$ (21.b)

4.2 FIC-FEM equations

We interpolate the transported variable ${\textstyle \phi }$ in the standard FEM fashion over a mesh of elements with ${\textstyle n}$ nodes [67,68] as

${\displaystyle \phi \simeq {\hat {\phi }}=\sum \limits _{i=1}^{n}N_{i}\phi _{i}}$

(22)

where ${\textstyle N_{i}}$ are the space shape functions and ${\textstyle \phi _{i}}$ the nodal variables.

Introducing Eq.22 into the FIC variational form 18 and using a Galerkin approach (${\textstyle W_{i}=N_{i}}$) gives the final system of discretized equations as

${\displaystyle \mathbf {M} {\dot {\boldsymbol {\phi }}}+[\mathbf {K} +\mathbf {C} +\mathbf {S} ]{\boldsymbol {\phi }}={\boldsymbol {f}}}$

(23)

where ${\textstyle {\boldsymbol {\phi }}=[\phi _{1},\phi _{2},\cdots ,\phi _{N}]^{T}}$ is the vector of nodal unknowns, with ${\textstyle N}$ being total number of nodes in the mesh.

The rest of matrices and vector ${\textstyle {\boldsymbol {f}}}$ in Eq.23 are obtained in the standard FEM fashion by assembling the element contribution given by

${\displaystyle \mathbf {M} _{ij}^{e}=\int _{\Omega ^{e}}\rho cN_{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 }$ (24) ${\displaystyle \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 }$ (25) ${\displaystyle \mathbf {C} _{ij}^{e}=\int _{\Omega ^{e}}\rho cN_{i}{\boldsymbol {v}}^{T}{\boldsymbol {\nabla }}N_{j}d\Omega }$ (26) ${\displaystyle \mathbf {S} _{ij}^{e}=\int _{\Omega ^{e}}sN_{i}N_{j}d\Omega +\int _{\Omega ^{e}}{\frac {1}{2}}({\boldsymbol {\nabla }}^{T}N_{i}){\boldsymbol {h}}_{v}sN_{j}d\Omega }$ (27) ${\displaystyle {\boldsymbol {f}}_{i}^{e}=\int _{\Omega ^{e}}\left(N_{i}+{\frac {1}{2}}({\boldsymbol {\nabla }}^{T}N_{i}){\boldsymbol {h}}\right)Qd\Omega -\oint _{\Gamma _{q}^{e}}N_{i}q_{n}^{p}d\Gamma }$ (28)

We note that the second integral in Eq.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 ${\textstyle {\boldsymbol {M}}}$. An exception is the problem of the evolution of a concentration field presented in Section 7.6 for which using the consistent form of matrix ${\textstyle {\boldsymbol {M}}}$ is essential to reduce excessive diffusion.

5. Transient solution scheme

5.1 Implicit time integration scheme

We discretize in time the system of equations (23) using a Generalized Trapezoidal rule [11,67]. The solution for the nodal values at a time instant is found using an incremental iterative strategy as

${\displaystyle {}^{i}{\boldsymbol {H}}^{n}\Delta {\boldsymbol {\phi }}=-{}^{i}{\boldsymbol {r}}_{t}^{n}}$

(29)

where ${\textstyle \Delta {\boldsymbol {\phi }}}$ is the increment of the nodal variables, ${\textstyle \theta }$ is a non dimensional time parameter (${\textstyle 0.5<\theta \leq 1}$ is required for the integration scheme to be unconditionally stable [11,67,68]), ${\textstyle (\cdot )^{n}}$ denotes values at time ${\textstyle t=t_{n}}$, ${\textstyle {}^{i}(\cdot )}$ denotes values at the ${\textstyle i}$th iteration and

${\displaystyle {}^{i}{\boldsymbol {H}}^{n}={\frac {1}{\theta \Delta t}}{\boldsymbol {M}}+{}^{i}{\boldsymbol {K}}^{n}+{\boldsymbol {C}}+{\boldsymbol {S}}}$

(30)

${\displaystyle {}^{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}}$

(31)

In Eq.31 we define ${\textstyle {\dot {\boldsymbol {\phi }}}={\frac {{\boldsymbol {\phi }}^{n+\theta }-{\boldsymbol {\phi }}^{n}}{\theta \Delta t}}}$.

From the value of ${\textstyle \Delta {\boldsymbol {\phi }}}$ obtained from solving Eq.29 we compute the value of ${\textstyle {\boldsymbol {\phi }}^{n+\theta }}$ at the ${\textstyle i+1}$ iteration as

${\displaystyle {}^{i+1}{\boldsymbol {\phi }}^{n+\theta }={}^{i}{\boldsymbol {\phi }}^{n+\theta }+\Delta {\boldsymbol {\phi }}}$

(32)

The iterative solution at ${\textstyle t_{n+1}}$, ${\textstyle {}^{i+1}{\boldsymbol {\phi }}^{n+1}}$, is obtained as

${\displaystyle {}^{i+1}{\boldsymbol {\phi }}^{n+1}={\frac {1}{\theta }}{}^{i+1}{\boldsymbol {\phi }}^{n+\theta }+\left(1-{\frac {1}{\theta }}\right){\boldsymbol {\phi }}^{n}}$

(33)

The non-linearity in the expression of matrix ${\textstyle {\boldsymbol {H}}}$ in Eq.30 is due to the dependence of the shock-capturing term in matrix ${\textstyle {\boldsymbol {K}}}$ via matrix ${\textstyle {\boldsymbol {D}}_{sc}}$ (see Eq.21.b).

The iterations proceed until convergence is achieved for both the unknown ${\textstyle \phi }$ and the residual measured in a ${\textstyle L_{2}}$ norm. In the transient problems solved in this work, convergence within each time step was typically achieved in 2–3 iterations.

5.2 Explicit time integration scheme

If we rewrite the residual in Eq.31 choosing ${\textstyle \theta =0}$ and defining now ${\textstyle {\dot {\boldsymbol {\phi }}}}$ as ${\textstyle {\frac {{\boldsymbol {\phi }}^{n+1}-{\boldsymbol {\phi }}^{n}}{\Delta t}}}$ leads to the following explicit solution scheme

${\displaystyle {\boldsymbol {\phi }}^{n+1}={\boldsymbol {\phi }}^{n}+\Delta t{\boldsymbol {M}}_{d}^{-1}\left[{\boldsymbol {f}}^{n}-({\boldsymbol {K}}^{n}+{\boldsymbol {C}}-{\boldsymbol {S}}){\boldsymbol {\phi }}^{n}\right]}$

(34)

where ${\textstyle {\boldsymbol {M}}_{d}}$ denotes the lumped diagonal form of matrix ${\textstyle {\boldsymbol {M}}}$.

In the expression of ${\textstyle {\boldsymbol {K}}}$ of Eq.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 ${\textstyle \alpha _{v}}$ and ${\textstyle \alpha _{r}}$ 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 ${\textstyle \alpha _{v}}$ and ${\textstyle \alpha _{r}}$ 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 [61].

The transient equation (1.b) can be written as follows

${\displaystyle \rho c{\boldsymbol {v}}^{T}{\boldsymbol {\nabla }}\phi -{\boldsymbol {\nabla }}^{T}({\boldsymbol {D}}{\boldsymbol {\nabla }}\phi )+{\bar {s}}\phi -Q=0}$

(35)

where

${\displaystyle {\bar {s}}=s+s_{t}\quad {\hbox{with }}s_{t}=\rho c{\frac {\dot {\phi }}{\phi }}}$

(36)

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

The non linear reaction term ${\textstyle s_{t}}$ can be approximated as

${\displaystyle s_{t}={\frac {\rho c}{\theta \Delta t}}f(\kappa )}$

(37)

with

${\displaystyle f(\kappa )\geq 0\quad {\hbox{and}}\quad \kappa ={\frac {\phi ^{n+\theta }-\phi ^{n}}{\phi ^{n+\theta }}}}$

(38)

The value of ${\textstyle s_{t}}$ in Eq.37 has been limited in this work to a ten percent of the value of the actual absorption parameter ${\textstyle s}$.

From Eq.(36) we can define an equivalent Damköhler number ${\textstyle {\bar {\sigma }}_{v}}$ as

${\displaystyle {\bar {\sigma }}=\sigma +\sigma _{t}\quad {\hbox{with }}\quad \sigma _{t}={\frac {1}{\theta C}}f(\kappa )}$

(39)

where ${\textstyle \sigma }$ is the Damköhler number defined in Eq.44 and ${\textstyle C={\frac {v\Delta t}{l^{e}}}}$ is the element Courant number where ${\textstyle l^{e}}$ is a characteristic element dimension. For 1D problems ${\textstyle l^{e}}$ is taken as the element length. For 2D problems we have chosen ${\textstyle l^{e}={\sqrt {2\Omega ^{e}}}}$ where ${\textstyle \Omega ^{e}}$ is the area of the element.

The following definition for ${\textstyle f(\kappa )}$ has been chosen in the examples solved in this work

${\displaystyle 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)}$

(40)

where ${\textstyle |a|_{\infty }^{e}}$ denotes the maximum value of ${\textstyle a}$ within an element and ${\textstyle \delta }$ is a positive number that controls the slope of the function ${\textstyle \tanh(\cdot )}$ that ranges from zero to one.

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

Good results were obtained in the transient problems solved using Eq.(40) with ${\textstyle \delta =1}$. The optimal choice of ${\textstyle \delta }$ 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 ${\textstyle \theta =1}$ in the definition of ${\textstyle s_{t}}$ and ${\textstyle \sigma _{t}}$.

The stabilization parameters ${\textstyle \alpha _{v}}$ and ${\textstyle \alpha _{r}}$ in Eqs.11 and 12.a for the transient convection-diffusion-radiation problem are defined as follows

${\displaystyle \alpha _{v}={\begin{cases}\displaystyle {\frac {2}{\bar {\sigma }}}\left(1-{\frac {{\bar {\sigma }}\tanh \gamma }{\xi -1}}\right)&,\quad {\bar {\sigma }}\geq 2^{-12}\\\displaystyle {\frac {\bar {\sigma }}{3}}+{\bar {\alpha }}_{v}\left(1-{\frac {\sigma }{\gamma }}\right)&,\quad {\bar {\sigma }}<2^{-12}\end{cases}}}$

(41)

with

${\displaystyle {\bar {\alpha }}_{v}=\coth \gamma -{\frac {1}{\gamma }}}$

(42)

and

${\displaystyle \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}}}}$

(43)

In the above expressions

${\displaystyle {\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 cs(l_{v})^{2}}{D}}\quad ,\quad \sigma ={\frac {sl_{v}}{|{\boldsymbol {v}}|}}={\frac {w}{2\gamma }}\end{array}}}$

(44)

The expressions of ${\textstyle \alpha _{v}}$ and ${\textstyle \alpha _{r}}$ in Eqs.41 and 43 are an extension of the values for the steady-state problem presented in [62].

In Eq.43 ${\textstyle \varphi }$ is a constant such that ${\textstyle 2\leq \varphi \leq 3}$. The “exact” expression of ${\textstyle \alpha _{r}}$ for 1D problems requires choosing ${\textstyle \varphi =3}$ [62]. In our computations for 2D problems we have obtained good results using ${\textstyle \varphi =2}$.

The characteristic length ${\textstyle l_{v}}$ is defined as the size of the element in the direction of the velocity ${\textstyle {\boldsymbol {v}}}$. Representing the element's edge joining nodes ${\textstyle a}$ and ${\textstyle b}$ with the vector ${\textstyle {\boldsymbol {l}}_{ab}}$, ${\textstyle l_{v}}$ is computed as

${\displaystyle l_{v}={\underset {\hbox{edges}}{\max }}\left\{{\boldsymbol {v}}^{T}{\boldsymbol {l}}_{ab}\right\}}$

(45)

The characteristic length ${\textstyle l_{sc}}$ is generally defined as

${\displaystyle l_{sc}={\sqrt {2\Omega ^{e}}}}$

(46)

where ${\textstyle \Omega ^{e}}$ is the element area. For elements belonging to the boundary in meshes of 3-noded triangles we have chosen

${\displaystyle l_{sc}=2{\sqrt {\Omega ^{e}}}}$

(47)

The expression of ${\textstyle l_{sc}}$ in Eq.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 ${\textstyle \rho }$ and the specific flux c are chosen such that ${\textstyle \rho c=1}$. We have also assumed an isotropic diffusion of value ${\textstyle k}$.

7.1 Transient advection-diffusion-absorption problem

The analysis domain ${\textstyle (x,y)=[0,8]\times [0,8]}$ is discretised into a regular mesh of ${\textstyle (2\times 8)\times (2\times 8)}$ 3-noded triangles of unit rectangular side (${\textstyle l_{v}=1}$) (fig:dac_scheme). The advection, diffusion and reaction coefficients are chosen as ${\textstyle v_{1}=8}$, ${\textstyle k=2}$ and ${\textstyle s=2}$. 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 ${\textstyle \gamma _{v}=2}$ and ${\textstyle \omega _{v}=1}$.

The Dirichlet boundary conditions ${\textstyle \phi (x=0)=3}$ and ${\textstyle \phi (x=8)=8}$ are employed. The initial solution is chosen to have a linear profile. The transient solution was obtained iteratively using the implicit scheme of Eqs.29–33 with ${\textstyle \theta =1.0}$ and a time step of ${\textstyle \Delta t=0.0625s}$ (fig:dac_results). This corresponds to an element Courant number ${\textstyle C=0.5}$. 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 3 show a perspective view of the quasi-steady-state solution for the temperature for ${\textstyle t=2s}$. The steady state solution matches the 1D result reported in [61].

Figure 1: Transient advection-absorption problem. Square domain with linear velocity and zero source.

Figure 2: Transient advection-absorption problem. Transient FIC-FEM solution obtained with a structured mesh of ${\displaystyle 2\times 8\times 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.

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 ${\textstyle (x,y)=[0,1]\times [0,1]}$. The problem data is: ${\textstyle {\boldsymbol {v}}=[-5\cdot {10}^{6},-9\cdot {10}^{6}]^{T}}$, ${\textstyle k=1}$, ${\textstyle \rho c=1}$, ${\textstyle s=0}$ and ${\textstyle Q=0}$. The transported variable is again the temperature. The boundary conditions are: ${\textstyle \phi =1}$ on (${\textstyle x=1,y>0.7)\cup (x<1,y=1)}$, ${\textstyle \phi =0.5}$ at (${\textstyle x=1,y=0.7}$) and ${\textstyle \phi =0}$ on the rest of the boundary (fig:SC_scheme). The transient solution was obtained using the implicit scheme here presented with ${\textstyle \theta =0.8}$ and a time step of ${\textstyle \Delta t=1}$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 6 shows a perspective view of the quasi steady-state solution for ${\textstyle t=2}$e-07s. The steady-state solution agrees with that reported in [51,55,62].

Figure 4: Square domain with non-uniform Dirichlet conditions, downwards diagonal velocity and zero source. Mesh of ${\displaystyle 2\times 20\times 20}$ three-noded triangles

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.

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: ${\textstyle {\boldsymbol {v}}=[1,0]^{T}}$, ${\textstyle k=10^{-8}}$, ${\textstyle \rho c=1}$, ${\textstyle s=0}$ and ${\textstyle Q=1}$. The homogeneous boundary condition ${\textstyle \phi =0}$ is imposed everywhere as shown in fig:Q_source_scheme. The solution develops an exponential layer at the outflow boundary ${\textstyle x=1}$ and parabolic layers at the boundaries ${\textstyle y=0}$ and ${\textstyle y=1}$. 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 ${\textstyle t=1.5s}$ is shown in Figure 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 ${\textstyle \theta =0.8}$ and a time step of ${\textstyle \Delta t=0.005s}$. Convergence within each time step was found in 2-3 iterations. The steady-state solution agrees with that reported in [62].

Figure 7: Square domain with linear velocity and constant source. Structured mesh of ${\displaystyle 2\times 20\times 20}$ three node triangles

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.

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 ${\textstyle \times }$ 2 units and the boundary conditions are shown in fig:rotational_scheme. The rotational velocity field is defined as ${\textstyle \mathbf {u} =10^{4}\times [y(1-x^{2}),-x(1-y^{2})]^{T}}$. A unit isotropic diffusion is assumed. The value of the temperature (${\textstyle \phi }$) is prescribed at part of the domain sides as shown in fig:rotational_scheme. In the rest of the sides, the Neumann condition (${\textstyle q_{n}^{p}=0}$) 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 ${\textstyle \phi =100}$ with a boundary layer at the right hand side where the solution drops towards the prescribed value of ${\textstyle \phi =0}$ and a circular sharp gradient region around the lower circular zone where the solution takes a zero value. 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 ${\textstyle t=1}$e-03s is shown in Figure 12. The transient solution was obtained using the implicit scheme presented in this work with ${\textstyle \theta =0.8}$ and a time step of ${\textstyle \Delta t=1}$e-05s. Convergence within each time step was found in 2-3 iterations. The steady-state solution agrees with that reported in [51].

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

Figure 11: FIC solution of the problem in Figure 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.

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 ${\textstyle (x,y)=[0,1]\times [0,1]}$. The problem data is: ${\textstyle {\boldsymbol {v}}=[1,0]^{T}}$, ${\textstyle k=0}$, ${\textstyle s=0}$, ${\textstyle \rho c=1}$ and ${\textstyle Q=0}$. The homogeneous boundary condition ${\textstyle \phi =0}$ 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 ${\textstyle \phi =1.0}$ from ${\textstyle x=0.1}$ to ${\textstyle x=0.2}$. The problem has been solved using a regular mesh of ${\textstyle 50\times 50}$ four-noded quadrilateral elements. Figures 14 and 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 ${\textstyle \Delta t=0.2s}$, giving an element Courant number of ${\textstyle C=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.

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.

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 [16].

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

The analytical solution for this 2D problem is [16]

${\displaystyle \phi (x,y,t)={\frac {\phi _{0}}{L_{3}4\pi Dt}}e^{-A}}$

(48.a)

with

${\displaystyle A=-{\frac {1}{4Dt}}\left\{[x_{1}-(x_{1}^{0}+v_{1}t)]^{2}+[x_{2}-(x_{2}^{0}+v_{2}t)]^{2}\right\}}$

(48.b)

In Eqs.7.6 ${\textstyle \phi }$ is the initial value of the concentration that is dropped at the point with coordinates ${\textstyle {\boldsymbol {x}}^{0}=[x_{1}^{0},x_{2}^{0}]^{T}}$, ${\textstyle L_{3}}$ is the vertical dimension of the analysis domain (in this problem we have taken ${\textstyle L_{3}=1}$), ${\textstyle v_{1}}$ and ${\textstyle v_{2}}$ are the horizontal and vertical components of the velocity vector, respectively and ${\textstyle D={\frac {K}{\rho c}}}$ where ${\textstyle K}$ is the average diffusion defined in Section 2.1.

The transient solution was obtained using an implicit scheme with ${\textstyle \theta =0.5}$ and a time step of ${\textstyle \Delta t=0.25s}$, giving a Courant number of ${\textstyle C=0.5}$. 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 ${\textstyle t=15s}$ can be seen in Figures 19 and 20, respectively.

Figure 16: Advection-diffusion of a concentration field. Rectangular domain with horizontal velocity. Structured mesh of 2 x 70 x 20 three-noded triangles

Figure 17: Advection-diffusion of a concentration field. FIC-FEM solution along line A-A' obtained with the mesh of Figure 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.

Figure 18: Advection-diffusion of a concentration field. Analytical solution along line A-A' plotted on the structured mesh of Figure 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.

Figure 19: Advection-diffusion of a concentration field. Contours of the FIC-FEM results at t = 15s.

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 [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 [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 ${\textstyle C=1}$. 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

[1] Baiocchi C, Brezzi F, Franca LP. Virtual bubbles and the GaLS. Computer Methods in Applied Mechanics and Engineering 1993; 105:125–141.

[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–129.

[3] Brezzi F, Russo A. Choosing bubbles for advection-diffusion problems. Mathematical Models and Methods in Applied Sciences 1994; 4:571–587.

[4] Brezzi F, Hauke G, Marini LD, Sangalli G. Link-cutting bubbles for the stabilization of convection–diffusion–reaction problems. Mathematical Models and Methods in Applied Sciences 2003; 13(3):445–461.

[5] Brooks AN, Hughes TJR. Streamline Upwind/Petrov–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–3):199–259.

[6] Burman E, Ern A. Nonlinear diffusion and discrete maximum principle for stabilized Galerkin approximations of the convection–diffusion–reaction equation, Computer Methods in Applied Mechanics and Engineering 2002, 191:3833-3855.

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

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

[9] Codina R. On stabilized finite element methods for linear systems of convection–diffusion–reaction equations, Computer Methods in Applied Mechanics and Engineering 2000a, 190:2681-2706.

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

[11] Donea J. A Taylor-Galerkin method for convective transport problems, International Journal for Numerical Methods in Engineering 1984, 20:101-119.

[12] J. Donea, A. Huerta, Finite element method for flow problems. J. Wiley & Sons, 2003.

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

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

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

[16] Fisher H, List J, Koh C, Imberger J, Brooks N. Mixing in Inland and Coastal Waters, Elsevier, 2013.

[17] Franca LP, Dutra do Carmo EG. The Galerkin gradient least-squares method, Computer Methods in Applied Mechanics and Engineering 1989, 74:41-54.

[18] Franca LP, Farhat C. Bubble functions prompt unusual stabilized finite element methods, Computer Methods in Applied Mechanics and Engineering 1995, 123:299-308.

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

[20] Galeo AC, Dutra do Carmo EG. A consistent approximate upwind Petrov–Galerkin method for the convection-dominated problems. Computer Methods in Applied Mechanics and Engineering 1988, 68(1):83-95.

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

[22] Hauke G, García-Olivares A. Variational subgrid scale formulations for the advection-diffusion–reaction equation, Computer Methods in Applied Mechanics and Engineering 2001, 190:6847-6865.

[23] Hauke G. A simple subgrid scale stabilized method for the advection-diffusion–reaction equation, Computer Methods in Applied Mechanics and Engineering 2002, 191:2925-2947.

[24] Hauke G, Sangalli G, Doweidar MH. Combining adjoint stabilized methods for the advection-diffusion–reaction problem, Mathematical Models and Methods in Applied Sciences 2007, 17(2):305-326.

[25] Hughes TJR, Brooks AN. A theoretical framework for Petrov–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.

[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–328.

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

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

[29] Idelsohn SR, Nigro N, Storti M, Buscaglia G. A Petrov–Galerkin formulation for advective-reaction-diffusion problems, Computer Methods in Applied Mechanics and Engineering 1996, 136:27-46.

[30] Idelsohn SR, Heinrich JC, Oñate E. Petrov–Galerkin methods for the transient advective-diffusive equation with sharp gradients, International Journal for Numerical Methods in Engineering 1996, 39:1455-1473.

[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–2215.

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

[33] Kikuchi F, Ushijima T. Theoretical analysis of some finite element schemes for convective diffusion equations, in: R. Gallagher, D. Norrie, J. Oden, O. C. Zienkiewicz (Eds.), Finite Elements in Fluids, Vol. IV, John Wiley and Sons Ltd, Chichester, 1982.

[34] Knobloch P. Improvements of the Mizukami-Hughes method for convection-diffusion equations, Computer Methods in Applied Mechanics and Engineering 2006, 196:579-594.

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

[36] Mizukami A, Hughes TJR. A Petrov–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.

[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–12):525-546.

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

[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–346.

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

[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–292.

[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–221.

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

[44] Oñate E, Sacco C, Idelsohn SR. A finite point method for incompressible flow problems. Computing and Visualization in Science 2000, 2:67–75.

[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–660.

[46] Oñate E, Perazzo F, Miquel J. A finite point method for elasticity problems. Computers and Structures 2001, 79:2151–2163.

[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–658.

[48] Oñate E. Possibilities of finite calculus in computational mechanics. International Journal for Numerical Methods in Engineering 2004, 60(1), 255–281.

[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–307.

[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–1500.

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

[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–36):3926–3946.

[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–455.

[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–3037.

[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–112.

[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–351.

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

[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–134.

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

[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–209

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

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

[63] J. Principe, R. 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–1402.

[64] Tezduyar TE, Park YJ. Discontinuity-capturing finite element formulations for nonlinear convection–diffusion–reaction equations. Computer Methods in Applied Mechanics and Engineering 1986; 59:307–325.

[65] Yu CC, Heinrich JC. Petrov–Galerkin methods for the time-dependent convective transport equations, International Journal for Numerical Methods in Engineering 1986, 23:883-901.

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

[67] Zienkiewicz OC, Taylor RL, Zhu JZ. The finite element method. The basis. 6th Ed., Elsevier, 2005.

[68] Zienkiewicz OC, Taylor RL, Nithiarasu P. The finite element method for fluid dynamics. 6th Ed., Elsevier, 2005.