Finite Element Method (FEM) Simulation
CS 482
Lecture, Dr. Lawlor
Our goal is to be able to simulate solid
objects, like wood, metal, and bone, as they deform, bend and
break. This is actually a problem from materials science,
which uses a thick pile of specialized terminology:
- Strain measures how much
the material has been deformed, usually measured from some
"rest" configuration. The symbol is usually e or
epsilon. Units are usually in inches of displacement per
inch of material. For Hooke's law on a one-inch spring,
strain is the end displacement x. Strain can be simulated
(for example, integrated from particle forces), or measured with
either accurate calipers or a strain gauge.
- Stress is the force the
material exerts to fight back against strain. The symbol
is usually the Greek letter sigma σ. The units for stress are in
pounds of force per square inch of material. In 1D
material, stress is a scalar. In 3D material, stress can be
different in different directions, but it can always be
represented by a 3x3 matrix, the "Cauchy stress tensor", which
maps each direction to the force the material exerts along that
direction.
Stress and strain can sometimes be
related via a stress-strain curve, with the strain on the X axis, and
stress on the Y axis. This is a somewhat misleading curve,
because stress also depends on the material's internal state, so
the stress-strain curve can become multi-valued. In general,
stress-strain behavior is usually modeled via a Constitutive Equation, such as:
- For a purely elastic material, the
stress-strain curve is linear, like a spring: push on it, and it
deforms (or deform it, and it pushes); let it go, and it returns
to its original shape and stops pushing. This linearity is
captured in Hooke's Law:
F=-kx in physics terminology; or σ=Ee (stress is stiffness time
strain) for engineers. The E there isYoung's
modulus, which measures the stiffness of an elastic
material; this varies from a few thousand psi for rubber to a
few million psi for metals. Of course, under enough strain
all real materials break down and go nonlinear at some
point--that is, no material completely follows Hooke's Law (it's
more a guideline, really.)
- By contrast, in a plastic material, like wet
clay, the stress-strain curve is highly nonlinear--enough strain
causes the material to yield and permanently deform. Many
materials are work hardening, displaying an increase in
stiffness upon plastic deformation.
- A brittle material, like
dried clay, fractures into pieces instead of deforming.
The stress-strain curve breaks down entirely. Most metals
are elastic for small strains, plastic for some distance outside
that, and then brittle.
For liquids, the terminology is a bit
different. Most liquids respond purely elastically to
stresses that are equal in all directions (called hydrostatic
stress), but flow in response to shear. A viscous material,
like honey, pushes back strongly in response to shear
stress. As many solids approach their melting point, they
often display a form of viscous flow calledcreep; like glaciers flow downhill. Non-Newtonian materials
respond differently depending on the rate at which the deformation
is applied. For example, ketchup is thixotropic, getting less viscous at higher strain
rates; this is why ketchup squirts easily, but takes a long time
to flow down the bottle. By contrast, cornstarch is shear
thickening--it gets stronger at higher strain rates, so in a large
vat of cornstarch, one person can float in it while
another person walks over it (quickly!).
In real life, material stresses often
depend somewhat on all the above factors--real materials are
elastoviscoplastic and subject to creep and fracture.
Luckily, most of the time only one or two of these phenomena are
important at once!
Practice--Finite
Elements
So that's what really happens--each
point in the material examines its local strain and internal
properties, and then figures out how much stress to apply to its
neighboring points. To discretize this, we usually just take
an element of finite, nonzero size, like a triangle or
tetrahedron, calculate the strain over that element, compute the
stress based on that strain, and apply that stress to element's
boundary points ("nodes"). Like any discrete computation,
this isn't exact, but it's close enough to use to design working
machines.
In code, a typical finite element
simulation looks a lot like our little spring system:
// zero out node forces
for (n=0;n<nodes.size();n++) nodes[n].F=vec3(0);
// add element forces to nodes
for (e=0;e<elements.size();e++) {
mat3 strain = compute_strain(e); /* looks at node positions */
mat3 stress = constitutive_model(strain,e); /* e.g., elastoplastic */
for (n=0;n<elements[e].numnodes;n++) /* apply our stress to each node */
nodes[n].F+=apply(stress,elements[e].nodes[n]);
}
// move nodes based on net force
for (n=0;n<nodes.size();n++)
{
vec3 A=nodes[n].F/nodes[n].m; /* F=mA */
nodes[n].V+=dt*A; /* dv/dt = A */
nodes[n].P+=dt*V+0.5*dt*dt*A; /* dp/dt = V */
}
All the good stuff above is hidden
inside the constitutive_model call.
Implementation
I've written several FEM implementations
in 2D and 3D, and they work, but they tend to freak out in one
situation or another, such as large rotations, large deformations,
etc. This is a common bug. For engineering type
simulations, you don't care what happens if your bridge gets tied
into a knot (it's not a bridge by that point anyway); but for
games, you really want something reasonable to happen.
I'm excited by Irving, Teran, and Fedkiw's 2004
paper, which proposes a simple but robust
matrix-based approach for keeping track of the rest configuration
and computing strains, from which you can then calculate stresses
according to any constitutive response you like. He's got a
bunch of amazing
simulation videos on his site. The one downside
is the dense thicket of matrix math required to disentangle object
deformation from simple rotation. This matrix math is both
difficult to write and slow to run: in 2004, around 1 second per
frame (!) for just 11K elements.
A paper from France in 2005
by Nesme, Payan, and Faure recommended a QR decomposition to
find the rotation portion of the matrix, and so achieve 30fps on
small meshes of a few thousand elements.