In this section, we present several finite element (FE) models, all based on triangular elements. The development of the FE method was contiguous to the development of computers in the early 60s and 70s. By the end of 70s, they were well established. They became particularly popular in engineering for the computation of stresses over structures, and somewhat popular in fluid mechanics and electricity. In all cases, they were and are still used because of the great flexibility they offer in term of geometrical representation, sometimes despite the cost or the lack of stability of the method. In structural mechanics, they have the decisive advantage of being able to follow the deformation of the mesh due to stress (Lagrangian time formulation), making the methods quite ``natural'' to this field. In electricity, the absence of nonlinear terms in most applications render the method reasonably successful. But in fluid mechanics, the method has always suffered from a lack of overall stability, from a lack of accuracy in the computation of the advective terms and from a very large cost, to the point that most commercial models used for engineering applications preferentially use the finite volume (FV) method along with near regular meshes (c.f., IDEAS, Star-CD, ...). The disadvantage of the FV methods is that they are usually of low order and that they require near regular meshes to ensure good performances. This means that they are not very suitable to model the ocean.
The problem of stability in FE methods in fluid mechanics was early analized by Ladyzhenskaya lady69, Babouska babuska71 and Brezzi brezzi74, who gave their names to the so-called Ladyzhenskaya, Babouska and Brezzi (LBB) stability condition. Their work focuses on the Stokes equations and they demonstrated the need, in fluid mechanics FE methods, for using different basis functions for the velocity and pressure. This was equivalent to staggering the grid in space, as was done for the FD methods. Furthermore, not any combination of basis functions satisfies the LBB condition [Fortin and Fortin, 1985,Pierre, 1988,Idelsohn et al. , 1995,Le Roux et al. , 1998]. Arnold et al. arnold84 and Fortin and Fortin fortin85 emphasized that one simple way to stabilize equal-order schemes is to add the so-called bubble function, and that this method does not lead necessarily to an additional cost, thanks to static condensation techniques (some easy manual Gaussian elimination before solving numerically the matrix problem). But, since, according to Pierre pierre88, these methods are equivalent to adding a penalty term in the fluid equations, they may be over-dissipative in the context of unsteady flows and the more general Navier-Stokes equations.
Mainly, the LBB condition comes down to increasing the order (or the
number of degrees of freedom) of the basis functions for the velocity
compared to the basis function for the pressure. However, one unresolved
issue related to the LBB condition is its relevance for the shallow water
equations. The three shallow water equations are similar enough that they
can be generalized to one vector equation:
The fact that the pressure basis functions have to be of lower order compared to the velocity basis functions means that the overall truncation order of the stabilized FE methods for the shallow water equations is probably lower than the one permitted by the velocity basis functions. And, since the other disadvantage of these stable FE methods is the cumbersome and time consuming solving of a large matrix problem (especially when all the variables are coupled), it is doubtful that these methods can compete with todays FD ocean models in terms of cost and accuracy.
We focus herein on four different FE models: the Lynch and Werner model lynch87,lynch91 (also called the Quoddy model and the only one of the four used for coastal oceanography), the Le Roux et al. (2000; hereafter LLS) model, the Hua and Thomasset hua84 model and the Peraire et al. (1986; hereafter PZM) model. Only one model, the LLS model, among the four satisfies the LBB condition of stabilty. The Quoddy and the PZM models use a non-staggered (i.e., equal order) formulation of the variables and therefore require some kind of stabilizing ``trick'' which we will present and discuss. Due to their equal-order formulation, these two models are the simplest, in some sense, of the four for the same reason that the A-grid FD formulation is simpler than the other staggering techniques. In the context of the finite elements, there are some additional technical advantages to using non-staggering formulations which stems from a lower number of matrices to define and to inverse. Also, it unifies the use of gradient or divergence operators. In general, equal-order models are fairly easy to implement from scratch. Hence, they can be more appealing than more complex LBB complying formulations.
Most FE methods are based on the Galerkin formulation. In these models,
the domain, ,
is broken up into a set of conformal elements
(conformal in the sense that all elements connect to neighboring elements
through common vertices). The form of the elements is rather unspecified
but triangles or quadrangles are usually recommended. We favor the use of
triangles (Fig. 2.5) because complex domains are more easily
divided into triangles than quadrangles. For each vertex of the mesh,
Mi, and in the context of linear finite elements, there is an
associated basis function, .
This basis function is piecewise
linear in each triangle to which Mi belongs and forms a ``hat'' on top
of Mi (Fig. 2.6). Over the rest of the domain, the
basis function is zero. Let us consider the equation
We now consider the issue of using irregularly spaced grid points in FD
and FE methods. In the FD method, an irregular spacing of the nodes leads
to a loss of order. Let us consider the equation
The same occurs for the FE method. Using linear
``hat'' functions to discretize this axis, the Galerkin discretization of
(2.31) with the basis functions at node 2 as the test function
The Quoddy model of Lynch and Werner lynch87,lynch91 is a full
3D baroclinic finite element model. This model was successfully applied
for coastal and tidal studies on the Scotian Shelf [Hannah et al. , 2000] and the
Vancouver Island area [Foreman et al. , 2000]. It was modified to model the
shallow water equations, retaining the main characteristic of the Quoddy
model, which are: equal order of approximation for velocity and elevation,
the divergence of the vertically integrated momentum equations can be
recast by using the mass balance equation and
eliminating the divergence of the vertically integrated mass
flux. This yields a prognostic equation for the
elevation of second order in time (a wave equation). Solving numerically
for the three equations (two momentum equations and one wave equation)
is easy and leads to a stable model, but, does not balance mass locally.
To guaranty a better local conservation of mass, a weighted mass equation
is added to the wave equation (the mass equation can also be viewed as a
PZM developed an
interesting model using equal-order interpolation and a two step explicit
time-integration. First, mean values are computed for each triangle
centroid from fluxes at vertices and then values at vertices are computed
from flux computed at triangle centroids. This model is part of a broader family
of Taylor-Galerkin formulations. We reproduce the demonstration about the
Taylor-Galerkin formulation of Priestley pries92. Starting
from the following prognostic equation in a conservation form
Hua and Thomasset hua84 developed a finite element model staggered in space, using discontinuous linear non conforming (P1NC) basis functions for the velocity (Fig. 2.8) and the usual linear basis functions (P1) for the pressure. This formulation leads to a diagonal mass matrix for velocity, which leads to a simplified matrix problem to solve for the elevation when semi-implicitly discretizing in time. They claim the model to be oscillation-free, although LLS demonstrated that the Hua and Thomasset model does not satisfy the LBB condition of stability for the Stokes flow problem. In the shallow water context, it shows some signs of instability (not shown). After some tests, we chosed to integrate the equations explicitly in time using a Runge-Kutta integration technique instead of the semi-implicit technique proposed by Hua and Thomasset because the instability problem was then less severe.
LLS proposed to use a semi-implicit semi-Lagrangian time integration along with a spatial FE discretization that satisfies the LBB stability condition. The particularity of their choice for the basis functions resides in using macro-elements. Each macro-elements is cut into four sub-triangular elements. The basis functions for the velocity are linear inside each sub-triangle and the basis functions for the pressure (or elevation) are constant (see Fig. 2.9). The equation for the elevation can be inverted locally. Hence, the solution of the coupled system of shallow water equations can be reduced to solving Helmholtz-like coupled equations for the velocity components. In order to interpolate the variables at the previous time step on an unstructured mesh, they also proposed a high order kriging method (see Trochu, 1993, for a review). Using this interpolation technique, they found that the model was performing very well for the purely advective problem. However, the application to a finite element shallow water model was somewhat disappointing. The high order method destabilizes the model (personal com.). Therefore, a low order kriging method had to be used, leading to potentially high artificial viscosity. The mass was not conserved, forcing LLS to add a mass corrector. Another disadvantage of the LLS formulation is the fact that the elevation basis functions are piecewise constant. This means that the truncation order of the model for the elevation is lower than that for the velocity and might lower the order of the velocity as well, since the equations for velocities and elevation are coupled in the shallow water equations.
Get the PS or PDF version here