Latest revision as of 09:30, 7 June 2019

Abstract

A stabilized finite element method (FEM) for the multidimensional steady state advection-diffusion-absorption equation is presented. The stabilized formulation is based on the modified governing differential equations derived via the Finite Calculus (FIC) method. For 1D problems the stabilization terms act as a nonlinear additional diffusion governed by a single stabilization parameter. It is shown that for multidimensional problems an orthotropic stabilizing diffusion must be added along the principal directions of  curvature of the solution. A simple iterative algorithm yielding a stable and accurate solution for all the range of physical parameters and boundary conditions is described. Numerical results for 1D and 2D problems with sharp gradients are presented showing the effectiveness and accuracy of the new stabilized formulation.

1 INTRODUCTION

Considerable effort has been spent in recent years to derive finite element methods (FEM) 1 for the solution of the advection-diffusion-reaction equation. In this work we will focus on the so called exponential regime originated by large absorptive (dissipative) reaction terms. Here the solutions are of the form of real exponential functions. Numerical schemes find difficulties to approximating the sharp gradients appearing in the neighborhood of boundary and internal layers due to high Peclet and/or Damköhler numbers. Non physical oscilaltory solution are found with the standard Galerkin FEM unless some stabilization procedure is used.

Stabilized methods to tackle this problem have been based on streamline-upwind/Petrov-Galerkin (SUPG) 2, Galerkin/least-squares 5, Subgrid Scale (SGS) 5 and Residual Free Bubbles 14 finite element methods. While a single stabilization parameter suffices to yield stabilized (and even nodally exact results) for the one-dimensional (1D) advection-diffusion and the diffusion-reaction equations (Vol. 3 in 1 and 8), this is not the case for the diffusion-advection-reaction equation. Here, in general, two stabilization parameters are needed in order to ensure a stabilized solution for all range of physical parameters and boundary conditions 4. As reported in 12 the SUPG, GLS and SGS methods with a single stabilization parameter fail to obtain a stabilized solution for some specific boundary conditions in the exponential regime with negative (absorption) terms when there is a negative streamwise gradient of the solution.

Oñate et al. [18] have recently presented a stabilized FEM for the advection-diffusion absorption equation based on the use of a single stabilization parameter which has the form of a diffusion term. In [18] the formulation is detailed for 1D problems and only a brief introduction to the multidimensional case is given. This paper extends the ideas presented in [18] and provides evidence of the effectiveness and accuracy of the new formulation to deal with multidimensional advection-diffusion-absorption problems with sharp gradients.

The stabilized formulation is based on the standard Galerkin FEM solution of the modified governing differential equations derived via the Finite Calculus (FIC) method [19–20]. The FIC equations are obtained by expressing the balance of fluxes in a domain of finite size. This introduces additional stabilizing terms in the differential equations of the infinitessimal theory which are a function of the balance domain dimensions. Although the FIC–FEM formulation here presented is general, we will restrict its application in this work to linear finite element approximations only.

The Galerkin FIC-FEM formulation described here introduces naturally an additional nonlinear dissipation term into the discretized equations which is governed by a single stabilization parameter. In the absence of the absorption term the formulation simplifies to the standard Petrov-Galerkin approach for the advection-diffusion problem For the diffusion-absorption case the diffusion-type stabilization term is identical to that recently obtained by Felippa and Oñate using a variational FIC approach 15. The general nonlinear form of the stabilization parameter is a function of the signs of the solution and its first and second derivatives. This introduces a non-linearity in the solution scheme and a simple iterative algorithm is described. A simpler constant expression of the stabilization parameter is also presented.

Details of the 1D formulation and its extension to deal with multidimensional problems are given. For the multidimensional case Oñate et al. 27 have recently shown that a general form of the stabilization parameters can be found by writting the FIC equations along the principal curvature directions of the solution. The resulting FIC-FEM formulation is equivalent in this case (for linear elements) to adding a stabilizing diffusion matrix to the standard infinitessimal equation. The stabilizing diffusion matrix depends on the signs of the solution and its derivatives and on the velocities along the principal curvature directions of the solution. This introduces a nonlinearity in the solution process. We present a simple iterative scheme based in assuming that the main principal curvature direction at each point is coincident with the gradient vector direction. In the last part of the paper we present a collection of 1D and 2D examples showing the effectiveness and accuracy of the new FIC-FEM formulation for different values of the advective and absorptive terms.

2 FIC FORMULATION OF THE 1D STATIONARY ADVECTION-DIFFUSION-ABSORPTION EQUATION

The governing equation for the 1D stationary advection-diffusion-absorption problem can be written in the FIC formulation as

where

In above equations ${\textstyle \phi }$ is the state variable, ${\textstyle x\in [0,L]}$ is the problem domain, ${\textstyle L}$ is the domain length, ${\textstyle u}$ is the velocity field, ${\textstyle k\geq 0}$ is the diffusion, ${\textstyle s\geq 0}$ is the absorption, dissipation or destruction source parameter, ${\textstyle Q}$ is a constant source term, ${\textstyle q^{p}}$ and ${\textstyle \phi ^{p}}$ are the prescribed values of the total flux and the unknown function at the Neumann and Dirichlet boundaries ${\textstyle \Gamma _{q}}$ and ${\textstyle \Gamma _{\phi }}$, respectively and ${\textstyle h}$ is a characteristic length which plays the role of a stabilization parameter. In the 1D problem ${\textstyle \Gamma _{\phi }}$ and ${\textstyle \Gamma _{q}}$ consist of four combinations at the end points of the problem domain.

Eqs.(1) and (2) are obtained by expressing the balance of fluxes in an arbitrary 1D domain of finite size within the problem domain and at the Neumann boundary, respectively. The variations of the transported variables within the balance domain are approximated by Taylor series expansions retaining one order higher terms than in the infinitessimal theory [19,20]. The underlined stabilizing terms in Eqs.(1) and (2) emanate from these higher order expansions. Note that as the characteristic length parameter ${\textstyle h}$ tends to zero the FIC differential equations gradually recover the standard infinitessimal form.

Successful applications of the FIC method to a variety of problems in computational mechanics can be found in [19–30,37].

3 FINITE ELEMENT FORMULATION

We will construct a standard finite element discretization ${\textstyle \left\{l^{e}\right\}}$ of the 1D analysis domain length ${\textstyle L}$ with index ${\textstyle e}$ ranging from 1 to the number of elements ${\textstyle N}$ 1. The state variable ${\textstyle \phi }$ is approximated by ${\textstyle {\bar {\phi }}}$ over the analysis domain. The approximated variable ${\textstyle {\bar {\phi }}}$ is interpolated within each element with ${\textstyle n}$ nodes in the standard manner, i.e.

with

where ${\textstyle N_{i}}$ are the element shape functions and ${\textstyle \phi _{i}}$ are nodal values of the approximate function ${\textstyle {\bar {\phi }}}$. Substituting Eq.(5a) into Eqs.(1) and (2) gives

where ${\textstyle {\bar {r}}=r({\bar {\phi }})}$ and ${\textstyle r_{\Omega },r_{q}}$ and ${\textstyle r_{\phi }}$ are the residuals of the approximate solution in the problem domain and on the Neumann and Dirichlet boundaries ${\textstyle \Gamma _{q}}$ and ${\textstyle \Gamma _{\phi }}$, respectively.

The weighted residual form of Eqs.(6)–(8) is written as

where ${\textstyle W_{i}(x)}$ and ${\textstyle {\hat {W}}_{i}}$ are test functions satisfying ${\textstyle W_{i}={\hat {W}}_{i}=0}$ on ${\textstyle \Gamma _{\phi }}$.

Assuming smooth enough solutions and integrating by parts the term involving ${\textstyle h}$ in the first integral gives for ${\textstyle {\hat {W}}_{i}=-W_{i}}$

The third term in Eq.(10) is computed as the sum of the integrals over the element interiors, therefore allowing for the space derivatives of ${\textstyle {\bar {r}}}$ to be discontinuous. Also in Eq.(10) ${\textstyle h}$ has been assumed to be constant within each element, (i.e. ${\textstyle \displaystyle {dh \over dx}=0}$ within ${\textstyle l^{e}}$).

The weak form is obtained by integrating by parts the advective and diffusive terms within ${\textstyle {\bar {r}}}$ in the first integral of Eq.(10). This gives

with

where a prime denotes the derivative with respect to the space coordinate.

Wee see clearly that the last term of Eq.(11) introduces within each element an additional diffusion of value ${\textstyle \beta k}$.

Substituting expression (5b) into (11) and choosing a Galerkin method with ${\textstyle W_{i}=N_{i}}$ within each element gives the discrete system of FE equations written in the standard matrix form as

where ${\textstyle {\bar {\boldsymbol {\phi }}}}$ is the vector of nodal unknowns and the element contributions to matrix K and vector ${\textstyle f}$ are

The amount of balancing diffusion in Eq.(14) clearly depends on the (nonlinear) stabilization parameter ${\textstyle \beta }$. The element and critical values of ${\textstyle \beta }$ are deduced in the next section for linear two node elements.

We note that the integral of the term ${\textstyle \displaystyle {h \over 2}\displaystyle {dN_{i} \over dx}Q}$ in Eq.(15) vanishes after asssembly when ${\textstyle h}$ and ${\textstyle Q}$ are uniform over a patch of linear elements.

4 COMPUTATION OF THE STABILIZATION PARAMETER FOR LINEAR ELEMENTS

The matrix ${\textstyle {K}^{e}}$ and the vector ${\textstyle {f}^{e}}$ for two node linear elements are (for constant values of ${\textstyle u,k}$, ${\textstyle s}$ and ${\textstyle Q}$)

In Eqs.(16) index ${\textstyle e}$denotes element values.

Assuming ${\textstyle Q=0}$, a typical stencil for elements of equal size ${\textstyle l}$ can be written as

where for simplicity a constant value of ${\textstyle \beta }$ across the mesh has been assumed. In Eq.(17) ${\textstyle \gamma ={ul \over 2k}}$ and ${\textstyle w={sl^{2} \over 4}}$ are the Peclet number and a velocity independent dimensionless number, respectively.

From Eq.(17) we deduce

In the vecinity of a sharp gradient zone we can take

where ${\textstyle {\bar {\phi }}_{max}}$ is the maximum value of the approximate function ${\textstyle {\bar {\phi }}}$ in the patch of elements adjacent to the sharp gradient zone and

where sign ${\textstyle {\bar {(\cdot )}}}$ denotes the sign of the magnitude within the brackets computed at the patch mid point.

Substituting Eq.(19) into (18) leads to the following expression of the stabilization parameter

The element stabilization parameter ${\textstyle \beta ^{e}}$ is now defined as

where ${\textstyle \beta }$ is given by Eq.(21) and the signs ${\textstyle S_{0}}$, ${\textstyle S_{1}}$ and ${\textstyle S_{2}}$ are computed now at the element mid-point.

It is clear from above that the computation of the stabilization parameter ${\textstyle \beta ^{e}}$ requires the knowledge of the sign of the numerical solution ${\textstyle {\bar {\phi }}}$ and that of the first and second derivatives of ${\textstyle {\bar {\phi }}}$ within each element. This necessarily leads to an iterative scheme. A simple algorithm which provides a stabilized and accurate solution in just two steps is presented below.

4.1 Critical stabilization parameter and unstability conditions

The following constant value of ${\textstyle \beta }$ over the mesh ensures a stabilized solution for all ranges of ${\textstyle \gamma }$ and ${\textstyle w}$

where ${\textstyle \beta _{c}}$ is the critical stabilization parameter. Note that ${\textstyle \beta _{c}}$ corresponds to the maximum value of ${\textstyle \beta }$ in Eq.(21) for ${\textstyle {S_{0} \over S_{2}}={S_{1} \over S_{2}}=1}$. A mathematical proof of Eq.(23) is given in [18].

Clearly the value of ${\textstyle \beta _{c}}$ of Eq.(23) is meaningful only if ${\textstyle \beta _{c}>0}$ and this can be taken as an indicator of an unstable solution. Conversely, a value of ${\textstyle \beta _{c}\leq 0}$ indicates that no stabilization is needed.

4.2 Iterative solution scheme

The following two steps iterative scheme is proposed in order to obtain a stabilized and accurate solution.

step step

Compute a first stabilized solution ${\textstyle {\bar {\boldsymbol {\phi }}}^{1}}$ using the critical value ${\textstyle \beta ^{e}=\beta _{c}}$ given by Eq.(23). This ensures a stabilized, although sometimes slightly overdiffusive, solution.

~ 

step

Compute the signs of the first and second derivatives of ${\textstyle {\bar {\boldsymbol {\phi }}}^{1}}$ within each element. The second derivative field is obtained as follows

where ${\textstyle ({\hat {\cdot }})_{i}^{e}}$ denotes averaged values of the first derivative at node ${\textstyle i}$ of element ${\textstyle e}$. At the boundary nodes the constant value of the derivative of ${\textstyle {\bar {\phi }}}$ within the element is taken in Eq.(24); i.e. ${\textstyle ({\hat {\cdot }})_{i}^{e}=\left({d{\bar {\phi }} \over dx}\right)^{(e)}={{\bar {\phi }}_{2}-{\bar {\phi }}_{1} \over l^{e}}}$.

Compute the enhanced stabilized solution ${\textstyle {\boldsymbol {\phi }}^{2}}$ using the element value of ${\textstyle \beta ^{e}}$ given by Eq.(22).

In all the 1D examples solved the above two steps have sufficed to obtain a converged stabilized and accurate solution [18]. The reason of this is that the signs of the first and second derivative fields typically do not change any further after the second step over the elements adjacent to high gradient zones.

5 EXTENSION TO MULTI-DIMENSIONAL PROBLEMS

Consider the general steady-state advection-diffusion-reaction equation written for the zero constant source case (${\textstyle Q=0}$) as

For 2D problems

are respectively the velocity vector, the gradient vector and the diffusivity matrix, respectively. For simplicity we have assumed an isotropic diffusion matrix.

The FIC form of Eq.(25) is written as

where ${\textstyle r}$ is the original infinitessimal differential equation as defined in Eq.(25).

The Dirichlet and boundary conditions of the FIC formulation are

where ${\textstyle n}$ is the normal vector to the boundary where the normal flux is prescribed. As usual index ${\textstyle p}$ denotes the prescribed values.

In Eqs.(27) and (29) ${\textstyle {h}=[h_{x},h_{y}]^{T}}$ is the characteristic vector of the 2D FIC formulation which components play the role of stabilization parameters. The underlined terms in Eqs.(27) and (29) introduce the necessary stability in the numerical solution [19,20,21].

If vector h is taken to be parallel to the velocity u the standard SUPG method is recovered [18–23]. The more general form of h allows to obtain stabilized finite element solutions in the presence of strong gradients of ${\textstyle \phi }$ near the boundaries (boundary layers) and within the analysis domain (internal layers) 27. The FIC formulation therefore reproduces the best features of the so called shock-capturing or transverse-dissipation schemes 2.

Let us write down the FIC balance equation in the principal curvature axes of the solution ${\textstyle \xi ,\eta }$ (Figure 1). The FIC balance equation is

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

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

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

or

We can see clearly from Eq.(33) that the FIC governing equations introduce orthotropic diffusion parameters of values ${\textstyle {\beta _{\xi }k}}$ and ${\textstyle {\beta _{\eta }k}}$ along the ${\textstyle \xi }$ and ${\textstyle \eta }$ axes, respectively with

Also note that the last term of Eq.(32b) will vanish after discretization for linear elements.

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

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

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

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

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

where

and ${\textstyle {\bar {\phi }}}$ is as usual the approximate solution.

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