EN222:
Mechanics of Solids
Division of Engineering
Brown University
3. Finite Element Solution of Boundary and Initial Value Problems
Problem Definition
We have collected together all the information we need to set up a boundary value or initial value problem in the mechanics of solids. To recapitulate, the standard problem in solid mechanics is:
Given:
1. A reference configuration for the solid (this is usually, but not always, taken to be the initial configuration)
2. A constitutive law for the solid (including mass density & as well as a stress-strain-temperature law)
3. The initial displacement field and velocity field
4. A history of loading and displacement boundary conditions, which may include a time varying body force field, as well as either boundary tractions on part of the boundary and boundary displacements on the remainder of the boundary ;
5. Constitutive laws for interfaces or contacting surfaces (we’ll address this later)
6. A way to determine the temperature distribution in the solid, if necessary. In complex problems (e.g. models of machining) the heat equation is coupled to the mechanical problem, and must be solved simultaneously.
Determine:
1. The displacement field as a function of time;
2. The stress field in the solid as a function of time.
The solution is found by solving the following equations:
together with the constitutive relation (this may include constitutive relations for interfaces or contacting surfaces)
Some problems are simpler, in that we are only interested in a steady-state or quasi-static solution. In this case, we need to be given all the information in items 1-6 above except item 2 – the initial conditions (we will need initial values of state variables in a history dependent constitutive law, however). The governing equations are then
A veritable arsenal of
mathematical and computational tools is available to solve these problems. But undoubtedly the most versatile tool is
the finite element method. In less than
a month, you can write yourself a finite element program that will solve
virtually any boundary or initial value problem you care to pose. In a few more months you can extend the
method to model contact; problems involving coupled field equations (e.g. heat
generation); and add fancy features like automatic mesh generation. Moreover, you can solve fully 3D problems
involving nonlinear geometry & material behavior on a cheap PC. This all sounds great – and it is – except
for two problems. One is that you can
solve such complex problems that you might spend years trying to figure out
what the solutions mean (or to choose the right problems to pose). The other is that a lot of the time the
bloody thing won’t converge. The big bucks
go to the people that can deal with these problems – that’s why there are so
many finite element consulting firms out there.
In this section of the course, we will present the essential features of the finite element method. Some subtle issues (like treatment of incompressible materials) will be left for future courses, and we will not discuss how to actually code the method efficiently.
We will start with the simplest FEM application and gradually add features. Specifically, we plan
1. To implement a 3D static linear elastic FEM code
2. To extend the code to solve dynamic problems
3. To extend the method to history dependent materials (e.g. linear viscoelasticity)
4. To extend the code to nonlinear materials and nonlinear geometry
5. To discuss implementations of interfaces and contact.
This section assumes some prior familiarity with FEM. If you are not at all familiar with the technique you should read the elementary introduction to the topic in the EN175 notes.
3.1 Finite Element Method for Static Linear Elasticity
A static linear elasticity problem is posed as follows.
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 acting on the solid (Note that in this section we will use b to denote force per unit volume rather than force per unit mass, to avoid having to write out the mass density all the time)
6. Boundary conditions, specifying displacements on a portion or tractions on a portion of the boundary of R
Calculate displacements, strains and stresses satisfying the governing equations of linear elastostatics
As we discussed in section 1.3, we can use the principle of virtual work as an alternative way to state the stress equilibrium equations. In addition, we can eliminate stress and strain from our governing equations by substituting the strain-displacement expression into the constitutive law, and then substituting the result into the principle of virtual work. As a result, our problem can be restated as
Calculate a displacement field satisfying
for all virtual velocity fields satisfying on . Once the displacement field has been computed, we can find the strain and stress using the governing equations.
To find an approximate solution, 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 , where the superscript a ranges from 1 to N. The unknown displacement vector at each nodal point will be denoted 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 are functions of position only, which must have the property that
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
We will discuss specific
interpolation schemes in more detail later.
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 are arbitrary nodal values of virtual velocity.
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 – specifically .
The virtual work equation must be satisfied for all possible sets of with for nodes a that lie on . At these nodes, the displacements must satisfy . Evidently, this requires
This is a system of n linear equations for the n nodal displacements.
Implementation of the finite element method
The basic principle behind the finite element method is evidently exceedingly simple. The key to its success is to devise an implementation that minimizes the computational cost associated with calculating and storing the stiffness matrix, and solving the resulting system of linear equations.
1-D example
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 bar with cross section and length L constrained on all its sides so that , subjected to body force , and boundary conditions specifying either the traction or displacement at x=0 and x=L. Assume furthermore that the bar is an isotropic, linear elastic solid with shear modulus and Poisson’s ratio
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 into a series of elements, as illustrated in the figure. Each element is bounded by two nodal points, and may also contain one or more interior nodes. Within each element, we interpolate the displacement field through the nodes on that element only. So, we would use a linear interpolation between the nodes on a 2-noded element, a quadratic interpolation between the nodes on a 3 noded element, and so on.
Generic linear and quadrilateral 1-D elements are illustrated below. 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 . The displacements within the element are then interpolated as
where denotes the number of nodes on the element, denotes the value of the displacement at each node, and the shape functions are
Linear 1-D element
Quadratic 1-D element
Of course, the actual nodal coordinates do not lie at –1, +1 and 0 for all the
elements. For a general element, we map
this special one to the region of interest.
A particularly convenient way to do this is to interpolate the spatial
coordinate within the element using the same scheme that we used for
displacements. Elements of this kind
are known as isoparametric elements.
Thus, in general we write
where denotes the coordinate of each node on the element.
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 L is the total number of elements, and and denote the coordinates of the ends of the lth element. We now notice an attractive feature of our interpolation scheme. The integral over the lth element depends only on the shape functions associated with the nodes on the lth element, since the displacement in this region is completely determined by its values at these nodes. We can therefore define element stiffness matrix, and element force matrix
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 . The global stiffness matrix is computed by summing all the element stiffness matrices
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 - thus
where
is the Jacobian associated with the mapping, which may be computed as
Note that the mapping also enables us to calculate the shape function derivatives in the element stiffness matrix as
Finally, recall that integrals may be approximated numerically using Gaussian quadrature as follows
where I=1…M denotes a set of integration points in the region [-1,+1], and is a set of integration weights, which are chosen so as to make the approximation as accurate as possible.
Specifically
1-D integration points and weights
M=1
M=2
M=3
Higher order integration schemes exist but are required only for higher order elements. For the linear 1-D element described earlier a single integration point is sufficient to evaluate the stiffness exactly. Similarly, for the quaratic element, two integration points will suffice.
To summarize, then, the finite element solution requires the following steps:
For each element, compute the element stiffness matrix as follows:
where
and the integration points are tabulated above, and shape functions were listed earlier.
Assemble the contribution from each element to the global stiffness
Similarly, if there is a non-zero body force then compute for each element
and assemble the global force vector
Add contributions to the force vector from prescribed traction boundary conditions at and
where the superscript denotes the node that lies at x=L.
Modify the stiffness matrix to enforce the constraints
Solve the system of linear equations
A simple example MAPLE code for this 1-D example is listed below. It is set up to solve for displacements for a bar of length 5, unit x-sect area, shear modulus 50, Poisson’s ratio 0.3, uniform body force magnitude 10, displacement u=0 at x=0, and traction t=2 at x=L.
You can download a copy of a slightly more general version of the code (which can use both linear and quadratic elements) and play games with it.
# EXAMPLE FEM CODE FOR 1D PROBLEM
> # Calculates displacement for bar constrained to displace in 1 direction only
> # subject to uniform body force, constrained so u=0 at x=0 and subjected to
> # traction at x=L
> # The example uses quadratic elements with 2 point Gauss integration
> # To run the program just keep hitting return to execute each chunk of the worksheet...
> #restart():
> with(linalg):
> #
> # Define values for parameters below
> #
> # Length and x-sect area of bar
> Length := 5.:
> A := 1.:
> # Material props
> mu := 50.:
> nu := 0.3:
> #
> const := 2*mu*A*(1-nu)/(1-2*nu):
> # Loading
> bodyforce := 10.:
> traction := 2.:
> # total no. elements and nodes; no. nodes on each element (3 for quadratic elemnts)
> L := 10:
> nnodes := 2*L+1:
> Ne := 3:
> #
> # Set up some data structures storing the mesh
> #
> coords := array(1..nnodes):
> for i from 1 to nnodes do
> coords[i] := Length*(i-1)/(nnodes-1):
> od:
> #
> # Element connectivity (speficies node numbers on each element)
> #
> connect := array(1..Ne,1..L):
> for lmn from 1 to L do
> connect[1,lmn] := 2*lmn-1:
> connect[3,lmn] := 2*lmn:
> connect[2,lmn] := 2*lmn+1:
> od:
> #
> # Integration points and weights for 2 point integration
> #
> npoints := 2:
> w := array(1..2,[1,1]):
> xi := array(1..2,[-0.5773502692,0.5773502692]):
> #
> # Assemble the global stiffness and force vector
> #
> K := array(sparse,1..nnodes,1..nnodes):
> F := array(sparse,1..nnodes):
> #
> for lmn from 1 to L do
> #
> # Extract the coords of each node on the current element
> #
> lmncoords := array(1..Ne):
> for a from 1 to Ne do
> lmncoords[a] := coords[connect[a,lmn]]:
> od:
> #
> # For the current element, loop over integration points and assemble element stiffness
> #
> kel := array(sparse,1..Ne,1..Ne):
> fel := array(sparse,1..Ne):
> #
> for II from 1 to npoints do
> #
> # Compute N and dN/dxi at the current integration point
> #
> N := array(1..Ne):
> dNdxi := array(1..Ne):
> N[1] := -0.5*xi[II]*(1.-xi[II]):
> N[2] := 0.5*xi[II]*(1.+xi[II]):
> N[3] := (1.-xi[II]^2):
> dNdxi[1] := -0.5+xi[II]:
> dNdxi[2] := 0.5+xi[II]:
> dNdxi[3] := -2.*xi[II]:
> #
> # Compute dx/dxi, J and dN/dx
> #
> dxdxi := 0.:
> for a from 1 to Ne do
> dxdxi := dxdxi + dNdxi[a]*lmncoords[a]:
> od:
> J := abs(dxdxi):
> dNdx := array(1..Ne):
> for a from 1 to Ne do
> dNdx[a] := dNdxi[a]/dxdxi
> od:
> #
> # Add contribution to element stiffness and force vector from current integration pt
> #
> for a from 1 to Ne do
> fel[a] := fel[a] + w[II]*bodyforce*J*N[a]:
> for b from 1 to Ne do
> kel[a,b] := kel[a,b] + const*w[II]*J*dNdx[a]*dNdx[b]:
> od:
> od:
> #
> od:
> #
> # Add the stiffness and residual from the current element into global matrices
> #
> for a from 1 to Ne do
> rw := connect[a,lmn]:
> F[rw] := F[rw] + fel[a]:
> for b from 1 to Ne do
> cl := connect[b,lmn]:
> K[rw,cl] := K[rw,cl] + kel[a,b]:
> od:
> od:
> od:
# Add the extra forcing term from the traction at x=L
> #
> F[nnodes] := F[nnodes] + traction:
>
> # Modify FEM equations to
enforce displacement boundary condition
> # To do this we simply replace the equation for the first node with u=0
> #
> for a from 1 to nnodes do
> K[1,a] := 0.:
> od:
> K[1,1] := 1.:
> F[1] := 0.:
>
> u := array(1..nnodes):
> u := linsolve(K,F):
> print(u):
[0, .06964285584, .1357142831, .1982142819, .2571428523,
.3124999942, .3642857073, .4124999925, .4571428490,
.4982142776, .5357142776, .5696428490, .5999999919,
.6267857067, .6499999926, .6696428501, .6857142791,
.6982142797, .7071428511, .7124999942, .7142857085]
The predicted displacement field is
plotted below
Of course in general we want to calculate more than just displacements – usually we want the stress field too.
We can calculate the stress field anywhere within an element by differentiating the displacements to calculate strains, and then substituting into the constitutive relation. This gives
This works well for a uniform body force with quadratic (3 noded elements) as the plot below 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 this
Notice that the stresses are most accurate near the center of each element (at
the integration point). This turns out
to be a general result. To get the best
accuracy for stresses, one should
integrate the stiffness matrix using the smallest number of integration
points required to evaluate the integrals exactly (2 for quadratic 1D elements,
1 for linear 1D elements), and then sample stresses at the integration
points. For this reason, FEM codes
generally output stress and strain data at integration points.
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 – most of the elements are zero. This allows the matrix to be stored in compact form – for very large matrices indexed storage (where only the nonzero elements together with their indices are stored) is the best approach; for smaller problems skyline storage or band storage (where only the central, mostly nonzero, band of the matrix is stored) may be preferable. In this case equation numbers need to be assigned to each degree of freedom so as to minimize the bandwidth of the stiffness matrix.
General 2D and 3D static linear elastic FEM code
It is straightforward to extend the 1-D case to more general problems. All the basic ideas remain the same.
In both 2D and 3D we divide up our solid of interest into a number of elements, shown schematically for a 2D region below
The displacement field within each element is interpolated between nodes that lie on the boundaries or within the element. There are many possible interpolation schemes. A few of the more common 2D and 3D schemes are summarized below. As for the 1-D problem, we develop an interpolation scheme for a simple geometry (cube or regular tetrahedron) and then map it to a more general region.
For both 2D and 3D elements, the general expression for the displacement field within an element is
where the shape functions must be chosen appropriately. A few of the more common interpolation schemes are listed below.
2D elements and interpolation functions
Linear and quadratic triangle elements (the numbers inside the elements indicate the face numbering scheme)
Linear and quadratic quadrilateral elements (numbers inside the elements indicate the face numbering scheme)
3D Element Interpolation Functions
Linear and quadratic tetrahedra
Element faces:
Face 1 has nodes 1,2,3
Face 2 has nodes 1,4,2
Face 3 has nodes 2,4,3
Face 4 has nodes 3,4,1
Linear and quadratic brick elements
Element faces:
Face 1 has nodes 1,2,3,4
Face 2 has nodes 5,8,7,6
Face 3 has nodes 1,5,6,3
Face 4 has nodes 2,6,7,3
Face 5 has nodes 3,7,8,4
Face 6 has nodes 4,8,5,1
These simple element geometries are mapped into a general shape by setting
within each element, where are the coordinates of the nodes.
We now wish to calculate the stiffness matrix and force vector for our interpolation scheme.
As for the 1-D example, we divide the volume and surface integrals up so as to integrate over each element in turn
As in our 1-D example, we notice that the volume integral over the l th element depend only on the shape functions for that element. We can again define an element stiffness matrix, and element force vector, which are functions only of the geometry of the element, properties of the material within that element, and the load applied to the element, and are independent of the rest of the solid
and form the global stiffness matrix and force vectors by adding together all the element stiffness matrices and force vectors
Finally, we need efficient techniques to evaluate the integrals for the element matrices. As before, it is convenient to map the region of integration to a convenient `unit’ area or volume – once again the space used to define the element shape functions is the best choice for the region of integration.
To map the region of integration we define the Jacobian matrix for the mapping as
where the shape function derivatives are easily computed. Note that for a 2D problem is a 2x2 matrix; for a 3D problem it is a 3x3 matrix. Then define
whereupon
We note in passing that the boundary integral in the traction 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 the preceding section (1D element) 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.
Finally, we need to devise a way to evaluate and an integration scheme. The former is straightforward – we write
where the derivatives with respect to are easy to compute (just differentiate the expressions given earlier…) and
For the latter, we once again adopt a quadrature scheme, so that
The integration points and weights depend on the element geometry, and are listed below for a few common element types
Triangular elements
1 point
3 point or
(the first scheme here is optimal, but has some disadvantages for quadratic elements because the integration points coincide with the midside nodes. The second scheme is less accurate but more robust).
4 point
Tetrahedral elements
1 point
4 point
where
Quadrilateral and hexahedral elements
For quadrilateral elements we can simply regard the integral over 2 spatial dimensions as successive 1-D integrals
which gives rise to the following 2D quadrature scheme: Let and for I=1…M denote 1-D quadrature points and weights as follows
M=1
M=2
M=3
then in 2D, an quadrature scheme can be generated as follows:
for J=1…M and K=1…M let
Similarly, in 3D, we generate an scheme as:
for J=1…M , K=1…M L=1…M let
With these definitions, then, we write the element stiffness matrix as
where
There is one final issue to discuss. How many integration points should be used to evaluate the integrals? There are two considerations. If too many integration points are used, time is wasted without gaining any accuracy (the solution is approximate anyway so there is no need to evaluate the integrals more accurately than we can interpolate the displacement field). If too few integration points are used, then one of two things will happen: the stiffness matrix may be singular, or convergence to the exact solution with mesh refinement will be reduced. The following schemes will avoid both
Linear triangle: 1 point
Quadratic triangle: 3 points
Linear quadrilateral: 4 points
Quadratic quadrilateral: 4 points
Linear brick: 8 points
Quadratic brick: 27 points
Sample linear elastostatic FEM code
Implementation of a finite element program is best understood by looking at a sample program. Here, we provide a general 2D or 3D linear elasticity code programmed in MAPLE. All the element types discussed here have been implemented.
The program reads an input file that looks like this:
No._material_props: 2
Shear_modulus: 10.
Poissons_ratio: 0.3
No._coords_per_node: 2
No._DOF_per_node: 2
No._nodes: 4
Nodal_coords:
0.0 0.0
1.0 0.0
1.0 1.0
0.0 1.0
No._elements: 2
Max_no._nodes_on_any_one_element: 3
element_identifier; no._nodes_on_element; connectivity:
1 3 1 2 4
1 3 2 3 4
No._nodes_with_prescribed_DOFs: 3
Node_#, DOF#, Value:
1 1 0.0
1 2 0.0
4 1 0.0
No._elements_with_prescribed_loads: 1
Element_#, Face_#, Traction_components
2 1 1.0 0.0
This input file sets up the boundary value problem illustrated below
The elements are linear elastic plane strain with . Here is a brief explanation of the data in the file
1. The first part of the input file specifies material properties. This should be self-explanatory!
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 – the identifier is provided to allow addition of more sophisticated element types such as reduced integration elements), then enter the nodes on each element following the convention shown earlier.
4. The fourth part of the code 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 code 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 the preceding section – note that you must be consistent in numbering nodes and faces on each element), and the components of traction acting on the element face, as a vector with 2 or 3 components.
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 – often during element stiffness assembly.
For the input file shown, the program produces an output file that looks like this
Nodal Displacements:
Node Coords u1 u2
1 0.0000 0.0000 0.0000 0.0000
2 1.0000 0.0000 .0350 .0000
3 1.0000 1.0000 .0350 -.0150
4 0.0000 1.0000 0.0000 -.0150
Strains and Stresses
Element: 1
int pt Coords e_11 e_22 e_12 s_11 s_22 s_12
1 .3333 .3333 .0350 -.0150 0.0000 1.0000 .0000 0.0000
Element: 2
int pt Coords e_11 e_22 e_12 s_11 s_22 s_12
1 .6667 .6667 .0350 -.0150 .0000 1.0000 .0000 .0000
This should be self-explanatory.
Here is a link to a MAPLE worksheet for the FEM code
Here are is a link to the sample input file
Additional demo files are available on the computing page