Why do we need the finite element method and how does it work?

To answer this question, consider the following example of a homogeneous earth dam as depicted in Fig. 1. The dam is located in a region of seismic activity and we are faced with the question of how this dam will behave under an earthquake load. How will it deform? Even if the problem seems quite simple at first glance (it is, after all, a homogeneous dam) the answer to this question is very difficult.

Figure 1. Homogeneous earth dam subjected to earthquake loading.

Although it was just stated, that it is difficult to answer such a question, it is not that complicated to describe the problem in a mathematical way, which can be done with the help of differential equations. In general, (partial) differential equations (P)DE describe nature. For the present case, let us consider the balance of linear momentum:

\[\begin{aligned} \frac{\partial \sigma_{ij}}{\partial x_j} + \rho b_i = \rho \ddot{u}_i ~~~\text{in}~~~ \Omega \end{aligned}\label{eq:linear-momentum}\]

Therein, \(\ddot{u}_i\) is the acceleration and \(\rho\) is the density of the material. Equation (\(\ref{eq:linear-momentum}\)) states that the material time derivative of the momentum \(\rho \ddot{u}_i\) is equal to the sum of external forces, namely the surface forces \(t_i\) (acting on the surface element \(ds\)) as well as the body forces \(\rho b_i\) (acting on the volume element \(dv\), e.g. the gravity in the present case). Very important: this equation is continuous, i.e. must be fulfilled in every point of the area \(\Omega\).

Note that in Eq. (\(\ref{eq:linear-momentum}\)) we used the Cauchy theorem:

\[ t_i = \sigma_{ij} n_j \]

where \(\sigma_{ij}\) is the Cauchy stress tensor and \(n_j\) is the unit normal vector at the surface together with the Gauss divergence theorem. A probably somewhat inaccurate, but hopefully understandable reason for the use of the Cauchy theorem is that it allows us to define the stress \(\sigma_{ij}\) not only at the external boundary of the body \(\Gamma\) but in every point of the body \(\Omega\). In addition, boundary conditions (BCs) can be formulated which constrain the solution of the PDE describing the problem:

\[ \begin{align} \hat{u}_i ~~~\mbox{on} ~~~ \Gamma_u \label{eq:bc-dirichlet}\\ \hat{t}_i ~~~\mbox{on} ~~~ \Gamma_t \label{eq:bc-neumann} \end{align} \]

The first type of BCs, Eq. (\(\ref{eq:bc-dirichlet}\)), is called essential or Dirichlet boundary condition. They prescribe the sought variable (or the unknown of the PDE) directly. The second type of BCs, Eq. (\(\ref{eq:bc-neumann}\)), is referred to as natural or von Neumann boundary conditions. To relate it to the question under consideration here: the Dirichlet BC \(\hat{u}\) corresponds zu the displacements at the bottom of the dam (\(\Gamma_u\)) arising from the earthquake loading displayed in Fig. 1. They can be determined by integrating the acceleration signal twice with respect to time. The von Neumann BCs are the forces acting on the surface \(\Gamma_t\) of the dam and, since no external forces are considered in the example, they are zero in the present case.

To this point, we know how to formulate a mathematical model of the physical problem. Unfortunately we are unable to solve it exactly. The reasons for this are numerous, e.g. that the geometry of the considered body is too complicated for an exact integration or that the variables which find input into the differential equation (e.g. the stiffness or hydraulic conductivity) are strongly nonlinear. This is where the Finite Element Method (FEM) comes into play. The Finite Element Method is a numerical method for solving differential equations (and/or integrals).

In the following sections important components of the FEM are introduced: