Latest revision as of 12:35, 27 May 2019

Published in Finite Element Methods: 1970's and Beyond, Franca L.P., Tezduyar T.E. and Masud A. (Eds.),CIMNE, Barcelona, Spain, 2004

Abstract

The expression “finite calculus” refers to the derivation of the governing differential equations in mechanics by invoking balance of fluxes, forces, etc. in a space-time domain of finite size. The governing equations resulting from this approach are different from those of infinitesimal calculus theory and they incorporate new terms which depend on the dimensions of the balance domain. The new governing equations allow to derive naturally stabilized numerical schemes using any discretization procedure. The paper discusses the possibilities of the finite calculus method for the finite element solution of convection-diffusion problems with sharp gradients and incompressible fluid flow.

Keywords Stabilization, finite calculus, finite element method.

1 INTRODUCTION

It is well known that standard numerical methods such as the central finite difference (FD) method, the Galerkin finite element (FE) method and the finite volume (FV) method, among others, lead to unstable numerical solutions when applied to problems involving different scales, multiple constraints and/or high gradients. Examples of these situations are typical in the solution of convection-diffusion problems, incompressible problems in fluid and solid mechanics and strain or strain rate localization problems in solids and compressible fluids using the standard Galerkin FE method or central scheme in FD and FV methods [1,2]. Similar instabilities are found in the application of meshless methods to those problems [3-5].

The sources of the numerical instabilities in FE, FD and FV methods, for instance, have been sought in the apparent unability of the Galerkin FE method and the analogous central difference scheme in FD and FV methods, to provide a numerical procedure able to capture the different scales appearing in the solution for all ranges of the physical parameters. Typical examples are the spurious numerical oscillations in convection-diffusion problems for high values of the convective terms. The same type of oscillations are found in regions next to sharp internal layers appearing in high speed compressible flows (shocks) or in strain localization problems (shear bands) in solids. A similar problem of different nature emerges in the solution of incompressible problems in fluid and solid mechanics. Here the difficulties in satisfying the incompressibility constraint limit the choices of the approximation for the velocity (or displacement) variables and the pressure [1].

The solution of above problems has been attempted in a number of ways. The underdiffusive character of the central difference scheme for treating advective-diffusive problems has been corrected in an ad-hoc manner by adding the so called “artificial diffusion” terms to the standard governing equation [2]. The same idea has been successfully applied to derive stabilized FV and FE methods for convection-diffusion and fluid-flow problems [1,2]. Other stabilized FD schemes are based on the “upwind” computation of the first derivatives appearing in the convective operator [2]. The counterpart of upwind techniques in the FEM are the so called Petrov-Galerkin methods [1,6,7]. Among the many methods of this kind we can name the SUPG method [8-10], the Galerkin Least Square (GLS) method [11,12] the Characteristic Galerkin method [13-15] the Characteristic Based Split (CBS) method [16,17] and the Subgrid Scale (SS) method [18-21].

The Finite Calculus (FIC) is a different route to derive stabilized numerical methods. The starting point are the modified governing differential equations of the problem derived by expressing the balance of fluxes (or equilibrium of forces) in a space-time domain of finite size [22]. This introduces naturally new terms in the classical differential equations of the infinitesimal theory which are a function of the balance domain dimensions. The merit of the modified equations via the FIC approach is that they lead to stabilized schemes using any numerical method. In addition, the different stabilized FD, FE and FV methods typically used in practice can be recovered using the FIC equations [22,23]. Moreover, these equations are the basis for deriving a procedure for computing the stabilization parameters [24,25].

Most stabilized FEM schemes can be framed within an extended Galerkin approach where the standard Galerkin expression is modified by adding adequate residual-based terms as

where [Galerkin] denotes the standard Galerkin expression, ${\textstyle {\boldsymbol {P}}(N_{j})}$ is a vector which terms depend on the shape functions ${\textstyle N_{j}}$ (and the physical parameters of the problem), ${\displaystyle r}$ is the vector of residuals of the finite element equations and ${\textstyle [\tau ^{e}]}$ is a matrix of stabilization parameters.

The computation of the stabilization parameters is still an open problem and much effort has been devoted to this topic [1,2,7-21,36].

Different stabilized methods can be derived from Eq.(1) by chosing expressions for matrices ${\displaystyle P}$ and ${\textstyle [\tau ^{e}]}$. In this paper we will focus in the stabilized finite element formulation using the FIC method. The equivalent form of Eq.(1) is written in the FIC formulation as [1-7]

where ${\textstyle n_{r}}$ is the number of residual equations, ${\textstyle r_{i}}$ is the ${\textstyle i}$th residual equation and ${\textstyle h_{j}}$ are characteristic length parameters which are typically of the order of the element dimensions.

The FIC formulation has been used in conjunction with the finite element formulation to solve a variety of problems in convection-diffusion [23-27] incompressible fluid dynamics involving free surfaces [28-32] and non linear solid mechanics problems allowing for large strains [28,41,42] using in all cases linear triangles and tetrahedra with equal interpolation for all variables.

The layout of the paper is the following. In the next section the main concepts of the FIC method are introduced. Applications of the FIC method to convection-diffusion problems with sharp gradients are detailed and some examples of application are given. Finally the possibilities of the FIC method in incompressible fluid mechanics are discussed and a finite element formulation is presented.

2 THE FINITE CALCULUS METHOD

We will consider a convection-diffusion problem in a 1D domain ${\textstyle \Omega }$ of length ${\textstyle L}$. The equation of balance of fluxes in a subdomain of size ${\textstyle d}$ belonging to ${\textstyle \Omega }$ (Figure 1) is written as

where ${\textstyle q_{A}}$ and ${\textstyle q_{B}}$ are the incoming and outgoing fluxes at points ${\textstyle A}$ and ${\textstyle B}$, respectively. The flux ${\textstyle q}$ includes both convective and diffusive terms; i.e. ${\textstyle q=v\phi -k{d\phi \over dx}}$, where ${\textstyle \phi }$ is the transported variable (i.e. the temperature in a thermal problem), ${\textstyle v}$ is the velocity and ${\textstyle k}$ is the diffusitivity of the material.

Figure 1: Equilibrium of fluxes in a space balance domain of finite size

Let us express now the fluxes ${\textstyle q_{A}}$ and ${\textstyle q_{B}}$ in terms of the flux at an arbitrary point ${\textstyle C}$ within the balance domain (Figure 1). Expanding ${\textstyle q_{A}}$ and ${\textstyle q_{B}}$ in Taylor series around point ${\textstyle C}$ up to second order terms gives

Substituting Equation (2) into Equation (1) gives after simplification

where ${\textstyle h=d_{1}-d_{2}}$ and all the derivatives are computed at the arbitrary point ${\textstyle C}$.

Standard calculus theory assumes that the domain ${\textstyle d}$ is of infinitesimal size and the resulting balance equation is simply ${\textstyle {dq \over dx}=0}$. We will relax this assumption and allow the space balance domain to have a finite size. The new balance equation (3) incorporates now the underlined term which introduces the characteristic length ${\textstyle h}$. Obviously, accounting for higher order terms in Equation (2) would lead to new terms in Equation (3) involving higher powers of ${\textstyle h}$.

Distance ${\textstyle h}$ in Equation (3) can be interpreted as a free parameter depending on the location of point ${\textstyle C}$ within the balance domain. Note that ${\textstyle -d\leq h\leq d}$ and, hence, ${\textstyle h}$ can take a negative value. At the discrete solution level the domain ${\textstyle d}$ should be replaced by the balance domain around a node. This gives for an equal size discretization ${\textstyle -l^{e}\leq h\leq l^{e}}$ where ${\textstyle l^{e}}$ is the element or cell dimension. The fact that Equation (3) is the exact balance equation (up to second order terms) for any 1D domain of finite size and that the position of point ${\textstyle C}$ is arbitrary, can be used to derive numerical schemes with enhanced properties simply by computing the characteristic length parameter from an adequate “optimality” rule leading to an smaller error in the numerical solution [22-24].

Consider, for instance, Equation (3) applied to the 1D convection-diffusion problem. Neglecting third order derivatives of ${\textstyle \phi }$, Equation (3) can be rewritten in terms of ${\textstyle \phi }$ as

We see clearly that the FIC method introduces naturally an additional diffusion term in the standard convection-diffusion equation. This is the basis of the popular “artificial diffusion” procedure [1,2,6] where the characteristic length ${\textstyle h}$ is typically expressed as a function of the cell or element dimension. The critical value of ${\textstyle h}$ can be computed by introducing an optimality condition, such as obtaining a physically meaningful solution. An interpretation of the FIC equations as a modified residual method is presented in [28].

Equation (3) can be extended to account for source terms. The modified governing equation can then be written in compact form as

with

where ${\textstyle Q}$ is the external source.

The essential (Dirichlet) boundary condition for eqn. (7) is the standard one (i.e. ${\textstyle \phi ={\bar {\phi }}}$ on ${\textstyle \Gamma _{\phi }}$ where ${\textstyle \Gamma _{\phi }}$ is the boundary whete the prescribed value ${\textstyle {\bar {\phi }}}$ is imposed). For consistency a stabilized Neumann boundary condition must be obtained.

Figure 2: Balance domain next to a Neumann boundary point B

The length of the balance segment ${\textstyle AB}$ next to a Neumann boundary is taken as one half of the characteristic length ${\textstyle h}$ for the interior domain (Figure 2). The balance equation, assuming a constant distribution for the source ${\textstyle Q}$ over ${\textstyle AB}$, is

where ${\textstyle {\bar {q}}}$ is the prescribed total flux at ${\textstyle x=L}$ and ${\textstyle x_{A}=x_{B}-{\frac {h}{2}}}$.

Using a second order expansion for the advective and diffusive fluxes at point ${\textstyle A}$ gives [22]

where ${\textstyle r}$ is given by eqn. (8).

Note that for ${\textstyle h=0}$ the infinitesimal form of the 1D Neumann boundary condition is obtained.

The underlined terms in Equations (7) and (10) introduce the necessary stabilization in the discrete solution using whatever numerical scheme.

The time dimension can be simply accounted for the FIC method by considering the balance equation in a space-time slab domain. Application of the FIC method to the transient convection-diffusion equations and to fluid flow problems can be found in [23,26,29-33]. Quite generally the FIC equation can be written for any problem in mechanics as [22]

where ${\textstyle r_{i}}$ is the ith standard differential equation of the infinitesimal theory, ${\textstyle h_{j}}$ are characteristic length parameters, ${\textstyle \delta }$ is a time stabilization parameter and ${\textstyle t}$ the time; ${\textstyle n_{b}}$ and ${\textstyle n_{d}}$ are respectively the number of balance equations and the number of space dimensions of the problem (i.e., ${\textstyle n_{d}=2}$ for 2D problems, etc.). Indeed for the transient case the initial boundary conditions must be specified. The usual sum convention for repeated indexes is used in the text unless otherwise specified.

For example, in the case of the convection-diffusion problem ${\textstyle n_{b}=1}$ and Equation (8) is particularized as

with

For a transient solid mechanics problems Equation (11) applies with ${\textstyle n_{b}=n_{d}}$ and

where ${\textstyle u_{i}}$ are the displacements, ${\textstyle \sigma _{ij}}$ are the stresses and ${\textstyle b_{i}}$ the external body forces.

The modified Neumann boundary conditions in the FIC formulation can be written in the general case as [22]

where ${\textstyle q_{ij}}$ are the generalized “fluxes” (such as the heat fluxes in a thermal problem or the stresses in solid or fluid mechanics), ${\textstyle {\bar {t}}_{i}}$ are the prescribed boundary fluxes and ${\textstyle n_{j}}$ are the components of the outward normal to the Neumann boundary ${\textstyle \Gamma _{q}}$.

In above equations we have underlined once more the terms introduced by the FIC approach which are essential for deriving stabilized numerical formulations.

3 WHAT DO THE CHARACTERISTIC PARAMETERS MEAN?

The characteristic parameters in the space and time dimension (${\textstyle h_{i}}$ and ${\textstyle \delta }$) can be interpreted as free intrinsic parameters giving the “exact” expression of the balance equations (up to first order terms) in a space-time domain of finite size.

Let us consider now the discretized solution of the modified governing equations. For simplicity we will focuss here in the simplest scalar 1D problem. The variable ${\textstyle \phi }$ is approximated as ${\textstyle \phi \simeq {\hat {\phi }}}$ where ${\textstyle {\hat {\phi }}}$ denotes the approximated solution. The values of ${\textstyle {\hat {\phi }}}$ are now expressed in terms of a finite set of parameters using any discretization procedure (finite elements, finite diferences, etc.). The discretized FIC governing equation would read now

where ${\textstyle {\hat {r}}:=r({\hat {\phi }})}$ and ${\textstyle r_{\bar {\Omega }}}$ is the residual of the discretized FIC equations.

The meaning of the characteristic parameters ${\textstyle h}$ and ${\textstyle \delta }$ in the discretized equation (16) changes (although the same simbols than in Eqs.(11) and (12) have been kept). The ${\textstyle h}$ and ${\textstyle \delta }$ parameters become now of the order of magnitude of the discrete domain where the balance laws are satisfied. In practice, this means that

where ${\textstyle l}$ is the grid dimension in the space discretization (i.e. the element size in the FEM or the cell size in the FDM) and ${\textstyle \Delta t}$ is the time step used to solve the transient problem.

The values of the characteristic parameters can be found now in order to obtain an a correct numerical solution. The meaning of “correct solution” must be obviously properly defined. Ideally, this would be a solution giving “exact” values at a discrete number of points (i.e. the nodes in a FE mesh). This is infortunatelly impossible for practical problems (with the exception of very simple 1D and 2D cases) and, in practice, ${\textstyle h}$ and ${\textstyle \delta }$ are computed making use of some chosen optimality rule, such as ensuring that the error of the numerical solution diminishes for appropriate values of the characteristic parameters [22-27]. This suffices in practice to obtain physically sound (stable) numerical results for any range of the physical parameters of the problem, always within the limitation of the discretization method chosen.

It is quite usual to accept that the characteristic parameters are constant within each element. This assumption is not justificable “a priori” and, in general, we should regard those parameters as functions of the space and time dimension and of the solution at each point of the analysis domain.

4 FINITE ELEMENT DISCRETIZATION OF THE FIC EQUATIONS FOR ADVECTIVE-DIFFUSIVE PROBLEMS

Let us consider the FIC governing equations for the steady-state advective-diffusive problem defined in vector form for the multidimensional case as

with

with

In above equations ${\textstyle {\boldsymbol {h}}}$ is the characteristic length vector, ${\textstyle {\boldsymbol {\nabla }}}$ is the gradient vector, ${\textstyle {\boldsymbol {D}}}$ is the diffusivity matrix, ${\textstyle {\boldsymbol {n}}}$ is the normal vector and ${\textstyle {\boldsymbol {v}}}$ is the velocity vector.

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

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

Application of the Galerkin FE method to Equations (12)-(13) gives, after integrating by parts the term ${\textstyle {\boldsymbol {\nabla }}r}$ (and neglecting the space derivatives of ${\textstyle {\boldsymbol {h}}}$)

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

Note that the residual terms have disappeared from the Neumann boundary ${\textstyle \Gamma _{q}}$. This is due to the consistency of the FIC terms in Equation (20b).

The last term in the third integral involving the derivatives of ${\textstyle {\boldsymbol {h}}}$ vanishes if ${\textstyle {\boldsymbol {h}}}$ is assumed to be constant wihtin each element.

4.1 Equivalence with SUPG form

The definition of vector ${\textstyle {\boldsymbol {h}}}$ is a crucial step as the quality of the stabilized solution depends on the module and direction of ${\textstyle {\boldsymbol {h}}}$.

We could, for instance, assume that vector ${\textstyle {\boldsymbol {h}}}$ is parallel to the velocity ${\textstyle {\boldsymbol {v}}}$, i.e. ${\textstyle {\boldsymbol {h}}=h{{\boldsymbol {v}} \over \vert {\boldsymbol {v}}\vert }}$ where ${\textstyle h}$ is a characteristic length. Under these conditions, Equation (23) reads (for ${\textstyle {\boldsymbol {h}}}$ assumed to be constant within each element)

Equation (24) coincides precisely with the so called Streamline-Upwind-Petrov-Galerkin (SUPG) method [1,6,7,10]. The ratio ${\textstyle \displaystyle {h \over 2\vert {\boldsymbol {v}}\vert }}$ has dimensions of time and it is termed element intrinsic time parameter ${\textstyle \tau }$.

It is important to note that the SUPG expression is a particular case of the more general FIC formulation. This explains the limitations of the SUPG method to provide stabilized numerical results in the vicinity of sharp gradients of the solution transverse to the flow direction. In general, the adequate direction of ${\textstyle {\boldsymbol {h}}}$ is not coincident with that of ${\textstyle {\boldsymbol {v}}}$ and the components of ${\textstyle {\boldsymbol {h}}}$ introduce the necessary stabilization along the streamlines and the transverse directions to the flow. In this manner, the FIC method reproduces the best-features of the so-called stabilized discontinuity-capturing schemes [1,34,35].

5 INTERPRETATION OF THE DISCRETE SOLUTION OF THE FIC EQUATIONS

Let us consider the solution of a physical problem in a space domain ${\textstyle \Omega }$, governed by a differential equation ${\textstyle r(\phi )=0}$ in ${\textstyle \Omega }$ with the corresponding boundary conditions. The “exact” (analytical) solution of the problem will be a function giving the sought distribution of ${\textstyle \phi }$ for any value of the geometrical and physical parameters of the problem. Obviously, since the analytical solution is difficult to find (practically impossible for real situations), an approximate numerical solution is found ${\textstyle \phi \simeq {\hat {\phi }}}$ by solving the problem ${\textstyle {\hat {r}}=0}$, with ${\textstyle {\hat {r}}=r({\hat {\phi }})}$, using a particular discretization method (such as the FEM). The distribution of ${\textstyle \phi }$ in ${\textstyle \Omega }$ is now obtained for specific values of the geometrical and physical parameters. The accuracy of the numerical solution depends on the discretization parameters, such as the number of elements and the approximating functions chosen in the FEM. Figure 1 shows a schematic representation of the distribution of ${\textstyle {\hat {\phi }}}$ along a line for different discretizations ${\textstyle M_{1},M_{2},\cdots ,M_{n}}$ where ${\textstyle M_{1}}$ and ${\textstyle M_{n}}$ are the coarser and finer meshes, respectively. Obviously for ${\textstyle n}$ being sufficiently large a good approximation of ${\textstyle \phi }$ will be obtained and for ${\textstyle M_{\infty }}$ the numerical solution ${\textstyle {\hat {\phi }}}$ will coincide with the “exact” (and probably unreachable) analytical solution ${\textstyle \phi }$ at all points. Indeed in some problems the ${\textstyle M_{\infty }}$ solution can be found by a “clever” choice of the discretization parameters.

An unstable solution will occur when for some (typically coarse) discretizations, the numerical solution provides non physical or very unaccurate values of ${\textstyle {\hat {\phi }}}$. A situation of this kind is represented by curves ${\textstyle M_{1}}$ and ${\textstyle M_{2}}$ of the left hand side of Figure 1. These unstabilities will disappear by an appropriate mesh refinement (curves ${\textstyle M_{3},M_{4}\cdots }$ in Figure 1) at the obvious increase of the computational cost.

In the FIC formulation the starting point are the modified differential equations of the problem as previously described. These equations are however not useful to find an analytical solution, ${\textstyle \phi (x)}$, for the physical problem. Nevertheless, the numerical solution of the FIC equation can be readily found. Moreover, by adequately choosing the values of the characteristic length parameter ${\textstyle h}$, the numerical solution of the FIC equations will be always stable (physically sound) for any discretization level chosen.

This process is schematically represented in Figure 3 where it is shown that the numerical oscillations for the coarser discretizations ${\textstyle M_{1}}$ and ${\textstyle M_{2}}$ dissapear when using the FIC procedure.

We can conclude the FIC approach allows us to obtain a better numerical solution for a given discretization. Indeed, as in the standard infinitesimal case, the choice of ${\textstyle M_{\infty }}$ will yield the (unreachable) exact analytical solution and this ensures the consistency of the method.

Figure 3: Schematic representation of the numerical solution of a physical problem using standard infinitesimal calculus and finite calculus.

6 COMPUTATION OF THE CHARACTERISTIC LENGTH VECTOR

The computation of the characteristic lengths is a crucial step as its value affect to the stability (and accuracy) of the numerical solution. This problem is common to all stabilized FE methods and different approaches to compute the stabilization parameters using typically extensions of the optimal values for simple 1D case (giving a nodally exact solution) have been proposed [1,2,6-21].

The computation the characteristic length values in 2D and 3D problems is of much higher complexity than in 1D problems. Indeed in 1D situations ${\textstyle h}$ is an escalar (either positive or negative) while it becomes a vector 2D/3D problems. Moreover, numerical experiments indicate that in 2D/3D situations the direction of vector ${\textstyle h}$ is crucial in order to obtain stabilized numerical solutions in problems where boundary layers and arbitrary internal sharp layers exist. It is well known that the SUPG assumption (${\textstyle h}$ being parallel to ${\textstyle u}$) generally does not suffice to give stable results in those cases [1,2,8,10,34,35]. Conversely, the numerical results in those cases are more insensitive to the module of ${\textstyle h}$ which can be taken of the order of a typical element dimension.