MECH ENG 4118 & 7059

Finite Element Analysis of

Structures
Course Notes
Semester 1, 2014

Lecturer: Dr John Codrington

These notes were prepared by Nick Kastelein, John Codrington and Carl Horward. The notes are
made available for study in MECH ENG 4118 & 7059, Finite Element Analysis of Structures at
The University of Adelaide. Any reproduction of these notes is not permitted.

i
ii
Contents

Preface iv
1 - Introduction 1
2 - Procedures for FEA 11
3 - Theory of FEA 39
4 - Meshing 63
5 - Heat transfer 75
6 - Non-linear analysis 83
7 - Dynamic Structural Analysis 103
8 - Buckling 121
9 - Explicit Dynamics 133
10 - Acoustics 153

iii
Preface
The course will equip students with the necessary knowledge to use finite
element analysis to solve problems related to solid mechanics, dynamics, heat-
transfer and acoustics. FEA is a design/research tool that is extensively used in
industry and research institutions. Students also will gain hands-on experience in
using the finite element analysis software ANSYS to solve realistic engineering
problems.

On completion of the course, students should:

1. Have a good understanding of the principles of current finite element
modelling techniques applied to solid mechanics, dynamics, heat transfer,
and acoustics;
2. Be equipped with basic understanding of the mathematical representation
of the various processes involving stresses in structures, vibration of
structures, transfer and conservation of heat;
3. Have a deep understanding of limitations and applications of current
techniques and codes to solve complex engineering problems;
4. Have developed analytical cognitive skills and improve problem solving
skills in these areas;
5. Undertake a self directed project to use FEA as a tool to solve an
engineering problem;
6. Be able to write a professional engineering report;
7. Be able to critically assess a finite element analysis for correctness;
8. Understand the need to undertake lifelong learning.

iv
Chapter 1
Introduction

1
Course Details
The Purpose of the Course
This course is about finite element analysis (FEA). FEA is used by engineers, scientists and
designers across a wide range of disciplines. It is a method of modelling and then analysing
physical systems. This particular course focuses on use of the finite element method (FEM) to
analyse physical structures for engineering purposes. The purpose of the course is:

to teach the fundamentals of the finite element method for the analysis of structures,
and
to equip students with the skills to do, or review, finite element analyses of structures.

The specific aims of the course are:

to provide an introduction to the theoretical basis of FEA

to provide the capabilities and understanding of all stages of the FEA process
to introduce students to a range of analysis types (e.g. structural, vibration, heat
transfer, etc)
to make students aware of the limitations of FEA
to provide a strong foundation for future learning in FEA.

Course Content
Part 1 Fundamentals

What and why do we use FEA?

Procedures and workflow of FEA
Methods of modelling real structures
Theory of the final element method

Part 2 Applications of FEA

Heat transfer
Vibration
Non-linear analyses & contact
Buckling
Acoustics
Dynamics

The course is designed to cover core analysis types, that students are likely to come across as a
professional engineer, plus a couple of advanced topics.

Assumed Knowledge
It is assumed that all the students already have a good understanding of the following topics:

Structural Statics, Dynamics

Stress Analysis
Heat Transfer and Thermodynamics
Basic Vibrations
Engineering Design

2
Recommended Reading
The following books are recommended reading:

Cook R, Finite Element Modeling for Stress Analysis, John Wiley & Sons, 1995.
Fagan M, Finite Element Analysis, theory and practice, Pearson Education Limited,
1992.
Moaveni S, Finite Element Analysis Theory and Application with ANSYS (3rd Edition),
Prentice Hall, 2007.
Adams V, Askenazi A, Building better products with finite element analysis, OnWord
Press, 1999.
Young W, Budynas R, Roark's formulas for stress and strain, McGraw-Hill, 2001.

It is not required to purchase any textbooks for this course. However, if you would like further
reading or a good reference book for your collection, then these are the recommended books.
Note that the books are listed in order of recommendation based on content only, without
taking into consideration the cost.

The following sources are also useful:

http://www.xansys.org/ An ANSYS user community

http://ansys.net/ Another ANSYS user community
http://www.eng-tips.com/ General engineering forum that discusses FEA
http://www.efunda.com/ General engineering reference website

Assessment
Assignments 20%

12 Weekly online quizzes. Due and closed 24 hours after distribution. 3 attempts
allowed.
Written Assignments. Due on specified due date (after approximately 2 weeks) at 5pm
in submission boxes (unless specified otherwise).

Project 20%

This project is used to reinforce the course material, and is a chance to demonstrate the
detailed application of FEA, ability to solve a complete engineering problem and ability
to present results professionally
FEA is used to answer any engineering question and it is strongly suggested for students
to relate their FEA project to their Final Year project
No group sizes larger than 3 people (ideally 2; no exceptions)
The software to be used for the project is preferably ANSYS Classic or ANSYS Workbench
The project must be challenging and solve an interesting real-world problem

Exam 60%

3 hours
Open Book

3
What and Why FEA?
Engineering Problems
Most engineering problems reduce to simple questions. For example: How did this component
fail? At what rate will thermal energy dissipate away from this component? or How much
will this component vibrate? In industry these questions are often driven by the broader
commercial implications, will this product/idea work? how long will this component last?
and how can we improve this product?

An engineer can use several approaches to find the answer(s) to such questions. These
approaches fall into one of two categories: experimental methods or theoretical methods.
Experimental methods ideally provide the most real-world results. However, they are often
impractical or too costly. Experimental methods may also be very difficult for isolating the
behaviour of individual parameters in the system. The alternative (and favourite of Engineers!)
is to use theoretical methods.

In most disciplines, governing mathematical descriptions of physical systems have been derived
in the form of Differential Equations (DEs). For example, the equations of motion or the heat
equation. To obtain the specific solution to a given problem the governing differential equations
need to be solved with boundary and initial conditions that represent the original system.

Theoretical methods can be further divided into analytical and numerical methods.

Analytical Methods: Provide a solution for the exact behaviour of a system at any point.
For some physical scenarios, the DEs reduce to mathematical problems that can be
solved directly using analytical methods.
Numerical Methods: Approximate exact solutions at a discrete number of points.
Solutions to DEs can be obtained by discretisation, which is dividing a mathematical
domain into a finite number of pieces in any dimension (time and/or space). The
solution is then approximated at each discretisation point or over each segment by
curve fitting (usually step or linear, but occasionally higher-order polynomials).

The aim of these approaches is to create a representation (model) of a physical system, which
can be used to predict the systems behaviour.

The engineering problem/system consists of the following:

The governing differential equations.

...applied over a control volume/geometry.
...with Boundary/initial conditions.

All the governing differential equations used in this course involve the continuity approximation
(i.e. Continuum Mechanics). This means treating the physical world like it is not made up of
atoms, but is instead described by continuous scalar and vector fields, such as Temperature,
Stress and Displacement. This approximation is generally valid above the 10-9 m length scale.

Finite Element Analysis

Of the three approaches mentioned above, FEAs rely specifically on numerical methods, though
in some instances incorporate analytical solutions (e.g. beam theory, vibration theory) and
experimentally obtained information (e.g. material properties). Complicated finite element
analyses also need to be verified and validated using analytical and experimental techniques, so
a finite Element Analysis is never stand-alone.

4
FEA is a widely used analysis tool because it is often the most practical way to solve an
engineering problem. FEA is

generally cheaper than creating a physical prototype,

is flexible enough to handle problems for which no analytical solutions are known, and
takes advantage of the opportunities created by modern computers, which have seen
increasing computational abilities over the last century.

FEA has been used in the analysis of Micro Electromechanical Machines, off-shore oil drilling
platforms, bones and joints, Cardiac heart valves, Welded joints, Satellite components, Aircraft
fuselages, Under-sea pipelines, Loudspeakers and Bridges, and the list goes on...

The Basic Steps of FEA

The finite element method uses the following procedure for creating a model:

The continuous geometry of the control volume is replaced with a finite number of
simpler units called elements, which are connected at points called nodes. This step is a
type of discretisation, called spatial discretisation (distinct from temporal discretisation,
which is dividing a period of time into discrete time-steps), or more commonly,
meshing. The array of elements and nodes are called a mesh.
A simplified continuous function is used for the solution within an individual element,
which gives a set of relationships between all the variables (scalar and/or vector field
values) at the nodes. This is called a shape function.
The individual element solutions are assembled to form the complete problem and
continuity is maintained at inter-element boundaries.
Boundary and initial conditions are applied so that the problem must have a unique
solution. The result is a set of simultaneous equations that are best solved using
computers. Sometimes they must be solved iteratively.
The unique solution is obtained, and the results are interpreted in relation to the
purpose of the analysis.

5
You can use the following steps to perform a thorough finite element analysis. Most of these
will be covered in much more detail during this course. You might need to redo any of these
steps several times and you should always have a healthy skepticism of your own work.

1. Start at the end. What results do you want from the analysis? This will determine the
analysis method. Do you even need to do a FEA? Know at the beginning what you want
to have at the end!
2. Choose your assumptions carefully. Determine: What boundary conditions, element
types, etc are applicable? How will you convert the real-object into a mathematical
model? The discretisation is usually not the only approximation. Physical
approximations are generally involved as well, for example assuming small
displacements, neglecting buckling, etc.
3. Do a simple hand-calculation first. Get text books for similar standard problems to
get an approximate answer.
4. Plan your Finite Element Model on a piece of paper. How many elements per
segment will you need? How will the model be divided to enable a simple mesh? How
will the loads/constraints be applied?
5. Create a solid-model.
6. Mesh the solid-model.
7. Apply the BCs.
8. Get someone to independently check the entire model (Quality Assurance).
9. Solve the model.
10. Review and post-process the results. The FEM section is over, you can now
extract the numerical answers.
11. Do a sanity check of the numerical answers.
12. Interpret the numerical answers.
13. Write your answers clearly to explain the results. Any analysis in engineering is
useless unless the results can be explained and interpreted.

Advantages and Disadvantages

Some advantages are that FEA:
Can be used for a variety of engineering problems
o Solid mechanics
o Heat transfer
o Fluid dynamics
Can readily handle very complex geometry
Can manage complex constraints
o E.g. indeterminate structures can be solved
Can allow for complex loading
o Nodal and element loads
o Time varying loads

6
Some disadvantages are that FEA:
Doesnt provide a closed-form solution
Only provides results for the given inputs
Only obtains approximate solutions
Has inherent errors
Mistakes by the user can be fatal
GIGO garbage in, garbage out!

As engineers, we are primarily interested in the A of FEA. That is, the analysis. Without it, you
merely have a mathematical equation solving method; the finite element method. Finite
element analysis can be powerful, but it cant be separated from an understanding of the
physics involved. To be a master at FEA, you need to have a solid understanding of the physics
and real-world behaviour of the problems you are trying to solve.

It is very easy to obtain a meaningless solution using finite element analysis. FEA analyses are
easier to use to get a qualitative understanding of a system; creating a model that is
quantitatively accurate is a difficult task (and usually not of interest to engineers!).

Difficulties in FEA
FEA is difficult because a physical system needs to be converted into a mathematical one, in the
form of a finite element (FE) model. This is where your engineering knowledge and experience
are crucial. There is a big difference between a CAD model and a good FE model. Modern FEA
packages often advertise ready conversion from one to the other and make it seem too easy.
The following advice is useful to remember:
Beware of the Wow factor. Pretty pictures, animations, etc are a distraction.
Check the underlying assumptions. Remember GIGO. Check Material properties,
Element types, BCs, Method of modelling, Density and distortion of the mesh
Are the results reasonable? Are there inconsistencies, are the numerical results near
expected, can you do physical testing to compare with FEA predictions?
Could the FEA have been done better and does it matter for your question? If you only
needed an answer within 30% of a particular value, there is no point in fine-tuning.

Historical Overview

1940s: Triangular elements for torsion problems (Courant)

1950s: Introduction of the displacement and stiffness methods for complex aerospace
structures
1960: Coining of the term Finite Elements by R Clough
1970s: FEM applications began appearing in a wide variety of engineering applications
1980s: Introduction of powerful computer graphics and better elements
1990s: FEM is the standard tool for structural analysis
2000s: FEM is the standard tool for engineering simulation

7
Final Words
FE computer programs have become widely available, easier to use, and can display results
with very attractive graphics. It can be hard to disbelieve FE results because of the polish of
their presentation, but remember; even an inept user can produce smooth and colourful
stress contour plots. The fancy results can hide incorrect boundary conditions, inappropriate
element types, poor mesh quality, and incorrect conceptualisation of the original problem.

A responsible user must understand the physical nature of the problem and the behaviour of
finite elements well enough to prepare a suitable model and evaluate the quality of the
results. Expertise in using FE for one physics (say, solid mechanics) does not imply
competence in using FE for another (e.g. heat transfer). Responsibility for results produced is
taken by the engineer who uses the software.

The above comments may sound pessimistic, but they do well to outline the purpose of this
course. The intention of this course is to provide a solid foundation in the key aspects of
finite element analysis for engineers. From this, and through time and experience, the user
will be able to use FEA to its full and powerful extent. Happy FEAing!

8
9
10
Chapter 2
Procedures for FEA

11
Every finite element analysis involves four main steps:

Conceptualisation
o Define the problem
o What is the question?
o Preliminary decisions
Pre-processing
o Create the model
o Define material properties
o Mesh the geometry
Solution
o Apply Boundary Conditions
o Solve
Post-processing
o Check the solution
o Interpret the results

Conceptualisation
Solving engineering problems requires converting physical systems into mathematical ones.
Through knowledge and experience of the field, the real-world problem is broken down into
simpler questions that can be answered within the limitations of modelling. A geometric domain
with boundary conditions and governing differential equations is then created (the FE model) to
match this conceptualisation.

Before attempting any finite element analysis it is imperative that you develop a sense for the
problem. You can save lots of time and money by spending a little time with a piece of paper to
develop your understanding of the problem and to plan your analysis. There are many questions
that should be asked before proceeding with the analysis. These preliminary decisions are a
compulsory step of ALL finite element analyses.

As part of the conceptualisation process, you will need to make many decisions about the:

Analysis type: what kind of simulation will be conducted e.g. static/transient,

linear/non-linear, or ?
Geometry: how much of the real geometry is included in the model, 2D or 3D? etc.
Element type: will 3D solid elements be required or can simplified versions be used?
Boundary conditions: how and where will the loads and constraints be applied?
Material properties: is the material homogenous and isotropic, or. ?

How you choose to create the model depends on what you want to get from the analysis. This
chapter details the different options available and the basic procedure that should be followed
when conducting a finite element analysis.

12
Which Analysis Type?
There are many types of analysis that can use Finite Element Modelling. FEA software packages
are commonly available in the following disciplines:

Structural: Deformation of solid bodies, motion of solid bodies, or contact of solid

bodies
Thermal: Applied heat, high/low temperatures, or changes in temperature
Electromagnetic: Devices subjected to electric currents (AC or DC), electromagnetic
waves, and voltage or charge excitation
Fluid: Motion of gases/fluids, or contained gases/fluids. Generally these analyses have
their own acronym, CFD (Computational Fluid Dynamics). Because of convection terms
and essential non-linearities which are absent from structural analyses, even basic CFD
requires an iterative approach.
Coupled-Field: Combinations of any of the above.

Within these physical disciplines, there are several types of analyses. For example:

Static
Dynamic (Transient, Periodic)
Linear
Non-linear

What to Model?
In an FEA analysis, a real and physical system is converted into a Finite Element model. This
involves changes to the geometry. The modeller has to ask:

What things are we concerned about with the model geometry?

How much detail is really necessary?
Are there any simplifications that can be made?
o Symmetry
o Plane stress/strain
o De-featuring
Other simplifications and idealisations?

The way you chose to model something depends on the question you are asking: What are the
local stresses at a geometric feature? What is the global structural displacement? or What is
the first vibration mode shape and resonance frequency? There are a variety of ways to model
a given problem.

Usually the whole geometry does not need to be modelled. The scope of the analysis is
determined, and often the model can also take advantage of symmetry.

The scope of the analysis: Depending on the desired outcome, an entire car could be modelled,
merely the chassis, or perhaps only part of the chassis. Whatever scope of analysis is chosen,
boundary conditions will be used to describe the relationship between the system and anything
outside of the system.

13
Planning the Structural Model
A typical structural model is defined by:
Volumes representing solid objects
Areas, which are faces of volumes, or represent planar or shell objects
Lines, which are edges of areas, or represent beam objects
Keypoints, which are locations in 3-D space and also represent vertices of objects

When creating an FE structural model, simplification is the key. The simpler the model, the
easier it is to create, mesh, solve and most importantly, understand. Put another way, the more
complex the model, the easier it is to make mistakes! However, correctly simplifying a real-
world structure requires experience and understanding of the underlying physical behaviour.

The main simplifications that can be made to the geometry when creating an FE model are:
Symmetry
De-featuring
The state of the stresses/strains

The model can use symmetry in different ways:

Reflective symmetry. Simplifies a scenario that is reflected about a plane.
Cyclic symmetry. Enforces the same solution on one surface to another identical
surface. This is useful when a 2D approximation is used for an infinitely deep problem,
or when a scenario includes a continuous pattern. For example, the array of tubes in a
heat exchanger.
Axisymmetry. Axisymmetry uses the same type of boundary conditions as a cyclic
symmetry scenario, it describes a situation that is revolute around a line axis.
Antisymmetry. Antisymmetry is essentially cyclic rotational symmetry.

14
Note that symmetry cannot be used when only the geometry is symmetrical. Everything in the
model, including the loads and constraints, must be symmetrical.

De-featuring of a model involves removing features of the real geometry. Features can be
removed or suppressed, when they are not going to affect (significantly) the analysis objectives.
For example, if the local stresses at a joint or connection are not of interest, then the joint can
be de-featured to a simpler geometry e.g by not modelling welds. As long as the region of
interest is sufficiently remote from the joint, then the de-featuring will have minimal influence.
Other common de-featuring includes removing fillets or radii to create straight edges.

The goal of de-featuring is to improve the mesh quality, thus increasing the accuracy of the
results, and to eliminate regions of unnecessarily high mesh density. Excessive mesh density can
significantly increase the computational time, reduce the overall mesh quality, and lead to
convergence issues. Complex features can also create difficulties when trying to generate the
mesh (especially when using auto meshers).

Another way to simplify the structural model is to take advantage of the state of the stresses
(or strains) developed within a particular geometry for a given loading configuration. You
should already be familiar with the concepts of plane stress and plane strain, as well as how
slender beams develop stresses differently compared to short stocky beams or compared to
plates. We can make use of these mechanical behaviours by geometrically modelling only the
material directions where significant stresses are developed. Further discussion on this is
found in the next section on element types.

15
It should be noted at this point, that the structural model plays no part in the final solution
stage. The model is used as a domain over which the mesh of nodes and elements is
generated. The structural model also provides an effective way to apply boundary conditions
(which are ultimately passed to the mesh).

There are two ways to create an FE mesh: 1) by meshing a structural model, and 2) by directly
creating the nodes and elements. A combination of each may also be used. Further discussion
on elements and meshing is provided in the next section and in Chapter 3.

What Element Type?

Elements are the foundation of any finite element model. They represent the physics (e.g. solid
mechanics, heat transfer) of the problem through mathematical equations. There are a large
variety of elements to choose from. Elements can vary in dimension, shape, number of nodes,
and function. The choice of element will determine how the geometry will be broken up or
discretised. It is important to select the elements that are best for the simulation being
undertaken. This is not always an easy task, and may have more than approach.

The choice of element is strongly dependent on

Geometry of the structure

Loading conditions
Type of analysis
Expected behaviour (stresses, strains, temperature)

Element selection includes

Physics: structural, thermal, electromagnetic, fluid, coupled-field

Dimension: beams (1D), plates/shells (2D), solids (3D)
2D Shape: quadrilateral, triangle
3D Shape: hexahedral, pentahedral, tetrahedral
Order: linear, quadratic, p-elements

16
The 1D and 2D element types (beams, shells, and planar solids) incorporate analytical solutions
to describe what is happening in the dimensions that are not directly modelled. This reduces the
computational requirements for solving the model.

First order elements have linear (first order polynomial) interpolation of field values, whereas
the second order elements have quadratic (second order polynomial) interpolation of field
values throughout the element. This will be further discussed in Chapter 4.

Some advice for choosing the right elements:

Use the simplest element type:

o If you can make a beam model ... do it
o if you can make a shell model ... do it
In general
o Quads are far better than triangles
o Bricks are better than wedges, which are better than tets
o Higher order elements are more accuracte than lower order elements
For structural analysis, dont use first order tri/tet elements (but they are ok for thermal
analysis)

There are also a range of special elements, which are useful for simplifying the model, such as
for describing idealised body interactions. Sometimes springs and dampeners can link the
modelled system to the surrounding environment, these are described as restraints as distinct
from constraints.

Degrees-of-Freedom (DOF) represent the allowable movement or freedom of each structural

node (or equivalent for thermal or simulations of other physics). In the case of a structural
analysis this means the orthogonal (i.e. independent) directions of translation and axes of
rotation. A 1D spring for example, will only have 1 DOF if the node can only move in the one
translation direction (see diagram below).

Each node has a certain number of degrees-of-freedom (each creating an equation in the
mathematical system for that element). The final model then has a total number of degrees-of-
freedom, which is the sum of all DOF from each node. This total DOF relates to the number of
equations in the final system and thus the computational power required to solve the model.

17
Beam Elements can be used when the generated stresses conform to beam theory assumptions.
This usually occurs for members where the cross-sectional dimensions are small relative to the
length. Beams are used to model trusses, frames and stiffeners (beams attached to wall-
structures). The nodes of beam elements can have up to 6 structural DOF, which are 2 in-plane
rotation DOF (i.e. bending), 3 translation DOF (i.e. x, y, z), and axial rotation (i.e. torsion about
the long axis). The cross-sectional geometry of the elements is not physically modelled; instead
it is included as a constant property of the elements (which is then incorporated into the
mathematical equations of the elements via beam theory).

Beam elements should be used when the generated stresses/strains in the real structure
conform to the assumptions of beam theory (that is, there are only axial normal stresses,
transverse shear and torsional shear). Typically this is for sections where the cross-section is
small compared to the length.

There are a range of beam elements based on different beam theories:

Slender beam theory (Euler-Bernoulli) Rule-of-thumb: L/h > 10
Stout beam theory (Timoshenko) Rule-of-thumb: 3 < L/h < 10
Torsion theory (St Venant & Vlasov)
Thin-walled theory

18
Shell Elements are 2D elements. They are used to model shells and thin-walls, and are
appropriate when the stresses and strains fit into shell-theory assumptions. The nodes of shell
elements can have up to 6 structural DOF (3 translation and 3 rotation), which is the same as for
beam elements. The thickness of the shell is included as a constant property of the elements
(which is then incorporated into the mathematical equations of the elements via shell theory).
The thickness can be made up of layers, such as in the case of a composite laminate.
Shell elements should be used when the generated stresses/strains in the real structure
conform to the assumptions of shell theory (that is, in-plane normal stresses, transverse shear
and torsional shear). Typically, this requires that the thickness of the shell is small compared to
the other dimensions.
Three shell theory formulations are:
Membrane theory - no bending or torsion
Thin shell theory (Kirchoff theory) Rule-of-thumb: L/t > 10
Thick shell theory (Mindlin theory) Rule-of-thumb: 3 < L/t < 10

Shell elements should also only be used when the radius of curvature of a wall is more than ten
times greater than the thickness of it. If 3D elements are to be used to model a thin wall, there
should be at least 3 to 4 elements through the thickness to achieve acceptable accuracy for the
stress distribution in bending, so a shell element representation can require far less nodes, and
take less time to compute. An axisymemtric version of shell elements is also available.

19
Planar solid Elements are 2D elements that are used to represent a 3D solid by only modelling a
single plane. These elements are used to model the situations of plane stress, plane strain and
axisymmetry. Planar solid elements only have 2 DOF, which are in-plane translation. They do not
allow rotation as it is not required for the theoretical formulation. All the elements must lie in
the same plane.

Solid Elements make no geometrical assumptions, so they look the same as what they model.
They do not have rotational degrees-of-freedom, which is appropriate for a 3D representation.
Solid elements have the most nodes of any element type. They are therefore the least
computationally efficient element and should only be used when absolutely necessary.

20
Pre-processing
Define the Material Model
Material models describe the properties of the materials from which the structure is made.
Every analysis requires some material property input, such as:

Youngs modulus for structural elements

Thermal conductivity for thermal elements
Viscosity for CFD
etc.
Many FEA packages (such as ANSYS) come with a library of material models. They will also allow
customised material input. Material models can be very varied depending on the dominant
material behaviour. Many material models are described later in these notes. Some examples of
materials and the type of mathematical model required for them are:

Metals: Steel / Aluminium / Cast Iron Linear Elastic / Isotropic

Rubber Hyper-elastic
Ceramics Brittle behaviour
Composites: Carbon fibre / Reinforced Concrete Orthotropic / Anisotropic
Even these standard material model choices require assumptions. Metals also exhibit plastic
behaviour above their yield stress, and rubbers may be linear elastic only over a very low stress
range. The modeller must ask:
Where in the stress-strain curve does the model have to operate?
Is there relevant time-dependent material behaviour (relaxation/creep/impact)?
Are there structural loads due to thermal influences?
Do the material properties change (decay / UV damage / corrosion)?
The element type must be compatible with the material model (some models assume small
displacements, and some cannot include thermal loads). Some FEA packages have special
element types for each material. Other packages have discrete element types and allow the
material model to be changed.

Create or Import the Model Geometry

The physical geometry must be transformed into a finite element model using the chosen
element types. Some FEA packages allow geometry to be imported from CAD models. At other
times, the geometry has to be created within the FEA package.

Some geometric detail is lost when a 3D object is described by 2D and 1D elements. A

mathematical description of the detail that is lost must be included as a property of these
simplified elements. These geometric properties are (as a minimum):

shell elements
o wall thickness
beam elements
o first moment of area about the two normal axis
o cross-sectional area
rods (1D elements with no bending)
o cross-sectional area

21
Mesh the geometry
Meshing is the process of filling the model with nodes and elements. The array of elements and
nodes are called the mesh. The solid model does not participate in the FE solution; it is simply a
means of effectively generating the mesh. Each node has several degrees-of-freedom, each of
which creates an unknown in the system that must be solved. The elements provide information
of how the nodes are connected together and thus are used to generate the final system of
equations. Element theory and formulation will be presented in Chapter 3.

There are many approaches to the meshing stage of the FEM. The structure can be filled
automatically, or it can be meshed more manually.

The size of the elements determines the precision of the solution. In general: the smaller the
elements, the more accurate the solution. In most analyses, the mesh is refined (made
increasingly fine) until the solution stops changing; this is called achieving mesh convergence.
Refinement is often not done uniformly throughout a domain; it can be focussed on a region of
interest to save time spent on unnecessary computation in the rest of the domain. The following
general procedures for testing mesh convergence should all be carried out in a finite element
analysis:

Iteratively double the mesh refinement and check the solution

View the error estimates
Compare to hand calculations

Note that the structure of the mesh influences the solution accuracy. Generally more structured
meshes (such as grids) are better, though they cannot easily handle complex geometries.
Patterns in meshes can also cause errors, where the material properties begin to behave
anisotropically (different in different directions).

Meshing is covered in detail in the Chapter 4.

Solution
The variables that are relevant to an analysis are divided into two categories, Degrees-of-
Freedom (DOF), such as displacement and temperature, and Loads, such as forces or heat flows.
Once a problem is discretised, a system of equations can be formulated that relate the loads to
the DOFs (using the assumed behaviour of the elements). The final step is to apply boundary
conditions (BCs) that represent the external loads and constraints on the structure.

22
Apply Boundary Conditions
There are two types of boundary conditions:

Loads are inputs to the system, such as Force, pressure, acceleration, convection,
radiation, heat flux etc.
o Concentrated Loads. Point loads, such as forces or heat flow rates.
o Surface Loads. Loads distributed over a surface, such as pressures or
convections.
o Body Loads. Volumetric or field loads, such as internal heat generation.
o Inertia Loads. Loads due to structural mass or inertia, such as gravity and
rotational velocity.
Constraints fix DOFs of the system, such as rotation, displacement, temperature etc.
o Constraints can be used to enforce symmetry or antisymmetry BCs.
o Constraints describe the relationship between features represented by the FEA
model and everything that lies outside the modeled system.

Loads and constraints are often distinguished based on what side of the governing equation
they are on.

Boundary conditions need to be applied with care to ensure that the problem has a unique
solution. The BCs will also determine how the structure responds (displaces, vibrates, heat-flow,
etc). It is VITAL that the boundary conditions accurately represent the real problem in order to
get relevant results.

Things to look out for:

Check the constraints. Is the model:

o over-constrained,
o under-constrained, or
o have redundant constraints?
Know which direction the loads are going both before and after deflection
Do the BCs represent the real-world problem?

Loads and constraints can represent a number of physical realities:

Bearing supports
Cylindrical support
Frictionless support
Thermal conditions
Pressure loading
Additional mass
Bolt pre-tension

The method by which the loads are applied must work with the finite element method. The
loads are discretised just like the geometry. The loads are applied at nodes or over elements if
the element type allows it. Generally a sufficiently fine mesh will mean that purely nodal loading
is sufficient.

23
To accurately model real boundary conditions, sometimes you need to be creative. Constraints
or loads may not exist for exactly what you want. A crude method to cope is to investigate the
limits or extreme cases (i.e. separately modeling simply supported and clamped constraints
when the reality is somewhere in between). Otherwise modeling can combine other BCs or use
special elements.

Sometimes loads can be applied in discrete load-steps. Several reasons for doing this are
outlined in Chapter 6 on non-linear problems. The main reason for deliberately using load-steps
is to account for geometric or material non-linearity, or specifically, to reduce the error
associated with small displacement assumptions.

Loads and constraints can be applied to the structural model or directly to the mesh. Structural
model loads are the easiest to apply as you deal with surfaces, edges, etc (i.e. geometric
features). They are also independent of the mesh, which means there is no need to reapply
them after re-meshing. Applying BCs directly to the mesh has advantages in certain situations
when structural features do not exist where you want to apply the loads or constraints.

Several issues to look out for when applying boundary conditions are:

Overconstraint - too many constraints are applied, which overly stiffen the model
compared to the real structure
Underconstraint - insufficient constraints are applied such that the model is not in static
equilibrium (which will cause a static analysis to fail)
Redundant constraints - multiple loads/constraints are applied that do the same thing,
thus creating additional and unnecessary calculations in the solve process

The step of applying boundary conditions is listed here as being part of the solution stage.
However, that the choice and application of boundary conditions is closely tied with the
conceptualisation stage (i.e. pre-processing). Careful selection of the BCs must be made that suit
the simplified model and requirements of the analysis of the real-world problem.

24
Solve the Model
This step is the computers job. However, it is good to have a basic understanding of what is
happening. Things can go wrong leading to the model not solving. There are also a plethora of
solver settings that can help with numerical convergence and speeding up the solution process.

The equations that the computer solves are usually expressed as a linear relationship between
the load vector, F, and DOF vector, U. The multiplier is the stiffness matrix, K. This gives:

The stiffness matrix is determined by the material properties, the shape functions, and the
element geometries.

The solution stage is more complicated when the problem is non-linear, i.e. the stiffness is not
constant, but is a function of other variables in the problem. Transient problems are also more
complicated as derivative terms are included, such as acceleration and velocity.

Computers use a wide variety of methods to solve linear matrix equations. Directly solving the
equations can be done using factorising techniques or matrix inversion. Row reduction, a high-
school maths method, is a well-known technique for solving systems of linear equations. The
direct solver in ANSYS uses LU factorisation to decompose the stiffness matrix into lower and
upper triangular matrices (the algorithm for which is similar to Gaussian elimination). The
solution is then found through forward and back substitution. If the matrix K is symmetric, then
Cholesky decomposition is used instead of LU decomposition. The direct solver is best for sparse
matrices (matrices primarily populated with zeroes) of medium sized systems. The direct solver
is also good for ill-condition matrices. Otherwise iterative techniques need to be used.

You would be familiar with iterative solution methods from previous studies, including the
Gauss-Siedel and Jacobi methods. ANSYS, and most FEA packages, offer a large range of iterative
solvers as alternatives to the direct solver. These iterative solvers in many cases can result in
less RAM usage, less total elapsed time, and better parallel processing performance. However,
in general, iterative solvers are not as robust as direct solvers. The iterative solvers can be better
suited for very large (DOFs) models. The main iterative solvers in ANSYS use conjugate gradient
methods. The default for many analysis types in ANSYS is the Preconditioned Conjugate
Gradient (PCG) solver.

Full mathematical procedures of the different solvers are very complicated and not part of this
course.

Post-Processing
Post-processing is the final step in FEA. You need to interpret your results with consideration of
the assumptions made during model creation and solution. You may be required to make design
decisions based on the results. It is therefore a good idea not only to review the results
carefully, but also to verify and validate of the model.

25
Getting accurate results depends on assumptions, analysis method, boundary conditions,
material properties, loading, and mesh quality. These all have to be an accurate representation
of the physical system to get accurate results.

Different elements can provide different outputs. For example, shell elements cannot provide
through-thickness stress distributions, so if this type of data was desired; a different element
type should have been used.

Review the Results

Firstly it is important to know what results to look at. For this it is important to know the physics
behind a problem. For example, when viewing the stress in an object, there are many
components of stress that can be considered: normal stresses, shear stresses, principal stresses,
or Von-Mises stresses, to name a few. The relevance of these depend on what is being asked of
the analysis.

The output gained from an FEA is ultimately to answer an engineering question. Usually this
answer is needed to help make an engineering decision. Good planning will put the results in a
useful form to help with this. Often graphs between two variables (between loads and DOFs or
geometry and DOFs, or geometry and loads) are useful, other times a contour plot or vector plot
of the domain provides a qualitative view of the DOFs and this is useful as well.

A single number output may be desired for answering the question, but it wont be as useful in
verifying or validating the model.

Some output from an FEA is derived. For example, a static structural analysis calculates the
forces and displacements only, so stress and strain are derived from these. Derived variables
(like stress, strain and thermal flux or temperature gradient) are not calculated at the nodes, but
at the integration points inside an element, as shown for a square shell element below:

Singularities
When reviewing the results one must be aware of singularities. A singularity is something that
can exist in a vector or scalar field where the field limits to an infinite value; for example: the
stress at a crack-tip or at a point-load. It is worth being aware that finite element modelling will
never fully capture a singularity by returning an infinite value.

The loads (internal and external) in an FE model are not fluxes, i.e. they are not per unit area.
Stresses and fluxes are interpreted outputs that take into consideration the geometry around a
node. Refining the mesh near a singularity will increase the interpreted value towards infinity.

26
Singularities do not occur in real problems as the material will yield and redistribute the stresses
(in the case of a structural singularity). In an FEA, singularities lead to really high local stresses,
which can trick the untrained user. Care must therefore be taken when interpreting the results.

The following chart indicates when a singularity occurs for structural loading:

Load or constraint Beam Shell Solid

applied to:
Point No Yes Yes
Line No No Yes
Area N/A No No
Singularities also occur at sharp re-entrant corners and at interfaces between different
materials. Analogous singularities occur for other physics such as thermal problems.

How should singularities be dealt with?

Generally, ignore them

However, if they are in a region of interest:
o Model geometry with fillets or other details which do not cause geometric
discontinuities (de-featuring)
o Apply loads and/or constraints spread over areas rather than point locations
o Use singular elements (mainly used for fracture mechanics)

Verify and Validate the Results

There are several ways to verify and validate an FEA solution.

Mesh independence is the first and essential check. The solution must be independent of the
mesh density, arrangement and element shape quality. The most common check for mesh
independence is to check the mesh convergence by increasing the mesh refinement (though not
necessarily changing the arrangement). Mesh convergence is important to confirm that the
mesh is fine enough to achieve the required accuracy. In most cases a mesh convergence test
will also test independence from the mesh arrangement. However, this is not always the case
and other checks should be made. Furthermore, the solution should always be checked at more
than one location and for the output values of interest. Degrees-of-freedom solutions (e.g.
displacement, temperature) will converge faster than derived quantities.

27
Other ways to check the mesh are through Error Estimates, which are built into most FEA
packages, or through the unaveraged plots of derived values (e.g. stress). Error estimates
analyse the deviation of the values taken by the DOFs between the nodes in adjacent elements.
A qualitative approach to check the mesh is to look at continuity in the unaveraged contour
plots for derived quantities. Error Estimates will be further detailed in Chapter 4.

The shape of the elements can also affect the accuracy of the solution. Elements are
mathematically best in certain optimal shapes. The quality of the element shape can be
measured in term of different geometric parameters e.g. the aspect ratio and skewness.
Element shape quality measures will be addressed in Chapter 4.

For transient problems the time-step must be small enough to achieve convergence, which will
be discussed in Chapter 7.

Sanity Checks are ways to confirm that the solution is sensible. These are very powerful and
easy checks to perform. Ask questions like:

Do the reaction forces balance the applied loads?

Is the deflection shape sensible?
Are the results expected?
Where is the maximum stress?

Finite element models should be tested for their sensitivity to all of the input parameters. This
includes the mesh, as noted above, as well as the material properties, the boundary conditions,
the geometric assumptions, and so forth. A sensitivity analysis involves varying the input
parameters to see how sensitive the final results are to these changes. This process can help to
confirm decisions made with the approximations and simplifications. i.e. by showing that a
particular parameter has little effect (or a huge effect) on the final result(s) of interest.
Obviously, it could be very tedious and time consuming to test every single parameter over a
large range of values. The user therefore needs to utilise their own intuition and understanding
of the problem to decide what parameters are critical and what are suitable ranges to test over.

The above checks of the FE model are what are known as verification. Verification checks that
the FE model accurately represents the conceptual description of the real-world problem, and
that the numerical implementation by the computer has been effective.

Verification asks have I done the model right?

Verification can also be achieved through hand calculations (usually of a simpler or stripped-
back FE model) or by exactly duplicating another already verified FE model (some FEA packages
come with libraries of verification models for a range of cases).

Validation is slightly different to verification. Validation is the process of confirming that the
analysis accurately represents the original real-world problem, and that the approximations and
simplifications were justified.

Validation asks have I done the right model?

28
Validation is done by comparison with experimental results, either new experiments or
previously published. Validation is a very important step (arguably the most-important!). The
entire purpose of conducting the FEA is to solve a real-world problem, so the model needs to be
representative of that problem.

Final Words
This is the procedure for FEA. Every analysis will use this same procedure regardless of the
complexity, physics or software package used.

Some common mistakes are:

Not understanding the basic physics of a problem

Wrong boundary conditions
Wrong element types (including contact)
Using a complex element without understanding it
Poor mesh quality
Forgetting to do a hand calculation (what answer do you expect?)
Using FEA when it is not needed

29
30
31
32
33
34
35
36
37
38
Chapter 3
Theory of FEA

39
Finite element analysis is all about discretising differential equations to create approximated
mathematical models of a real system. To be truly good at FEA, you need to understand how the
FE method works in order to create these models.

FEA uses very generic methods that have been developed to handle the modelling of any
physical system. These methods use elements, nodes, properties and shape functions to
describe a problem, and solvers to find the solution to it. In this chapter the theory of FEA will
be presented.

The Finite Element Method

The equations that are solved in the finite element method are usually expressed as a linear
system in the form:

where F is the load vector, U is the DOFs vector, and K is the stiffness matrix.

You should already be familiar with this form of equation, for example the 1D spring equation.
In this case, F would be the force applied to the spring (or reaction due to an applied
displacement), U the spring displacement, and K the spring stiffness.

In more general finite element models, such as 3D solid mechanics problems, the above still
applies. However, there are more degrees-of-freedom, including displacements in the other two
dimensions, and the stiffness matrix incorporates the material constitutive behaviour for all
geometric dimensions, e.g. Hookes law.

The F and U vectors consist respectively of all of the loads, and all of the DOFs for all of the
nodes. Some values of these vectors are set by the boundary conditions, the remainder form

loads and reaction loads, + , to distinguish between inputs and outputs of the model:
the output of the analysis. It is easiest to express the load vector as a combination of applied

The stiffness matrix, K, relates the loads on each node to the degrees-of-freedom of that same
node and any neighbouring nodes (nodes with a common element). Whilst the loads and DOFs
are defined by each node, the stiffness matrix is determined by looking at each element.

This Chapter is all about the development of the stiffness matrix, K. This will provide an insight
to how differential equations can be transformed into a discretised problem and a system of
simultaneous equations. There are several approaches for formulating finite element problems
and the stiffness matrix, including Direct Formulation, Minimum Total Potential Energy
Formulation, and Weighted Residual Formulations. Each will be discussed in this chapter.

The main focus in this chapter, for the sake of simplicity, is on static, linear-elastic structural
analyses. The terminology and basic form of the above equations still applies for other physics,
such as thermal or electromagnetic.

A Background on numerical methods section is included at the end of this chapter to assist in
understanding, but it is not for assessment.

40
Finite Element Formulations
Finite Element Formulations are the mathematical methods that turn a continuous problem into
a discretised one. They transform continuous differential equations into a system of
simultaneous equations. The best way to begin the theory of FEA is with some simple worked
examples.

Direct Formulation
Consider the problem shown below consisting of two different width members joined together
and loaded in tension by a force P. Lets assume that for this configuration only the uniform axial
stresses and strains are significant. The geometry can then be discretised into two elements
based on the 1D spring equation; with one element for each member.

Each element has its own nodes with corresponding displacement at each node. These are the
degrees-of=freedom of the problem. Each element also has a spring stiffness, k.

The equations for a single element in tension are dealt with first. The equations will be written
in generalised notation for the j-th element, which links the i-th and i+1-th nodes. The stiffness
matrix of the element uses the spring stiffness, , which represents the tensile force per unit
extension (in N/m). The stiffness is given by:

41
The spring equation can be written for the force at each node (i and i+1) using the spring
stiffness and the nodal displacements and . Note that we have defined the forces and
displacements as being positive downwards (corresponding to tensile stress as being positive).

A linear system incorporating the stiffness matrix, , is then formed for this element using the
spring equations for each node. This system relates the two loads and the two displacements:

= =

Now each element has its own stiffness matrix, . Each individual element stiffness matrix can
be re-expressed as a larger global stiffness matrix due to the j-th element, , seen in the
equation below. The final global stiffness matrix , is created by superimposing these
matrices. For this formulation this yields a symmetric matrix, and so a well-conditioned
problem:

0 0 0 0 0
= = 0" + 0 "= + "
0 0 0 0 0

For problems with more elements and more degrees of freedom per node, larger stiffness
matrices are developed, like:

0 0 0 0 0
% 0 0 .
+ ' ' 0 0
$ 0 0 -
$ 0 ' '+ ( ( 0 -
$ 0 0 ( ( + ) ) 0 0 -
$ 0 0 0 ) )+ * * 0 -
$ 0 0 0 0 * *+ + +-
# 0 0 0 0 0 + + ,

The direct formulation worked well for this simple example, but it relied on the known solution
for a 1D spring. Assumptions about the stress distribution were also implicit in the use of the
spring stiffness. This is the disadvantage of direct formulation; it appeals to known solutions,

law, and with consideration to Newtons Third Law ( =

and not first principles. For example, the spring stiffness used above was derived from Hookes
).

/= 0 1D Hookes law

= 1 2,

= 1 2
56 76
86
Newtons Third Law

1 1
=
1 1

42
For more complex geometry and boundary conditions we need more versatile methods to
formulate the stiffness matrix. The formulation should allow the generation of a stiffness matrix
from any fundamental differential equations. These methods are the various weighted residual
formulations and the minimum total potential energy formulation. The latter is of most interest
as it is used for all solid mechanics problems. However, we will first briefly look at weighted
residual formulations as they are commonly used for heat transfer problems.

Weighted residual Formulations

Weighted residual methods are based on assuming an approximate solution for the governing
differential equation. Because the solution is not exact, substitution of the solution back into the
differential equation will lead to some residuals or errors. A variety of weighted residual
formulations exist, but they all use an approach to make the residual or error vanish over some
selected range or at certain points.

Here is an example of this method. Consider the differential equation, wchi describes the
problem of a tensile loaded tapered member:

1:2 > =0
;<
;=

Lets assume the solution:

1:2 = ? : + ?' : ' + ?( : (

Substituting the assumed solution into the differential equation provides the following, where @
is the error function:

1:2 1? + 2?' : + 3?( : ' 2 > = @

If the error is minimised, then the solution will be more accurate. In order to minimise the error
function for a discretised problem, several methods exist, including:

Collocation method
Subdomain method
Least-squares method
Galerkin method

The Galerkin method is most common, which uses weight functions, , and enforces:

D E @ F: = 0

43
Minimum Total Potential Energy Formulation
The minimum total potential energy (MTPE) formulation is used for almost all solid mechanics
problems. It is based on the principle that a structure or body will deform or displace to a
position that minimises the total potential energy. This position is the stable configuration for
equilibrium of the system. The MTPE formulation stems from the second law of
thermodynamics.

When a load is applied to a body, it deforms.

In a static structural analysis, all the work that is done on the body is by the external loads (W).

potential energy, , is then defined as:

The internal energy of deformation (strain) is the internal strain energy created (). The

= W

where:

= D JD /F0K FL

M = D NF

For a static structural analysis with no energy losses (e.g. heat, sound, hysteresis, fracture) the
potential energy will be zero (conservation of energy). If approximating functions are used for
the internal stresses & strains then a solution that minimises this expression will be the most
accurate. In other words, by minimising the total potential energy equation the approximate
solution will best match the exact solution. The potential energy is minimised when:

O
=0
O

This principle can be used to determine the stiffness matrix, K and the load vector, F.

44
We return to the example of the two different width members connected together and loaded
in tension (used for the direct formulation). For an individual linear spring the strain energy is:

1
P = 1 2'
2

and the work done by the nodal force is:

WP = 1 2

The potential for the spring is therefore:

1
P = P WP = 1 2' 1 2
2

The potential energy is minimised when:

OP
=0= 1 2
O1 2

and therefore:

= 1 2

which is simply the spring equation as expected. To consider more advanced problems, the
concept of shape functions needs to be introduced.

Shape Functions
In an analytical problem, degrees-of-freedom correspond to continuous scalar and vector fields,
such as temperature or displacement. In a discretised problem, these values are defined only at
each node. However, to determine the stiffness matrix, K, a method must be used that can
interpolate the values taken by these fields between nodes, throughout an element.

The FE method uses Shape Functions, or interpolation functions, to do this. The shape
functions are not arbitrarily selected. A shape function represents assumed behavior within a
given element. How well each assumed element shape function matches the true behavior
directly affects the accuracy of the solution. In an FEA, shape functions are used to not only
approximate the solution field, but also to describe the element (thus component) geometry in
different coordinate systems.

The simplest form of shape function is a linear one. Consider a single 1D link element of length L
(as shown over the page). The element has two nodes i and j (one at each end), where the
nodes are located at Xi and Xj. The element may be stretched or squashed with one DOF at each
node, ui and uj, respectively. The displacement field, u, within the element may be
approximated using a linear function (which in this simple case is actually the exact solution).

45
The linear displacement field within the element is given by:
Q
1R2 = ? + ?' R

where c1 and c2 are constants. By using the elements end conditions, which are:
Q
1R 2 =
Q
1R 2 =

we can arrive at:

Q 1R2
R R
= + R
R R R R

or written another way as:

Q 1R2
R R RR
=J K +J K

For the FE method, we re-write this in matrix form:

Q 1R2
= ST 1R2 T 1R2U

where S are called shape functions. In this case, the shape functions are:

R R RR
T = , T =
R R R R

This corresponds to simple linear interpolation and as such these shape functions are referred
to as being linear, or first order. Additional nodes can be added to increase the order and
improve the approximation of the curve to the displacement field (i.e quadratic and cubic
elements). Shape functions can also be derived for 2D and 3D elements as well as for beam and
shell elements.

The above derivation uses shape functions to describe the solutions field (in this case
displacement field) in terms of the nodal displacement values. Shape functions are also used to
transform between local element and global coordinate systems. First it is necessary to
introduce the concept of global, local and natural coordinates.

46
Finite element modelling involves the use of several different coordinate systems. Global
coordinates are used to represent the location of each node, orientation of each element, and
for the application of boundary elements. Local and natural coordinates systems are also used,
because they offer mathematical advantages when formulating the element geometry and
stiffness matrix. For the 1D element considered previously, the local and natural coordinates are
given as:

Local coordinates: V = R R

Natural coordinates: W = Y 1
'X
6 ZY[

range from W = 1 at node i to W = 1 at node j. Natural coordinates are very useful for
Local coordinates start at node i while natural coordinates start in the centre of the element and

numerical integration, which is needed to generate the stiffness matrix. The shape functions for
a 1D linear element in natural coordinates are:
1 1
T = 11 W2, T = 11 + W2
2 2

Using these shape functions yields the same results in the displacement field equation as for the
global shape functions given above.

Shape functions can be used to represent the solution field (such as displacement or
temperature) as well as for coordinate transformation between natural and global systems.
Elements that have the same shape functions for the solution field and geometric coordinate
transformation are called isoparametric. Solid 3D and planar solid elements are isoparametric.
Beam and shell elements are not isoparametric. This is because beam and shell elements use
higher order shape functions to describe the solution field (well come back to this later in the
chapter). Further discussion on the reason for using shape functions and their implementation is
also provided later in this chapter.

functions. Generally \ order elements will have \ 1 mid-side nodes. Higher order elements
The order of an element describes the order of the polynomial of the (geometry) shape

better provide a better approximation of the geometric domain as mid-side nodes allow the
element sides to be curved (e.g. a quadratic function).

47
Most FE software uses linear or quadratic elements (though in some cases, P-elements can be
employed). Higher order elements provide better accuracy for fitting the geometry as well as
the solution field. However, they significantly increase the number of degrees-of-freedom and
thus the solve time. This may particularly be an issue for non-linear (iterative) solutions.

There are two classes of elements that are distinguished by their use of shape functions:
H-elements have their shape functions fixed during analysis, and solution convergence
occurs due to mesh refinement. These
o Require less calculations to ensure convergence
o Have decreased computational time
P-elements can increase in order. These automatically iterate until convergence or until
reaching the maximum order of polynomial, which for ANSYS is 8. These elements are
useful for high stress-gradients. These elements
o Have more degrees-of-freedom
o Are more tolerant of element shape distortion (skewing)
o Require lower mesh density
Using P-elements will significantly increase the solution time as it is required to keep iterating
and increasing the order of the elements until a convergence criterion is reached. This can
create difficulties for models with singularities. The best approach is to use quadratic elements
with sufficient mesh density when high stress gradients in the solution field are an issue.

48
Shape Functions and the Minimum Total Potential Energy Formulation
Shape functions will now be used to solve the two member tension problem from before. The
displacement within a given element, assuming linear shape functions, is:

= ]T T' ^

:' : ::
T = , T' =

This provides the following longitudinal strain equation:

F FT FT' 1 1
0= = = =
F: F: F:

Hence, the strain energy of an element _ is:

8
0'
P = D JD /F0K FL = D JD` 0aF0K FL = D . F:
c 2

'
= 0 =
'
d e
2 2

1 2'
=
2

The work done at a node f is:

M =

The Total Potential Energy for the whole domain with N elements and n nodes is given by:
g h

=M = P M

For the current problem this gives:

1 2' 1 '2
'
= +
' (
2 2 ' ' ( (

Minimising the Total Potential Energy requires equating the partial derivatives (with respect to
the DOFs) to zero. This provides the following system of algebraic equations:

O
= 1 2 =0
O '

O
= 1 2 1 '2 ' =0
O ' ' (

O
= 1 '2 ( =0
O ( (

49
this is vectorised to form;

% 0 .
$ - 0
$ -
$ + - ' " i ' j = 0"
$ - ( ( 0
$ 0 -
# ,

Because these shape functions represent the exact distribution (which is linear), this
representation has yielded the exact results and is the same formulation as the direct
formulation. For problems where the shape functions represent an approximation of the actual

the best discretisation of the problem. In other words, even when does not turn out to be the
solution field (such as for a tapered bar in tension) this method will still be capable of providing

true zero value, this method will at least locate the smallest . This is why MTPE is more
versatile than a direct formulation.

Integration Points
Derived variables like stress and strain are determined at integration points in an element. The
integration points are also used to evaluate the strain energy integral equation when
formulating the stiffness matrix. In the previous example, the integration was a simple analytical
one. For more complex and arbitrary element shapes a numerical approach must be used.

The integral of a function, f(x), can be determined by a weighted sum of the values that the
function takes at specific points in the element:

D 1V2FV = M 1Vk 2

where W is the weighting factor at each summation point. You should have encountered this
concept in previous years when looking at numerical integration methods, such as the midpoint
or trapezoidal rules.

There are several ways to determine the points, Vk , that will work in this equation. These points
are called the integration points. For computational efficiency, it is aimed to minimise the
number of integration points required while maximising the accuracy. The Gauss-Legendre
scheme provides such a compromise.

50
The integration points are also good locations for approximating derived variables like strain and
stress, even though these variables actually involve a derivative not an integral. If these values
were derived at the nodes, then different values would be derived depending on which element
was used because derived variables have a discontinuity at the nodes. The integration points are
mathematically good points to use. They lie inside the elements, away from the nodes.

Derived values that are calculated at the integration points are usually extrapolated to the nodal
points (or in some cases element centroid) using basic shape functions. The shape functions
used for this extrapolation will not necessarily match (generally are lower) the order of the
shape functions used for the element formulations.

General Shape Functions and the Stiffness Matrix

Shape functions are typically derived in dimensionless form (i.e. natural coordinates), for ease of
implementation with elements of any shape, orientation, and location. Lists and tables of shape
functions for elements of different dimension and shape can be found in many FEA textbooks as
well as in theory reference guides for the FEA packages.

The general form of the equation for the stiffness matrix is:

= D lm n lFL

constitutive law (o = n p). The superscript t refers to the matrix transpose, and the integral is
where B is a matrix of shape function derivative terms and D is a matrix for the material

evaluated over the volume of the element V.

The matrix B of shape function derivatives defines the strain-displacement relationship:

p=l

51
Recall from Solid Mechanics that the differential form of these equations is:

F
0XX =
FV
Fq
0== =
F:

Fs
0rr =
Ft

for the normal strains and is:

F Fq
uX= = +
F: FV

Fq Fs
u=r = +
Ft F:

F Fs
uXr = +
Ft FV

for the shear strains.

The tables below provide just a few examples of 2D planar shape functions. (Source: Moaveni S,
Finite Element Analysis Theory and Application with ANSYS (3rd Edition), Prentice Hall, 2007).

52
Beam and Shell Elements
Beam and shell elements are slightly different in formulation compared to solid 3D and planar
elements. We saw in Chapter 2 that beam and shell elements have rotational DOF in addition to
the translational DOF. These extra DOF are required because the solution field shape functions
for beam and shell elements are third order (i.e. cubic). (Recall: beam and shell elements are not
isoparametric). The higher order shape functions are required to correctly model the bending
behaviour. A derivation of beam elements for pure bending is provided below. Shell elements
are formulated in a similar manner. Beam and shell elements have less nodes than solid
elements. The idea behind these elements is to reduce the node count (which reduces
computational effort) and to increase the accuracy of the solution. Accuracy is increased by
using specific theory for the type of mechanical behaviour that beams and shells undergo (i.e.
in-plane normal stress, transvers shear and torsional shear). Unnecessary mechanical behaviour
is removed, which simplifies the governing differential equations. The trade-off is that beam and
shell elements are limited to particular ranges of geometry and loading (as discussed in Chapter
2), i.e. to modelling beam and shell type structures.

Beam Elements with Pure Bending

The simplest form of beam element is one that has pure bending. Consider a 2D beam element
as shown below. We will assume that the beam is slender and thus will Euler-Bernoulli beam
theory (which you should all be familiar with from the Stress Analysis and Design course). Each
node has two degrees-of-freedom: transverse translation (U1) and in-plane rotation (U2).

We assume a displacement solution field in the form of a cubic equation:

q = ? + ?' V + ?( V ' + ?) V (

53
Next we use the known conditions at each node for the displacement (v) and slope (dv/dx) to
solve for the four unknown constant terms (c1to c4). This leads to the following matrix equation
for the displacement field in terms of the shape functions:

q = T v + T 'v ' + T v + T 'v '

where:

3V ' 2V (
T =1 '
+ (

2V ' V(
T' = V + '

3V ' 2V (
T = '
(

V' V(
T' = + '

and:

v q
%v . x
$ ' - = wq y
$v -
#v ' , x

The strain energy from bending is given by:

/0 0'
Q = D FL = D FL
2 2
z z

Using the curvature-strain relationship from beam theory leads to:

F' q 8 '
= D d ' e FV D : ' F
Q
2 c FV
7

The last integral term is simply the moment of inertia of the beam cross-section. The strain
energy equation then becomes:

{ 8 F' q
'
Q = D d e FV
2 c FV '

The strain energy equation can be put into matrix form using the shape functions. First, the
differential term becomes:

v
F q %v .
F T F T' F T F T'
' ' ' ' '
= $ '- = n
FV '
FV ' FV ' FV ' FV ' $v -
#v ' ,

54
Note that here D is a derivate vector not the a matrix for the material constitutive law (as
described in a previous section). Therefore the strain energy for the element in matrix form is:

{ 8
| = D
Q }
n} n FV
2 c

where use has been made of the relationship:

F' q
'
d ' e = 1n 21n 2 = }
n} 1n 2
FV

The stiffness matrix is then obtained by differentiating the strain energy equation above with
respect to each DOF.

F|Q 8
= { D n} nFV =
Fv~ c

12 6 12 6
{
6 4 ' 6 2 '

= (w y
12 6 12 6
6 2 ' 6 4 '

The load matrix is given by:

T
% .

= $T -
$ -
# ,

where here S are the shear transverse forces applied at each node and M are the moments
applied at each node.

Beam elements can also have axial deformation and torsion. These are normally decoupled from
each other and the bending, and derived in a similar manner using the appropriate theory. The
axial deformation is just the solution weve looked at previously in this chapter for a 1D spring or
link element.

In a real model, beam elements will be orientated in all directions and not just horizontal as was
used in this derivation. In 2D and 3D, the local coordinate stiffness matrix for each element is
transformed to the global system (to account for the elements orientation) using transformation
matrices that pre- and post-multiple with the element stiffness matrix (Note: this will be
required for the theory assignment you will need to do some extra reading to learn how to do
this).

55
NOTES

56
Background on Numerical Methods
Numerical methods describes a broad range of problem solutions and methods, that all have
one approach in common: Discretisation. To provide a background on numerical methods,
examples are included here of the use of numerical methods to solve three types of problems:

Integration
Time-integration
Solving an equation

Integration
Consider the following simple integration problem:
'
D 12V V ' 2FV
c

The analytical solution is simple:

V( 8 4
' '
D 12V V ' 2FV = V ' =4 =
c 3 c
3 3

A Trapezoidal approximation and the Simpson rule both use numerical methods. The
problem is discretised into six segments, the curve is approximated as linear over each section,
and the function value is evaluated at the boundaries between each segment.

The trapezoidal rule uses a linear interpolation for each discrete interval, summarised by the
following formulae:
h
' 1V V Z 2 1 + 2 +2 +2 + + 2 +
D 12V V ' 2FV
= V
Z c ' hZ h

c 2 2

= 1 2 Weight, 12
0 0 0 1 0
1 0.33 0.555 2 1.111
2 0.66 0.889 2 1.778
3 1 1 2 2
4 1.33 0.889 2 1.778
5 1.66 0.555 2 1.111
6 2 0 1 0
TOTAL 7.778

'
V
D 12V V ' 2FV 7.778 = 1.2963
c 2

57
The Simpson rule uses a quadratic approximation across each pair of intervals, according to:
'
+4 +2 +4 +2 + + 4 +
D 12V V ' 2FV V
c ' ( ) hZ h

c 3

= 1 2 Weight, 12
0 0 0 1 0
1 0.33 0.555 4 2.222
2 0.66 0.889 2 1.778
3 1 1 4 4
4 1.33 0.889 2 1.778
5 1.66 0.555 4 2.222
6 2 0 1 0
TOTAL 12

'
V
D 12V V ' 2FV 7.778 = 1.333
c 3

Because the Simpson rule uses quadratic interpolation, and this function is a quadratic, it
achieves the exact solution.

These methods differ in how a function is interpolated from data known at discrete points.
These types of interpolation are both used in FEA; first order H elements use linear
interpolation, whereas higher order H-elements, and some P-elements use higher-order
polynomials and include mid-side nodes.

The complication above what is shown here is that FEA requires linear interpolation for multi-
dimensional functions. The interpolation methods are all generalised by the shape function

between points V and V' is:

method discussed in detail previously. A simple linear interpolation for a 1-dimensional element

V' V VV
1V2 = 1V 2 + 1V' 2
V' V V' V

The shape functions used here were:

V' V VV
T = , T' =
V' V V' V

The integral of this function between V and V' is the same expression seen in the summation
for the trapezoidal integration above. Working this out is left as an exercise.

X 1V' V 2` 1V 2 + 1V' 2a
D 1V2FV =
X 2

58
Differential Equations
Consider the following differential equation:

d:
= :
d

The analytical solution is simple:

1= m
D d: = D d
= : c

\1:2 = + \1:c 2

: = S:c U ~m

The solution to this problem requires an initial condition, the value of :c .

The euler scheme involves extrapolating a function from one time-step to another, assuming a
constant slope over each interval. This method is first-order accurate, which means halving the
time-step will double the accuracy. The iteration formula is:

F:
: =: + 1 2
F

The results of doing this for several time-steps are shown below (for k = 1, f0 = 1).

150
t = 1
t = 0.5
t = 0.25
exact

100
y

50

0
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
t

The Runge-Kutta scheme is much more complicated (there are several intermediately
complicated schemes in between Euler and Runge-Kutta). This method is fourth-order accurate,
meaning halving the time-step will make it 16 (2) ) times more accurate. This method is more

derivative is represented as ;m = 1, :2 is:

computationally expensive. The iteration method for a general first-order DE where the
;=

1
:h = :h + 1 +2 +2 + ) 2, = h h
6 ' (

59
where:

= 1h , :h 2

= Jh + , :h + K
'
2 2

= Jh + , :h + K
'
(
2 2

) = 1h + , :h + (2

For the problem being considered, this yields:

150

100
t = 1
t = 0.5
t = 0.25
y

exact

50

0
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
t

It is difficult to see the differing accuracy from this graph, because the method is so good that all
of the answers appear essentially exact. The values of the function at t = 5 are presented below,
showing that accuracy still increases as delta t decreases.

Interval, 12
1 145.7171
0.5 148.1579
0.25 148.3935
Exact 148.4132

This type of analysis is very relevant to transient FEA problems. A form of Euler iteration, despite
being much less accurate than higher-order techniques, is the most common method for FEA
time integration. When it is implemented, it is used in conjunction with error-controlling
methods, such as limiting courant-numbers and tracking conservation of energy.

The implementation of the Euler method for a multi-dimensional problem can be described as
implicit or explicit. These two broad approaches to dynamic problems will be discussed and
compared in the later Chapter on explicit dynamics.

60
Solving Equations
Consider the function,

1V2 = V (

We intend to find the value of x that solves the equation,

1V2 = 8

V( = 8

V( 8 = 0

Analytical analysis tells us that an equation of this kind has three solutions, though two of these
solutions may be imaginary (fundamental theorem of algebra). These solutions are found by
factorising the equation to yield:

1V 221V ' + 2V + 42 = 0

V = 2, 1 3f

Newton Raphson Iteration is a simple iterative technique that will converge gradually towards

using the x-intercept of this line as the new estimate. For solving 1V2 = , this method uses:
the result. The method requires creating a tangent to the curve at the current estimate and

1 2
= +

1 2

This iterative technique converges very rapidly for the function being considered. For a starting
estimate of x = 5:

0 5
1 3.44
2 2.5187
3 2.0995
4 2.0046
5 2.000011

This technique and variants of it are used in FEA, especially for non-linear analysis. For example,
when the analysis must locate a point on the stress-strain curve for a material, outside of the
linear elastic region, this technique may be employed. Again, there are multi-dimensional forms
of this algorithm.

The Newton-Raphson iteration technique for FEM will be covered in greater detail in Chapter 6,
which looks at non-linear FE problems.

61
62
Chapter 4
Meshing

63
A critical step in the FEA procedure that was outlined in Chapter 2, is meshing. Meshing is worth
separate consideration because:

Meshing takes time and thought to implement well.

Determining mesh convergence takes patience.
Poor meshes can cause undetected errors.
There are many useful tricks and tips that can assist in meshing an object well,
specifically for connections, and complex geometries such as in large systems with
multiple element types.

The finite element mesh mathematically represents the geometry of the problem, and the
solution field. The accuracy of the nodal DOF solution is directly dependent on the mesh quality.

Mesh Attributes
Mesh Density and Refinement
Generally, the smaller the elements are, the greater the accuracy of the solution. Mesh
convergence is when the solution reaches the desired precision. Traditional definition of mesh
convergence says the ideal model has just enough nodes to capture the required level of detail
without over-meshing. A better definition says that an ideal mesh enables you to solve the
problem with the desired level of accuracy in the time allotted. Note that for structural
problems, the coarser the mesh, the stiffer the structure; this means that a coarse mesh will
tend to underestimate deflections and overestimate natural frequencies.

Nodal DOFs, such as displacement and temperature will converge before other derived
variables, including stress and heat flux. Mesh convergence is also dependent on loads,
constraints and geometry. When any of these change, mesh convergence has to be tested again.

There is an opposing phenomenon to mesh convergence, called truncation error. Computers

only store a finite number of significant figures. In the last digit of a number, rounding off
occurs. This loss of precision is very small, but can accumulate for problems with increasing
numbers of related computations. This cumulative error is called truncation error. As the mesh

64
becomes more fine (or for transient problems, as the time-step becomes more small),
truncation error becomes more significant.

Highly refined meshes can also introduce further errors if the elements are poorly shaped.

Local mesh refinement is increasing the mesh density in a region of interest, to save
computation time. As a rule of thumb, the ratio between the largest and smallest element
length should not be greater than 100, and the ratio between the largest and smallest area
should not be greater than 10 000.

Regions of high stress (or heat flux, etc) gradient need a finer mesh to capture their true
behaviour. In structural analyses, the following sections are often good candidates for mesh
refinement:

65
Mesh Structure
How the mesh is structured will affect the solution. The mesh should be arranged to best fit the
expected solution profiles, which can be very difficult for complex (3D) geometry and loading.
There are two categories of mesh:

Free Meshing: Mesh is generated by an automated algorithm that progresses through

the geometry.
Mapped/structured: Size controls are used on all edges (or surfaces) to force the mesh
to a particular pattern. Mapped meshing:
o requires setting all of the size controls
o usually requires subdividing the geometry into mappable shapes
o is restricted to square type elements for most software
o is closely related to Sweep meshing
o can be very time-consuming

Generally structured meshes are better, though periodic patterns in the modelling can
sometimes bias the solution, and decrease the accuracy, as demonstrated in the image below.
For the mesh on the left, the asymmetry of the mesh has caused asymmetry in the solution,
even though the physical scenario is symmetrical.

66
Element Shape Quality
Remember that shape functions represent an approximation of the distribution of a scalar or
vector field inside an element. These shape functions are most accurate when the element has a
specific shape. For triangular and tetrahedral elements, the ideal shapes are equilateral. For
rectangular and rectangular prism elements, the ideal elements have right-angles and are also
equilateral.

Element shape quality can be best understood in terms of the purpose of elements. Elements
are fitting the solution field with the approximated function (linear or quadratic). Thus, the
larger the element and the more 'non-ideal' the shape (e.g. skewed, stretched) then the larger
the area that is being solution fitted or the more distorted the area is that is being solution
fitted. This means that the approximation will not fit the actual solution over the element area
very well. This is also why the elements should try to follow the predicted flow of the solution
field. Errors with the numerical integration (used for finding the stiffness matrix) of badly
shaped elements can further add to this.

There are a number of variables used to describe the quality of a mesh. Each of these relate to
the shape of the elements and their alignment with respect to neighbouring elements. They are:

Skewness
Aspect Ratio
Warping Factor
Maximum Corner Angle
Jacobian Angle

Some mesh quality variables that are particularly important for CFD are:

Orthogonality
Volume difference between neighbouring cells.

Mesh quality matters both before and after loading. For some problems, deflection during the
simulation can deform the mesh and decrease the element quality. Shell and Beam elements
should experience no more than 30 degrees of flexure:

Shell and beam elements also have recommended limits on the initial curvature to thickness
ratio for the geometry that they represent. High curvature across an individual element will lead
to inaccurate solutions. Refining the mesh will help, however it may be necessary to use solid
elements. This is due to the stresses developed by the high curvature to thickness ratio not
being properly approximated using the shell or beam formulations.

67
Skewness is one measure of how much the shape of an element is from its optimal shape.
Elements with poor skewness look sheared. There are two methods for determining
Skewness:

The first method is for triangular and tetrahedral elements. For triangular elements, a
circle is constructed that intersects all the nodes of the element. An equilateral triangle
is then created that also just fits inside the circle. This equilateral triangle represents the
optimal cell area. The skewness is:

A A
A

For a tetrahedral element, the same method is used, but with a sphere, and using the
volume.
The second method is based on the optimal internal angle, . The maximum and
minimum internal angles of a cell are considered, and the skewness is given as:

max ,
180

The mesh quality relates to the skewness according to the following chart:

0.00 0.25 0.50 0.80 0.95 0.98 1.00

Excellent Very Good Good Acceptable Bad Unacceptable

Aspect Ratio, for a 2D element, is the ratio of the shortest to the longest side. As a general rule,
this should be less than 5.

Averaged vs Unaveraged solution fields

In the theory chapter we saw that the DOF values are calculated at the nodes. However, derived
values like stress and strain are calculated inside the element at the integration points. How
come FE software still presents the results at the nodes? This is because the derived values are
extrapolated out to the nodal points using linear shape functions. If a node is connected to
multiple elements then each element will lead to an extrapolated value at the node. From this
there are two ways of presenting results for derived values. The first is to average the calculated
values at each node weighted based on the area/volume of the adjoining elements. This is often
referred to as the nodal or averaged solution. The second way to present the results is to plot

68
(or list) them separately for each element. This leads to discontinuities in the solution field at
element boundaries and nodal points. These solutions are referred to as the element or
unaveraged solution.

A comparison of the averaged and unaveraged solutions can provide an indication of mesh
independence. Large discontinuities in the element solution means that greater refinement or a
change in mesh shape/structure is needed in those regions. Note that physical singularities will
always have large discontinuities in the solution field that cannot be removed.

Mesh Error Estimates

FEA software packages, such as ANSYS, usually contain in-built methods to estimate the error
produced in the results by the mesh. In ANSYS, these are known as the Element Error Estimates.
These estimates are based usually on differences in derived values when using averaged (nodal)
and unaveraged (element) solutions.
ANSYS has several different methods available. The main methods are:
Nodal-Component value deviation: Difference between the averaged stress (or heat
flux) at a node and the unaveraged value. (SDSG for structural and TDSG for thermal
analyses in ANSYS)
Element Energy Error: Error in the element strain energy (or thermal energy), calculated
using the difference between the averaged stress (or heat flux) at a node and the
unaveraged value. (SERR for structural and TERR for thermal analyses in ANSYS)
Percent error in the energy norm: Percentage of the element energy error normalised
to the total strain (or thermal) energy. (SEPC for structural and TEPC for thermal
analyses in ANSYS)

69
Mesh Independence
All of the aforementioned mesh attributes play a key role in determining the independence of
the obtained results from the mesh. The simplest approach is to double the mesh refinement
and compare the solution, although other checks are needed to test the mesh structure and
element shape quality. Mesh independence is affected by the geometry as well as the location
of the applied loads and constraints. What works for one model (or load case) might not work
for another.

How important is mesh independence? The question that you should ask is how important are
accurate results? The accuracy and level of detail that you require to get out of the model sets
the level of effort for creating the model. That is, if only a quick and nasty solution is needed,
then there is no need spending weeks checking every possible mesh configuration to determine
the most optimal. Furthermore, the final results depend not only on the mesh, but also on the
assumptions and simplifications, chosen analysis method, accuracy of the boundary conditions,
material properties, etc. All aspects of the model need to be accurate representations of the
physical system to obtain accurate results. The uncertainty, , has many sources:

!" " #$ % & "# % ' ( ) #$ * ) + % , $+ + -" .+ % &"+/

An ideal mesh enables you to solve the problem with the desired level of accuracy in the time
allotted for the project.

Mesh Transitions
In complex models it is often necessary to join

Elements of different type

Elements of different order

These mesh transitions present a complication to modelling. The nodes that bound the element
may have incompatible degrees of freedom, or the edges that bound elements can have
incompatible numbers of nodes.

Mesh connecting can be done by either using extra elements to tie the nodes together, (such as
rigid links or bars, extra plate elements) or using Constraint Equations to link the degrees of
freedom together between the nodes. It is important, however not to add extra stiffness or
mass to the model which might alter the response.

70
The following examples are an incorrect approach because continuity is not maintained along
the element boundary.

Complicated methods can be used to include extra elements and make a mesh connection. A
few examples are shown in the following images:

Multi-Point Constraint equations (MPCs) are often required. These are expressed in the
following form of equation, where the 0& degree of freedom is unconstrained or unconnected
and so is made dependent on the other DOFs, 01 :

2& 0& % 3 21 01 0

How to Create a Good Mesh

There are two ways to create an FE mesh: 1) by meshing a structural model, and 2) by directly
creating the nodes and elements (without a structural mode). Meshing a structural model is the
most commonly used approach. Most FEA software packages have advanced auto-mehsers built
into them. These meshers progress through the geometry starting at particular locations, which
are referred to as seed points. Depending on the algorithm, multiple seed points may be used.

Auto meshing can be done completely automatically without any user input. In this case, the
software will choose the mesh size (usually, one size fits all; default values). This is easiest, but
generally results in the poorest quality mesh (meshing software is continually improving, but it
has its limits). It is best practice to manually define size controls on volumes, areas and/or lines

71
to help control the mesh. Modern FEA packages have a range of special meshing tools that help
with auto-meshing. These tools can sweep and sub-divide regions to help create uniform
meshes with good quality elements. Caution must always be exercised when using these auto
features as their success can be a bit hit-or-miss.

A technique that can help with mesh generation is to manually slice or combine regions of the
model before meshing. Try to create uniform and prismatic regions that match the shape of the
elements. The best approach is to plan ahead and create your geometry with meshing and
loading in mind.

Another technique to help with meshing is to concatenate lines and areas (called virtual
topology in ASNSYS Workbench). Concatenating lines and areas tricks the meshing program into
ignoring the divisions, that is, it treats concatenated objects as one. Auto-meshers usually
control the node placement and edge sizing based on point, line and area geometry. Using
concatenation allows features to be ignored (in terms of setting size controls), which can lead to
smoother more a continuous mesh.

Nodes and elements can be created directly without any solid modelling. This approach can
create high quality meshes, but is extremely painful and time consuming. Badly shaped
elements of an auto-generated mesh can also be manually modified to improve their quality.

72
Final Words
Creating an acceptable mesh is often a process of trial and error (the amount of iterations
decreases with increased user experience). Think ahead when creating your model to help make
the meshing process easier.

Run mesh independence checks, use the built in mesh quality measures, review the solution
fields, and always verify the model!

73
74
Chapter 5
Heat Transfer

75
Heat transfer in real designs, such as heat exchangers and combustion chambers, is often
analysed using FE methods. The governing principles of heat transfer can be summarised as:

The temperature of an object is determined by the thermal energy in it. The relationship
between the change in the amount of thermal energy in an object and the temperature
of that object is defined by the specific heat, .
Energy flows from one object to another when they are at a different temperature.
Thermal energy can flow using three different mechanisms:
o Conduction: through a solid
o Convection: from a solid to a liquid
o Radiation: from one object to another via electromagnetic waves

A thermal or heat transfer analysis requires thermal elements. These elements have one DOF at
each node (Temperature), and the loading is expressed in terms of energy flow rates, or Powers.
Thermal elements can be 1D, 2D and 3D and can have mid-side nodes, just like structural
elements.

In this chapter, the theory of Heat transfer will be discussed, as well as the methods for
implementing it in FEA analyses.

Heat Transfer Theory

Conduction
The flow of heat inside a domain is governed by conduction. For a continuous problem,
conduction is proportional to the temperature gradient at a point and the conductivity. This is
called Fouriers Law of conduction:

In this equation, is the scalar field temperature, is the flow of heat per unit area in the
direction (heat flux) with units J/m2-s, and is the conductivity in the direction, . If is
isotropic, this reduces to the following vector representation, where is a vector in the
opposite direction to the temperature gradient:

The energy balance of a volume works according to the following simple equation:

Energy in Energy out + Energy generated = Energy Stored

By considering a infinitesimally small volume, the energy balance equation can be used along
with Fouriers law to derive the governing equation for heat transfer in 3D (details omitted
here). The resulting differential equation is known as the heat equation:

!+ ## !+ %% ! + & = '(
" " $ $ )

This equation includes the effect of conduction in three dimensions through a solid as well as
the materials capacity to store heat energy. Internal heat generation is also included through

76
the term &, which represents energy generated in the domain and is a scalar field with units
J/m3-s. Flow of thermal energy in a fluid works differently because the temperature can advect
away from a control volume, and because the flow of energy can induce convection currents:
flow of the fluid. Problems like this are dealt with using CFD, which is outside the scope of this
course.

If is isotropic and the object is solid, the governing differential equation simplifies to:

* + & = '(
)

Convection
Convection is energy flow between a surface and a fluid. It is described by the following
equation, where + is the surface temperature, , is the bulk temperature of the fluid, . is the
convective film coefficient and is the flow rate into the fluid per unit surface area in J/m2-s:

= . / + ,0

Convection can occur due to the forced flow of a higher or lower temperature fluid over a
surface. It can also occur naturally due to buoyancy forces developed in the surrounding (initially
still) fluid due changes in density of the fluid with temperature.

Radiation
Radiation is the transfer of heat energy via electromagnetic waves. It is the direct result of the
random movements of atoms and molecules in matter.

Net radiation flow from one surface (1) to another surface (2) is governed by the relationship:

12* = 3452 62* / 7

2 7
*0

In this equation:

4 is the Steffan-Boltsman Constant, which describes the radiation power of an ideal

black body
3 is the emissivity, which represents the radiation power of the surface as a proportion
of the radiation power of an ideal black-body
6 is the form factor, which represents the alignment and proximity of the two
interacting surfaces
1 is the flow rate of heat in J/s

Radiation is related to the fourth-power of temperature, which creates an essential non-linearity

in the problem. For radiation to the environment from a hot object, the form factor would be
set to 1.

77
Finite Element Models of Heat Transfer
When the differential heat equation is discretised, each node has one degree of freedom:
Temperature. The differential equation is formulated into a finite element model represented
by the following matrix equation:
8 = 9:; + <:
Here : is a vector of nodal temperatures, and 8 is a vector of the energy flow at each node in
J/s or W. The C matrix (thermal damping matrix) contains the specific heats, and K (thermal
stiffness matrix) contains the conductivities.

Heat transfer problems can be considered analogous to structural problems, as shown below:
Structural Thermal
= = > + >; + ?>@ 8 = :; + :
Displacement (DOF) (m) Temperature (DOF) (K)
Force (N) Heat flow rates (J/s, W)
Pressure (MPa, N/m2) Heat Flux (applied) (J/s/m2, W/m2)
Strain (m/m) Thermal Gradient (K/m)
Stress (MPa, N/m2) Heat Flux (calculated) (W/m2)
Reaction Forces (N) Reaction Heat Flow (N)
Elastic Modulus (MPa, N/m2) Conductivity (W/m/K)
Stiffness (N/m) Heat Transfer Coefficient (W/m2/K)
Damping Coefficient (N-s/m) Specific Heat (J/kg/K)

The procedure for a thermal analysis is much the same as for a structural one. The geometry
needs to be modelled with appropriate material properties, and then loads (e.g. heat
flows/fluxes) and constraints (temperature) are applied. One difference between structural
problems is that thermal analyses commonly use convection and radiation boundary conditions.
To continue the analogy, these BCs can be thought of as linear and non-linear, respectively,
spring foundations. Lastly, the model is meshed and solved with all the same considerations as
for a structural analysis. However, this time the physics is heat transfer, and thus the user must
be experienced with this type of problem to best conduct the analysis.

Like structural models, thermal analyses are subject to many simple modelling errors. Most of
which are analogous to the structural situation. For example, if a problem has one heat source
(such as an applied nodal heat flow or internal heat generation), but has no heat outlet (such as
prescribed nodal temperatures or convection loads); then the domain will continue to fill up
with energy, so the steady state solution will be an infinite temperature. This is the structural
equivalent of rigid-body motion in an unconstrained problem, which leads to infinite
displacement.
There are also thermal gradient/flux singularities caused by such features as:
Point heat sources
Re-entrant corners and cracks in the mesh
Poorly shaped elements
Again, these are similar to situations with stress analysis, where point loading and re-entrant
corners in the geometry cause singularities.

78
Non-linear Problems
A heat transfer problem becomes non-linear when the:

Specific Heat Matrix

Thermal Conductivity Matrix, or
Equivalent Nodal Heat Flow Vector
...are functions of temperature, or mathematically:

8/:0 = /:0:; + /:0:

These equilibrium equations must be solved iteratively. Non-linear cases help model systems
that include:
temperature-dependent material properties
temperature-dependent film coefficients
use of radiation elements
temperature-dependent heat sources (flow of flux)

Solving Transient Problems

The thermal analysis is transient when it:

is concerned with the behaviour before steady state is reached

has no steady state solution
has boundary conditions that vary with time
Only in transient analyses are the temporal derivatives (AA)) and specific heats relevant.
Solving the discretised transient problem involves dividing the time into discrete time-steps,
separated by ).

For this transient analysis, an implicit Euler time-integration scheme is presented. This scheme is
equivalent to a general trapezoidal rule where the derivative used over the time-step is a
weighted average of the derivative at the beginning and end of the time-step:

: D2 = : + /1 F0):;G + F):;GDH

At each time-step, the general equation still applies:

8 = :; D2 + : D2

This yields the equation for solving transient problems, which is unconditionally (always) stable
for 0.5 < F < 1:
1 1 1F
+ !: D2 =8+ : + :; !
F) F) F

Selecting an appropriate time-step, ), is important. If the time-step is too low, then non-
physical oscillations can occur, called thermal undershoot. If the time-step is too large then
thermal gradients will not be captured adequately. A good method for setting the time-step is
achieved by setting the Fourier Modulus equal to 1 ( relates to element size):
4 )
6M = =1
'(/ 0*

79
Coupled Field Analysis
Solids will expand as their temperature increases. This can influence the stresses in a solid if:

the solid is constrained to prevent expansion,

the thermal expansion coefficient is not uniform throughout the solid, such as when
two materials are joined together, or
there is a temperature gradient. If an object is at a uniform temperature, then the
expansion will also be homogeneous. However, when there is a temperature gradient,
the non-uniform strain will cause stresses to be induced in the solid.

ANSYS allows for two methods to deal with these types of analysis, the Sequential Coupled
Field (weak) method and the Direct Coupled Field (strong) method, which each have their
advantages:

The Sequential coupled field method is an older method. It involves two steps with different
element types. This method is efficient for running many thermal transient time points but few
structural time points, and with ANSYS is easily automated with the use of input files.

Step 1: Do the heat transfer analysis. Model it with thermal elements and apply thermal
loading.
Step 2: Perform a structural analysis.
o Use structural elements
o Define the material properties including the thermal expansion coefficient
o Apply the structural loading, including temperatures from thermal analysis
o Solve and review the results

For some coupling situations, the sequential method can be more efficient and flexible, because
you perform the two analyses independently of each other. In a sequential thermal-stress
analysis, a nonlinear transient thermal analysis can be followed by a linear static stress analysis,

80
where the nodal temperatures from any load step or any time-point in the thermal analysis can
be used as loads for the stress analysis.

The Direct coupled field method uses one element type to solve both physics problems and
allows true coupling between thermal and structural phenomena. This method may carry
unnecessary overhead for some analyses. For this analysis the temperature represents an
additional degree of freedom for all the elements. The two aspects of the system, structural and
thermal, are combined to form the following discretised system:

? P >@ P >;S + O UT > =

O QR S + O TQ R QO Q = R S
P P :@ TU
:; P T : 8

where Ctu and Kut are coupling terms between the different physics.

Direct coupling is advantageous when the coupled field analysis

is highly nonlinear
is best solved in a single solution using a coupled formulation such as for the following
examples:
o piezoelectric analysis
o conjugate heat transfer with fluid flow
o circuit-electromagnetic analysis

A harmonic thermal-structural analysis (harmonic analysis is finding the steady-state periodic

response to periodic loading, see Chapter 7) can only be conducted using the direct-coupled
method. Notice that this coupled problem has non-symmetrical matrices and so is not well
conditioned.

81
82
Chapter 6
Non-linear Analysis

83
Many problems in FEA include non-linearities. There are three sources of non-linearity in a
structural problem. These are:

Material non-linearities: Non-linear stress-strain constitutive behaviour

Geometric non-linearities: Change in shape of the structure during loading changes the
stiffness
Changing status: Abrupt changes in stiffness due to contact, element birth/death

The technology and methods for solving systems of linear equations are well advanced, and are
very suitable for computers. Hence, the main challenge in a non-linear FEA is to use linear
solution methods to solve the non-linear problem.

This chapter will discuss situations where non-linearities arise, and methods that can be used to
solve them. The examples that are referred to will again be based on solid mechanics problems.

Non-linearity
A FEA problem is linear when it can be expressed with a system of linear equations, i.e. an
equation like

= ,

when K is constant. This equation indicates that the outputs are proportional to the degrees of
freedom (DOFs) i.e. doubling the load will double the displacement, or for a thermal problem,
doubling the temperatures will double the thermal energy flows.

A problem is non-linear when it cannot be expressed this way. Non-linear problems can be
thought of as having variable stiffness, such as in the equation,

= ,

where K is not known without solving for .

Solving Non-Linear Finite Element Models

Consider the load-deflection curve for a uniaxial tensile test of a non-linear material.

Load, F

Deflection, u

The load-deflection curve can be described as a non-linear function of u:

84
How do we turn the non-linear equation into a linear equation?

There are main two methods for handling non-linearities in the stiffness matrix, i.e. where the
stiffness is a function of the deflection or strain. The first approach is the Purely Incremental
Method. This method is less used in practice; however, it is instructive to briefly look at how it
works. The second and superior approach is the Newton-Raphson Method. This method is
implemented in most FEA packages.

Method 1: Purely Incremental Method

Consider a small increment of force, dF, and displacement, du, for the non-linear load-deflection
curved described above. If at the current load level, F, the displacement, u is known; the internal
force f(u + du) for the new external force F + dF can be approximated by the linear function:

d
+d = + d
d

This gives:

d
+ d = +d
d

From which we can assume:

d
d = d
d

or:

u d =d

For a 1D problem, i.e. a typical load-displacement curve, the tangent stiffness, KT , is simply the
changing slope of the load-displacement curve.

In matrix (FE) formulation:

or:

where KT is the tangent stiffness matrix (sometimes called the Jacobian stiffness matrix) and the
subscript i refers to each incremental step. The above form is now a linearised system that by
solving incrementally represents the original non-linear one.

The applied load is divided into increments. Over each interval, the stiffness matrix is assumed
constant and calculated using the nodal displacements for the current configuration. The
displacements for the next step, Ui+1, are then calculated and used to update the stiffness
matrix.

85
The Problems with this method are that errors accumulate with each load increment and the
final results are out of equilibrium. If the step size is sufficiently small then the calculated
response will match the actual non-linear behaviour. Overall the purely incremental method is
robust and has no convergence issues (as there is no convergence checks). Though as stated,
the increment sizes are required to be very small to keep the error low.

This approach is an explicit one and is similar Eulers method for solving ordinary differential
equations. We need an improved approach that can correct for the error at each step.

Method 2: Newton-Raphson Method

The Newton-Raphson method uses a series of equilibrium iterations to solve an equation (or
system of equations in the case of a finite element model). Recall the introduction to the
method in the Background on Numerical Methods section at the end of Chapter 3. To solve
for = , we can use:

= +

The equation can be rearranged to:

or:

In matrix form the above equation can be written as:

!=

This equation is very similar to that of the purely incremental method. Except now equilibrium is
enforced at each increment to correct the build-up of errors. The Newton-Raphson method is
therefore an implicit approach. The equilibrium correction is applied through the nodal restoring
force vector f(Ui). The restoring forces are calculated from the internal stresses at the current
displacement configuration. Due to the linearisation and assumption of a constant stiffness

86
matrix for each increment, the internal restoring forces will not match the applied forces, F.
That is, the system will be out of equilibrium. The nodal restoring forces are determined by:

"
= = # $% & '(

The tangent stiffness matrix is obtained in various ways, depending on the type of non-linearity.
This is further discussed in the next sections.

The Newton-Raphson method is implemented by first calculating the tangent stiffness matrix,
"
) , and restoring force vector, ) , for the current displacements, Ui. Using the applied
forces vector, F, the linear system can be solved to give U1. The procedure is repeated (with the
stiffness matrix updated) until the solution converges.

If the full applied forces are used at the first step, then algorithm looks like the image below (in
the case of a 1D load-deflection curve).

The applied loading minus the restorative loading, * , is called the residual and is a
measure of how much the forces are out of balance in a structure. Iterations continue until this
residual is below a tolerance value. A plot of the residuals can be used to debug a failed
solution.

There are several variations of the Newton-Raphson procedure. The most commonly used in
FEA packages is an Incremental Newton-Raphson Method whereby only part of the load, F, is
applied at each step. This helps to improve convergence and is a must for path-dependent non-
linearities (such as plasticity).

At this point it is necessary to introduce the concept of load steps, sub-steps, and equilibrium
iterations (terminology as used in ANSYS similar in other software). Load steps are changes in
the externally applied loading on the model; they are defined by the user. Equilibrium iterations
are each time the Newton-Raphson procedure iterates until the solution is converged. A single
load step with four equilibrium iterations is depicted in the previous figure. Sub-steps are the
subdivisions of a given load step, as used in the Incremental Newton-Raphson Method. They can
be automatically set by the software or directly defined by the user. The use of sub-steps helps
to improve convergence.

87
Convergence is not guaranteed in all cases, sub-step algorithm will only converge if the starting
point is within the radius of convergence. Most software packages use automatic stepping
algorithms that predict and controls the step size for all sub-steps in a load step. This includes
bisection controls that automatically decrease the sub-step size if a given sub-step fails to
converge. Though, convergence is still not guaranteed.

Each load-step and sub-step are associated with a value of time. In most static analyses time is
simply used as a counter and does not mean actual, chronological time.

A variation of the Incremental Newton-Raphson Method is the Modified Newton-Raphson

Method. This approach updates the stiffness matrix less frequently, such as only at the start of
each sub-step. The solution converges more slowly (i.e. requires more iterations) than the
standard method, but has less matrix formulations, which are computationally expensive.

Geometric Non-linearities
Geometric non-linearities occur in a structure or component due to the changing geometry as it
deflects. The stiffness changes as the shape changes and/or the material rotates. Previous
discussion in this course has been based on small rotation and small strain theories. Under
these assumptions strains are calculated based on the undeformed shape (i.e. engineering
strains) and original body location (i.e. no rigid body motion). The benefit is that the governing

88
equations can be simplified and linearised, and DOFs can be de-coupled. Small strain and
displacement theories are suitable for many engineering materials and structures as the
deflections/strains are not large enough to cause significant changes in the stiffness. There are
however, many situations when it is necessary to consider geometric non-linearities.

There are three main types of geometric non-linearity (though there are others):

Large strain assumes that the strains are no longer infinitesimal. Shape changes (e.g.
area, thickness) are accounted for. Deflections and rotations may be arbitrarily large.
Large rotation assumes that the rotations are large, but the mechanical strains (those
that cause stresses) are small. The structure is assumed not to change shape except for
rigid body motions.
Stress stiffening assumes that both strains and rotations are small, but internal stresses
can contribute to the stiffness.

Geometric non-linearities arise due to inclusion of non-linear terms in the strain-displacement

relationships (the matrix B).

In most cases the tangent stiffness matrix is given by:

,& = ++ &

where K is the usual stiffness matrix and S is the stress stiffness matrix. The stiffness matrix and
stress stiffness matrix are calculated using the appropriate strain-displacement laws, through
matrix B, and the material law, through matrix D. The stress stiffness matrix may or may not be
used, depending on the analysis.

Deflections, strains, and stresses are determined using the theories appropriate for each type of
non-linearity (not part of this course see any standard text on non-linear structural
mechanics). The outputs for large strain analyses will be in the form of true strain and true
stress (not engineering stress and strain as for all other analyses).

Stress Stiffening
An interesting behaviour that arises with geometric non-linearities, is stress stiffening. Stress
stiffening is the stiffening (or weakening) of a structure due to its stress state. The stiffening
effect normally needs to be considered for thin structures with bending stiffness very small
compared to axial stiffness, such as cables, thin beams, and shells. Stress stiffening couples the
in-plane and transverse displacements.

89
To understand the coupling effect, consider a thin beam (with low bending stiffness)
constrained at both ends and loaded transversally at the mid-point. Under linear theory, the
shear loading and the bending moments are associated with the lateral deflection only, and the
axial loading is associated with the axial deflection only. These two systems are assumed to be
de-coupled. However, this assumption is not always valid. As the beam deflects laterally, it will
also try to contract axially. This generates membrane stresses (if the contraction is constrained),
which resist the lateral deflection (bending) and increase the overall bending stiffness. De-
coupling the different directions is the limiting case i.e. as the lateral and axial deflections
approach zero. The more the beam deflect (and the less the bending stiffness compared to the
axial stiffness), the less valid these approximations become.

The stiffness of a member is increased by tensile membrane stresses and decreased by

compressive stresses. Sufficiently large compressive membrane stress reduces the bending
stiffness to zero; that is, the structure buckles. Buckling will be dealt with in Chapter 8.

A problem need not be non-linear for stress stiffening to play a role. The effect is essential in a
linear buckling analysis and can strongly affect the vibration frequencies of a tensioned string for
example. In these cases, the stress stiffening can be considered as a prestress.

Material Non-Linearities
Material non-linearities occur because of non-linear relationships between stress and strain.
That is, the stresses developed are a non-linear function of the strains in the body. The most
well known material non-linearity is plastic yielding. There are however, a range of different
material non-linear behaviours, including many variations of how a material yields and
plastically deforms (e.g. strain rate dependency, strain hardening). An additional factor with
material non-linearities is that most are path-dependent. This means that the stress also
depends on the strain history and not just the current strain values. The exceptions to this are
non-linear elasticity and hyperelasticity, which are path-independent.

90
There are six main types of material non-linearity (though there are others):

Rate-independent plasticity is characterised by the irreversible instantaneous straining

that occurs in a material.
Rate-dependent plasticity allows the plastic strains to develop over a time interval
(sometimes termed viscoplasticity).
Creep is an irreversible straining that occurs in a material and is rate-dependent so that
strains develop over time. The time frame for creep is usually much larger than that for
rate-dependent plasticity and the applied loads are typically much lower. Creep is
significantly more severe at elevated temperatures.
Non-linear elasticity allows non-linear stress-strain relationships, but all straining is
reversible.
Hyperelasticity is defined by a strain-energy density potential that characterises
elastomeric and foam type materials. All straining is reversible.
Viscoelasticity is rate-dependent material characterisation that includes a viscous
contribution to the elastic straining.

Material non-linearities are accounted for through the stress-strain relationship matrix D. In
most cases the tangent stiffness matrix is given by:

, = # $- . , $'(

where . is a matrix of the jacobian of the material constitutive law:

d& ,
. , =
d,

Plasticity theory provides a mathematical relationship that characterises the elastoplastic

response of materials (not part of this course see any standard text on non-linear structural
mechanics). There are three key components in modelling plasticity: the yield criterion,
hardening rule, and the flow rule.

The yield criterion determines the stress level at which yielding is initiated, for example
the Tresca criterion or Von Mises criterion. The Von Mises criterion is the most
commonly used.
The hardening rule describes the changing yield surface with progressive yielding (e.g.
Kinematic hardening or isotropic hardening). It can be thought of as the way in which
the resistance to plastic flow increases/decreases with continued loading. The hardening
rule also incorporates how the yielding will change with each excursion from tensile to
compressive yielding (Bauschinger effect).
The flow rule determines the direction and amount of the plastic straining and thus
provides a way of calculating the plastic strains.

All of these calculations are carried out at each Newton-Raphson iteration step. They may
require their own iteration procedure.

91
Changing Status Non-Linearity
Changing status non-linearities are situations where a particular feature of the model changes
due to the loading and results in an immediate re-formulation of the stiffness matrix.

Contact is one of the most frequently encountered non-linearities, and is also one of the most
difficult non-linearities to process. It is a changing status non-linearity because the stiffness can
change suddenly depending on objects being in contact or not. Only once two objects touch the
contact region between them needs to be modelled to enforce impenetrability in a stress
problem, or to allow conduction in a thermal analysis.

Physical surfaces that are in contact do not interpenetrate (known as contact compatibility).
They can transmit compressive normal forces and tangential friction forces but do not transmit
tensile normal forces. They are therefore free to separate and move away from each other.

A contact changing status non-linearity actually involves three possible statuses, with different
stiffnesses: Open, closed and sliding, or closed and sticking. ANSYS allows for four statuses,
because the open status has both near-field and far-field options (for computational reasons).

The stiffness changes abruptly as two parts come into or out of contact. That is, when the
contact is open, there is no load transfer and thus zero stiffness. Once in contact, load can be
transferred and there is a stiffness. This is why contact is a changing status non-linearity.

92
Factors that make contact analyses complicated are:

The region of contact is typically unknown at the start of the analysis (and may change
during the analysis)
Many contact problems include friction
Parts might be unconstrained except for contact with other parts (unconstrained free
bodies have zero overall stiffness)

There are five types of contact in ANSYS. Bonded and No Separation contact are. Frictionless,
Rough and Frictional contact are non-linear and require multiple iterations.

Contact Type Iterations Normal Behaviour Tangential Behaviour

Bonded 1 No Gaps No Sliding
No Separation 1 No Gaps Sliding Allowed
Frictionless Multiple Gaps Allowed Sliding Allowed
Rough Multiple Gaps Allowed No Sliding
Frictional Multiple Gaps Allowed Sliding Allowed

Frictional contact uses Coulombs law:

-/0120- /3 5 06"7/3

where is the coefficient of static friction. Once the tangential force exceeds the above limit;
sliding will occur.

Contact Algorithms
There are four solution methods for contact problems used in ANSYS. These are:

Pure Penalty
Augmented Lagrange
Normal Lagrange
Multi-Point Constraint

The Pure Penalty method enforces contact compatibility by using a spring between the
contact surfaces. A simple linear spring relationship is formulated for the spring when the
surfaces have the closed contact status:

06"7/3 = 806"7/3 9:202-"/- 60

93
where Fnormal is the contact normal force, x is the penetration depth, and knormal is the normal
contact stiffness (also known as the penalty parameter).

Under this formulation, some interpenetration of the surfaces will occur (i.e. compatibility is not
truly satisfied). Ideally for an infinite contact stiffness, there would be no penetration. However,
high values for the contact stiffness create convergence issues. As long as the penetration is
small, the results will be suitably accurate.

A contact force is required mathematically to generate a contact force at the interface. Thus, the
penetration must be greater zero than zero for equilibrium. Therefore, for best accuracy, the
goal is to minimize the amount of penetration that occurs at the contact interface.

The contact stiffness is the most important parameter in penalty-based methods. The ideal
choice of contact stiffness, k, is dependent on the problem and some trial and error may be
required. When using the pure penalty method, a higher stiffness increases accuracy, but is can
negatively affect convergence, as high stiffness can also cause oscillations.

Usually, the contact stiffness is calculated based on the stiffness of the underlying elements, as
shown in the equation bellow. In ANSYS, for bonded and no-separation cases, FKN is given a
value of 10. For other cases in general, FKN = 1.0 is used. For particularly flexible scenarios, FKN
= 0.1 is a good choice.

8;60-/;- = < =8 0>2"3? 01 @

Pure penalty method is always used for the tangential contact that is bonded or rough, as:

-/0120- /3 = 8-/0120- /3 9A3 > 01

The Augmented Lagrange method is also a penalty based method. It is less sensitive to the
magnitude of the contact stiffness than pure penalty because it incorporates an additional
constant in the equation. However, the method can require more iterations than pure penalty.

= 8 9:202-"/- 60 + B

The Normal Lagrange method includes the normal contact force directly as another DOF that
must be solved for. This method has two advantages:

It enforces zero/nearly zero penetration

It does not require normal contact stiffness

94
However, this method:

Can be more computationally expensive

Can have convergence difficulties related to over constraining

Chattering is an effect that happens in normal Lagrange methods, when the status alternates
rapidly between closed and open. This does not happen in penalty-based methods because
there is not a discontinuity in the contact force function.

Multi-point constraint (MPC) methods add constraint equations that tie nodes together. This
method is not penalty or lagrange multiplier based, and is only applicable for bonded and no-
separation contact. Large deformations are also supported.

For bonded contact ANSYS uses a Pure Penalty formulation with large contact stiffness by
default. MPC formulation is also a good alternative for bonded contact because of its many nice
features.

For frictionless or frictional contact, consider using either the Augmented Lagrange or Normal
Lagrange methods. The Augmented Lagrange method is recommended because of its attractive
features and versatility. The Normal Lagrange method can be used if the user does not want to
deal with contact stiffness values and wants zero penetration. However, for the Normal
Lagrange formulation, only the direct solver can be used, which may limit the size (number of
nodes) of the models that can be efficiently solved.

Contact Pairs
For some of the above solution methods (the ones that apply an opposing force) a contact
surface and a target surface must be defined.

95
The target surface is treated as continuous geometry and is passive in enforcing contact
compatibility. Penetration through this target surface is calculated at detection points on the
contact surface, and forces are applied to the contact surface to oppose this penetration. The
detection points can be the nodes. The detection points can also be the integration points,
which means more detection points and greater accuracy. In ANSYS:

Pure Penalty and Augmented Lagrange Formulations use integration point detection
Normal Lagrange and MPC Formulation use nodal detection (normal direction from
Target)

Nodes (or integration points) of the contact surface cannot penetrate into the target surface (or
at least, ANSYS tries to minimise the penetration).

There are two ways to handle a pair of surfaces that may contact: Asymmetric contact, and
symmetric contact.

Asymmetric contact. For asymmetric contact, one of the surfaces is the contact surface, and the
other is the target surface. Guidelines for choosing which surface should be the contact, and
which should be the target are:

The surface with the coarse mesh should be the Target surface
Stiffer surfaces should be the Target surface
Larger surfaces should be the Target surface
The lower order surface should be the Target surface
If one surface is convex surface and the other is flat (or concave surface), the flat (or
concave) surface should be the Target surface.

96
Generally asymmetric contact is the most efficient way to model surface-to-surface contact (as
there are less required calculations and status checks). However, sometimes it will not perform
satisfactorily, because no clear distinction exists between target and contact surfaces or perhaps
because both surfaces have coarse meshes.

For asymmetric contact in ANSYS:

ANSYS Workbench can automatically perform the surface designation (Auto-

Asymmetric) or,
The user can designate the appropriate surface(s) for contact and target manually
Selection of inappropriate Contact vs.Target may affect results
Reviewing results is easy and straightforward

Symmetric contact. For symmetric contact, two relationships are created: in one relationship
the first surface is the contact and the second surface the target, in the other relationship, the
second surface is the contact and the first surface is the target.

For symmetric contact in ANSYS:

It is easier to set up (Default in ANSYS Workbench)

It is more computationally expensive
Interpreting data such as actual contact pressure can be more difficult

There are important things to note:

Only Pure Penalty and Augmented Lagrange formulations allow Symmetric behaviour
Normal Lagrange and MPC require Asymmetric behaviour
For Symmetric contact, results (in ANSYS) are reported separately for both Contact and
Target surfaces (and the true result is an average of both)
For Asymmetric contact, results are only available on Contact surfaces

Rigid Body Contact

In some cases contact problems can be simplified by considering two general classes:

Rigid-to-flexible
Flexible-to-flexible

In a rigid-to-flexible problem, one or more contacting surfaces are treated as rigid, for example,
metal forming problems. Stresses within the rigid body are not calculated, and the rigid surface
is always the target. This significantly decreases the required calculations.

97
In a flexible-to-flexible problem, both or all contacting bodies are deformable. This is used when
all surfaces have similar stiffnesses. For example, a splined shaft interference fit, where both
parts are equally flexible.

Contact Controls
Most FE packages have a range of extra features to help control the contact analysis:

Convergence
Near/far detection
Interface treatment (e.g. gaps)

Pinball Region is a feature in ANSYS Workbench to detect surfaces that are near/far open. In
other words, it detects contact elements that are likely to contact and those that are not. This
provides computational efficiency in contact calculations. The Pinball also has other uses: In
bonded contact, any gap smaller than the pinball radius will be treated as bonded, and in the
multi-point constraint method it affects how many nodes will be included in the MPC equations.

For frictional or frictionless behaviour, the initial gaps between parts are not automatically
ignored because they may represent the geometry. In some cases if contact is not established,
this may cause rigid body motion (un-constrained geometry), which is not allowed for in static
structural analyses. Interference treatment is designed to solve this problem: the contact
surfaces can be internally offset by a specified amount. This will essentially change the geometry
at the contact interface. This is intended for applications where the required adjustment is small
(e.g. CAD).

Convergence of contact problems can be a major issue. Most software packages proved tools to
help check the model to find the source of the issues. One approach is to use plots of the
Newton-Raphson residuals to see which elements are not converging. From this the mesh may
be refined, contact stiffness adjusted, or step size controls modified.

98
Applying Loads in a Non-Linear Analysis
In non-linear analyses, the deformation of components can be quite large. It is therefore very
important to know which way your loads are pointing at the start of an analysis through to the
end of the analysis. For example, consider a point force applied perpendicular to the end of a
cantilever beam. As the amount of deflection increases the surface normal will change direction
and no longer be aligned with the original line of force. The real structure must be assessed to
ensure that the applied loads behave in the correct way as the body deforms. Several examples
of how the force direction can change are shown in the figure over the page.

99
100
101
102
Chapter 7
Dynamic Structural Analysis

103
Dynamic problems are problems that vary with time. The dynamic differential equations that
govern a problem are more complicated than static differential equations (but they will reduce
to the static DEs if the time derivatives are set to zero). In dynamic problems inertial forces and
damping can play a significant role.

All dynamic problems can be solved using discretisation in the time domain. That is, by dividing
the simulation time into a number of discrete time-steps. The problem with discretisation in the
time domain is that the amount of required computation drastically increases. For periodic
problems (problems where the loading repeats in a cyclic pattern), a different approach can be
taken. These analyses are very important to learn, because vibrations and cyclic behaviour are
common-place in real-world structures.

This chapter will progress through the different types of dynamic analysis and their associated
advantages and limitations. The different dynamic analyses include:

Modal analysis
Harmonic response
Response spectrum analysis
Random vibration analysis
Linear and non-linear transient analysis (implicit methods only in this chapter)

Dynamic Equation of Motion

Dynamic problems have variables that change with time, so both and are functions of time.
Because of this, the time-dependent terms in the system equation are relevant, and the
behaviour becomes dependent on the derivatives of the degrees-of-freedom. The differential
equation of motion for a structural system is second-order and given by:

= + +

In the previous equation, C is the damping matrix, M is the mass matrix, is the velocity vector,
and is the acceleration vector. This equation is easiest to consider for the situation when ,
and are all constant, making it a linear differential equation. In this case, still describes
the conservative elastic reaction to deformation in the structure, the vector now describes
damping, which is an internal non-conservative force that removes energy from the system, and

but rather the other side of Newtons second law equation, =

represents the inertial force in the object (unlike the other terms, this is not really a force
). This representation will
be used throughout this chapter. However, it is important to note that there are other
formulations of dynamics problems (e.g. Lagrangian mechanics).

Different Types of Dynamic Analysis

A full transient analysis will solve the dynamic equation of motion with all of the terms included.
This is mathematically intensive and can lead to long solution times. Various simplifications can
be made to the governing equation in certain situations, in order to reduce the complexity.

104
The simplifications made in each analysis type are:

Modal: F(t) = 0, and usually Damping, C, is ignored.

Harmonic: F(t) and U(t) are assumed sinusoidal.
Response Spectrum: Input is a frequency domain spectrum and the solution is obtained
from modal superposition (i.e. weighted summation of the mode shapes).
Random Vibration: Similar to the response spectrum analysis except that the input is a
probabilistic spectrum and the output is a statistical distribution.
Transient: No assumptions are made. The full equation of motion is solved.

The table below summarises the required input and the output obtained for the different types
of dynamic analysis. It is useful to review this table to help decide what type of analysis is
required for your problem at hand.

Analysis Type Input Output

Modal No load inputs Natural frequencies and mode
Damping can be included, but shapes
is not usually Stress/strain, etc profiles (not
Behaviour must be linear absolute values)
Harmonic Sinusoidally varying loads Sinusoidally varying response
Behaviour must be linear Absolute values for stress/strain,
displacement, etc
Response over a range
frequencies for the input
Response Spectrum Base excitation frequency Maximum response for the
spectrum (displacement/ given input spectrum
velocity/acceleration)
Behaviour must be linear
Random Vibration Base excitation power spectral 1 contour results for
density (PSD) function stress/strain, displacement, etc
Behaviour must be linear Response PSD at single points

Transient (implicit) Time varying loads Full time varying response

Behaviour can be non-linear

Modal Analysis
Every structure has the tendency to vibrate at certain frequencies, called natural frequencies.
Each natural frequency is associated with a certain shape, called a mode shape, which the
structure tends to assume when vibrating at that frequency. The goal of modal analysis is to
determine these modes, i.e. natural frequency and corresponding mode shape, of the object or
structure during free undamped vibration.

There are theoretically an infinite number of modes for a continuous structure. In practical
vibration problems, however, only a few lowest modes are important in finding the vibration

105
response to excitation. Further, the higher modes will have low energy (for the lower frequency
excitation) and are damped out very quickly.

The lowest frequency mode is called the fundamental frequency, or the first mode. The modes
are numbered in order of lowest to highest frequencies, and as the frequency increases, the
mode shapes become more complicated and may involve coupling between different mode
shapes. Note that in FEA software, orthogonal modes at the same frequency are numbered in
sequence and thus the numbering will not match the true mode number.

A discretised model will also have natural frequencies and mode-shapes. These can be used to
estimate the vibration response of the real continuous structure. The discretised model will
have only a finite number of modes equal to the number of degrees of freedom considered in
the model. The lower-frequency modes will be most accurate compared to the continuous
structure, because the model (mesh) is less able to capture the real mode shapes as they
become more complicated.

Modal analyses can provide information about:

mode shapes,
natural frequencies, , and
mode participation factors, which are a measure of the amount of mass moving in
each direction for each mode.

These results are all useful in modelling a structures response to a broad variety of transient
loading scenarios. Modal analysis is the most fundamental dynamic analysis.

A modal analysis is useful when:

designing a structure to avoid resonant vibrations

designing a structure to make use of resonant vibrations
getting an indication of how a structure will respond to dynamic loads
calculating solution controls (time steps, etc.) for other dynamic analyses

A modal analysis should therefore ALWAYS be conducted before any other dynamic analysis.

A modal analysis requires some important assumptions. Modal analyses are always linear
(constant K and M), and independent of damping (no C matrix). Although damping can be
included if desired, by using a complex eigenvalue problem incorporating imaginary numbers.
Modal analyses are conducted without time varying applied loads, which makes it a free
vibration problem. Modal analyses can include constant applied loads known as pre-stress
loads. These loads provide stress stiffening to the structure and change the modal behaviour.

Because excitation and damping arent considered in a modal analysis, the magnitude of the
calculated displacements (or other outputs) of a particular component are meaningless. The
information is valid in the sense that it tells how much some areas will move relative to other
areas on the same part. However, the information does not reveal the actual displacement
values (this is what harmonic or transient analyses are for).

106
Lastly, the results from a modal analysis can be used to directly conduct other dynamic
simulations. This is achieved through the method of mode superposition. It is based on the
concept that any loading can be represented as a superposition of the individual mode shapes.
Mode superposition can be used for harmonic, response spectrum or transient analyses. For
harmonic and transient analyses, the mode superposition method may, in some cases, provide a
faster solution time.

Theory of Modal Analysis

Because free vibration and no damping are assumed, the governing equation of motion for a
modal problem is given by:

+ =

The modal model involves sinusoidal responses (with respect to time), so is assumed that
displacement and its derivates are:
= sin

= cos

= sin

Here is the mode shape (a vector representing the relative amplitude of the vibration of each
degree of freedom), and is the natural frequency of each node. By combining the above
equations, we get:
=

This is an Eigenvalue problem. The trivial solution is when is the zero vector and so there is no
vibration. Otherwise, non-trivial solutions can be obtained whenever

det = 0.

The roots of this equation are the eigenvalues, which are the natural frequencies of the
discretised system. A modal analysis can return as many natural frequencies as DOFs in the
model, though at best only the first half can be expected to be accurate outside of the
discretised model.

The solutions for the natural frequencies are substituted back into

= .

This problem is now unconstrained. Because of the zero-determinant conditions required to find
any non-trivial solutions, there are now an infinite number of solutions. Just like in all eigenvalue
problems, the eigenvectors are obtained by applying another constraint through normalizing .
The mode shapes can be normalised in several ways. ANSYS offers two methods:

Normalised by the mass matrix:

%
& =1

Normalised by unity, which means set the largest value to 1.

107
In most FEA packages, the eigenvalues and eigenvectors are found using a Lanczos Algorithm,
commonly the Block Lanczos method.

Mode participation factors are given by the following equation, where D is a unit displacement
spectrum in each direction and rotation.

= %
&(

If the mode shape has been mass-normalised according to the previous equation, then he
square of the mode participation factor gives the effective modal mass for that mode. The sum
of the effective modal masses of all the modes should equal the mass of the whole structure.

Pre-Stress Effects
The vibrational behaviour of a structure is changed when there are static loads in addition to the
vibrational loading. This is due to stress stiffening (which was discussed in Chapter 6). The static
pre-load changes the stress state in the material, which in turn changes the stiffness response. A
classic example of this is a guitar string. As the string tension is increased, the frequency of
vibration also increases (which translates as an increased pitch of the sound produced).

A pre-stress modal analysis works by first conducting a linear static analysis (with only the static
loads applied). From this the stresses in the model are obtained. These stresses are used to
calculate the stress stiffness matrix S. The problem is then dealt with by adding the stress
stiffness matrix with the structural stiffness matrix:

+) & = .

Harmonic Analysis
A harmonic analysis is more specific than a modal analysis. A harmonic analysis seeks the
steady-state solution to known harmonic loading. Harmonic analyses are useful for determining
the damped resonance frequencies of a structure i.e. so that they can be avoided during
operation. They are also useful for checking the displacements and stresses (or whatever is
desired) of a structure for a given sinusoidal loading case (e.g. an engine running at different
speeds). The outputs from a harmonic analysis are true values, and are not normalised as is the
case for a modal analysis. Damping can also be included into the analysis.

There are two approaches to a harmonic analysis. The first is to use complex phasor notation to
represent a vibration, and to solve the harmonic problem directly. The second method is to use
the results of a modal analysis via mode superposition.

In a harmonic analysis it is assumed that the stiffness K, damping C, and mass M matrices are all
constant (i.e. linear mechanical behaviour). All loads and displacement are assumed to vary at
the same known frequency, although the phases do not have to be the same.

108
Complex Phasor Notation
In this method for a harmonic analysis, a complex phasor is employed, because (unlike the use
of the sine function in the modal analysis) with complex analysis techniques the whole DE can
be factorised without neglecting damping and loading.

The complex phasor is given by:

* +
= cos, - + .sin, -

Both the applied loading and the vibration displacement response are assumed to have a
harmonic nature with frequency, . In phasor notation this gives:

= /01 cos,2 + - + sin,2 + - = /01 *

45
* 65
= 7 + * 65

= /01 cos,2 + - + sin,2 + - = /01 *

45
* 65
= 7 + * 65

displacement, and is 2 the phase lag of the driving force displacement output.
where Fmax and Umax are vectors of the magnitude of the input applied force and output

109
The time derivatives of the displacement are therefore:

= 8 7 + * 65

= 8 7 + * 65

Using the phasor representations, the overall dynamic equation becomes:

= + +

7 + =, + - 7 +

This is a simple linear system of equations except it incorporates complex numbers throughout.
The system can be solved in two ways. This first method is to use complex arithmetic and invert
the system matrix, just as would be done in a normal static analysis. In this case the system
equation can be thought of as being simply:

: = : :

The second approach is to use mode superposition, which is discussed in the next section. This
approach can be quicker and less computationally expensive than using the complex solver.

The results from a harmonic analysis provide the absolute values for the displacement and
derived parameters (e.g. stress and strain). The results will be complex if damping is included or
the applied loads are complex (i.e. different loads are applied with a phase lag between them).
A complex solution means that not all points on the structure are in phase with the applied load
or each other. Results can be viewed in different ways including real and imaginary parts or
amplitude and phase angle.

Mode Superposition
Mode superposition is based on the idea that any dynamic response can be represented as the
superposition of the individual mode shapes by using the appropriate scaling factors (see figure
below for an example). The mode superposition method is suitable for harmonic analysis,
general transient loading, and the response spectrum method.

110
The first step of the method is to change from nodal coordinates, the usual system for FEA, to
generalised coordinates. The relationship between the nodal and generalised coordinates is
defined as:

=; <
=7

where n is the number of modes to be included and Yi are the generalised coordinates (in this
case, modal coordinates). The more modes used, the more accurate the calculation will be.
Substitution into the equation of motion and pre-multiplication by 5 leads to:

5
> ; <+ 5
> ; <+ 5
> ; < = 5
>
=7 =7 =7

Using the orthogonal property of the eigenvectors and the previous normalisation of the mode
shapes (and a few other assumptions), the following equation can be derived:

<> + 2 > @> <> + > <> =A

where @> is the damping ratio for each mode (j=1 to n) (damping is discussed at the end of this
5
chapter) and fi is equal to > . This gives a system of n equations in the n unknown modal
coordinates (in the case of a transient analysis). This is far less equations then compared to the
number of DOFs, which is the number of equations for the full solution method. For a harmonic
analysis, the loads and modal coordinates are replaced with phasor notation, similar to as was
done previously. Since the system is a linear combination of modes, only linear behaviour is
allowed. A further disadvantage of the mode superposition method is that non-zero
displacement loads cannot be applied.

Response Spectrum Analysis

Buildings, ships, and other structures are often required to withstand transient shock and
impact loadings, such as seismic events or blasts. A response spectrum analysis is a fast
alternative to a full transient analysis. A response spectrum analysis provides an indication of
the maximum response of a structure to a transient loading. The output is a maximum because
the phase separation of the various frequency inputs is not taken into account.

The analysis must be preceded by a modal analysis because it uses mode superposition; the
maximum response is computed for the each mode due to the time history, and then the
individual responses are combined together to give the total response.

The input to the mode is a base excitation frequency spectrum, and may be expressed in terms
of displacement, velocity, or acceleration. ANSYS allows for two types of response spectrum
analysis: the single-point response (where a single spectrum excites all specified points in the
model) and the multi-point response (different spectrums excite different points in the model).

At the end of the modal analysis, the mode participation factors are calculated with ( as a
vector that represents the direction of excitation (this is different to previous definition):

= %
&(

111
Once the modal analysis has been conducted, the degree to which each mode is excited by the
specific spectrum is determined. The basic process is best outlined by the table below, which
progresses from left to right.

Mode Frequency Mode Participation Mode Spectrum Mode Response

Shape Factor Damping Value Coefficient

1 @7 )7 B7 C7
Ratio
7 7 7
2 @ ) B C
3 E E E @E )E BE CE

TOTAL C%

The values for ) are obtained from the response spectrum. Note that a response spectrum is
not the same as the input base excitation spectrum. Instead, ), - represents the degree to
which a single DOF oscillator with frequency would respond to the entire input spectrum. In a
response-spectrum analysis, a multiple-DOF system is considered to be a collection of single-
DOF oscillators, whose responses are super-imposed.

Mode coefficients , Ai, are then determined from the participation factors, natural frequencies
and spectrum values. The next step is to calculate the response amplitudes, R, for each mode
by scaling the mode shapes with the mode coefficients. The equations for both of these steps
are given below (depending on the type of spectrum input):

B =) , H=B , displacement spectrum input

)
B = , H= B , velocity spectrum input

)
B = , H= B , acceleration spectrum input
7

The final step in the process is to superimpose each scaled mode shape response to give an
overall response for the structure. The overall response represents the maximum that would

this process, Q > , though the will not be discussed here.

occur for the given input spectrum. There are several ways to calculate the weighting factors in

112
The basic summation equation is:

7/
S S

C = R; ; Q > H H> T
=7 >=7

where N is the number of modes included in the summation.

Random Vibration Analysis

An extension of the response spectrum method is a random vibration analysis. It is a spectrum
analysis technique based on probability and statistics. Random vibration analyses are meant for
loads such as acceleration loads in a rocket launch that produce different time histories during
every launch. Transient analysis is not an option since the time history is not deterministic
(sample is not repeatable). Instead, using statistics the sample time histories are converted to
Power Spectral Density function (PSD), a statistical representation of the load time history.

A Random Vibration analysis computes the probability distribution of different results, such as
displacement or stress, due to some random excitation. The analysis follows a modal analysis
and uses mode superposition. An internal combination is done to compute the combined effect
from each mode and their interactions.

Random vibration analyses are commonly used for airborne electronics, acoustic loading of
Airframe parts, jitter in alignment of optical equipment, etc.

Damping
Damping forces are present in any real-system. Damping results in energy-loss, which leads to a
decay in amplitude of the motion. Without damping, a vibrating structure would theoretically
oscillate forever. Physical damping behaviour is very complicated, though through appropriate
approximations, the mathematical implementation can be simplified.

113
There are three types of damping:

Viscous damping: Damping applied when a body moves through a fluid. E.g. shock
absorber.
Structural damping: This is also called hysteresis or solid damping and is due to the
internal material behaviour. E.g. Internal friction that turns kinetic energy to heat.
Material damping is not very well understood.
Coulomb damping: This is also called dry-friction damping, and is due to external
friction. It generally is modelled as a fixed force that is independent of velocity. E.g.
sliding friction.
For viscous damping, the damping force is proportional to velocity by the damping coefficient, V,
such that W = V X. If the damping is too large, the system will not oscillate. The critical
damping, V:Y , is the threshold between oscillatory and non-oscillatory behaviour, and the
damping ratio, @, is the ratio of the damping in a system to the critical value.
V
@=
V:Y

Note that this linear damping, which is characteristic of viscous or solid damping, causes a
difference between the natural frequencies of a system and the resonance frequencies.
Coulomb damping has no effect on frequency. The relationship between the damped
frequency, W , and the natural frequency, , for a 1 DOF system is given by:

W = Z1 @

In the frequency domain, viscous damping is linearly increasing. That is, the amount of damping
increases with frequency. On the other hand, structural or solid damping is generally considered
as being constant with respect to frequency. Coulomb damping is not dependent on the
displacement or velocity, but rather on the normal force and the friction coefficient.
Damping can be specified in many forms. The table below provides a summary of the conversion
between different measures of damping (note that U is the strain energy).

114
Rayleigh Damping
Mathematically, damping can be quite complicated to include into an analysis. There are many
damping models and approaches, some simpler than others. The most common approach in FEA
is to use Rayleigh damping coefficients, and . This approach says that the damping matrix is a
linear combination of the mass and stiffness matrices. The coefficient is the mass proportional
damping coefficient and is the stiffness proportional damping coefficient. This approach
decouples the problem, so that the linear form of the equation of motion can still be used. The
Rayleigh coefficients have no real direct physical meaning, they are simply the result of the
maths required to decouple damping in the problem. The coefficient acts to damp out low
frequencies while damps out higher frequencies.
The difficulty with Rayleigh damping is choosing appropriate values for and , such that the
damping characteristics of the real system are best modelled. Viscous damping is linear with
frequency, so it is well represented by . Structural or solid damping is ideally constant with
frequency, so it usually requires a combination of both terms.
ANSYS (and similarly in most FEA packages) allows damping to be applied to a model in multiple
ways, depending on the type of analysis.
The damping matrix for a transient analysis (and damped modal analysis) is a linear combination
of the mass and stiffness matrices through the Rayleigh damping coefficients:

=[ +\

Note that in ANSYS there are additional terms in the above equation, which are for damping
applied directly through the elements, and for and damping applied through the material
properties. These extra terms simply add to the total damping. If both methods are defined,
then the total damping will be additive. The moral here, is to be aware of where and how you
are adding damping terms into the system.

damping ratio, @, with frequency. The justification for this is that the response frequency is not
In a transient analysis (and damped modal) it is not possible to directly define a constant

known in advance, and thus the constant damping cannot be defined (as it requires prior
definition of the response frequency in the damping matrix). In this case, the Rayleigh damping
coefficients need to be calculated to best fit the desired damping over the expected response
frequency range.

In the absence of damping, can be related to the damping ratio through:

[ = 2@

and in the absence of damping, can be related to the damping ratio through:

2@
\=

In many practical problems, damping (or mass damping) may be ignored ( = 0). To represent
a constant damping ratio with frequency, a combination of Rayleigh damping coefficients can be
determined by:

115
[ \
+ =@
2 2

range of frequencies. This provides simultaneous equations from a given value of @:

One approach is to assume that the sum of and damping is nearly constant over a small

2@
[ = 2@ and \ =
7
7 + 7+

The damping matrix for a harmonic analysis is given by:

2
=[ + ]\ + @^

equal to the driving frequency, . As before, there are additional terms in the damping matrix
where the constant damping ratio is now included as the response frequency is known and is

for material and element defined damping that are not shown here.
Coulomb damping is included by defining contact surface elements with friction capability.
Coulomb damping is non-linear and thus requires a non-linear analysis.
In mode superposition analyses, the damping matrix is not explicitly computed, but rather the
damping is defined directly in terms of a damping ratio.

Transient Structural Analysis

The dynamic behaviour of a system can be difficult to model. In certain cases of loading and
response, simplifications can be made (as outlined in the previous sections). However, a fully
transient simulation is often required. All terms in the dynamic equation of motion are included:

& +_ + =

Dynamic simulations require discretisation in the time-dimension. The duration of the

simulation is divided into a number of discrete time-steps. The simulations are then gradually
solved from the beginning of the simulation period to the end. This process is called time-
integration.

116
There are a number of approaches to time-integration. A couple of these, called the Euler and
Rung-Kutta schemes were presented in the background of numerical methods, at the end of
Chapter 3. Other higher-order and intermediate methods also exist.

The different time-integration schemes are grouped into two categories: implicit and explicit.
For the explicit approach, only the known derivatives from the current time-step are used to
solve for the next time-step. For the implicit approach, the derivatives from the next time-step
are used. The implicit method is generally more accurate and allows for a longer time-step, but
each step is more computationally expensive, because it requires matrix inversion. The explicit
method requires the use of a small time-step, but each time-step can be solved more quickly.

The implementation of implicit time integration is detailed in this chapter. Explicit dynamics is
presented later in Chapter 9.

Implicit Time Integration - Newmark Scheme

The most common implicit method of solving this is the Newmark scheme. In the Newmark
scheme, the derivatives used are determined by linear interpolation:

`7 = + a,1 b- +b `7 c

`7 = + + a,0.5 [- +[ `7 c

This is unconditionally stable if b 0.5, [ 0.25,0.5 + b- and 0.5 + b + [ > 0. When this is
applied to the structural differential equation, to enforce equilibrium at `7 , the following
equations result:

1 b
h &+ _+ i
[ [ `7

1 1 1 2[
= +&h + + i
[ [ 2[
b b[ b 2[
+_h + + i
[ [ 2 [

The advantages of implicit time integration are that it is unconditionally stable, and that the
time step will vary only to satisfy accuracy requirements and non-linear iteration. To determine
an appropriate time-step it is often worth conducting a modal analysis first. The time-step must
be small enough to capture the time-varying loads and the time-varying response, but large
enough to reduce solution time. The time-step will affect the convergence behaviour of non-
linear solutions.

A good initial guess for a suitable time step size is:

1
=
20Ajkl

where fmax is the maximum frequency of interest in the structural response.

117
The disadvantages of implicit time integration are:

that non-linear solutions are obtained using a series of linear approximations (Newton-
Raphson method), and each time step may have many equilibrium iterations
that the solution requires inversion of the equivalent stiffness matrix at each time step
that small, iterative time steps may be required to achieve convergence of non-linear
problems

Transient analyses offer a wide range of options for modelling real world systems. Both non-
linear and time-dependent behaviour can be included. Because of this increased level of
complexity, the accuracy and interpretation of the results becomes exponentially more difficult.
In ANSYS, convergence tools are provided, but convergence is not guaranteed for highly non-
linear problems.

Final Words
When conducting dynamic analyses (particularly non-linear transient ones), be aware of the
limitations of your conceptualisation of the real world problem, and how it is implemented in
the FE software. Verification and validation are extremely important. Transient analyses can be
very sensitive to the initial and boundary conditions as well as to the model geometry and
material constitutive behaviour. The more there is to model (and variability in the model), the
more there is to get wrong.

Summary of Structural Analyses

If inertial and damping effects are small, use a linear or non-linear static analysis
If the loading is purely sinusiodal and the system is linear, use a harmonic analysis
If only the maximum response is required for an input spectrum, use a response
spectrum analysis
If the bodies can be assumed to be rigid then the solution time can be significantly
reduced
In all other cases, a full transient analysis can be used as it solves the full system
equation with no simplifications

118
119
120
Chapter 8
Buckling

121
Load-carrying structures may fail in a variety of ways:

Material failure plastic yielding, fracture, fatigue, creep, corrosion

Structural failure excessive deformation, loss of stiffness
Combination of the above

Buckling is an elastic instability that accompanies compressive stresses, and is a type of

structural failure. Failure by buckling can occur at loads far less than those expected for material
failure (e.g. compressive yielding). At the onset of buckling the structural member will
experience a very large change in displacement with very little change in load. Buckling can lead
to collapse of a structure and is usually accompanied by plasticity.

Buckling is a geometric non-linearity. In a similar manner to dynamic analyses, there are several
types of buckling analysis. A linear buckling or eigenvalue analysis is analogous to a modal
analysis, except that the buckling modes are obtained. For the true buckling response, a non-
linear static analysis is required.

Buckling occurs by compressive and/or shear strains and is therefore particularly important for:

Columns and beams

Compression members
Thin-walled vessels
Stiffened skin structures

Buckling Theory
Buckling is an elastic instability. Stability is a measure of how well a system can handle
disturbances. When a structure is stable, it can handle disturbance deflections. However, if the
load on a structure exceeds the critical buckling load, then it becomes unstable so that a tiny
disturbance deflection that slightly changes the load position, will radically upset the
equilibrium.

Around the critical buckling load, the structure undergoes a loss of stiffness. That is, the
structure can undergo large deflections corresponding to small changes in load. This translates
as a horizontal or flat region in the load-displacement curve.

122
Even though it looks the same here, buckling is not the same as yielding. It is possible for a
structure to buckle elastically and return to its former shape after the load is released. Hence
the phrase, elastic instability.

Fixed-end Column under Axial Compression

The typical buckling scenario that can be used to demonstrate this is the fixed-end column
under axial compressive loading. When the structure is in compression, the compressive loads
can cause bending moments if some deflection occurs perpendicular to the load.

When the load becomes too large, the bending moment due to the compressive load will cause
the component to bend and deflect more, which will increase the eccentricity of the load (the
lever arm) and will make the bending moment even larger, which will cause the component to
bend even more, and so on.

For small loads, the material may be able to bend a finite amount and bear the load. However, if
the bending moment exceeds the bending strength of the material, then the component will
yield and potentially collapse.

Buckling can be understood best through the load-deflection curve (in this case, axial load vs
transverse deflection). The application of an axial compressive load will cause only axial
compression until the bifurcation point is reached. A bifurcation point is when there is more

123
than one possible outcome. At the bifurcation point in the fixed-end column, with continued
loading the structure could keep deflecting purely axially, deflect to the right, or deflect to the
left (for a 2D problem). For a load below the bifurcation point, the structure will be in stable
equilibrium. At the bifurcation point, the structure is in neutral equilibrium. Above the
bifurcation point, the structure will be in unstable equilibrium.

How the column buckles will depend on how it is constrained. You may recall from the Stress
Analysis and Design course that columns with different end constraints have quite different
buckling loads and shapes (as shown below).

A structure will also have more than one buckling load (or bifurcation point). This leads to
multiple critical buckling loads each with a corresponding mode shape.

124
In the idealised perfect structure the column will indefinitely compress only in the axial
direction. So what makes a real structure buckle?

A structure or component will buckle due to real-world imperfections:

Geometric Imperfections
Material defects
Eccentric axial loading
Perturbation Loads: small transverse loads can cause load eccentricity below the critical
load

These imperfections will decide the direction that a component buckles and will lower the
actual buckling load below the theoretical bifurcation point. The true buckling problem is non-
linear.

Modelling Buckling
There are several methods used to model buckling, that are suitable for FEA:

Linear (eigenvalue) Buckling

Non-linear Buckling

The linear buckling analysis is inexpensive, and gives an indication of the buckling loads and
mode shapes, though it provides overestimates. The mode shapes from a linear buckling
analysis can be included as a geometric imperfection in a non-linear buckling analysis, which will
provide more accurate results.

For each method, it is important to get the constraints and the loads correct (especially
important is the relationship between the loading and the deflection of the structure).

Different types of structure will buckle in different ways. For example, some structures will
collapse and have a negative stiffness after the onset of buckling. Others may tend towards a
peak value or even start to increase in stiffness. Several examples are shown in the figure below.
The different types of behaviour are important when deciding on what type of buckling analysis
will be appropriate.

125
Linear Eigenvalue Buckling
A linear or eigenvalue buckling analysis predicts the bifurcation points and the corresponding
mode shapes. The analysis obtains the same solutions as the analytical Euler method. It
considers an idealised model where there are two main assumptions:

1. The stiffness is based on the original geometry configuration

2. Pre-buckling behaviour is a linear function of any initial applied load; such that if the
applied load is doubled, then the stress-stiffness will also double

The derivation starts with the incremental form of the equation of motion as derived in the
Chapter for non-linear problems. The tangent stiffness matrix, KT, is given by the sum of the
regular stiffness matrix, K, and the stress-stiffness matrix, S.

A linear buckling analysis has two steps. Firstly, a pre-load is applied (usually a one-unit load)
and the stress stiffness is determined based on that. The magnitude of the load is arbitrary,
although careful choice can make interpretation of the results easier. The stress stiffness matrix
is then calculated using the results:
=

The stress stiffness matrix is then substituted into the incremental system equation. Recall that
it has been assumed that the stress stiffness matrix at an applied load F = F0, is a linear function
of the stress stiffness for the initial load, = . This gives:

+ =

126
At the onset of buckling, there will be a large change in deflection for a small change in load.
This can be written as F 0, which gives the system equation to be:

+ =

A zero change in load can always cause a zero change in deflection, so the trivial solution is
= . The non-trivial solution will find the bifurcation points. This is an eigenvalue problem
just like for a modal analysis. The roots of the characteristic equation will give the eigenvalues:

| + |=0

The eigenvalues, , are scale factors for the critical buckling loads of each mode:

Note here that the critical buckling load is a vector NOT a single value. This means that if more
than one load is applied to the model (e.g. perturbation load), all loads will be scaled and
included in the analysis. The eigenvectors, , are the corresponding buckling mode shapes.
Normalisation of the eigenvectors is achieved by making the largest value equal to 1.

The results displacements and stresses, etc are all relative values. The magnitudes have no
meaning, though the relative values provide an indication of the nature of the distributions.

The assumptions used in the linear buckling analysis are idealised. The critical loads that are
calculated will be over-estimates. In reality there will be imperfections in the structure or
perturbation loads that lead to early onset of instability. Furthermore, the stiffness matrix, , is
not constant, but changes due to the change in geometry. Similarly, the stress-stiffness matrix is
not only a function of load, but of the deflected geometry. A more accurate prediction of
buckling loads requires a non-linear analysis.

Linear buckling analyses do have some advantages:

Relatively inexpensive (fast) analysis

Gives an indication of the expected buckling loads and mode shapes
Provides information for a non-linear analysis
o Setting the mesh density
o Mode shapes can be used as a geometric imperfections

Non-linear Buckling Analysis

Real buckling is a non-linear problem. There are two main stages of buckling behaviour:
Pre-buckling
Post-buckling

There are several approaches for solving a non-linear buckling problem. These include:
Static analysis using Load Control
Static analysis using Displacement Control
Dynamic solution (load or displacement control)

127
Additional instability analysis tools that will be discussed are the Arc-Length Method (alternative
to Newton-Raphson) and Stabilisation (similar to weak-springs).

It is not necessary to run a linear buckling analysis before a non-linear one; however it is
recommended. The linear analysis can provide an indication of expected mode shapes, expected
critical load values, and a verification check of the non-linear solution. A non-linear buckling
analysis involves incrementally increasing the load/displacement and monitoring the load-
deflection curve for the point of instability (and beyond if desired).

Unlike linear buckling, a non-linear buckling analysis can include

Initial geometric imperfections

Large-deformation
Contact
Plastic yielding behaviour
Other nonlinear behaviour

Non-linear buckling is therefore more accurate than eigenvalue buckling.

Depending on the problem not all non-linear methods will be sufficient to capture post-buckling
behaviour. This will be looked at in the following sections.

Non-linear analysis requires initiation of buckling. Because of this the non-linear analysis usually
requires a slightly different model/BCs to the linear analysis. Without initiation of buckling, the
load-displacement can follow the unstable equilibrium path. In some cases, the original applied
loading is such that it will initiate the buckling, but this must be considered carefully. The
methods to initiate buckling are:

Geometric imperfections
o Defects
o Extra curvature
o Superimposed mode shape
Perturbation loads
o Out-of-plane forces
o Eccentric axial loads
o Moments

Buckling initiation will remove the sharp discontinuity in the load-displacement curve. However,
the magnitude of the imperfection or perturbation loads will influence the results. The value
should be small relative to the overall dimensions and actual loads of the structure and should
best match the real structures expected behaviour. Manufacturing tolerances can be used to
estimate the magnitude of the imperfection, and the mode shapes from a linear analysis can be
used as the form of the geometric imperfection.

Load Control Method

For the load control method, the solution is obtained by incrementally applying forces. The load-
displacement curve is then obtained from which the critical buckling load can be extracted. A
disadvantage of this approach is that it cannot always provide information beyond the critical

128
buckling load, because the Newton-Raphson method cannot handle zero-stiffness (as the KT
matrix becomes singular).

An example of this is snap-through buckling of a simply supported shallow arch. At the critical
buckling load, the load-displacement curve has a local maximum, the tangent stiffness (slope) is
zero, and the Newton-Raphson method collapses.

The failure of the solution procedure can be used as an advantage for determining the buckling
load. That is, the load at the point when the solver stops is the buckling load. Automatic
bisection features of most FE solver software can be used to help narrow down the buckling
load to a reasonable level of accuracy. Auto-bisection works by the solver shrinking the load
increment step if it fails to converge. This means it will automatically hone in on the buckling
load. The chosen minimum step size will therefore affect the resolution of the buckling load.

Lack of convergence does not necessarily indicate that the critical buckling load has been
reached. It is important to always check the load-displacement curve, and that the tangent
stiffness is near zero.

129
Displacement Control Method
The first alternative option to use the displacement control method. With displacement control,
a displacement is imposed on the structure, and the load is the reaction load at the imposed
displacement. This method remains stable beyond the critical buckling load (usually), and
provides full pre- and post-buckling behaviour.

A downside to the displacement control approach is that you need to infer the load control
response from the displacement-control solution. This will only work when you know what
displacements to impose!

An additional drawback of the displacement (and load) control method is that it does not allow
for dynamic behaviour. The figure shown above is not quite correct. This is because despite the
centre point being constrained, other parts of the structure will still snap-through. The snap-
through behaviour is dynamic and will induce inertial forces. This causes problems for the static
solver. There are several ways around the problem. One approach is to use a transient analysis,
the other is to use special features built into FE software to cope with these types of analyses.

Transient Analysis Method

The final approach is the transient analysis method. A period of time is simulated with either
load control or displacement control. From this the load-displacement curve for the dynamic
response can be obtained.

Time-stepping is very important for a transient buckling analysis. The minimum time-step
influences precision, and large time-steps can allow the structure to jump the instability
region. Automatic time-stepping can be used to locate the critical buckling load.

130
A disadvantage of a transient simulation is that the softening response will not be calculated if
the structure snaps-through (It will be a dynamic snap). However, that said this dynamic
behaviour is what real structure would follow. Care must be taken to ensure the chosen analysis
method gives the information required for the analysis.

Dynamic simulations also require the inclusion of dynamic properties such as density and
damping. It can often be difficult to damp out unwanted dynamic effects such as ringing.

Instability Analysis Tools

The arc-length method is a modification of the Newton-Rhapson method, which is designed to
handle instabilities and negative tangent stiffness, so it can evaluate the load-displacement
curve into the post-buckling region. Instead of finding the intersection of the tangent to a
horizontal line at the applied load, it finds the intersection of the tangent to a circular path
around the initial point on the load-displacement curve.

As a result, the arc-length method has a variable applied load per sub-step. An additional
unknown is introduced into the solution, ( 1 ! ! 1). This provides the equation:
" %&
= #$$

The load factor is given by the basic formula for a sphere of radius '. However, the DOF
dimensions of the sphere are scaled by (, and the load dimensions are effectively scaled by
#$$ :

' = )( * " + * = )( * ,* + * 1D

131
The arc-length method can be used for static problems with ramped loading/displacement. It
cannot be used with rate-dependent materials (e.g. visco-elasticity, creep), and is best suited to
solve problems where the response is smooth without sudden bifurcation points

Stabilisation is another method used for non-linear buckling. Artificial dampers are added to the
outside of the structure. This is good for problems with snap-through and collapse behaviour.
However, this approach can over-dampen the response and increase the buckling load. It is
best to only apply stabilisation for post-buckling behaviour. In ANSYS, a scale factor for the
applied damping coefficient can be set by the user. The maximum allowable stabilisation force
can also be set.

132
Chapter 9
Explicit Dynamics

133
The dynamic equation of motion can be solved using implicit or explicit time integration
schemes. Chapter 7 outline the implicit approach. Explicit dynamics offers an alternative for
dynamic problems.

Problems with high strain-rates require small time-steps because the system is changing rapidly.
Standard implicit methods, despite being more accurate, are computationally expensive for
each time step, so implicit methods are used instead. In addition, non-linear implicit transient
analyses use the Newton-Raphson method, which iterates until convergence at each time step.
Non-linear problems may take a very long time to solve using standard implicit methods. In
some cases, they may not solve at all due to convergence issues.

Explicit time integration is more efficient for short duration and highly non-linear simulations,
for example impact, penetration and shock wave propagation problems. It is also easier to
incorporate non-linearities with the explicit method, such as problems involving contact,
buckling, plasticity, and material failure.

Time-Integration Scheme
There are several time-integration schemes that can be used in explicit dynamics. The central-
differencing scheme is generally used because it is more accurate than forward or backward
differencing, but is still simple to implement. In the central differencing scheme, the first
derivative is evaluated at the half time-steps, so that:

/ /

Because this is an explicit method, when it is used to solve the equation of motion the stiffness
matrix goes to the right hand side and doesnt require inverting.

134
Explicit Dynamics Analysis
Physical Formulation
The formulation of the stiffness matrix was discussed in Chapter 3, and could be based on
direct formulations or on minimising total potential energy.

Explicit methods generally use a different formulation to the structural models that have been
considered in earlier chapters. Explicit dynamics is formulated by using the momentum and
mass equations shown below. With this method, and because small time-steps are used, it is
easier to incorporate non-linear behaviour into an explicit dynamics model than an implicit one.

An accurate solution should also satisfy conservation of energy. For an explicit dynamics
simulation, conservation of energy is not enforced. Instead the energy is monitored to provide a
measure of the solution accuracy.

1
2 2 2 "

Solution Procedure
In the explicit dynamics scheme, the accelerations are solved for first, and then the velocity and
displacements are derived using the central-difference scheme. At the end of each sequence of
integrations, the structural properties are derived from the change in geometry, so that the
change in internal forces is determined before the accelerations are calculated for the following
time-step. This is shown in the following diagram.

135
In explicit dynamics, the element equivalent has nodes and a centre point, which is required to
trace and apply body accelerations and body forces. Positions, velocities and accelerations are
all defined at the nodal points. Material quantities (e.g. mass, temperature, stress) are defined
at the element centres. This makes explicit dynamics analysis not very good for investigating
stress distributions, for example. However, explicit dynamics is aimed at problems with large
deformations and considerable plastic flow. As such, the displacements are of primary interest.

A detailed summary of the explicit dynamics mathematical procedure for 2D plane strain
elements is included at the end of this chapter.

Time-stepping
The time-step in an explicit dynamics model, must be bound in order to keep the central-
differencing scheme accurate. A limitation is applied to the time-step size so that a stress-wave
shouldnt be able to travel further than the smallest dimension of an element in a single time-
step. This is represented by the Courant Condition:

$& )
( *+
where is the characteristic length of an element, ( is the local stress wave-speed, $ is a
stability factor (usually 0.9), and is the time-step. The characteristic length of an element is
calculated according to the following table:

Hexahedral/Pentahedral The volume of the element divided by the square

of the longest diagonal and scaled by 2/3
Tetrahedral The minimum distance of any element to its
opposing element face
Quad Shell The square root of the shell area
Tri Shell The minimum distance of any elements node to
its opposing element edge
Beam The length of the element

Most explicit dynamics software automatically sets the time-step in order to maintain the
Courant condition. Unfortunately that can mean that one or two very small elements will reduce
the time-step for the whole domain. It may be possible to ensure that this isnt the case for the
original mesh, but it may still occur over time as the mesh deforms.

Mass scaling can be used to prevent small elements from changing the time-step. As the time-
step is inversely proportional to the speed of sound, it is related to the mass of each element:

1 1
1
( 0++ 0++
/

If the mass of very small elements is artificially increased, then it can change the stress wave-
speed in that cell, and will stop that element from unnecessarily limiting the time-step. This is
appropriate when there are only a few small elements, but it is important to realise that it is
artificial mass scaling; it is an inaccuracy and changes the inertia of the structure. Mass scaling is

136
only permissible because it applies to small elements that have very small inertia anyway, it is
not always a good solution to time-stepping problems.

Corrective forces
In the explicit dynamics formulation described above, the strain rates in an element are
expressed in terms of diagonally opposite nodes. For example, the x-direction strain is (see
slides at end of chapter for derivation):

1
45 7 85 98 5 9 85 7 8:
3 3

Certain combinations of values can lead to the resultant strain be calculated as zero even
though the nodal strain and displacement values are non-zero. These particular deformation
shapes are called hourglass deformations (shown in the image below). If no strain is calculated
then no stress will be either, and thus no reaction forces to oppose the deformation and
prevent it from increasing. Hourglass deformation can occur for 2D and 3D shapes. In ANSYS
Explicit Dynamics (and most other packages), corrective nodal forces called hourglass damping
are introduced to the elements to try and prevent this type of deformation happening.

A second type of corrective force that can be applied in explicit dynamics is anti-tangle forces.
These forces help to prevent nodes in an element passing back through themselves. In an
explicit dynamics simulation, the mesh moves with the object, which then often undergo large
and complicated plastic deformations. In such simulations, there is a risk of mesh
entanglement, which is when re-entrant corners in a deformed element cross back over
themselves. This can be dealt with by the application of nodal anti-tangle forces.

137
Static Modelling using Explicit Dynamics
Explicit dynamics is aimed at solving highly non-linear dynamic problems. However, it can also
be used static equilibrium problems. This is achieved through applying damping on the velocity
calculations, such that:

; 51 2=>? 8; 51 =>? 8;

where Rd is a damping factor.

Elements and Meshing

Meshing considerations for explicit dynamics are generally similar to other FEA. There are
however a couple of key differences. The first major difference is the element formulation,
which was seen in the previous sections of this chapter. The other major difference is to do with
generating the mesh and its refinement.

Explicit dynamics has a similar range of element configurations compared to implicit methods,
including solid elements, shell elements, and beam elements. Beam and shell elements behave
in much the same way with the cross-sectional properties being entered as parameters. These
parameters will not be used to control the time-step size, as this would create unacceptably
small time-steps.

When meshing an explicit dynamics model the element sizing must be carefully set by the user.
The requirements for the mesh refinement however, are much different to that of an implicit
analysis. In an implicit analysis, the regions of stress concentration are known in advance and
generally stay the same throughout the whole solve time. In an explicit analysis the regions of
high stress continually change due to the large amounts of elastic and plastic deformation. It is
therefore not recommended to locally refine the mesh at initial stress concentration points. This
will simply increase the solve time due to the small step-size forced by the small elements. Fine

138
geometric details should also be avoided in explicit dynamics models as they require small
element sizes (which again, will decrease the step-time size). Because of this, explicit dynamics
meshers tend to behave differently to implicit ones. They are set to automatically defeature the
geometry prioritising larger element size over correctly capturing the geometric detail. This has
little effect on the accuracy of the results, as explicit dynamics modelling is all about simulating
the overall deformation behaviour (i.e. it is not typically used to obtain stress values). The level
of error in the results due to defeaturing the geometry is usually very small compared to the
overall deformation in the model. Furthermore, the aim of explicit dynamics is to obtain
deformation results and not stress results (which are single valued for the whole element).

Material Models
In general, materials have a complex response to dynamic loading. The following phenomena
may need to be modelled

Strain hardening (rate independent/dependent)

Orthotropic behaviour (e.g. composites)
Compaction (porous materials)
Crushing damage (e.g. ceramics, glass, rock, concrete)
Failure
Thermal softening
Chemical energy deposition (e.g. explosives)

No single material model incorporates all of these effects.

The following are commonly used for structural analyses, and are worth understanding:

Elasticity is used to describe purely linear behaviour below the yield stress. The standard
isotropic elastic model has been used throughout these notes. Another elastic model is
orthotropic. This means that a different set of elastic variables are associated with three
perpendicular directions. This is convenient for some composites, and possibly some other
solids, such as metals with crystals elongated in one direction.

Plasticity is the behaviour of a material past its yield point. Below the yield point, there are no
material non-linearities, but above it, a material will require much less additional stress to
deform. Elastic deformation is recoverable, i.e. the deformed solid will spring back. Plastic
deformation is permanent. The Von-Mises yield criterion is used for all plasticity models.

5 8 5 78 5 7 8 2A

which can be written in terms of purely deviatoric stresses:

2A
B B B7
3

This criterion is like a multi-dimensional surface. Isotropic hardening assumes that this surface
dilates as hardening ensues, whereas kinematic hardening assumes that this surface translates
as hardening ensues.

139
Various models exist for the yield strength function. The first models are based on constructing
a non-linear stress-strain curve.

The bi-linear (isotropic or kinematic) model assumes two linear functions can be used,
the elastic linear function which passes through the origin, and the plastic linear
function which does not.
A A 3 C

The multi-linear (isotropic or kinematic) model breaks the stress-strain curve into a
greater number of linear segments. As the loading is applied, elements can be promoted
to the next linear function when they exceed the bounding deformation.

The remaining models allow for temperature and strain-rate dependence

Johnson-Cook strength

A D3 E C FD1 log C F41 KL* :

Cowper-Symonds strength

R9
CO
A D3 E C F M1 N Q S
P

and others including Steinberg Guinan Strength, and Zerilli Armstrong Strength

Plasticity is not the only form of material strength behaviour. There are also models for brittle
materials (e.g. ceramics), granular materials (e.g. soil), and porous materials (e.g. foam).

Visco-elastic behaviour allows for creep and relaxation. Relaxation is when the load in a strained
solid gradually approaches zero. Creep is when the strain of a loaded solid gradually increases.

Hyper-elastic behaviour occurs when a solid undergoes very large deformation, which is all
elastic. Different models are appropriate for different strain ranges:

Model Strain Range

Neo-Hookean Up to 30%
Mooney-Rivlin 30% to 200%
Polynomial
Yeoh Up to 300%
Ogden Up to 700%

Failure Modelling
Modelling material failure can also be included into the analysis. The failure has two
components:
Material failure is the first phase, when specified criteria are met within a material,
which initiates post failure response.
Post failure response occurs. The mechanical behaviour change is dependent on the
failure model. In instantaneous failure, the element can only transfer compressive
hydrostatic pressure. In gradual failure (damage), stresses are limited by a damage
evolution law.

140
Various models include:
Plastic strain failure, failure occurs when strain exceeds a critical value
Principal stress/strain failure accommodates for either a max tensile stress/principal
strain, or max shear stress/shear strain.
Stochastic failure introduces distribution to the failure values.
Crack softening failure is a damage model often used for concrete. Note it doesnt
actually simulate a growing crack.
Erosion is a numerical mechanism for the automatic removal of elements. There are three
options to initiate erosion: on geometric strain (eff > critical value), on material failure and on
minimum time step (t < tcrit ). This method retains element inertia by replacing elements with
mobile point masses. It can be used to keep element mass in the model once eroded.

Body Interactions
Body interactions are used to define relationships between different bodies in a model. There
are three main ways of treating body interaction (i.e. contact) in explicit dynamics:
Trajectory based contact
Proximity based contact
Bonded bodies and rigid links

Trajectory Based Contact

For trajectory based contact, the trajectory of nodes and faces are tracked as they deform. If the
trajectories intersect, then contact is detected. Forces are added to repel the surfaces back
away from each other. A problem with this technique is that the time step is not adjusted as the
contact occurs. This means that significant penetration can occur over a given time step.

141
Proximity Based Contact
For the proximity based contact method, a detection zone is added around each body (or each
body for which contact is activated). If nodes enter the zone, a penalty force is added. This
approach is similar to the penalty method for implicit analysis. The time step then reduces so
that no penetration will occur.

Bonded Bodies and Rigid Links

Fixed contact between bodies can be modelled in two ways. Bonded bodies are where External
nodes/faces of bodies are tied together at initialisation. The connections can be breakable. The
alternative is to use rigid links. This approach allows remote points of separate bodies to be
joined together as desired. Additional rigid links are automatically created by the software to
ensure all DOF are transferred correctly. These connections can also be breakable.

Final Words
Explicit dynamics can be used to produce very spectacular and entertaining simulations. With an
understanding in the theory and implementation of explicit dynamics, it is also possible to
produce meaningful results for highly non-linear problems. As always, verify and validate your
models. There are so many modelling settings and adjustments within FE software, especially for
explicit dynamics.

Lastly, the possibilities from non-linear transient analyses are endless. Have some fun with it!

142
143
144
145
146
147
148
149
150
151
152
Chapter 10
Acoustics

153
Acoustics is the study of sound, including its origin and propagation. Sound travels as
mechanical waves through structures and fluids. Acoustics is therefore directly related to the
field of vibrations.

Acoustic analyses are preformed to gain an understanding of how waves are going to travel
through a system. In particular, acoustics studies are usually interested in how the pressure
waves travel through a fluid. The field of acoustics is not just about the sound we can hear; it is
also often concerned with infra-sound and ultra-sound.

The two difficulties required for many acoustic analyses are

the need to model fluids (for an acoustic analysis, the fluid can often be modelled
without using CFD because several assumptions are made), and
the need to couple the fluid to a structure (so that waves can move from a solid into a
liquid and vice versa).

This chapter will cover these two topics.

Theory of Acoustics in a Static Fluid

The motion of a fluid is governed by several physical principles, which can be described in three
equations: the mass continuity equation, the momentum equation (forces in the fluid), and the
energy equation (heat transfer and the relationship between temperature and density). These
form the fundamentals of CFD analyses.

For the acoustic problems considered in this course, several assumptions are made to simplify
things:

The fluid is compressible, and density changes with variations in pressure.

The fluid is inviscid, so there is no energy dissipation due to viscosity.
There is no mean flow of the fluid.
The mean density and pressure are uniform throughout the fluid.
Thermal effects can be neglected.

Implementing these simplified assumptions does not require CFD. Following these assumptions
leads to the acoustic wave equation, also referred to as the lossless wave equation, for sound in
fluids:

1
= , =

Where is pressure, is time, is the speed of sound in the fluid, is the bulk modulus of the
fluid, and is the mean fluid density.

154
FEM of Acoustics
In an acoustics problem, the structural and fluids physics need to be considered simultaneously
to capture the fluid-structure interaction. Direct coupling (not sequential coupling) must be
used.

At the fluid-structure interface, the acoustic pressure exerts a force applied to the structure and
the structural motions produce an effective fluid load.

Two approaches are taken to discretise the static fluid:

Pressure formulated elements

Displacement formulated elements

Pressure Formulated Elements

The pressure formulated elements have one degree of freedom at each node, which is the
scalar pressure field.

The structural system has the equation:

+ + = +

The fluid system has the equation:

+ + = +

The coupling loads are only applied to the elements at the interface, and are given by the
following formula. There is an extra displacement degree of freedom at the interface.

= , =

The combined system can be represented by the following matrix equation:

0 0
+ 0 + 0 !=

The fluid loads are defined using a flow load, according to the following equation. Flow loads
defined internally instead of at the structure interface represent point sound sources.

Flow Load = Particle Acceleration Effective Surface Area Density

Constraints are applied through the pressure DOF. Zero pressure represents a free surface,
otherwise a boundary in the fluid domain represents an ideal rigid wall.

Pressure-formulated elements are advantageous because they require less DOFs and because
linear shape functions are used for the pressure fields. There is direct coupling with the
structure and direct application and output of pressure. The disadvantages of this technique are
that the system matrices are un-symmetric (hence poorly conditioned), that there are no
displacement degrees of freedom for elements away from the interface, and that energy losses
due to viscosity are difficult to incorporate.

155
Displacement Formulated Elements
For displacement formulated elements, the same element types are used for the fluid and the
structure, but the material model is changed to reflect the fluid behavior. The normal stiffness is
determined by the bulk modulus of the fluid, and the shear stiffness terms are set to near-zero.

Conceptually, the fluid-structure boundary is easier to understand than the pressure-formulated

elements. However, there is still a complication. The fluid must be allowed to move tangentially
to the structure surface, so there are no shared nodes between the solid structural domain and
the fluid domain. Along the boundary, coincident nodes are created, that are constrained using
multi-point constraint equations.

Displacement formulated elements are being phased out of ANSYS. This kind of analysis would
normally be conducted using CFD software.

156