Finite Element Discretisation

Galerkin himself, could only take a few functions to approximate the solution. So the quality of the approximation depended on how good the estimate of the actual solution was. The basic idea of the FEM is to divide the domain \(\Omega\) into a finite number of sub domains \(\Omega_e\) - finite elements - with simple geometries. This is also called the finite element discretization of the model area. The elements are interconnected at nodes at which the solution (physical field which is sought) is unknown. As an example, the finite element mesh of the homogeneous earth dam is given in Figure 1.

Figure 1. Finite element mesh of the homogeneous earth dam depicted in Figure 1.

Instead of defining the trial functions \(N_I(x)\) on the entire domain \(\Omega\), the trial functions are defined on the element domain \(\Omega_e\), i.e. they approximate the solution only in a sub domain rather than the solution over the entire domain \(\Omega\) (see right-hand side of Figure 1).

As introduced in the previous section, we require a weighted integral of the residual \(R\) to zero in an average sense over the domain. By dividing the domain \(\Omega\) in a many subdomains \(\Omega_e\), the solution of the whole system becomes equal to the sum of the solutions of all subdomains. This allows us to express the Integral over the weighted residual \(R_I(u^*(x),x)\) as follows:

\[ \begin{align} R_I(u^*(x),x) &= \int_\Omega N_I(x) D(u^*(x),x) \\ &\approx \sum_{Ie}^{Ne} \int_{\Omega_e} N_I(x) D(u^*(x),x) d \Omega_e ~~~\text{with:}~ I = \{1,2,...N\} \label{eq:weighted-integral2} \end{align} \]

A crucial feature of the FEM is that the geometry of the finite elements is selected to be simple (e.g. triangular or rectangular), so that the trial functions can be generated systematically and after a certain pattern. The choice and properties of the test functions in the FEM will be discussed in more detail below. Compared to the (original) Galerkin approach, we are talking about thousands of very simple trial functions instead of only a few but very complicated ones. At this point it becomes obvious that the development of the FEM would not have been possible without the development of the computer. For the discretisation (sub division) of the domain different element shapes can be used. The choice of element geometry depends mostly on the dimension (1D, 2D or 3D) and the complexity of the geometry of the body to be discretized. A set of the most commonly used element shapes is shown in Figure 2.

Figure 2. Different element shapes for the discretisation of the physical domain.

To calculate Equation (\(\ref{eq:weighted-integral2}\)), the individual terms must be integrated over the volume \(\Omega_e\) of each element. The evaluation of these integrals is however difficult and their analytical solution if often impossible. Therefore, numerical integration is usually used. Relevant information is given in section Numerical Integration. Applying Eq. (\(\ref{eq:weighted-integral2}\)) to the balance of linear momentum (see the Introduction to this section), we obtain:

\[ \underbrace{\sum_{Ie}^{Ne} \int_{\Omega_e} \rho N_I(x) \ddot{u}^*_i d \Omega_e}_{\mathbf{M}\ddot{\mathbf{u}}^*} - \underbrace{\sum_{Ie}^{Ne} \int_{\Omega_e} \frac{\partial N_I(x)}{\partial x_j}\sigma_{ij}(u^*_i) d \Omega_e}_{\mathbf{F}^{int}} - \underbrace{\sum_{Ie}^{Ne} \int_{\Omega_e} N_I(x) \rho b_i d \Omega_e}_{\mathbf{F}^{ext}} \approx R_I(u^*(x),x). \label{eq:weighted-integral3} \]

Therein, \(\mathbf{M}\ddot{\mathbf{u}}^*\) are the inertia forces, \(\mathbf{F}^{int}(\mathbf{u}^*)\) the internal and \(\mathbf{F}^{ext}\) the external forces, respectively. In fact, Eq. (\(\ref{eq:weighted-integral3}\)) is a system of equations, with one equation per degree of freedom/unknown. In the present example, the unknowns are the displacements in x- and y-direction. Assuming the embankment is discretised using 500 nodes, the total number of unknowns would amount to 1000 resulting in a system of 1000 equations. The objective is now to find \(\mathbf{u}^*\) such that the system is in equilibrium, i.e. \(\mathbf{R}(\mathbf{u}^*)\) vanished. How appropriate \(\mathbf{u}^*\) are found is discussed in section Newton-Raphson method.