![]() |
|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|
|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Chapter 8
Theory and Implementation of the Finite Element Method
The
derivation and implementation of the finite element method outlined in the
previous chapter is simple and easy to follow, but it gives the misleading
impression that the finite element method relies on the principle of minimum
potential energy, and so is applicable only to linear elastic solids. This is not the case, of course
This chapter contains
In addition, a set of sample finite element codes (implemented in MAPLE and MATLAB) are provided to illustrate the how the various finite element procedures are implemented in practice.
8.1 Generalized FEM for static linear elasticity
This section gives a more general derivation and implementation of the finite element method for static linear elastic solids than the energy-based derivation given in Chapter 7.
8.1.1 Review of the principle of virtual work
Governing equations: We begin by summarizing the usual governing equations of linear elasticity, which must be solved by the FEA code. Â
Given: 1. The shape of the solid in its unloaded condition 2. The initial stress field in the solid (we will take this to be zero in setting up our FEM code) 3. The elastic constants for the solid 4. The thermal expansion coefficients for the solid, and temperature distribution (we will take this to be zero for our FEM code, for simplicity) 5. A body force distribution 6. Boundary conditions, specifying displacements
Calculate
displacements, strains and stresses 1. The strain-displacement equation 2. The elastic stress-strain law 3. The equation of static equilibrium for stresses 4. The boundary conditions on displacement and stress
The principle of virtual work: As we discussed in section 2.4, the principle of virtual work can be used to replace the stress equilibrium equations. Â
To
express the principle, we define a kinematically admissible virtual
displacement field
In addition, we define an associated virtual strain field
The
principle
of virtual work states that if the stress field
for
all possible virtual displacement
fields and corresponding virtual strains, it will automatically satisfy the
equation of stress equilibrium
8.1.2 Integral (weak) form of the governing equations of linear elasticity
The
principle of virtual work can be used to write the governing equation for the
displacement field in a linear elastic solid in an integral form (called the
`weak form’). Instead of solving the
governing equations listed in the preceding section, the displacements,
strains and stresses 1. Find a displacement field
for all
virtual velocity fields 2. Compute the strains from the definition 3. Compute the stresses from the stress-strain law The stress will automatically satisfy the equilibrium equation and boundary conditions, so all the field equations and boundary conditions will be satisfied.
The significance of this result is that it replaces the derivatives in the partial differential equations of equilibrium with an equivalent integral, which is easier to handle numerically. It is essentially equivalent to replacing the equilibrium equation with the principle of minimum potential energy, but the procedure based on the principle of virtual work is very easily extended to dynamic problems, other stress-strain laws, and even to problems involving large shape changes.
Derivation: start with the virtual work equation
Recall that
Finally,
recall that
Substituting into the virtual work equation gives the result we need.
8.1.3 Interpolating the displacement field and the virtual velocity field
To
solve the integral form of the elasticity equations given in 8.1.2, we discretize
the displacement field. That is to
say, we choose to calculate the displacement field at a set of n
discrete points in the solid (called `nodes’ in finite element
terminology). We will denote the
coordinates of these special points by
The displacement field at an arbitrary point within the solid will be specified by interpolating between nodal values in some convenient way. An efficient and robust implementation of the finite element method requires a careful choice of interpolation scheme, but for now we will denote the interpolation in a general way as
Here,
x denotes the coordinates of an arbitrary point in the solid. The interpolation functions
for all b=1…N. (This is to make sure that the displacement field has the correct value at each node). Recently developed meshless finite element methods use very complex interpolation functions, but the more traditional approach is to choose them so that
The simple constant strain triangle elements introduced in 7.1 are one example of this type of interpolation scheme. We will define more complicated interpolation functions shortly.
We can obviously interpolate the virtual velocity field in exactly the same way (since the principle of virtual work must be satisfied for all virtual velocities, it must certainly be satisfied for an interpolated velocity field…) so that
where
8.1.4 Finite element equations
Substituting the interpolated fields into the virtual work equation, we find that
where summation on a and b is implied, in addition to the usual summation on i,j,k,l.Â
Note that the interpolation functions are known functions of position. We can therefore re-write the virtual work equation in matrix form as
where
Here
K is known as the `stiffness matrix’ and f is known as
the force vector. K is a function
only of the elastic properties of the solid, its geometry, and the
interpolation functions and nodal positions. It is therefore a known
matrix. Similarly, f is a
function only of the known boundary loading and body force field, and the
interpolation scheme and nodal positions.Â
Observe that the symmetry of the elasticity tensor implies that
K also has some symmetry
The
virtual work equation must be satisfied for all possible sets of
This is a system of n linear equations for the n nodal displacements.
8.1.5 Simple 1D Implementation of the finite element method
Before describing a fully general 3D implementation of the finite element method, we will illustrate all the key ideas using a simple 1-D example. Consider a long linear elastic bar, as shown in the picture. Assume 1. The bar has shear modulus 2. The bar has cross section 3. It is constrained on all its sides so that 4. The bar is subjected to body force 5. The bar is either loaded or constrained at its ends,
so that the boundary conditions are eitherÂ
For this 1-D example, then, the finite element equations reduce to
where
We could obviously choose any interpolation scheme, evaluate the necessary integrals and solve the resulting system of equations to compute the solution. It turns out to be particularly convenient, however, to use a piecewise-Lagrangian interpolation scheme, and to evaluate the integrals numerically using a Gaussian quadrature scheme.
To
implement the Lagrangian interpolation scheme, we sub-divide the region
Generic
linear and quadrilateral 1-D elements are illustrated in the table. The local nodes on the element are numbered
1 and 2 for the linear element, and 1,2,3 for the quadratic element as
shown. We suppose that the element
lies in the region
where                             Â
where
Next, we need to devise a way to do the integrals in the expressions for the stiffness matrix and force vector. We can evidently divide up the integral so as to integrate over each element in turn
where
for
each element, which depend on the geometry, interpolation functions and
material properties of the element.Â
The first and last elements have additional contributions to the
element force vector from the boundary terms
Finally
we need to devise a way to compute the integrals for each element stiffness
matrix. It is convenient to map the
domain of integration to [-1,+1] and integrate with respect to the normalized
coordinate
where
Note that the mapping also enables us to calculate the shape function derivatives in the element stiffness matrix as
Finally, note that integrals may be computed numerically using a quadrature formula, as follows
whereÂ
8.1.6 Summary of the 1D finite element procedure
To summarize, then, the finite element solution requires the following steps: 1. For each element, compute the element stiffness matrix as follows:
where
and the integration points 2. Assemble the contribution from each element to the global stiffness
3. Similarly, if there is a non-zero body force, then compute for each element
and assemble the global force vector
4. Add contributions to the force vector from
prescribed traction boundary conditions at
where
the 5. Modify the stiffness matrix to enforce the constraints
6. Solve the system of linear equations
for
the unknown displacements
8.1.7 Example FEM Code and solution
A simple example MAPLE code for this 1-D example can be found in the file FEM_1D_Static.mws
It is set up to solve for displacements for a bar with the following parameters: 1. Length L=5, unit x-sect area, 2. Shear modulus 50, Poisson’s ratio 0.3, 3. Uniform body force magnitude 10, 4. Displacement u=0 at x=0 5. Traction t=2 at x=L.
The code computes the (1D) displacement distribution in the bar. The predicted displacement field is plotted on the right.
Of
course in general we want to calculate more than just displacements
This works well for a uniform body force with quadratic (3 noded elements) as the plot on the right shows.
However, if we switch to linear elements, the stress results are not so good (displacements are still calculated exactly). In this case, the stress must be uniform in each element (because strains are constant for linear displacement field), so the stress plot looks like the figure to the right.
It is interesting also to examine the stiffness matrix (shown below for 3 linear elements, before addition of the u=0 constraint for the first node )
Notice
that stiffness is symmetric, as expected, and also banded. A large FEM matrix is sparse
8.1.8 Extending the 1D finite element method to 2 and 3 dimensions
It is straightforward to extend the 1-D case to more general problems. All the basic ideas remain the same.  Specifically 1. In both 2D and 3D we divide up our solid of interest into a number of elements, shown schematically for a 2D region in the picture on the right. 2. We define interpolation functions
where 3. We introduce an element stiffness matrix for each element by defining
where 4. The volume integrals over each element are
calculated by expressing the volume or surface integral in terms of the
dimensionless coordinates
Here, 5. The global stiffness matrix
is then computed by summing the contribution from each element as
6. The stiffness matrix is modified to enforce any prescribed displacements 7. The system of equations
is solved for the unknown nodal displacements. 8. The stresses and strains within each element are then deduced.
To implement this procedure, we must (a) Define the element interpolation functions; (b) Express the integrals for the element stiffness matrices and force vectors in terms of normalized coordinates; (c) Formulate a numerical integration scheme to evaluate the element stiffness matrices and force vectors. These details are addressed in the sections to follow.
8.1.9 Interpolation functions for 2D elements
The 2D interpolation functions listed below are defined for the region
The numbers shown inside the element show the convention used to number the element faces.
                                          Â
8.1.10 Interpolation Functions for 3D elements
The 3D interpolation functions
listed below are defined for the region                                        Â
    The element faces are numbered as follows.
8.1.11 Volume integrals for stiffness and force in terms of normalized coordinates
In
this section we outline the procedure that is used to re-write the integrals
for the element stiffness and force in terms of the normalized coordinates
To evaluate them, we need to 1. Find a way to calculate the derivatives of the shape
functions in terms of 2. Map the volume (or area) integral to the region
Calculating the shape function derivatives. The shape function derivatives can be evaluated by writing
where
the derivatives
where
Note that
Mapping the volume integral: To map the region of integration we define
where
the matrix
We note in passing that the boundary integral in the element force vector can be regarded as a 1-D line integral for 2D elements and a 2D surface integral for 3D elements. So the procedures we developed in 8.1.5 (1D elements) can be used to evaluate the surface integral for a 2D element. Similarly, the procedures we develop to integrate stiffness matrices for 2D elements can be used to evaluate the surface integral for a 3D element.
8.1.12 Numerical integration schemes for 2D and 3D elements
Finally, to evaluate the integrals, we once again adopt a quadrature scheme, so that
The
integration points
Choosing the number of integration points: There are two considerations. If too many integration points are used, time is wasted without gaining any accuracy). If too few integration points are used, the stiffness matrix may be singular, or else the rate of convergence to the exact solution with mesh refinement will be reduced. The following schemes will avoid both
There are situations where it is preferable to use fewer integration points and purposely make the stiffness singular.  These are discussed in more detail in Section 8.5.
8.1.13 Summary of formulas for element stiffness and force matrices
With these definitions, then, we write the element stiffness matrix as
where
8.1.14 Sample 2D/3D linear elastostatic FEM code
You can find a MAPLE implementation of a simple 2D/3D static linear elasticity code in the file   FEM_2Dor3D_linelast_standard.mws
The code reads an input file. Several examples are provided : 1. Linear_elastic_triangles.txt: Simple 2D plane strain problem with two triangular elements 2. Linear_elastic_quad4.txt: Simple 2D plane strain problem with eight 4 noded quadrilateral elements 3. Linear_elastic_quad8.txt: Simple 2D plane strain problem with two 8 noded quadrilateral elements 4. Linear_elastic_brick4.txt: Simple 3D problem with 8 noded brick elements
As
an example, we show how to run the program with the first input file. The
file sets up the problem illustrated in the figure above. The elements are linear elastic plane strain with
The program input file is listed on the right. Here is a brief explanation of the data in the file 1. The first part of the input file specifies material properties. A number ‘1’ on the Plane strain/stress line indicates a plane strain analysis; a number ‘0’ indicates plane stress. 2. The second part specifies properties and coordinates of the nodes. For a 2D problem each node has 2 coordinates and 2 DOF; for a 3D problem each node has 3 coordinates and 3DOF. Then enter nodal coordinates for each node. 3. The third part lists the element
properties. Here, you must specify the
number of elements, and the maximum number of nodes on any one element (you
can mix element types if you like).Â
Then you must specify the nodes connected to each element (known as
element connectivity).  For each
element, you must specify the number of nodes attached to the element; an
identifier that specifies the element type (you can enter any number in this
version of the code 4. The fourth part of the file specifies boundary constraints. For any constrained displacements, enter the node number, the displacement component to be prescribed, and its value. 5. The last part of the file specifies
distributed loading acting on the element faces. The loading is assumed to be uniform. For each loaded boundary, you should
specify the element number, the face of the element (the face numbering
convention was described in section 7.2.9 and 7.2.10
Note that the program performs absolutely no error checking
on the input file. If you put in a typo,
you will get some bizzarre error message from MAPLE
For the input file shown, the program produces an output file that looks like this The code prints the displacements at each node in the mesh, and also the strains and stresses at each integration point (where these quantities are most accurate)Â for each element.
To run the code, you must complete the following steps 1. Open the maple executable file; 1. Edit the code to insert the full path for the input file in the line near the top of the code that reads > # Change the name of the file below to point to your input file > infile :=fopen(`D:/fullpathoffile/Linear_elastic_triangles.txt`,READ): 2. Scroll down near the bottom to the line that reads > #     Print nodal displacements, element strains and stresses to a file > # > outfile := fopen(`path/Linear_elastic_triangles.out`,WRITE): and enter a name for the output file. 3. Return to the top of the file, and press <enter> to execute each MAPLE block. You will see the code plot the undeformed and deformed finite element mesh at the end. The stresses and strains in the elements are printed to the output file.
|
|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
(c) A.F. Bower, 2008 |