Subsections

• 5.1.1 Mesh Aspects
• 5.1.2 Shape Function
• 5.1.3 Weighted Residual Method

5.1 Basics

The basic aim of the finite element method is to solve a PDE, or a system of coupled PDE s, numerically. Instead of finding the analytic solution of the PDE, which is usually a function of the coordinates, it is tried to determine this function values for discrete coordinates on grid points. For this purpose the continuum is discretized with a number of so-called finite elements which results in a mesh with grid nodes. If the finite elements are appropriately small, the solution of the PDE can be approximated with a simple function, the so-called shape function, in each element, which acts as a contribution to the approximation of the global solution of the PDE. All together a linear or non-linear equation system must be obtained. The real advantage is that such (non-)linear equation systems can be quite easily solved today by computer programs.

5.1.1 Mesh Aspects

For solving a PDE numerically, in the first step the simulation domain $\Omega$ is divided into a number of $M$ (as possible geometrically simple) elements $E_i$ (like tetrahedrons or cubes). The quality of the approximated solution depends on the size and so on the number of finite elements in the discretized domain. The elements must be located in such a manner, that there are no empty spaces between them and that they do not overlap:

$$\Omega=\bigcup_{i=1}^M E_i \qquad \mathrm{and} \qquad \forall i \neq j: Int(E_i) \cap Int(E_j)=0.$$ (5.1)

Here $Int(E_i)$ is the set of all points in the element $E_i$, except those which are located on the surface. Furthermore, surface conformity of neighbor elements is demanded. This means that on each surface of an element inside the domain exact only one neighbor element is bordered.

If $P$ is the set of all grid nodes of the discretized domain $\Omega$, than each grid node $p_k \in P$ has a unique global index $k=1,\dots,N$. $N$ is the number of all nodes in the whole mesh. A node $p_k$ can also have several local indices.

5.1.2 Shape Function

From the mathematical point of view the shape function shall interpolate the discrete solution function values between the grid nodes. If a PDE is written in the form

$$\mathcal{D}[u](x,y,z)+g(x,y,z) = 0,$$ (5.2)

where $\mathcal{D}$ is a second order differential operator, the desired solution of the PDE is a function $u(x,y,z)$, which is approximated by [85]

$$\tilde{u}(x,y,z)=\sum_{j=1}^N u_j N_j(x,y,z) = \{u^T\}\{N\}$$ (5.3)

The summation is performed over all grid nodes. Thereby, $u_j$ is the respective value of the function at a grid node with the number $j$ and $N_j(x,y,z)$ is the shape function.

The shape function can be chosen quite freely, but it must be appropriate for interpolation. In addition every shape function $N_j(x,y,z)$ must have a value of 1 on the grid node j (with coordinates $\vec{p_j}=\{x_j,y_j,z_j\}^T$) and values of 0 on all other grid nodes:

$$N_j(x_i,y_i,z_i)=\left\{ \begin{array}{ll} 1 & \textrm{for } \mathrm{i= j} 0 & \textrm{for } \mathrm{i\neq j} \end{array} \right.$$ (5.4)

In practical applications linear shape functions or polynomials of low order are used. The numbering of the grid nodes should be carried out in a way, that nodes at Dirichlet-boundaries, where the solution of $u(x,y,z)$ is already known, are ranked behind the others. The nodes where the value $u_j$ must be calculated get the indices $j=1,\dots,N_A$, the other $N_B$ nodes at the Dirichlet-boundaries are indexed with $j=N_A+1,\dots,N_A+N_B$.

In order to get the solution of the PDE, ``only'' the unknown coefficients $u_j$ in (5.3) must be obtained. For determining the values $u_j$ there are two available ways, the method of Ritz or the method of weighted residuals [86]. The method of Ritz seeks a stationary point of the variational functional [87]. This variational approach leads to an integral formulation of the PDE. The Ritz method can only be applied, if for the boundary problem an equivalent variational formulation exists. The more universal method is based on the weighted residuals and therefore it is used in this work.

5.1.3 Weighted Residual Method

When the approximated solution $\tilde{u}(x,y,z)$ is reinserted in the PDE (5.2), there exists a residual

$$\mathcal{R} = \mathcal{D}[\tilde{u}(x,y,z)]+g(x,y,z).$$ (5.5)

If the residual would disappear ( $\mathcal{R}=0$), the exact solution $\tilde{u}(x,y,z)= u(x,y,z)$ would have been found [88]. Since the function space of the approximation solutions $\tilde{U}$ is a subset of the function space of the exact solutions ${U}$, such a solution does not exist for the general case. Therefore, it is tried to fulfill the residual condition not exactly, but with $N$ weighted or averaged linearly independent weight functions (also called test functions) $W_i(x,y,z)$, so that [89]

$$\int_{\Omega}[W_i(x,y,z)\mathcal{D}[\tilde{u}(x,y,z)]+W_i(x,y,z)g(x,y,z)] d\Omega=0, \qquad \mathrm{for} \quad i=1,2,3,\dots N.$$ (5.6)

The weight functions $W_i$ must be chosen in a suitable way, because the quality of the solution depends on them. If the weight functions are identical with the shape functions ( $W_i(x,y,z)= N_i(x,y,z)$), this approach is called Galerkin's method [85].

With the approximated solution (5.3) the weighted residual (5.7) can be rewritten in the form

$$\sum_{j=1}^N u_j \int_{\Omega} W_i(x,y,z)\mathcal{D}[N_j(x,y,z)] d\Omega + \int_{\Omega}W_i(x,y,z)g(x,y,z) d\Omega=0.$$ (5.7)

It is possible to determine the $N$ unknown function values $u_j$ on the grid nodes with this system of $N$ equations.