### Louis Moresi

Professor of Geophysics Australian National University Canberra, Australia +61 4 0333 1413

Email Twitter Facebook Google+ Github Youtube Researchgate

# Computational Geodynamics: Finite Elements 2

$\newcommand{\dGamma}{\mathbf{d}\boldsymbol{\Gamma}} \newcommand{\erfc}{\mbox{\rm erfc}} \newcommand{\curly}{\sf } \newcommand{\Red }[1]{{ #1}} \newcommand{\Green }[1]{{ #1}} \newcommand{\Blue }[1]{{ #1}} \newcommand{\Emerald }[1]{{ #1}}$

### Weak Forms of Some Useful Equations

As we discovered earlier, the Galerkin form of weighted residual approximate method can be made equivalent to a full variational problem if integration by parts is employed in the correct manner. This then allows the development of fully automatic variational methods for arbitrary strong forms of the governing equations. From here on we follow the notation of Hughes (The Finite Element Method) which is reasonably clear.

#### 1D Heat Conduction

The differential or strong form of the equation we have encountered a number of times is $\begin{equation} \begin{split} \frac{d^2 u}{d x^2} - {\curly f} &= 0 \;\;\; \text{on} \;\;\; \Omega u(1) & = {\curly g} -\frac{d u}{d x} (0) &= {\curly h} \end{split} \end{equation }$ where the boundary conditions are supplied at either end of a domain of unit length.

To find the weak form of the equation we need to have a set of functions $\curly S$ to use for the trial functions which satisfy the boundary condition $\curly g$ on $u$ at $x=1$. We also need a set of weighting functions which are zero at $x=1$ which we call $\curly V$. The weighting functions play the role of the variations.

The statement of the weak form of the problem (as distinct from its solution which will wait till later) is then to find $u \in {\curly S}$ such that for all $w \in {\curly V}$ $\begin{equation} {\int_0^1 \frac{d w }{d x} \frac{d u}{d x} dx - \int_0^1 w f dx + w(0) h = 0} \end{equation}$

This looks a lot like the variational solutions we found before. How do we get to it ? The procedure is relatively general (except that in higher dimensions it becomes more time consuming). First assume we have found the solution to the strong form, $u$. This must satisfy $\begin{equation} 0 = \int_0^1 w( \frac{d^2 u}{d x^2} - {\curly f} )dx \end{equation}$ by the definition of the problem, as long as $w$ is well behaved, which is ensured by a sensible choice of weighting functions. This starts to look a lot like the Galerkin approach although we are not yet seeking an approximate solution. Now integrate by parts to give $\begin{equation} 0 = \int_0^1 \frac{d w }{d x} \frac{d u}{d x} dx - \int_0^1 w f dx - \left[ w \frac{d u}{d x} \right]_0^1 \end{equation}$ the boundary conditions on $du / dx$ are substituted to give the weak form as above. We write the equation in terms of the symmetric operators: $\begin{equation} \begin{split} a(w,u) & \equiv \int_0^1 \frac{d w }{d x} \frac{d u}{d x} dx (w,{\curly f}) & \equiv \int_0^1 w f dx \end{split} \end{equation }$ to obtain an abstract form which can be used for all the finite element formulations we derive here $$$a(w,u) = (w,{\curly f}) + w(0){\curly h} \label{eq:FEabst}$$$

Note, that the partial differentiation step produces a symmetric operator from the non symmetric initial form. This approach also works in 2D as follows.

#### 2D Heat Conduction

Define a heat flux (vector), $\mathbf{q}$ related to temperature, $u$ as follows $\begin{equation} q_i = -K_{ij} \frac{\partial u}{\partial x_j} \end{equation}$ where $\mathbf{K}$ is the symmetric conductivity tensor. (This is the generalized version of Fourier’s law). A volumetric heating rate of $\curly f$ gives the following strong form $\begin{equation} \begin{split} \nabla q - {\curly f} &= 0 \;\;\; \text{on} \;\;\; \Omega u & = {\curly g} \;\;\; \text{on} \;\;\; \Gamma_{\curly g} -q_i n_i &= {\curly h} \;\;\; \text{on} \;\;\; \Gamma_{\curly h} \end{split} \end{equation }$ where $\Gamma_ {\curly g}$ and $\Gamma_{\curly h}$ are the regions of the boundary with fixed temperatures and fixed fluxes respectively, $\mathbf{n}$ is the unit normal over $\Gamma_{\curly h}$.

The corresponding weak form (found in the usual way) is $\begin{equation} -\int_\Omega \frac{\partial w}{\partial x_j} q_i d\Omega = \int_\Omega w {\curly f} d\Omega + \int_{\Gamma_{\curly h}} w{\curly h} d\Gamma \end{equation}$ The symmetry of the operator is not obvious here until we substitute the constitutive law to obtain $\begin{equation} \int_\Omega \frac{\partial w}{\partial x_j} q_i d\Omega = \int_\Omega \frac{\partial w}{\partial x_j} K_{ij} \frac{\partial u}{\partial x_j} d\Omega = \int_\Omega (\nabla w)^T \mathbf{K} (\nabla u) d\Omega \end{equation}$

#### 2D/3D Fluid Flow

Now we finally come around to the problem we really need to solve. The strong form of the Stokes’ flow problem is identical to that of linear elasticity (which is what all the finite element textbooks deal with). The general constitutive law is $\begin{equation} \sigma_{ij} = c_{ijkl} \epsilon_{kl} \end{equation}$ which reduces to $\sigma_{ij} = \eta \epsilon_{ij}$ for homogeneous fluid. $\epsilon_{ij}$ is defined by $\begin{equation} \epsilon_{ij} = \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \end{equation}$ where $\mathbf{u}$ now represents the fluid velocity.

The strong form of the equation is $\begin{equation} \begin{split} \frac{\partial \sigma_{ij}}{\partial x_j} + {\curly f}_i = 0 \;\;\; \text{on} \;\;\; \Omega u_i & = {\curly g}_i \;\;\; \text{on} \;\;\; \Gamma_{ {\curly g}_i} \sigma_{ij} n_i &= {\curly h}_i \;\;\; \text{on} \;\;\; \Gamma_{ {\curly h}_i} \end{split} \end{equation }$ with a corresponding weak form $\begin{equation} \int_\Omega \left( \frac{\partial w_i}{\partial x_j} + \frac{\partial w_j}{\partial x_i} \right) \sigma_{ij} d \Omega = \int_\Omega w_i {\curly f}_i d\Omega + \sum_{d=1}^{n_{\rm dim}} \left[ \int_{\Gamma_{ {\curly h}_i}} w_i h_i d\Gamma \right] \end{equation}$ Which is symmetric in the trial functions and the unknown velocities since $\begin{equation} \int_\Omega \left( \frac{\partial w_i}{\partial x_j} + \frac{\partial w_j}{\partial x_i} \right) \sigma_{ij} d \Omega = \int_\Omega \left( \frac{\partial w_i}{\partial x_j} + \frac{\partial w_j}{\partial x_i} \right) c_{ijkl} \left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right) d \Omega \end{equation}$

### Galerkin Approximate Weak Form

We have one last task before we can develop a real numerical method. How do we get away from the highly abstract notion of trial functions, weighting functions, and weak forms to produce a compact system of equations which can be written out in a nice tidy matrix form. The answer lies in the Galerkin approximate solution method. Let us go back to the 1D problem for a moment. We have to express our solution in terms of a sum over basis functions as before $\begin{equation} w^h = \sum_{A=1}^{n} c_A N_A \end{equation}$ where the $h$ superscript indicates that we have moved from an infinite choice of functions to a finite one based on the discrete functions $N_A$ which we will define in a minute, although for now we need to ensure that they are zero where the boundary conditions apply in the same way as $\curly V$. We define a corresponding expansion for $u$ $\begin{equation} \begin{split} u^h &= \sum_{A=1}^{n} d_A N_A + {\curly g}N_{n+1} {\curly g}^h &= {\curly g} N_{n+1} \end{split} \end{equation }$ which includes an extra basis function to allow the boundary condition ($u(1)={\curly g}$) to be satisfied.

Substituting into (\ref{eq:FEabst}) gives $\begin{equation} a\left( \sum_{A=1}^n c_A N_A , \sum_{B=1}^n d_B N_B \right) = \left( \sum_{A=1}^n c_A N_A , {\curly f} \right) + \left[ \sum_{A=1}^n c_A N_A(0) \right] {\curly h} -a\left( \sum_{A=1}^n c_A N_A , {\curly g} N_{n+1} \right) \end{equation}$ Our operators are symmetric and wholly linear, so the order of summation and application of the operator can be interchanged freely to give: $\begin{equation} 0 = \sum_{A=1}^n c_A G_A \end{equation}$ where $\begin{equation} G_A = \sum_{B=1}^n a(N_A,N_B)d_B - (N_A,{\curly f}) - N_A(0){\curly h} + a(N_A,A_{n+1}){\curly g} \end{equation}$ As we have done many times before, we appeal to the argument that the particular choice of variation must be totally arbitrary: this equation must hold no matter what $w^h \in {\curly S}^h$ we choose, and hence no matter what the combination of $c_A$ may be. This then means that $G_A$ must be identically zero independent of $c_A$, i.e. $$$\sum_{B=1}^n a(N_A,N_B)d_B = (N_A,{\curly f} ) - N_A(0){\curly h} + a(N_A,A_{n+1}) {\curly g} \label{eq:gal2}$$$

As in the standard variational method, we have eliminated all references to the actual variation and left unknowns which only related to the physical variables (i.e. the coefficients of the expansion for $u$).

If we simply write $\begin{split} K_{AB} &= a(N_A,N_B) F_A &= (N_A,{\curly f}) + N_A(0)h - a(N_A,N_n+1){\curly g} \end{split} \label{eq:FEstdform}$ then we have a matrix formulation immediately since (\ref{eq:gal2}) now becomes $\begin{equation} \sum_{B=1}^n K_{AB} d_B = F_A \;\;\;\;\;\; A=1,2,3,\ldots \end{equation}$ or $\begin{equation} \mathbf{K d} = \mathbf{F} \end{equation}$ The matrix $\mathbf{K}$ is known as the stiffness matrix — the association with the matrix of stiffnesses from our discrete set of springs being obvious.

### Generalization

The same argument can be applied to higher dimensions, and to problems with different constitutive laws, and vector unknowns, however, the identification of the components of the stiffness matrix with the operator $a(\cdot,\cdot)$ acting on every possible combination of the galerkin basis functions still holds. When the unknowns are vectors, the basis functions are used in each respective direction which complicates the notation even more than before. Essentially, though it means that the entries to $\mathbf{K}$ from each $A$ and $B$ are actually matrices which are $n_{\rm dim} \times n_{\rm dim}$.

For example, in the constant viscosity Stokes’ flow problem, we now have $\left. K_{AB}\right|{ij} = \int\Omega \left( \frac{\partial N_A}{\partial x_j} + \frac{\partial N_A}{\partial x_i} \right) \eta \left( \frac{\partial N_B}{\partial x_j} + \frac{\partial N_B}{\partial x_i} \right) d \Omega \label{eq:festokes}$

### Discretization & Shape Functions

Although we have now nominally made our problem finite by expressing everything in terms of Galerkin basis functions, we have yet to decide what those functions ought to be. Whatever our choice, it should make the problem easy to compute numerically.

The choice of basis functions is very broad but we narrow it down by tying it very closely to the way we choose to split up the domain of our problem. With an integral method, the domain decomposition is simple since $$$\int_\Omega d\Omega = \int_{\Omega_1} d {\Omega_1} + \int_{\Omega_2} d {\Omega_2} + \int_{\Omega_3} d {\Omega_3} + \ldots \label{eq:decomp}$$$ which means that the decomposition shown above is as accurate as the representation of the boundary allows.

The functions which interpolate the node points at the corners of subdomains are a suitable set for use as an approximate basis as we require for the Galerkin formulation. This is easy to see if we consider a one dimensional problem.

In the one dimensional case, the choice of subdomains is limited to breaking up the line into a number of segments, though not necessarily of equal length. If the approximation to the continuum field, $u$ is made by linear interpolation to give $u^h$, then $u^h$ can also be expressed as a sum on the local triangular functions $N_A$

$\begin{equation} u^h = \sum_{A} u(x_A) N_A(x) \end{equation}$ This is an exact representation of the interpolation provided the $N_A$ takes the value one at node $A$ and zero at all other nodes, and varies in a linear manner in between.

Note that the functions all take the same form apart from scaling. Also, because the basis functions are localized around a node and its neighbours, the direct interaction between a node and its immediate neighbours is non-zero but the interaction with more distant nodes is zero. This makes the stiffness matrix banded — more importantly it is sparse and therefore the potentially enormous number of interactions (the square of the number of unknowns) is contained. This procedure can be extended to higher dimensions and to higher order interpolations as shown in [the Figure below][#shape-functions]

#### Naming of Things and What They are Like}

The interpolation functions are known as shape functions. The subdomains are known as elements. The elements correspond to the individual springs of our discrete example. The shape functions are pure interpolation functions of the required order within the element and can be differentiated the appropriate number of times. This means that the order of interpolation must match the order of derivatives in the FE operator. Crossing the element boundaries, the shape functions have discontinuous derivatives (as do the piecewise interpolations). Continuous derivative functions are possible (e.g. splines) but add significant complexity. Shape functions in 2D and 3D can be formed from the product of 1D shape functions in each direction.

#### Element Matrices

The domain decomposition illustrated above produces a minature version of the matrix problem on the individual elements themselves. $\begin{equation} {k^e}_{ij} {d^e}_j = {f^e}_i \end{equation}$ The interpretation of the local problem is similar to that of the individual springs in the discrete example we derived earlier. By assembling the individual under-constrained element problems into the full problem we obtain a soluble system with the correct number of boundary conditions etc. The local equilibrium problem is identical with that of the global problem. That is, we use the same variational form and the same operators to build the local matrices. In general, finite element methods are designed around the local element frame of reference and large problems assembled from smaller ones. The book-keeping involved is to track which degrees of freedom in the local problem are associated with which others in the global problem. Clearly the interaction coefficient between one degree of freedom and another is derived from the local equilibrium relationships of a number of elements because the shape function for a given nodal point is spread across a number of elements.

### Numerical Integration

The real power of finite elements comes from its ability to handle remarkably distorted geometries with ease. (Obviously this only applies in two or more dimensions). This is possible because simple mappings can be derived to convert a distorted element into a regular one (square, cubic, regular-tetrahedral). This takes the form of a change of variables which can be done {\em fully automatically} if we make minor restrictions on the choice of distortions we allow the elements to have. This restriction takes the form of ensuring that the mappings which transform the element shape can be described by the shape functions of the element. This means that only linear mappings can be applied to linear elements. If we wish to map a curved boundary exactly, then elements which have the appropriate order shape functions have to be used even if they are not needed for the representation of the differential operator. (e.g. quadratic elements for a circular boundary). This concept produces “isoparametric” elements.

The mapping shown in is achieved by a change of variables $\begin{equation} \int \int_{\Omega^e} \phi dx dy = \int \int_{\rm square} \phi {\curly j} d \xi d \eta \end{equation}$ where ${\curly j}$ is the jacobian of the transformation defined by $\begin{equation} {\curly j} = {\rm det} \frac{\partial \mathbf{x}}{\partial \boldsymbol{\xi}} = {\rm det} \left[ \begin{array}{cc} \frac{\partial x_1}{\partial \xi_1} & \frac{\partial x_1}{\partial \xi_2} \frac{\partial x_2}{\partial \xi_1} & \frac{\partial x_2}{\partial \xi_2} \end{array} \right] \end{equation }$ In the isoparametric concept, the components of the jacobian can be written as, for example,
$\begin{equation} \frac{\partial x_1}{\partial \xi_1} (\xi_1,\xi_2) = \sum_{a=1}^{n_en} \frac{\partial N_a}{\partial \xi_1} {x_1}_a \end{equation}$

Now we have the integral of a number of things defined over a regular, square domain (or cube etc). A number of ways to estimate these integrals is available. For example, in 1D the familiar trapezium rule can integrate linearly interpolated functions exactly. Given that our approximate method has already reduced the degree of polynomial within the element to something manageable, it is possible to find integration schemes which produce the {\em exact} integral for our approximation — in other words introducing no further error.

One possibility which is commonly used is Gaussian quadrature. This is a textbook method which can be made exact for any order of interpolation function desired. For linear interpolation in 2D, the following rule applies (assuming that we have transformed to a square element) $\begin{equation} \int_{-1}^{1} !! \int_{-1}^{1} \phi(\xi,\eta) d\xi d\eta \ \ \ \cong \ \ \ \sum_{l=1}^{n_{\rm int}} \phi( \tilde{ \xi_l}, \tilde{ \eta_l}) W_l \end{equation}$ In which $n_{\rm int}$ is the number of points in the quadrature rule with co-ordinates $(\tilde{ \xi_l}, \tilde{ \eta_l})$. Each point has a weight associated with it of $W_l$. For the four point rule:

### Stresses

Note Derivatives of the shape functions are commonly not defined at the node points and consequently it is not possible to compute stresses etc at the nodal points. They are defined on the interiors of the elements and must be extrapolated to the nodes using some form of best fit procedure.