A basic practical introduction to the Finite Element Method

1
A basic (practical) introduction to the Finite
Element Method (FEM)
2
What is the Finite element method?
The finite element method (or as it is commonly
abreviated to as FEM) is a numerical technique in
which the problem domain is discretionalized into
a series of interconnected regions called
elements which are joined together to form a mesh
which represents the computational volume. In the
finite element method it is assumed that there is
some continuous function (which represents the
solution) that exists over the entire
computational region which we wish to
solve. Originally the finite element method was
derived by accident to solve structaul problems
before it was realised that it could be used to
solve differential and partial differentail
equations of any system. FEM is now widely
accepted and used in may fields to solve complex
systems and there are a plethera of specialised
finite elements that have been developed over the
years to tackle various numerical obsticals.
3
Why use the Finite element method?
1) Meshing of a region using isoparametric
elements, a non fixed grid and adaptive
meshing 2) Although more complex to impliment
than finite difference methods the extension of
finite elements to higher order terms (as many as
you want) is straightforward. 3) Specialised
elements can be constructed and may be talior
made to fit your specific problem. General
numerical instabilities can be overcome through
the use of specialised elements. 4) Even though
using a talior series expansion of simple linear
finite difference and finite elements for a
generic system would suggest that there is
little difference between the two schemes (over a
uniformly structured mesh) except the
distribution of the source terms is somewhat
different however the extension to higher order,
arbitary meshing and the incorporation of
specialised functions/boundaries and continuity
equations gives finite elements wide scale
versatility and greater accuracy.
4
FEM aspects
There are four main aspects to finite elements on
the whole 1) The formulation and method used to
construct the system of equations that you
desire to be solved 2) The type of finite
elements that are used generic, specialised
elements such as edge elements/serendipty
elements, elements with singular functions, high
order C1, C2 interpolation functions.... to deal
with higher order derivative terms
etc.... 3) Numerical pitfalls or limitations of
the chosen elements/formulation always be
aware of the limitations of the
elements/formulation/system you
use..... 4) Meshing the chosen meshing scheme
that you use to be it an algebraic, manual or an
automatic one. Is the mesh going to be a
refineable mesh (either more elements are added
or a higher order function is used until some
tolerance is reached) or an adaptive one (a new
better mesh is constructed based on the old one
until some tolerance is reached).
5
FEM aspects continued...
The three most popular methods used to derive a
FEM formulation are 1) The Galerkin wieghted
residual method 2) The Variational method 3) The
Rayleigh-Ritz method In electromagnetic
calculations broadly speaking you will come
across the following specialized finite elements
that have been used to overcome certain numerical
problems serendipity elements, edge elements,
singular basis elements, upwinded elements. In
electromagnetic cavity problems one of the very
basic problems that is numerically encountered is
spurious solutions which occur form not enforcing
a divergence free condition on the field of
interest properly.
6
Spurious solution example
Traditional nodal elements can give rise to
spurious solutions if not properly formulated
(with say a penalty factor formulation or some
other such technique) edge elements on the
other hand are specially constructed in such a
way that the divergence free condition is
mathmatically guaranteed however edge elements
have their own pitfalls and are not the panacia
of all ills.....
Vector plot showing the correct (physical)
solution to an eigen modal problem suggested by
R.Crowley on the left. On the right a nodal based
finite element solution in which the divergence
of the field has not been enforced in which
spurious (unphysical) solutions are present.
7
An introductory 1D FEM problem that can be done
by hand
We wish to solve, as an introductory example, the
1D equation of
In the range of
Subject to the boundary conditions of
The analytical solution to the above problem is
8
Deriving the weak formulation
We will use the Galerkin method of weighted
residual to obtain a FEM formulation, in which we
will use FEM elements that have C0 continuity.
First we multiply both sides of the equation by a
test function w, which represents the sort
after numerical solution to the problem. Then we
rearrange the equation such that all terms are on
one side
We then integrate the above equation of the range
of interest and use integration by parts to
simplify the problem into derivative terms no
larger than the order of 1, so that C0 continuity
FEM elements can be used.
The above equation is the weak FEM form of the
problem. The terms on the LHS represent the
system matrix terms which are a representation of
the equations to be solved and the terms on the
RHS are the system vector which represent the
boundary conditions of the problem.
9
The FEM test function
We will use a linear 1D FEM element, which will
have the following form for the test function u
in terms of the spatial coordinate x
In which the terms a and b are constants.
Let us consider a typical FEM element in terms of
the above test function u. The 1D FEM element
will consist of two nodes, in which the solution
u within an element can be graphically
represented as
In the above diagram xi and xi1 are two adjacent
nodes used in the elements construction which
have the solutions of ui and ui1 respectively.
10
The FEM test function continued
We will use a linear 1D FEM element, which will
have the following form for the test function u
in terms of the spatial coordinate x
The individual solutions ui and ui1 can
therefore be written in terms of the spatial
coordinate and the chosen FEM test function as
Using these two equations the constants a and
b can be derived in terms of the solutions ui
and ui1 as
11
The FEM test function continued..
A general expression for the solution u of at
any node for any element may be derived in terms
of the individual solutions ui and ui1 by
substituting the equations for a and b into
the general expression for the solution u (i.e.
uaxb) as
In the above equation H1 and H2 are
where hi is defined as
The derivative of u can be derived in terms of
the derivatives of H1 and H2. The derivatives of
H1 and H2 with respect to x are therefore
derived as
12
Writing the matrix FEM formulation
Given that the weak formulation of the problem
has been derived as
In which the solution u (across an element in
terms of the nodal solutions) has been defined in
terms of a linear C0 FEM element as
Since the trial function w is assumed to be the
approximate solution to the problem (i.e. wu),
then setting wu and substituting the equation
above produces a linear matrix form of the
problem i.e.
13
Expanding the matrix
Now that we have the matrix form we must expand
and fully integrate it to form the linear matrix
problem, which we can solve. As an example of how
to do this lets expand the first term on the LHS
i.e.
14
The matrix for a single element
Expanding and integrating all the terms on the
LHS of the equation results in the matrix for a
single element as
It is now possible to extend the above analysis
for a single element to multiple elements
15
Adding more elements
Let us assume that the range of interest is split
up into two adjoining elements of Element1 (in
the range of 0ltxlt0.5) and Element2 (in the
range of 0.5ltxlt1) as
Performing the same process as before we can
construct matrices for each element. Note that in
this case this will simply mean applying the
previously determined elemental matrix we derived
in terms of the variable hi to each element,
which results in
16
The matrices for each element
Performing the same process as before we can
construct matrices for each element. Note that in
this case this will simply mean applying the
previously determined elemental matrix we derived
in terms of the variable hi to each element,
which results in
The element matrix for Element1 as
and the element matrix for Element2
17
The system matrix and system vector
Directly adding the expressions for the first and
second elements results in what is known as the
system matrix. Combining the natural boundary
conditions, which have arisen from the derivation
of the weak form, into the vector on the LHS of
the overall equation results in what is called
the system vector.
Substituting the known Boundary conditions, which
in this case happen to be Dirichlet boundary
conditions, results in the following matrix
system
18
Solving the system
The system can be written more compactly as
In which K is the system matrix, F is the
system vector and u is the sort after solution.
This linear algebra system can then be solved
using whatever linear algebra solver technique
you desire (provided the system meets the
requirements for the chosen solver technique)
be it a direct method or an iterative
one. Regardless of what ever method you choose
the effective solution which you will obtain for
the above system is
19
The Solution
A comparison between the solved two element
numerical FEM and the analytical answer are
The difference between the numerical and the
analytical answer is approximately 0.61 for a
two element system. If we used an infinite number
of elements the numerical answer would exactly
equal the analytical answer effectively as the
number of elements increases the numerical
solution will approach that of the true
analytical answer.
20
Isoparametric elements
The previous example used analytically calculated
trial functions this approach is rather
cumbersome and limited in scope. A better and
widely used approach is to use isoparametric
elements, in which the trial functions can be
calculated numerically using the Legendre-Gauss
integration. Isoparametric elements are
integrated in a normalised coordinate system
(ususally between the dimensionless points of -1
and 1), which has been derived from a physical
coordinate system.

• Isoparametric elements have a number of distinct
  advantages
• They have great versatility representing curved
  surfaces.
• They can be deformed to some extent (rule of
  thumb never squash an elements Jacobian more
  than a factor of two) and still give very
  accurate answers.
• It is realtivily straight forward to derive and
  use higher order polynomial test functions for
  use in these type of elements.

21
The isoparametric version of the introductory 1D
FEM
Let revisit the introductory test problem of
In the range of
Subject to the boundary conditions of
The analytical solution to the above problem is
22
The isoparametric FEM test function
This time we will use linear isoparametric
elements to solve the problem. Let us consider a
1D FEM element which has a natural coordinate
system described by the variable e which can be
mapped (translated) into a physical coordinate
system describe by the actual nodal coordinate
values of x.
Once again the subscript of i represents the
nodal position/number and the actual solution to
the problem is represented by the variable u.
23
The isoparametric FEM test function
Any point between e-1 and e1 in the actual
coordinate system can be mapped into the physical
coordinate system using the following linear
shape function
The physical coordinate system can then be
described in terms of the natural coordinate
system e as
24
The isoparametric FEM test function continued
The shape functions are also used to
interpolate/describe the solution u within an
element in much the same way as before.
25
The 1D problem revisited using isoparametric
elements
Lets examine our previous example of
With the boundary conditions of
Using the Galerkin weighted residuals method (for
elements with C0 continuity) in which the problem
is multiplied by a test function w and then
integrated using integration by parts to derive
the weak form (given as before) as
26
The Jacobian
The integration in the weak form is expressed in
terms of the physical coordinate x, but we want
to express this equation in terms of the
isoparametric natural coordinate system e
i.e. To do this we multiply the above equation by
a Jacobian J, which is a coordinate transform
described as
So the weak form written in terms of the natural
coordinate system becomes
27
Writing the isoparametric matrix
Since the trial function w is assumed to be the
approximate solution to the problem (i.e. wu),
then setting wu and substituting the equation
above produces a linear matrix form of the
problem i.e.
In which the solution u (across an element in
terms of the nodal solutions) has been defined in
terms of a linear C0 FEM element as
and
28
Computing the Jacobian
The derivative of the shape function in terms of
the physical coordinate system i.e. Are
derived by applying the inverse of the Jacobian
to the derivatives of the terms in tehr natural
coordinate system i.e.
The Jacobian is computed as
The inverse of the Jacobian is therefore
29
Expanding the isoparametric matrix
Putting all the terms into the matrix form of the
weak formulation (presented above), integrating
and expanding in terms of e results in the
following matrix for the LHS of the equation i.e.
the system matrix component
30
The isoparametric system matrix
Performing the integration from -1 to 1 results
in the following system matrix in terms of the
physical coordinate system
Which is exactly the same result as we derived
before by manually doing the integration in terms
of the physical coordinate system. Using
isoparametric elements allows one stress free FEM
programming, in which the results can easily be
extended to higher order test functions and
multiple dimensions.
31
Writing a FEM Code
Following the aforementioned logic the
construction of an FEM code (for any dimensional
system) is realtively straightforward. The
process can generally be split into three
distinct phases
1) The generation of a mesh and the construction
of the elements within the mesh. 2) Numerical
integration loops over all dimensions to
integrate the weak form of the system of
equations that are being solved. 3) Solution of
the final matrix system.
The above process is somewhat over simplified (to
the extreme) as there are of course rules and
regulations which must be enforced for certain
aspects of the numerical scheme..... as well as
rule of thumb aspects aswell....
32
A straightforward example
Shown below is the output from a simple 2D
axisymmetrical FEM code (written in Octave) which
solves a Skin depth problem using cubic
serendipity elements.
33
Extension of the simple example and meshing
aspects
As most of us here in the wakefield interest
group work with cavity problems, which are in
effect specialized pillboxes, it would of course
be nice to write a simple pillbox example. The
eigen modal pillbox problem is but a mere
extension of the weak form of the skin depth
problem and it makes a good benchmarking
exercise. The pillbox problem can be derived by
reworking the formulation (of Maxwell's
equations), the meshing and the solution of the
system of equations .