Fin Elem

SEACAS Formulation of Nonlinear Constitutive Finite Element SPH Contact
Library Nonlinear Continuum Models Formulations Formulation Surfaces

Problems Mechanics
Finite Element Formulations

Theory
Manuals Introduction
< Go Back Weak Form Revisited
Blue text
Discretization
i indicates
a link to more Quasistatics
information.
Dynamics
Nonlinear Equation Solving
Basics of Element Design
Advanced Element Design Issues
Element Examples
Eight-Node Uniform Strain Element
Four-Node Corotational Shell
Theory Manuals (2/19/99) Finite Element Formulations - Main Menu 1

Introduction Weak Form Discretization Quasistatics Dynamics Nonlinear Basics of Advanced Element
SEACAS Equation Element Element Examples
Revisited
Library Solving Design Design Issues
Introduction
Theory Introduction
Manuals
Finite Element
Formulations
< Go Back
Blue text
i indicates
a link to more
information.
Theory Manuals (2/19/99) Finite Element Formulations - Introduction 1

Introduction
SEACAS
Library In this chapter we explore the finite element techniques utilized in the description of large
deformation problems in solid mechanics. Beginning with the notational framework and
problem description discussed in Formulation of Nonlinear Problems and utilizing the
Theory nonlinear continuum mechanics and material modeling issues discussed in Nonlinear
Manuals Continuum Mechanics and Constitutive Models, we discuss in this chapter how discrete
approximations to the governing nonlinear field equations are generated and solved.
The discussion will take place in three general stages. The first stage, consisting of the first
Finite Element
Formulations five sections of the report, emphasizes the global formulation of the finite element
method and treats aspects best understood by considering the discretized system in its
entirety. Topics to be discussed in this way include a brief presentation of weak forms
appropriate for large deformation problems, in Weak Form Revisited, a general
Introduction discussion of Discretization, time independent and dependent problems (i.e.,
Quasistatics and Dynamics), and Nonlinear Equation Solving. These sections will
< Go Back emphasize the derivation of discrete system equations from the underlying variational
principle, the form of these system equations in matrix form, and the iterative solution of
these equations that is required for nonlinear applications.
The next stage treats element technology, presenting the fundamentals necessary to
formulate and implement the basic building block of the finite element method: the finite
element. Indeed, the most basic advantage of the finite element method over other more
classical variational methods is its modularity. That is to say that the method of
discretization is tailored to small systematically generated subdomains of the problem of
Theory Manuals (2/19/99) Finite Element Formulations - Introduction - Introduction 2
interest (i.e., elements) making the method applicable to a myriad of geometrical
situations. Importantly much of finite element technology is sufficiently generic so that
SEACAS many aspects of element formulation are virtually unchanged from application to
Library application. We will discuss these aspects in two sections. The first, Basics of Element
Design, will cover the most essential features of element design including requirements
for global convergence, shape function definition, and numerical integration to produce
Theory local contributions to the global equations. The second section, Advanced Element
Manuals
Design Issues, deals with concerns more specific to large deformation solid mechanics
with the primary concern being near incompressibility of materials and the effect
numerical treatment of such phenomena.
Finite Element
Formulations The third stage of our presentation will consist of specific element examples,
summarizing some formulations that are in particularly prevalent use in computational
solid mechanics. In Eight-Node Uniform Strain Element, we present some of the
implementational details associated with an element widely used for the description of
Introduction three-dimensional continuous media, particularly in explicit dynamic and matrix-free
quasistatic applications. In Four-Node Corotational Shell, a common structural element
< Go Back
is discussed in some detail.
Theory Manuals (2/19/99) Finite Element Formulations - Introduction - Introduction 3

SEACAS Revisited Equation Element Element Examples
Weak Form Revisited

Theory Weak Form Revisited
Manuals
Finite Element
Formulations
< Go Back
i Blue text
indicates
a link to more
information.
Theory Manuals (2/19/99) Finite Element Formulations - Weak Form Revisited 1

Weak Form Revisited
SEACAS
Library We begin by providing a brief review of the field equations to be considered. The problem
to be solved is as shown schematically in Figure 1.7, where the finite deformation
response of a body, denoted Ω in its reference configuration, is to be computed. Assuming
Theory that this time-dependent configuration mapping is denoted by ϕ t , the following problem
Manuals
is to be solved for each time, t , in the time interval of interest:
∇ ⋅ T + f = ρa on ϕ t ( Ω ) , (3.1)
Finite Element
Formulations
ϕ t = ϕ t on ϕ t ( Γ u ) , (3.2)
and
Weak Form
Revisited t = t on ϕ t ( Γ σ ) , (3.3)
< Go Back
where all notations are as discussed in Notational Framework. In particular a is the
material acceleration expressed in spatial coordinates, f is the body force per unit
(spatial) volume, and T is the Cauchy stress tensor. The vector t is the Cauchy traction
vector, obtained via t = Tn , where n is the outward unit normal to the spatial surface
ϕt( Γσ ) .
The problem is also subject to initial conditions of the form
Theory Manuals (2/19/99) Finite Element Formulations - Weak Form Revisited - Weak Form Revisited 2
ϕ ( X, 0 ) = ϕ 0 ( X ) on Ω , (3.4)
SEACAS
and
Library
∂ϕ
( X, 0 ) = V 0 ( x ) on Ω . (3.5)
∂t
Theory Recall that although Eqs. (3.1)-(3.3) are written in the so-called spatial configuration, we
Manuals
still consider ourselves to be working in a Lagrangian framework, where all quantities are
ultimately indexed to material points through the mapping x = ϕ t ( X ) (see Lagrangian
and Eulerian Descriptions).
Finite Element
Formulations
A prerequisite of the finite element method is that a weak, or variational, form of the above
field equations be available for discretization. This can be obtained, following the general
procedure outlined for linear problems in Weak Forms, by considering weighting
Weak Form functions ϕ∗ , defined over Ω , which satisfy the following condition:
Revisited
< Go Back ϕ∗ = 0 on Γ u (3.6)
(cf. (1.82)), where we also assume that all ϕ∗ are sufficiently smooth so that any desired
partial derivatives can be computed. In treating large deformation problems, it is useful to
consider spatial forms of the functions ϕ∗ , obtained by composition with the (unknown)
–1
mapping ϕ t . We denote these spatial variations in the sequel by w , and note that they
may be obtained via
w ( x ) = ϕ∗ ( ϕ t ( x ) )
–1
(3.7)
SEACAS for any x ∈ ϕ t ( X ) . This causes the condition

Library
w = 0 on ϕ t ( Γ u ) (3.8)
Theory to be satisfied, and provided the configuration mapping ϕ t is smooth (which we assume to
Manuals
be the case), all required partial derivatives of w can be computed.
With these definitions in hand, the development in Weak Forms can be reproduced in the
Finite Element current context to provide the following spatial representation of the variational form for
Formulations
large deformations:
Given the boundary conditions t on ϕ t ( Γ σ ) , ϕ t on ϕ t ( Γ u ) , the initial conditions ϕ 0
Weak Form and V 0 on Ω , and the distributed body force f on ϕ t ( Ω ) , find ϕ t ∈ S t for each time
Revisited
t ∈ ( 0, T ) such that:
< Go Back
∫ ρw ⋅ a dv + ∫ ( ∇w ):T dv
ϕt( Ω ) ϕt( Ω )
(3.9)
 
= 
ϕ
∫ w ⋅ f dv + ∫ w ⋅ t da

t( Ω ) ϕt( Γσ )
for all admissible w , where S t is as defined as
S t = { ϕ t ϕ t = ϕ ( t ) on Γ u , ϕ t is smooth } (3.10)
SEACAS and where admissible w are related in a one-to-one manner via (3.7) to the material
variations ϕ∗ ∈ W with the definition of W being
Library
W = { ϕ∗ ϕ∗ = 0 on Γ u , ϕ∗ is smooth } . (3.11)
Theory
Manuals Note that in contrast to previous development, the constitutive relation governing T is left
unspecified: it can, in general, be subject to both geometric and material nonlinearities.
The notation a for the acceleration is to be understood as the material acceleration, as
Finite Element defined by (2.27) in Material and Spatial Velocity and Acceleration.
Formulations
In addition, the solution ϕ is subject to the following conditions at t = 0 :
Weak Form
∫ ϕ∗ ⋅ ( ϕ t = 0 – ϕ0 ) dΩ = 0 (3.12)
Revisited
Ω
< Go Back and
 ∂ϕ 
∗
∫ ϕ ⋅  ∂ t – V 0 dΩ = 0 , (3.13)
Ω t=0 
both of which must hold for all ϕ∗ ∈ W .
Revisited
Discretization
Theory Introduction
Manuals
Galerkin Finite Element Methods
Generation of Matrix Equations
Finite Element
Formulations
Localization and Assembly
< Go Back
i Blue text
indicates
a link to more
information.
Theory Manuals (2/19/99) Finite Element Formulations - Discretization 1

Introduction
SEACAS
Library The process of numerically approximating a continuous problem is generically called
discretization. In the finite element method, the entity discretized is a weak form
(alternatively, variational equation). In the current context the variational form to be
Theory considered is that described in Weak Form Revisited. We now refer the reader to Figure
Manuals 3.1, which gives the general notation to be used in description of the discretization process
for the problem at hand.
Finite Element
Formulations
e
Ω
Discretization
< Go Back























ℜ
3
Ω
Figure 3.1 General notation for finite element discretization of the reference domain.
As referred to in Figure 3.1, the reference domain Ω is subdivided into a number of

e
element subdomains, Ω , where the superscript e is an index to the specific element in
Theory Manuals (2/19/99) Finite Element Formulations - Discretization - Introduction 2

question, running between 1 and n el , where n el is the total number of elements required
for the discretization. We assume in the figure and throughout the ensuing discussion that
SEACAS 3
Library Ω is a subset of ℜ , with the two-dimensional case readily obtained as a special case of
the theory we will discuss.
Note also from Figure 3.1 that a number of nodal points are indicated by the dots. We
Theory shall assume that all degrees of freedom in the discrete system to be proposed will be
Manuals
associated with these nodes. As one might also notice, these nodes may lay at corners,
edges, and in interiors of the elements with which they are associated. A key feature of the
finite element method will be that a specific element can be completely characterized by
Finite Element
Formulations the coordinates and degrees of freedom associated with the nodes attached to it. In the
following we will index the nodes with uppercase letters A, B, C, etc., with such indices
running between 1 and n np , the total number of nodal points in the problem.
Discretization
< Go Back
Theory Manuals (2/19/99) Finite Element Formulations - Discretization - Introduction 3

Galerkin Finite Element Methods
SEACAS
Library The essence of any finite element method lies in the discretization of a weak or variational
form. This discretization process involves two important approximations: approximation
of a typical member of the solution space S t , and approximation of the weighting space
Theory W . These approximations are typically expressed as an expansion in terms of prescribed
Manuals
shape or interpolation functions, usually associated with specific nodal points in the mesh.
Since the number of nodal points is obviously finite, the expansion is likewise finite,
giving rise to the concept of a finite-dimensional approximation of a space.
Finite Element
Formulations
Roughly speaking, the idea of discretization is as follows. We know from earlier chapters
that if the variational equation is enforced considering all ϕ t ∈ S t and ϕ∗ ∈ W as
mandated by its definition, then the solution of the weak form is completely equivalent to
Discretization that of the strong form (i.e., the governing partial differential equation with boundary/
initial conditions). This fact results because of the arbitrary nature of the ϕ∗ and because
< Go Back
of the very general definitions for S t and W . If we restrict our attention only to some
subset of the above spaces, we now make an error with the solution of our approximated
weak form no longer being identical to the solution of the strong form. If our choice for
the type of shape functions to be used is reasonable, however, we can represent the full
solution and weighting spaces with arbitrarily closeness by increasing the number of nodal
points and/or the degree of polynomial approximation utilized in the interpolation
functions. In the limit of such refinement, we should expect recovery of the exact solution
(i.e., convergence).
Theory Manuals (2/19/99) Finite Element Formulations - Discretization - Galerkin Finite Element Methods 4
Let us represent the shape function associated with node A as N A , and assume it to be as
follows:
SEACAS
N A :Ω → ℜ .
Library
(3.14)
h
Given a time, t , the finite-dimensional counterpart of ϕ t will be denoted as ϕ t and is
Theory
Manuals expressed in terms of the shape functions as
nnp
∑
h
ϕt= N Bd B( t ) , (3.15)
Finite Element
Formulations B=1
where d B ( t ) is a 3-vector containing the (in general unknown) coordinates of nodal point
B at time t. Given a prescribed set of nodal shape functions { N B } , B = 1, …, n np , the
Discretization h h
finite dimensional solution space S t is defined as the collection of all such ϕ t :
< Go Back
nnp
h  h 
∑ B B
h
St=  ϕt= N d ( t ) ϕ t ≈ ϕ t ( X ) for all X ∈ Γ u . (3.16)
 B=1 
In other words, we require members of the discrete solution space to (approximately)

satisfy the displacement boundary condition on Γ u . The approximation comes about
h
because, in general, we only force ϕ t to interpolate the nodal values of ϕ t on Γ u with the
N B serving as the interpolation functions. We might also note that Γ u itself is typically
geometrically approximated by the finite element discretization, contributing also to the
SEACAS approximation.
Library
h
This defines the discretization procedure for ϕ t , at least notationally. It still remains,
however, to approximate the weighting space. The (Bubnov-) Galerkin finite element
Theory
Manuals method is characterized by utilizing the same shape functions to approximate W as were
used to approximate S t . Accordingly, we define a member of this space, ϕ∗ , via
h
Finite Element nnp

Formulations
ϕ∗ = ∑
h
NAcA , (3.17)
A=1
where the c A are 3-vectors of nodal constants. We can then express the finite dimensional
Discretization
h
weighting space W via
< Go Back
nnp
h  
W =  ϕ∗ = ∗
∑ N A c A ϕ = 0 for all X ∈ Γ u  .
h h
(3.18)
 A=1 
h
Analogous to the situation for S t , Eq. (3.18) features a discrete version of the boundary
h
condition on Γ u . In other words, W consists of all functions of the form (3.17) resulting
in satisfaction of this condition. Note that the only restriction on the c A is that they result
in satisfaction of the homogeneous boundary condition on Γ u ; they are otherwise
arbitrary.
SEACAS
Library With these ideas in hand, the approximate Galerkin solution to the initial/boundary value
problem takes the form described below.
h h
Theory
Given the boundary conditions t on ϕ t ( Γ σ ) , ϕ t on ϕ t ( Γ u ) , the initial conditions
Manuals
h h h
ϕ 0 and V 0 on Ω , and the distributed body force f on ϕ t ( Ω ) , find ϕ t ∈ S t for each
time t ∈ ( 0, T ) such that:
Finite Element
h h h h
∫ ρw ⋅ a dv + ∫ ( ∇w ):T dv
Formulations
h h
ϕt( Ω ) ϕt( Ω )
(3.19)
h h
Discretization
= ∫ w ⋅ f dv + ∫ w ⋅ t da
h h
ϕt( Ω ) ϕt( Γσ )
< Go Back
h h
for all admissible w , where S t is as defined in (3.16) and where admissible w are
related to the material variations ϕ∗ ∈ W via

h h
–1
w ( x ) = ϕ∗ ( ( ϕ t ) ( x ) ) .
h h h
(3.20)
h h
In Eq. (3.19) T refers to the Cauchy stress field computed from the discrete mapping ϕ t
h
SEACAS through the constitutive relations, whereas a is the discrete material acceleration.
Library
The initial conditions are ordinarily simplified in the discrete case to simply read:
d B( 0 ) = ϕ0 ( X B ) (3.21)
Theory
Manuals
and
d˙B ( 0 ) = V 0 ( X B ) , (3.22)
Finite Element
Formulations
both of which must hold for all nodes B = 1, …, n np , where X B are the reference
coordinates of the node in question.
Discretization
< Go Back
Generation of Matrix Equations
SEACAS
Library We are now in a position to summarize the discrete equations that will result from Eq.
(3.19). Before doing so, let us develop one more notational necessity. We can reexpress the
nodal vectors c A and d B in terms of their components via:
Theory
Manuals c A = { c iA }, i = 1, 2, 3 (3.23)
and
d B = { d jB }, j = 1, 2, 3 .
Finite Element
Formulations (3.24)
Note that indices i and j are spatial indices, in general. It is useful in generating matrix
equations to have indices referring not to nodes A and B or spatial directions i and j, but
Discretization rather to degree of freedom numbers in the problem. Toward this end we define for
notational convenience the concept of an ID array that is set up as follows:
< Go Back
ID ( i, A ) = P (global degree of freedom number) . (3.25)
In other words, the ID array takes the spatial direction index and nodal number as
arguments and assigns a global degree of freedom number to the corresponding unknown.
In general, the number of degrees of freedom is n dof , given by
n dof = 3 × n np . (3.26)
With this notation in hand, the equation numbers P and Q are defined as follows:
Theory Manuals (2/19/99) Finite Element Formulations - Discretization - Generation of Matrix Equations 9
P = ID ( i, A ) (3.27)
and
SEACAS
Library
Q = ID ( j, B ) . (3.28)
We now generate the discrete equations by substitution of Eqs. (3.15) and (3.17) into
Theory
(3.19), causing the variational equation to read:
Manuals
n n
 np –1   np –1 
∫  ∑ A t
ρ  N ( ϕ ( x ) )c A

⋅ 

∑ B t
N ( ϕ ( x ) )ḋ ˙
B ( t ) dv

h
ϕt( Ω ) A = 1 B=1
Finite Element n
Formulations  np –1  h
+ ∫  ∑ ∇N A ( ϕ t ( x ) ) ⊗ cA :T dv
h A = 1  , (3.29)
ϕt( Ω )
n n
 np –1   np –1 
= ∫  ∑ N A ( ϕ t ( x ) )cA ⋅ f dv + ∫  ∑ N A ( ϕ t ( x ) )cA ⋅ t da
h A = 1  h A = 1 
Discretization ϕt( Ω ) ϕt( Γσ )
< Go Back
nnp
∑
h h
where we note in particular that T is a function of ϕ t = N B d B ( t ) through the
B=1
strain-displacement relations (nonlinear, in general) and the constitutive law (as yet
unspecified and perhaps likewise nonlinear).
Proceeding now to examine (3.29) term-by-term, the inertial term can be expanded as
follows
n n
 np –1   np –1 
∫ ρ  ∑ N A ( ϕ t ( x ) )cA ⋅  ∑ N B ( ϕ t ( x ) )ḋB
˙( t ) dv
h A = 1  B = 1 
ϕt( Ω )
SEACAS
nnp 3 n
Library  –1  np –1  
= ∑ ∑ ∫  ρN A ( ϕ t ( x ) )ciA  ∑ N B ( ϕ t ( x ) )ḋiB
˙  dv
A = 1 i = 1 ϕh( Ω )  B = 1  
t
, (3.30)
nnp 3 nnp 3
–1 –1
Theory
Manuals
= ∑ ∑ ciA ∑ ∑ ∫ ρN A ( ϕ t ( x ) )δ ij N B ( ϕ t ( x ) ) dvḋjB
˙
A = 1i = 1 B = 1 j = 1 ϕh( Ω )
t
ndof n
 dof 
= ∑ cP  ∑ M PQ ḋQ
˙
Finite Element P = 1 Q = 1 
Formulations
where M PQ is defined as follows:
–1 –1
Discretization
M PQ = ∫ ρN A ( ϕ t ( x ) )δ ij N B ( ϕ t ( x ) ) dv . (3.31)
h
ϕt( Ω )
< Go Back The second term of (3.29) can be simplified via
n
 np –1  h
∫  ∑ ∇N A ( ϕ t ( x ) ) ⊗ cA :T dv
h A = 1 
ϕt( Ω )
SEACAS
n
Library  np 3 3
–1 h
= ∫  ∑ ∑ ∑ N A, j ( ϕ t ( x ) )ciA Tij  dv , (3.32)
h A = 1 i = 1 j = 1 
ϕt( Ω )
ndof
int
Theory = ∑ cP FP
Manuals P=1
where
Finite Element 3
int –1 h
Formulations FP = ∫ ∑ N A, j ( ϕ t ( x ) )Tij dv . (3.33)
h
ϕt( Ω ) j = 1
Finally, the last two terms of (3.29) can be treated as

Discretization n n n
 np –1   np –1  dof
ext
∫  ∑ N A ( ϕ t ( x ) )cA ⋅ f dv + ∫  ∑ NA ( ϕt ( x ) )cA ⋅ tda = ∑ cP FP , (3.34)
< Go Back h A = 1  h
ϕt( Ω ) ϕt( Γσ ) A = 1 P=1
where
ext –1 –1
FP = ∫ N A ( ϕ t ( x ) )fi dv + ∫ N A ( ϕ t ( x ) ) ⋅ ti da . (3.35)
h h
ϕt( Ω ) ϕt( Γσ )
We now define the following vectors and matrices of global variables, all with dimension
of n dof :
c = { cp}
d ( t ) = { d Q( t ) }
SEACAS
Library int int
F ( d ( t ) ) = { FP } . (3.36)
ext ext
F = { FP }
Theory M = [ M PQ ]
Manuals
The results of Eqs. (3.30)-(3.35) can now be summarized as follows:
T int ext
Finite Element c [ Mḋ ˙
(t) + F (d(t)) – F ] = 0, (3.37)
Formulations
which must hold for all n dof -vectors c that result in satisfaction of the homogeneous
boundary condition imposed on W (i.e., Eq. (3.18)).
Discretization
Finally, we make the observation that not all of the members of d ( t ) are unknown; for
< Go Back those nodes lying on Γ u , these degrees of freedom are prescribed. Furthermore, the
corresponding entries of c at these nodes are typically taken to be zero, so that the
h
aforementioned condition on W is obeyed. Since the remainder of the vector c is
arbitrary, it must be the case that the elements of the bracketed term in (3.37)
corresponding to unprescribed degrees of freedom must be identically zero, so that (3.37)
will hold for arbitrary combinations of the c P . Thus we can write the following nonlinear
equation that expresses the discrete equations of motion:
int ext
(t) + F
Mḋ ˙ (d(t)) = F . (3.38)
SEACAS Here we employ a slight abuse of notation because we have asserted in (3.36) that all
Library vectors and matrices have dimension n dof ; and yet we only enforce Eq. (3.38) for
unprescribed degrees of freedom. Denoting the number of unprescribed degrees of
freedom as n eq , one can account for this difference in practice by calculating the vector
Theory
Manuals and matrix entries for all degrees of freedom and then by merely disregarding the
n dof – n eq equations corresponding to the prescribed degrees of freedom. The members
of d ( t ) that are prescribed do need to be retained in its definition, however, since they
Finite Element enter into both terms on the left-hand side of (3.38). It should simply be remembered that
Formulations
only n eq members of d ( t ) are, in fact, unknown.
Discretization
< Go Back
Localization and Assembly
SEACAS
Library The development to this point is mostly a matter of mathematical manipulation with little
insight gained into the character of the interpolation functions, N A . In fact, the basic nature
of these interpolation functions distinguishes the finite element method from other
Theory variational solution techniques.
Manuals
The details of shape function construction will be discussed in Basics of Element Design
in the context of element programming. However, it is useful to discuss now the basic
Finite Element character of finite element approximation functions to give general insight into the
Formulations
structure of the method. We refer then to Figure 3.2, which depicts a node, A, in Ω and
some generic elements attached to it. A basic starting point for the development of a finite
element method is as follows: the shape function associated with Node A, N A , is only
Discretization nonzero in that subportion of Ω encompassed by the elements associated with Node A and
is zero everywhere else in Ω .
< Go Back
This property of the shape functions is crucial to the modular character of the finite
element method. Shape functions N A having this property are said to possess local
support.
Theory Manuals (2/19/99) Finite Element Formulations - Discretization - Localization and Assembly 15
NA
SEACAS
Library
Node A
Theory
Manuals
Finite Element
Formulations
Figure 3.2 Local support of finite element interpolation functions. Region of support for NA
shown as shaded.
To gain insight into the effect of this property, let us examine the expression given in Eq.
Discretization (3.31) for an element of the mass matrix M PQ . We note in particular that the integrand of
< Go Back (3.31) will only be nonzero if both Nodes A and B share a common element in the mesh;
otherwise, M PQ must be zero. If we fix our attention on a given Node A in the mesh, we
can, therefore, conclude that very few Nodes B will produce nonzero entries in M . This
matrix is, therefore, sparse; and it would be a tremendous waste of time to try to compute
M by looping over all the possible combinations of node numbers and spatial indices
without regard to elements and the node numbers attached to them.
Instead the global matrices and vectors needed in the solution of (3.38) are more typically
computed using two important concepts: localization and assembly. Still considering the
matrix M as an example, we note that by the elementary properties of integration, we can
write:
SEACAS
–1 –1
∫ ρN A ( ϕ t ( x ) )δ ij N B ( ϕ t ( x ) ) dv
Library
M PQ =
h
ϕt( Ω )
nel
Theory
∑
–1 –1
Manuals = ∫ ρN A ( ϕ t ( x ) )δ ij N B ( ϕ t ( x ) ) dv , (3.39)
e = 1 ϕh( Ωe )
t
nel
∑
e
Finite Element = M PQ
Formulations
e=1
where
e –1 –1
Discretization M PQ = ∫ ρN A ( ϕ t ( x ) )δ ij N B ( ϕ t ( x ) ) dv . (3.40)
h e
ϕt( Ω )
< Go Back
Thus, the global mass matrix can be computed as the sum of a number of element mass
e
matrices. This fact in itself is not especially useful because each of the M is extremely
e
sparse, even more so than M . In fact, the only entries of M that will be nonzero will be
those for which both P and Q are degrees of freedom associated with element e .
e
This fact can be exploited by defining another local element matrix m containing only
degrees of freedom associated with that element. We introduce element degree of freedom
SEACAS
Library indices p and q , as indicated generically in Figure 3.3. Assuming that p and q can take
on values between 1 and n edf , where n edf is the number of degrees of freedom
e
associated with the element, an n edf × n edf matrix m is constructed as follows:
Theory
Manuals
e e
m = [ m pq ] . (3.41)
p=8
Finite Element
p=6
Formulations
p=7 a=4 a=3 p=5
Discretization
< Go Back a=1 a=2

p=2 p=4
p=1 p=3
Figure 3.3 Element (local) degrees of freedom for a sample finite element.
e
The m pq can be specified by introducing the concept of a local node number a or b , as
also shown in Figure 3.3. With these definitions we can write
e –1 –1
m pq = ∫ ρN a ( ϕ t ( x ) )δ ij N b ( ϕ t ( x ) ) dv , (3.42)
h e
ϕt( Ω )
SEACAS
Library
where a sample relationship between indices i , a , and p appropriate for the element at
hand might be
Theory p = (a – 1) × 2 + i (3.43)
Manuals
(similarly for j , b and q ). The notation N a simply refers to the shape function associated
with local Node a . By definition it is the restriction of the global interpolation function
Finite Element N A to the element domain.
Formulations
e
Calculation of local element entities, such as m , turns out to be a highly modular
procedure whose form remains essentially unchanged for any element in a mesh. Detailed
Discretization discussion of this calculation is deferred until Basics of Element Design. Let us suppose
for a moment, however, that we have a procedure in hand for calculating this matrix. We
< Go Back
might then propose the following procedure for calculating the global mass matrix M and
int
internal force vector F :
int
Step 1: Zero out M , F .
Step 2: For each element e , e = 1, …, n el :
e
• a) Prepare local data necessary for element calculations – e.g., X ( n edf -vector of
e
SEACAS element nodal coordinates), d ( n edf -vector of element nodal configuration
Library
mappings), etc.
int e  int e 
Theory
• b) Calculate element internal force vector f = f p  and element mass
Manuals  
e e
matrix m = [ m pq ] via
Finite Element 3
int e
∑
–1 h
∫ N a, j ( ϕ t ( x ) )T ij dv
Formulations
f p = (3.44)
h e
ϕt( Ω ) j = 1
and Eq. (3.42).

Discretization
• c) Assemble the element internal force vector and element mass matrix into their global
< Go Back counterparts by performing the following calculations for all local degrees of freedom
p and q :
e
M PQ = M PQ + m pq (3.45)
and
int int int e

FP = FP +f p, (3.46)
where local degrees of freedom are related to global degrees of freedom via the LM
array, defined so that
SEACAS
Library P = LM ( p, e ) (3.47)
and
Theory Q = LM ( q, e ) . (3.48)
Manuals
Step 2a) above is referred to as localization; given a particular element, e , it extracts the
local information from the global arrays necessary for element level calculations. Step 2b)
Finite Element
consists of element level calculations; these computations will be discussed in detail in
Formulations Basics of Element Design. Step 2c) is the process known as assembly and takes the data
produced by the element level calculations and places them in the proper locations of the
global arrays.
Discretization We can thus now summarize the effect of localization and assembly in a finite element
architecture. They act as pre- and post-processors to the element-level calculations,
< Go Back enabling the entities needed for global equilibrium calculations to be computed in a
modular manner as a summation of element contributions. Of course, the effectiveness of
this procedure, as well as the convergence behavior of the numerical method in general,
depends crucially on the interpolation functions chosen and their definitions in terms of
elements. We defer this topic for now and concentrate in the coming sections on the
classes of problems and global equation-solving strategies to be utilized.
Revisited
Quasistatics
Theory Introduction
Manuals
Internal Force Vector
External Force Vector
Finite Element
Formulations
Incremental Load Approach
< Go Back
i Blue text
indicates
a link to more
information.
Theory Manuals (2/19/99) Finite Element Formulations - Quasistatics 1

Introduction
SEACAS
Library As discussed previously in the context of Linear Elastic IBVP in The Quasistatic
Approximation, the quasistatic approximation is appropriate when inertial forces are
negligible compared to the internal and applied forces in a system. As discussed in
Theory Discretization, the quasistatic system of equations is obtained by omission of the inertial
Manuals term in the discrete equations of motion. Thus in this section we discuss solution of
problems of the form:
int ext
Finite Element F (d(t)) = F (t) (3.49)
Formulations
subject to only one initial condition of the form
d ( 0 ) = d0 . (3.50)
Quasistatics
Note that the time variable t may correspond to real time (e.g., if rate-dependent material
< Go Back response is considered) but need not have physical meaning for rate-independent behavior.
For example, it is common for t to be taken as a generic parameterization for the applied
loading on the system, as discussed below in Incremental Load Approach.
It could also be noted that if the initial condition were taken as the same as the reference
configuration of the body, then
d0 = XA . (3.51)
A
Theory Manuals (2/19/99) Finite Element Formulations - Quasistatics - Introduction 2

Internal Force Vector
SEACAS
Library int
The quantity F ( d ( t ) ) is known as the internal force vector and consists of that set of
forces that are variationally consistent with the internal stresses in the body undergoing
analysis. The generic expression for an element in this vector is
Theory
Manuals
 3 h 
 ∑ N A, j ( ϕ t ( x ) )T ij dv ,
int –1
FP = ∫ j = 1 
(3.52)
h
ϕt( Ω )
Finite Element
Formulations as given in Generation of Matrix Equations. This vector-valued operator is, in general, a
nonlinear function of the unknown solution vector d ( t ) due to the possible Material
Nonlinearity and/or Geometric Nonlinearity inherent in the definition of the Cauchy
h
Quasistatics stress T ij in (3.52). As implied by our notation, we assume the solution vector d to be
< Go Back smoothly parameterized by t , which may represent time or some other loading parameter.
Theory Manuals (2/19/99) Finite Element Formulations - Quasistatics - Internal Force Vector 3
External Force Vector
SEACAS
Library ext
The external load vector F ( t ) must equilibrate the internal force vector, as is clear
from Eq. (3.49). As first presented in Generation of Matrix Equations, the expression
ext ext
Theory
for an element F P of F ( t ) is as follows
Manuals
–1
∫ N A ( ϕ t ( x ) )f i ( t ) dv
h
ext ϕt( Ω )
Finite Element FP = , (3.53)
–1
∫ NA ( ϕt ( x ) ) ⋅ t i ( t ) da
Formulations
+
h
ϕt( Γσ )
Quasistatics
where the explicit dependence of f i and t i upon t has been indicated and where
P = ID ( i, a ) , as given in (3.27). In other words, we assume that the prescribed internal
< Go Back
force loadings f i and prescribed surface tractions t i are given functions of t .
Equation (3.53) as written implies no dependence of either t i or f i upon ϕ t ( x ) (and

thus d ). Provided no such dependence exists, the external force vector is completely
parameterized by t , and the sole dependence of the equilibrium equations upon d occurs
int
through F . However, it is important to realize that some important loading cases are
precluded by this assumption, with perhaps the most important being the case of pressure
Theory Manuals (2/19/99) Finite Element Formulations - Quasistatics - External Force Vector 4
loading, where the direction of applied traction is opposite to the outward surface normal,
which in large deformation problems depends upon ϕ t ( x ) . Such a load is sometimes
SEACAS called a follower force and will, in general, contribute additional nonlinearity to the
Library
problem. Such complications are readily handled but are not encompassed by our current
notational framework for the sake of simplicity.
Theory
Manuals
Finite Element
Formulations
Quasistatics
< Go Back
Theory Manuals (2/19/99) Finite Element Formulations - Quasistatics - External Force Vector 5
Incremental Load Approach
SEACAS
Library We may now summarize the global solution strategy most commonly applied to
quasistatic nonlinear solid mechanics applications. We assume that we are interested in the
solution d ( t ) over some interval of interest for t :
Theory
Manuals t ∈ [ 0, T ] . (3.54)
We subdivide this interval of interest into a set of subintervals via
Finite Element N–1

Formulations [ 0, T ] = ∪ [ t n, t n + 1 ] , (3.55)
n=0
where n is an index on the time steps or intervals, and N is the total number of such
Quasistatics
increments. We assume that t 0 = 0 and that t N = T , but we do not, in general, assume
that all time intervals [ t n, t n + 1 ] have the same width.
< Go Back
With this notation in hand, the incremental load approach attempts to solve the following
problem successively in each time interval [ t n, t n + 1 ] :
Given the solution d n corresponding to time level t n , find d n + 1 corresponding to

t n + 1 satisfying:
int ext
F ( dn + 1 ) = F ( tn + 1 ) . (3.56)
Theory Manuals (2/19/99) Finite Element Formulations - Quasistatics - Incremental Load Approach 6
This governing equation is also often expressed by introducing the concept of a residual
vector R ( d n + 1 ) :
SEACAS
Library ext int
R ( dn + 1 ) = F ( tn + 1 ) – F ( dn + 1 ) . (3.57)
Solution of (3.56), therefore, amounts to finding the root of the equation

Theory
Manuals R ( dn + 1 ) = 0 . (3.58)
The physical meaning of this approach is depicted graphically in Figure 3.4. Starting with
Finite Element
an initial equilibrium state at t n , so that R ( d n ) = 0 , we introduce a prescribed load
Formulations
ext ext ext
increment ∆F = F ( tn + 1 ) – F ( t n ) and attempt to find that displacement
increment, d n + 1 – d n , that will restore equilibrium (i.e., result in satisfaction of (3.58)).
Quasistatics
This will require a nonlinear equation solving technique for determination of d n + 1 , a
topic that will be discussed further in Nonlinear Equation Solving.
< Go Back
int
SEACAS F (d)
Library
 F ext ( t
 n + 1)
∆F ext
 ext
Theory
 F ( tn)
Manuals
d
dn dn + 1
Finite Element
Formulations Figure 3.4 Simple illustration of the incremental load approach to quasistatic problems.
Quasistatics
< Go Back
Revisited
Dynamics
Theory Introduction
Manuals
The Semidiscrete Approach
Time-Stepping Procedures
Finite Element
Formulations
Explicit Finite Element Methods
Implicit Finite Element Methods
< Go Back
Blue text
i indicates
a link to more
information.
Theory Manuals (2/19/99) Finite Element Formulations - Dynamics 1

Introduction
SEACAS
Library We now restore the inertial terms to the discrete equation system and examine prospective
techniques for solution. To recap the key result of Generation of Matrix Equations, the
problem we consider now takes the form
Theory
int ext
Manuals Mḋ˙( t ) + F (d(t)) = F , (3.59)
to be solved for t ∈ [ 0, T ] , subject to the initial conditions

Finite Element
Formulations
d ( 0 ) = d0 (3.60)
and
ḋ ( 0 ) = v 0 . (3.61)
Dynamics
< Go Back
Theory Manuals (2/19/99) Finite Element Formulations - Dynamics - Introduction 2

The Semidiscrete Approach
SEACAS
Library It might be noted from Eq. (3.59) that time remains continuous in our formulation at this
point, whereas the spatial discretization has already been achieved by the finite element
interpolations summarized in Discretization. This type of finite element approach to
Theory transient problems is sometimes referred to as the semidiscrete finite element method,
Manuals since the approximation in space is performed first, leaving a set of equations discrete in
space but still continuous in time. To complete the approximation, a finite differencing
procedure is generally applied in time, as discussed below.
Finite Element
Formulations
Dynamics
< Go Back
Theory Manuals (2/19/99) Finite Element Formulations - Dynamics - The Semidiscrete Approach 3
Time-Stepping Procedures
SEACAS
Library As discussed in Quasistatics, we subdivide the time interval of interest [ 0, T ] via
N–1
[ 0, T ] = ∪ [ t n, t n + 1 ] (3.62)
Theory n=0
Manuals
and consider the following generic problem. Given algorithmic approximations for the
solution vector ( d n ), velocity ( v n ), and acceleration ( a n ) at time t n , find approximations
Finite Element d n + 1 , v n + 1 and a n + 1 for these quantities at time t n + 1 . Note that in contrast to the
Formulations
quasistatic problem, the variable t here does have the interpretation of actual time.
Several time-stepping algorithms have been proposed for this incremental problem we
have posed. So that we might have a template with which to work, we will consider
Dynamics
perhaps the most pervasive of these schemes: the Newmark family of temporal integrators
< Go Back ([Newmark, N.M., 1959]). This algorithm can be summarized in a time step [ t n, t n + 1 ]
as follows:
int ext
Ma n + 1 + F ( dn + 1 ) = F ( tn + 1 )
2
∆t
dn + 1 = d n + ∆tv n + ---------- [ ( 1 – 2β )a n + 2βa n + 1 ] , (3.63)
2
vn + 1 = v n + ∆t [ ( 1 – γ )a n + γ a n + 1 ]
Theory Manuals (2/19/99) Finite Element Formulations - Dynamics - Time-Stepping Procedures 4

where β and γ are algorithmic parameters that define the stability and accuracy
characteristics of the method.
SEACAS
Library Although, obviously, a wide range of algorithms exist corresponding to the different
available choices of β and γ , two algorithms in particular are prevalent in common use:
1
Theory 1. Central differences ( β = 0 , γ = --- ). This integrator is second order accurate and
Manuals
2
only conditionally stable, meaning that linearized stability is only retained when ∆t
is less than some critical value. This algorithm is an example of an explicit finite
element integrator, to be discussed in Explicit Finite Element Methods.
Finite Element
Formulations
1 1
2. Trapezoidal rule ( β = --- , γ = --- ). This integrator is second-order accurate and
4 2
unconditionally stable for linear problems, meaning that the spectral radii of the
Dynamics integrator remain less than 1 in modulus for any time step ∆t (in linear problems).
This algorithm is an example of an implicit finite element integrator, to be discussed
< Go Back in Implicit Finite Element Methods.
Theory Manuals (2/19/99) Finite Element Formulations - Dynamics - Time-Stepping Procedures 5

Explicit Finite Element Methods
SEACAS
Library 1
Examining the central differences algorithm as an example, let us take β = 0 , γ = --- and
2
substitute into Eq. (3.63). Upon doing so, we obtain the following algorithm:
Theory
–1 ext int
Manuals
an + 1 = M ( F ( tn + 1 ) – F ( dn + 1 ) )
2
∆t
dn + 1 = d n + ∆tv n + ---------- a n , (3.64)
Finite Element
2
Formulations
∆t
vn + 1 = v n + ------- [ a n + a n + 1 ]
2
where the first equation has been written as solved for a n + 1 .

Dynamics
Equation (3.64) can be used to explain why this formulation is termed explicit. Consider
< Go Back the case where M is a diagonal matrix. This is not, in general, the case if we strictly follow
the variational formulation; reference to Eq. (3.31) will verify that unless two shape
functions N A and N B are mutually orthogonal, the mass matrix will not, in general, be
diagonal. However, it is common practice, as will be discussed in Basics of Element
Design, to diagonalize the mass matrix. In the event that this is done, Eq. (3.64) shows that
given the three vectors { a n, v n, d n } , the data at t n + 1 , { a n + 1, v n + 1, d n + 1 } can be
computed explicitly (i.e., without the need for solution of coupled equations).
Theory Manuals (2/19/99) Finite Element Formulations - Dynamics - Explicit Finite Element Methods 6
Although this form of the central difference formulation is readily obtained from the
Newmark formulae, it does not give insight into the source of the “central difference”
SEACAS terminology and, in fact, does not represent the manner in which the integrator is
Library ordinarily implemented. To see the usual form, let us define the following auxiliary
algorithmic velocity vector:
1
Theory v 1 = v n + --- ∆ta n , (3.65)
Manuals n + --- 2
2
which also implies a corresponding relation for the previous time step:
Finite Element
1
Formulations v 1 = v n – 1 + --- ∆ta n – 1 . (3.66)
n – --- 2
2
Subtracting Eq. (3.66) from (3.65) gives

Dynamics
1
v 1–v 1 = v n – v n – 1 + --- ∆t ( a n – a n – 1 ) . (3.67)
< Go Back n + --- n – --- 2
2 2
However, evaluation of (3.64) during the time step [ t n – 1, t n ] reveals that
1
v n – v n – 1 = --- ∆t ( a n + a n – 1 ) , (3.68)
2
so that upon substitution into (3.67) we find
v 1 –v 1 = ∆ta n . (3.69)
n + --- n – ---
2 2
SEACAS
Library Furthermore, substitution of (3.65) into the second equation of (3.64) gives
d n + 1 = d n + ∆tv 1. (3.70)
n + ---
2
Theory
Manuals Thus by collecting these latest two results, together with the equilibrium equation
evaluated at t n , we can reexpress the algorithm completely equivalently as
–1 ext int
Finite Element an = M ( F ( tn) – F ( dn) )
Formulations
v 1 = v 1 + ∆ta n
n + ---
2
n – ---
2
. (3.71)
d n + 1 = d n + ∆tv 1
Dynamics n + ---
2
< Go Back Note that the velocity and displacement updates emanate from centered difference
approximations to the acceleration a n and velocity v 1 , respectively, giving the
n + ---
2
algorithm its name. The velocity measures that are utilized by the algorithm are shifted by
a half step from the time values at which the acceleration and configuration are monitored.
As mentioned above, explicit finite element schemes are only conditionally stable,
meaning that they only remain stable when the time increment ∆t is less than some
critical limit. This limit, sometimes called the Courant stability limit, can be shown to be
as follows
SEACAS
2
Library
∆t ≤ ---- , (3.72)
ω
where ω is the highest modal natural frequency in the mesh. It can also be shown that this
Theory frequency can be conservatively estimated via
Manuals
ω ≈ 2  ---
c
, (3.73)
 h max
Finite Element
Formulations
where c and h are the sound speed and characteristic mesh size, respectively, associated
with the element in the mesh having the largest ratio of these two quantities. Combining
Eqs. (3.72) and (3.73) we find that
Dynamics
∆t ≤  ---
h
. (3.74)
 c min
< Go Back
In other words, the time step may be no larger than the amount of time required for a
sound wave to traverse the element in the mesh having the smallest transit time. This fact
tells us immediately that explicit finite element methods are most appropriate for those
problems featuring very high frequency response or wave-like phenomena. For problems
featuring low frequency response, literally thousands of time steps may be required to
resolve even a single period of vibration due to the stringent stability limit posed by (3.74).
For such problems an unconditionally stable algorithm is highly desirable, albeit at the
cost of explicit updates in each increment.
SEACAS
Library
Theory
Manuals
Finite Element
Formulations
Dynamics
< Go Back
Implicit Finite Element Methods
SEACAS
Library To introduce the concept of an implicit finite element method, we examine the trapezoidal
1
rule, which is simply that member of the Newmark family obtained by setting β = --- and
4
1
γ = --- . Substitution of these values into Eq. (3.63) yields
Theory
Manuals
2
int ext
Ma n + 1 + F ( dn + 1 ) = F ( tn + 1 )
Finite Element
2
Formulations ∆t
dn + 1 = d n + ∆tv n + ---------- [ a n + a n + 1 ] . (3.75)
4
∆t
vn + 1 = v n + ------- [ a n + a n + 1 ]
2
Dynamics
Insight into this method can be obtained by combining the first two equations in (3.75) and
< Go Back
solving for d n + 1 . Doing so gives
Theory Manuals (2/19/99) Finite Element Formulations - Dynamics - Implicit Finite Element Methods 11
4 F
ext
( tn + 1 )
---------- Md n + 1
2
∆t =
+ M  a n + ∆tv n + ---------2- d n
SEACAS 4
Library int
( dn + 1 )  
+F ∆t
. (3.76)
= ---------2- ( d n + 1 – d n )  – ------- v n – a n
4 4
an + 1
 ∆t 
Theory ∆t
Manuals
∆t
vn + 1 = v n + ------- [ a n + a n + 1 ]
2
Finite Element
Clearly solving the first equation in (3.76) is the most expensive procedure involved in
Formulations updating the solution from t n to t n + 1 . This equation is not only fully coupled, but also
is highly nonlinear, in general, due to the internal force vector. In fact, we could write the
first equation of (3.76) in terms of a dynamic incremental residual R n + 1 via
Dynamics
( t n + 1 ) + M  a n + ∆tv n + ---------2- d n
ext 4
F
< Go Back  
∆t
R ( dn + 1 ) =
–  ---------2- Md n + 1 + F ( d n + 1 )
4 int . (3.77)
 
∆t
= 0
This system has the same form as (3.57), which suggests that the same sort of nonlinear
solution strategies are needed for implicit dynamic calculations as in Quasistatics. Some
common equation-solving alternatives are discussed in Nonlinear Equation Solving.
Theory Manuals (2/19/99) Finite Element Formulations - Dynamics - Implicit Finite Element Methods 12
Revisited
Nonlinear Equation Solving

Theory Introduction
Manuals
Newton Raphson Framework
Line Search
Finite Element
Formulations
Quasi-Newton Methods
Conjugate Gradient Methods
< Go Back Preconditioning
i Blue text
indicates
a link to more
information.
Theory Manuals (2/19/99) Finite Element Formulations - Nonlinear Equation Solving 13

Introduction
SEACAS
Library In this section we explore some of the alternatives available for solving the nonlinear
discrete equations associated with computation of an unknown state at t n + 1 , in either the
context of a quasistatic problem (i.e., Eq. (3.57)) or an implicit dynamic formulation (Eq.
Theory (3.77)). In either case, the equation to be solved takes the form
Manuals
R ( dn + 1 ) = 0 , (3.78)
Finite Element where the residual R ( d n + 1 ) is considered to be a nonlinear function of the solution vector
Formulations
dn + 1 .
Nonlinear
Equation
< Go Back
Theory Manuals (2/19/99) Finite Element Formulations - Nonlinear Equation Solving - Introduction 14
Newton Raphson Framework
SEACAS
Library We now return to the general concept of a Newton-Raphson iterative solution technique,
as discussed in the one-dimensional context in Material Nonlinearity. To review, a
Newton-Raphson solution technique for (3.78) is defined in iteration i by
Theory
Manuals i ∂R
R ( dn + 1 ) + ∆d = 0 , (3.79)
∂d
i
dn + 1
followed by the update

Finite Element
Formulations
i+1 i
d n + 1 = d n + 1 + ∆d . (3.80)
i
Iterations on i typically continue until the Euclidean norm R ( d n + 1 ) is less than some
Nonlinear
i
tolerance; ∆d is smaller than some tolerance, the quantity R ( d n + 1 ) ⋅ ∆d is small, or
Equation
< Go Back
some combination of these three conditions.
It is instructive to examine the form taken by Eq. (3.79) for the quasistatic and implicit
i
dynamic cases. For the quasistatic case R ( d n + 1 ) takes the form
i ext int i
R ( dn + 1 ) = F ( tn + 1 ) – F ( dn + 1 ) , (3.81)
so that Eq. (3.79) can be rewritten as
Theory Manuals (2/19/99) Finite Element Formulations - Nonlinear Equation Solving - Newton Raphson Framework 15
i i
K ( d n + 1 )∆d = R ( d n + 1 ) , (3.82)
i
SEACAS
Library where the incremental stiffness matrix K ( d n + 1 ) is given by
int
i ∂R ∂F
Kd n + 1 = – = i . (3.83)
∂d ∂d
i
dn + 1 dn + 1
Theory
Manuals
Thus application of the Newton-Raphson method to quasistatic problems amounts to
solution of successive linear problems, as defined by (3.82).
Finite Element In the implicit dynamics case, let us consider the trapezoidal rule as a template. In this
Formulations
case, the residual is of the form
( t n + 1 ) + M  a n + ∆tv n + --------- 
ext 4
F d
 2  -
n
i ∆t
Nonlinear R ( dn + 1 ) = = 0. (3.84)
–  --------- 
Equation 4 i int i
Md + F ( d )
< Go Back  2 n -
+ 1 n + 1 
∆t
This causes (3.81) to take the form
4 i i
M +
---------
2
K ( d
- n+1 ) ∆d = R ( d n + 1) , (3.85)
∆t
i
where the stiffness matrix K ( d n + 1 ) is as given in (3.83). In either case solution of the
global incremental equations will require the assembly of the coefficient matrix on the
left-hand side. Following the same assembly procedures outlined in Localization and
i
k  d n + 1 ,
e e
Assembly, this matrix is given as an assembly of element stiffness matrices
SEACAS  
Library
each of which can be expressed generally as
i i
k  d n + 1 k pq  d n + 1
e e e e
= , (3.86)
Theory    
Manuals
where
int
i ∂f p i
k pq  d n + 1  de  ,
Finite Element e e
= (3.87)
Formulations
  e  n + 1
∂dq
int
where f p is as given in Eq. (3.44).
Nonlinear
Equation
We can, therefore, conclude that for a Newton-Raphson treatment of either a quasistatic or
< Go Back implicit dynamic system, an important function of the element subroutine is to return an
element stiffness in addition to the internal force vector and mass matrix that may also be
required. We will discuss in detail the mechanics of this operation in Basics of Element
Design.
Line Search
SEACAS
Library It is noteworthy that the Newton-Raphson method is only guaranteed to be convergent in
an asymptotic sense, subject also to some smoothness and differentiability conditions.
This means that solution updates may not be effective if one is excessively far from the
Theory solution or if significant nonsmoothnesses are present in the equation system. Indeed, in
Manuals many problems the early displacement updates in a given load increment given by (3.82)
or (3.85) may actually be counterproductive in that they take one farther from the solution
rather than closer. It is, therefore, imperative to have a technique that controls the manner
Finite Element in which the solution is sought such that bad displacement updates, as predicted by the
Formulations
linearized kernel, are not allowed to carry one too far from the desired solution.
The concept of line search, pervasive in nonlinear equation solving, is employed for this
purpose. To motivate the concept, we consider the case of a so-called quadratic system,
Nonlinear where the total system energy Π ( d ) can be expressed as a quadratic function of the
Equation
solution vector d via
< Go Back
T
1 T ext
Π ( d ) = --d- Kd – F d, (3.88)
2
ext
where for simplicity we assume K and F to be constant. We seek the minimizer of
Π ( d ) , which, of course, can be equivalently expressed as the solution of
ext
Kd = F . (3.89)
Theory Manuals (2/19/99) Finite Element Formulations - Nonlinear Equation Solving - Line Search 18
In fact, this problem statement is a finite-dimensional analogue of that discussed in Weak
Forms, where it was asserted that the linearly elastic boundary problem can be solved by
SEACAS finding the displacement field, minimizing the total potential energy of the system. We
Library could, therefore, think of the problem we have posed as a finite element discretization of
such a system. Although the system we consider here is quasistatic, the technique we
motivate is utilized for solution of nonlinear dynamic systems as well.
Theory
Manuals Neglecting the fact that this problem could be solved via Gaussian elimination, we
i
consider a generic iterative procedure for solving it.Suppose we have a current iterate d
as well as a proposed displacement increment ∆d . In a Newton-Raphson method ∆d
Finite Element would be computed in a given iteration by solving (3.82), which would produce the exact
Formulations
solution (to machine precision) after the displacement update (3.80). If ∆d is not such a
good choice for the displacement increment, however, we would like a method for
detecting this fact and for controlling growth in the residual. In this discussion we take ∆d
Nonlinear as a search direction in n eq -space and look for solutions in this direction reducing (or at
Equation
< Go Back least controlling growth in) the residual.

i
Given d and ∆d , then, we introduce a search parameter s and consider an update of the
form:
i+1 i
d ( s ) = d + s∆d . (3.90)
The line search parameter is chosen such that the update produced by (3.90) is in some
sense optimal. In this spirit we choose s as the minimizer of:
1 i T i
--(
- d + s∆d ) K ( d + s∆d )
Π̂ ( s ) = 2 , (3.91)
SEACAS i T ext
Library – ( d + s∆d ) F
∂
which can be found by finding the solution of ∂ s Π̂ ( s ) = 0 . If we assume that K is
Theory
Manuals
symmetric, positive definite, one finds in taking this derivative that s is given as the
solution of
T i ext
∆d ( K ( d + s∆d ) – F ) = 0. (3.92)
Finite Element
Formulations
Two forms of this equation are useful under various circumstances. First, in the linear
system we now consider, (3.92) is readily solved to explicitly yield s :
T i
Nonlinear ∆d R
Equation s = ---------, (3.93)
T
∆d K∆d
< Go Back
i ext i
where R = F – Kd . This form of the line search is actually used in some
implementations but depends, strictly speaking, on linearity to be effective. Thus a more
generally used form is generated by reexpressing (3.92) as the following problem:
Find s such that
T i
∆d R ( d + s∆d ) = 0 . (3.94)
i ext i
In (3.94), R ( d + s∆d ) = F – K ( d + s∆d ) in the linear case. The advantage of Eq.
(3.94), however, is that it admits more general representations of the residual; for a
SEACAS i
Library nonlinear quasistatic problem, we can use (3.94) with R ( d + s∆d ) given by
i ext int i
R ( d + s∆d ) = F –F ( d + s∆d ) . (3.95)
Theory
Manuals
Similar generalizations for the dynamic case are, of course, also possible with the dynamic
residual given for the trapezoidal rule in (3.84).
From (3.95) we conclude that a line search procedure looks for updated iterates where the
Finite Element search direction is orthogonal to the residual. This is equivalent to an energy minimization
Formulations
in the linear case, whereas the interpretation in the nonlinear case is not quite so
straightforward. Furthermore, in the nonlinear case it is not efficient or even necessary to
find the root of (3.94) to machine precision. More commonly one uses some sort of root
Nonlinear finder to find an s between 0 and 1 that satisfies (3.94) to some tolerance. Making the
Equation definition
< Go Back
T i
G ( s ) = ∆d R ( d + s∆d ) , (3.96)
a typical algorithm to find s could be outlined as:

i
Given d and ∆d
• IF ( G ( 1 ) > TOL × G ( 0 ) or G ( 1 ) × G ( 0 ) < 0 ) THEN
Iterate for s ∈ ( 0 , 1 ] such that G ( s ) < ( TOL ) G ( 0 ) (3.97)

• ELSE
s = 1. (3.98)
SEACAS
Library • ENDIF
The check in the IF statement amounts to checking whether a full step (with s = 1 ) leads
Theory
to an unreasonably large increase in G and whether a root might reasonably be expected in
Manuals the interval ( 0 , 1 ] .
Finite Element
Formulations
Nonlinear
Equation
< Go Back
Quasi-Newton Methods
SEACAS
Library One can establish that Newton-Raphson iteration is quadratically convergent
asymptotically, meaning that the error associated with a given iteration tends to the square
of the previous iteration’s error as iterations proceed. One can, in fact, roughly see the
Theory reason for this fact from Eq. (1.21), which states that the Newton-Raphson update is
Manuals motivated by a first order Taylor series expansion of the residual about the current solution
vector iterate. We might, therefore, expect that the error incurred from this approximate
update should tend toward the square of the displacement update as iterations proceed,
Finite Element which is, indeed, the case. A much more rigorous derivation of this property can be found
Formulations
in [Kelley, C.T., 1995]. This quadratically convergent behavior is a highly desirable
property and makes the Newton-Raphson method much more rapidly convergent than
many other equation-solving alternatives.
Nonlinear However, it is also true to say that a full Newton-Raphson method can be tremendously
Equation
expensive due to the necessity of solving the successive fully coupled linear systems
< Go Back implied by (3.82) and (3.85). If we choose to use a direct equation-solving technique for
solving these systems, such as Gaussian elimination, the cost of solving each linearized
problem will vary as the cube of the number of equations ( n eq ) . For very large problems
this cost can become prohibitive.
In response to this situation, a number of methods, known collectively as quasi-Newton
(alternatively, secant) methods, have been developed. These methods replace most of the
Newton-Raphson iterations with a cheaper update to the solution vector, sacrificing
convergence performance but making the average equilibrium iteration much less
Theory Manuals (2/19/99) Finite Element Formulations - Nonlinear Equation Solving - Quasi-Newton Methods 23
expensive. The interested reader should consult [Dennis, J.E. and Schnabel, R.B., 1996]
or [Kelley, C.T., 1995] for excellent overviews of these methods from a generic, nonlinear
SEACAS equation-solving viewpoint.
Library
Before discussing specific quasi-Newton methods, let us motivate them through
consideration of the term “secant method”. Suppose we have a scalar-valued, nonlinear
equation
Theory
Manuals
R(d) = 0 (3.99)
and wish to employ an iterative method to obtain the root. If we are currently performing
Finite Element iteration i and wish to obtain the next iterate i + 1 , a secant method will do this by
Formulations
d
replacing the tangent to the curve, K i = – [ R ( d i ) ] , with the secant
dd
R ( di) – R ( di – 1 )
K̃ i = ----------------- in the Newton-Raphson updating scheme (see Figure 3.5).
Nonlinear di – di – 1
Equation
Thus, the next iterate is obtained via
< Go Back
R ( di) – R ( di – 1 )
di + 1 R ( di) .
= d i – ------------------------ (3.100)
di – di – 1
R(d)
SEACAS
Library
K̃ i
di di – 1 d
Theory
Manuals
Figure 3.5 One-dimensional illustration of quasi-Newton (secant) iteration.
In generalizing this concept to multiple directions, we seek to approximate so-called

consistent tangents needed by Newton-Raphson updates (see (3.83)) with stiffnesses that
Finite Element
Formulations will be cheaper to compute and invert. In a secant method, using the one-dimensional
example as motivation, we demand that these approximate tangents obey the so-called
quasi-Newton equation:
Nonlinear K̃ i ( d i – d i – 1 ) = R ( d i ) – R ( d i – 1 ) , (3.101)
Equation
< Go Back or in terms of the inverse:
–1
d i – d i – 1 = K̃ i ( R ( d i ) – R ( d i – 1 ) ) . (3.102)
In the one-dimensional case either expression implies uniquely the secant method already
discussed but with multiple unknowns, expressions (3.101) or (3.102) place a less
stringent restriction. There are, therefore, a multitude of such methods, collectively termed
quasi-Newton or secant methods, whose defining feature is the satisfaction of (3.101) and
(3.102). Here we concentrate on one particular method, the BFGS (Broyden-Fletcher-
Goldfarb-Shanno) method, proposed originally and most coherently for specific use in
finite element calculations by [Matthies, H. and Strang, G., 1979].
SEACAS
Library
In the BFGS method, one typically starts with an assembled Newton-Raphson tangent
given, for example, by (3.87). One performs one iteration with this tangent (probably
including a line search). Rather than repeating this procedure for subsequent iterations, the
BFGS method takes the tangent from the Newton-Raphson iteration and updates it in a
Theory
Manuals manner consistent with (3.102) and uses this tangent to compute the next iterate for the
solution (also probably including a line search in the update).
The BFGS method is effective in many circumstances because the update to the tangent
Finite Element matrix is inexpensive and is actually done to a previously determined inverse so that no
Formulations
matrix inversion is necessary in most equilibrium iterations. To be more specific, we
suppose that the last tangent utilized in an iteration process is K̃ i – 1 . The BFGS update is
defined as
Nonlinear
Equation –1 T –1 T
K̃ i = ( I + v i w i )K̃ i – 1 ( I + w i v i ) , (3.103)
< Go Back
where
∆d i – 1
v i = ----------------------------, (3.104)
∆d i – 1 ⋅ ( R ( d i ) – R ( d i – 1 ) )
w i = – ( R ( d i ) – R ( d i – 1 ) ) + αiR ( d i – 1 ) , (3.105)
where
1
– s i – 1 ( R ( d i ) – R ( d i – 1 ) ) ⋅ ∆d i – 1 --
2
α i = ------------------------------------. (3.106)
SEACAS R ( d i – 1 ) ⋅ ∆d i – 1
Library
Throughout the above equations we have indexed the search directions and line search
parameters such that
Theory
Manuals
d i = d i – 1 + s i – 1 ∆d i – 1 . (3.107)
The next search direction ∆d i is then computed via
Finite Element –1
Formulations ∆d i = K̃ i R ( d i )
. (3.108)
T –1 T
= (I + v i w i )K̃ i – 1 ( I + w i v i )R ( d i )
Nonlinear
One may discern from (3.108) why a BFGS iteration is so much cheaper than a Newton-
Equation –1
Raphson iteration. Keeping in mind that K̃ i – 1 is typically stored in practice as a
< Go Back
factorized stiffness matrix, one can compute ∆d i efficiently by proceeding right to left in
the second line of (3.108). Thinking in this manner, the update consists only of dot
products, scalar vector multiplies, and a backsolve procedure.
The BFGS method is, like other quasi-Newton methods, superlinear in convergence rate,
meaning that the error decreases in a manner faster than linear but not as fast as the
quadratic rate displayed by Newton-Raphson. Thus it is most effectively used in large
problems, where this disadvantage in convergence behavior is offset by its great savings in
the iteration process. It should also be noted that the success of BFGS depends critically
upon the incorporation of line search. Since the iterations are based less directly on the
SEACAS underlying mechanics of the system than they would be in Newton-Raphson, it is
Library particularly important that the line search prevent excessive excursions away from the
solution in the case of bad search directions. Typically, BFGS solvers also contain
provisions to compute new Newton-Raphson tangents to restart the iteration process in the
Theory event the BFGS iterations are ineffective.
Manuals
Finite Element
Formulations
Nonlinear
Equation
< Go Back
Conjugate Gradient Methods
SEACAS
Library The desire to solve very large problems has recently led researchers to consider so-called
indirect iterative, or matrix-free strategies, where the set of nonlinear equations is
iteratively solved without any need to compute, store, or invert a tangent matrix at any
Theory stage of the iteration process. Perhaps the most celebrated iterative technique for solving
Manuals linear equations in the last two to three decades has been the conjugate gradient method. In
this section we briefly consider its extension to matrix-free, nonlinear equation solving. To
do so, however, it is useful to examine the linear case first, from which the nonlinear
Finite Element algorithms are readily derived.
Formulations
We begin then by considering the same linear system and potential energy function Π
given in (3.88) and consider that this potential energy is to be minimized. We note that at
any prospective solution point d , the steepest descent direction of the objective function
Nonlinear
Equation Π is given by the negative of the gradient:
< Go Back
∂ ext
– (Π )= F – Kd = R ( d ) . (3.109)
∂d
In other words, the steepest descent direction at any prospective solution d is merely given
by the residual R ( d ) . One of the most elementary methods from nonlinear equation
solving/optimization, the steepest descent method, utilizes at each iteration the current
steepest descent direction as the search direction, with a subsequent line search defining
the update to the solution vector. This method, while intuitive, is not as effective as other
alternatives including, in particular, the conjugate gradient strategy we now discuss.
Theory Manuals (2/19/99) Finite Element Formulations - Nonlinear Equation Solving - Conjugate Gradient Methods 29
The difficulty with steepest descent is that in many cases, successive search directions
“zig-zag,” meaning that a given search direction will contain significant components in the
SEACAS direction of previous steps. It is readily imagined that such repetition is wasteful, since
Library presumably these earlier iterations have already eliminated, or at least markedly reduced,
the error in their respective search directions.
This discussion is made more quantitative by recalling the formula for a generic line
Theory
Manuals search, Eq. (3.94). In the case of a linear problem, this condition takes the form
T i T ext i
∆d R ( d + s∆d ) = ∆d ( F – K ( d + s∆d ) )
. (3.110)
T i
Finite Element = ∆d K ( d – ( d + s∆d ) ) = 0
Formulations
i
In (3.110) the term d – ( d + s∆d ) is recognized as the error associated with the next
i
iterate d + s∆d . Thus we can see from (3.110) that the line search criterion causes the
Nonlinear error associated with the next iterate for the solution vector to be K-orthogonal to the
Equation
search direction. Thus if we consider the solution space to be described by a sequence of
< Go Back
vector spaces of increasing dimension, where search directions form the basis and the
stiffness matrix K serves as a metric, then each iteration with line search removes all error
in that direction. It is, therefore, wasteful to have subsequent search directions that have
non-zero components in previous directions. The aim of the conjugate gradient method is
to orthogonalize this direction as iterations proceed.
To begin derivation of the method, let us use the notations p for a search direction and α
for a line search parameter, as opposed to the ∆d and s used previously (in large part to
conform to usage in the literature). The restriction we will place on the search directions is
that they will be K-orthogonal:
SEACAS
Library T
p i Kp j = 0 , (3.111)
where i and j are indices for two different iterations. In each iteration a line search will
Theory be performed eliminating all error in the current search direction via
Manuals
d i + 1 = di + αipi , (3.112)
where α i is given by (cf. (3.93)):

Finite Element
Formulations
T
p iR ( d i )
α i = ---------. (3.113)
T
p i Kp i
Nonlinear
Equation Since the method operates by eliminating all error in each successive search direction,
< Go Back which is orthogonal to all previous directions, this method will yield the exact solution in
n eq iterations in perfect arithmetic. In this sense conjugate gradients can be viewed as a
direct method, although in practice iterations are terminated far before this point, making
the method approximate.
Yet to be discussed is the generation of the orthogonal search directions. This is done by
applying a Gram-Schmidt orthogonalization procedure to a set of linearly independent
vectors, taken in the case of conjugate gradients as residual vectors. The process is begun
by taking the original search direction as the residual, which would be the same direction
taken by steepest descent:
SEACAS
Library p0 = R ( d0 ) . (3.114)
Subsequent search directions are defined via
Theory i–1
Manuals
pi = R ( di) + ∑ β ik p k , (3.115)
k=0
where the coefficients β ik must be found to ensure K-orthogonality of p i with previous

Finite Element
Formulations search directions. This calculation can be done by taking the K-inner product of (3.115)
with p j :
T
Nonlinear 0 = p i Kp j
Equation
i–1
∑
T T
< Go Back = R ( d i ) Kp j + β ik p k Kp j . (3.116)
k=0
T T
= R ( d i ) Kp j + β ij p j Kp j
Eq. (3.116) then allows us to write a formula for β ij
T
R ( d i ) Kp j
β ij = – ------------. (3.117)
T
SEACAS
p j Kp j
Library
At first glance it would appear that to find search direction p i , all β ij would need to be
calculated for j < i . However, we can apply further reasoning to see that, in fact, only one
Theory of the coefficients required by (3.115) is nonzero. To begin, since the error in a given
Manuals
iteration is K-orthogonal to earlier search directions, we can write
T
p j Ke i = 0, j = 0, …, i – 1 , (3.118)
Finite Element
where e i = d i – d . Since Ke i = – R ( d i ) , it follows that
Formulations
T
p j R ( d i ) = 0, j = 0, …, i – 1 . (3.119)
Nonlinear
Equation We can further conclude that since the search direction in a given iteration is a linear
< Go Back combination of all the residuals that preceded it,
T
R ( d j ) R ( d i ) = 0, j = 0, …, i – 1 . (3.120)
We can write the residual at the j + 1 iteration as
R ( d j + 1 ) = – Ke j + 1 = – K ( e j + α j p j )
. (3.121)
= ( R ( d j ) – α j Kp j )
Taking the inner product of this equation with R ( d i ) , we find
T T T
SEACAS α j R ( d i ) Kp j = R ( d i ) R ( d j ) – R ( d i ) R ( d j + 1 ) . (3.122)
Library
Using now the orthogonality property (3.120), we find
Theory  1 T
Manuals T  – ------------
R ( d i ) -R ( d i ), j = i – 1
R ( d i ) Kp j =  α i – 1 . (3.123)

 0, j = 0, …, i – 2
Finite Element
Formulations Thus examining (3.115) and (3.117), we find that only one of the needed coefficients is
non-zero – namely, β i, i – 1 . Naming this quantity β i subsequently we can rewrite
(3.115) as
Nonlinear
Equation
p i = R ( d i) + βip i – 1 , (3.124)
< Go Back where
T
1 R ( di) R ( di)
β i = ------------ . (3.125)
α i – 1 --------------
p
T -
Kp
i–1 i–1
Eq (3.125) can be written even more simply by noting that
T T
α i – 1 p i – 1 Kp i – 1 = p i – 1 K ( d i – d i – 1 )
T
SEACAS = pi – 1( R ( di – 1 ) – R ( di) )
Library
T
= pi – 1R ( di – 1 ) . (3.126)
T
= ( R ( d i – 1 ) + βp i – 2 ) R ( d i – 1 )
Theory
Manuals T
= R ( di – 1 ) R ( di – 1 )
Thus,
Finite Element
Formulations T
R ( di) R ( di)
β i = -------------------. (3.127)
T
R ( di – 1 ) R ( di – 1 )
Nonlinear Eq. (3.127) is the Fletcher-Reeves version of the orthogonalization [Fletcher, R. and
Equation Reeves, C.M., 1964]. Another commonly used form, due to Polak and Ribiere, is trivially
< Go Back obtained from the orthogonality property as
T
R ( di) ( R ( di) – R ( di – 1 ) )
β i = --------------------------- . (3.128)
T
R ( di – 1 ) R ( di – 1 )
Collecting all of these results, the conjugate gradient algorithm for linear problems can be
summarized as
p0 = R ( d o )
T T
pi R ( d i ) R ( di ) R ( di )
SEACAS α i = ---------------------
= ---------------------------------
T T - -
Library pi Kpi pi Kpi
d i + 1 = d i + α i pi
R ( di + 1 ) = R ( di ) – α i Kpi
. (3.129)
Theory  T
R ( di + 1 ) R ( di + 1 )
Manuals  (Fletcher-Reeves)
 -----------------------------------------
T
R ( di ) R ( di )

βi + 1 = 
 R(d T
) ( R ( di + 1 ) – R ( di ) )
 i + 1
(Polak-Ribiere)
Finite Element  ---------------------------------------------
T
Formulations  R ( di ) R ( di )
pi + 1 = R ( d i + 1 ) + β i + 1 pi
With our derivation of the linear conjugate gradient method now complete, we return to
Nonlinear the actual problem of interest here, which is nonlinear conjugate gradients. At this point
Equation
we should distinguish between two alternatives. One could still adopt a Newton-Raphson
< Go Back nonlinear equation-solving strategy and employ a conjugate gradient-type algorithm to
solve the linear system of equations. This type of algorithm is sometimes termed a
Newton-iterative method (see [Kelley, C.T., 1995]). One should note, however, that use of
this method still requires formation of the stiffness matrix, so that even if equation-solving
savings are realized by using CG over a direct solver, the memory requirements of the
method will tend to be extensive. Also if determination of the global Jacobian matrix is
difficult, expensive, or impossible, then this method will be limited just as its more
traditional ancestors are.
Thus one is led to consider the possibility of using conjugate gradient iterations
themselves as the nonlinear solution iterates. The algorithms of this type commonly used
SEACAS look remarkably like the linear algorithm summarized by (3.129). We summarize here
Library perhaps the most successful of these, termed the Polak-Ribiere algorithm after the form of
β , used in the (now approximate) orthogonalization:
ext int
Theory p0 = R ( do) = F –F ( do)
Manuals
α i obtained from nonlinear line search
d i + 1 = d i + αip i
Finite Element
ext int
Formulations R ( di + 1 ) = F –F ( di + 1 ) . (3.130)
T
R ( di + 1 ) ( R ( di + 1 ) – R ( di) )
βi + 1 = ----------------------------------
T
R ( di) R ( di)
Nonlinear
Equation
p i + 1 = R ( d i + 1 ) + βi + 1 p i
< Go Back
Really only one difference from the linear algorithm is obvious. The line search in (3.130)
must now be given by an expression appropriate for nonlinear problems. Although several
alternatives might be devised, one of the simplest (attributed to [Bartels, R. and Daniel,
J.W., 1973]) is generated by considering only the first Newton-Raphson iterate (with
initial guess α i = 0 ) to the line search equation
T
p i( R ( d i + αip i ) ) = 0 , (3.131)
which would lead one to consider
T
SEACAS p iR ( d i )
Library α i = --------- , (3.132)
T
p i Kp i
where K is in practice a diagonalized estimate of the tangent (perhaps a secant) evaluated

Theory at the last iteration.
Manuals
Finite Element
Formulations
Nonlinear
Equation
< Go Back
Preconditioning
SEACAS
Library It is widely recognized that conjugate gradient methods are best behaved when the
condition number of the underlying stiffness matrix is small (i.e., when the eigenvalues are
clustered together). Due to this fact conjugate gradient algorithms are almost never applied
Theory without preconditioning, a term referring to the act of converting an equation system to
Manuals one having the same solution while possessing a tighter clustering of eigenvalues.
Thinking somewhat simplistically about this idea, we might conceive that the ultimate in a
well-conditioned coefficient matrix would be the identity matrix, which has all
Finite Element
Formulations eigenvalues equal to one. Relying once more on a linear system to motivate our ideas, let
us consider that our generic linear system
ext
Kd = F (3.133)
Nonlinear
Equation is ill-conditioned and that we wish to employ a well-behaved CG algorithm to solve it. If
< Go Back we devise a matrix, M , that we feel to be a good approximation to the inverse of K , we
might choose to iteratively solve the equation
ext
MKd = MF , (3.134)
since the matrix MK should have a tight cluster of eigenvalues about one (if M is, indeed, a
–1
good approximation of K ).
Theory Manuals (2/19/99) Finite Element Formulations - Nonlinear Equation Solving - Preconditioning 39
However, CG is only applicable to symmetric systems, and MK is not necessarily
symmetric. If the M we select is symmetric positive-definite, we can find a matrix, E , such
SEACAS that
Library
T –1
EE = M (3.135)
Theory
and consider solution of the system
Manuals
)
–1 –T – 1 ext
E KE d = E F , (3.136)
)
T
Finite Element where d = E d . It is to be noted that this system also will be satisfied by the solution of
Formulations (3.133), while remaining symmetric and positive-definite. If we straightforwardly apply
the linear CG algorithm (3.129) to (3.135), assuming the Polak-Ribiere form, we obtain
Nonlinear
Equation
< Go Back
)
)
)
– 1 ext –1 –T
p0 = R ( do) = E F – E KE d 0
)
T
)
SEACAS R ( d i) R ( d i)
Library α i = -----------------------------------
T –1 –T
)
p i E KE p i
)
d i + 1 = d i + αi p i
Theory . (3.137)
)
)
)
)
T
Manuals
R ( d i + 1 ) = R ( d i ) – α i EKE p i
)
T
)
R ( d i + 1 ) ( R ( d i + 1 ) – R ( d i) )
βi + 1 = -----------------------------------
)
T
)
Finite Element
Formulations R ( d i) R ( d i)
)
)
)
p i + 1 = R ( d i + 1 ) + βi + 1 p i
Algorithm (3.137) can be written in a form not explicitly involving the matrix E or the
Nonlinear
)
)
)
)
–1
)
Equation
supplementary vectors R , p , and d by noting that R ( d i ) = E R ( d i ) ,
< Go Back
–T –1
)
T
M = E E , and by taking p i = E p i . We therefore write the Preconditioned
Conjugate Gradient algorithm for linear systems as:
ext
R ( do) = F – Kd 0
SEACAS
p 0 = MR ( d o )
Library
T
R ( d i ) MR ( d i )
α i = -------------------------------
T
p i Kp i
Theory
Manuals
d i + 1 = di + αipi . (3.138)
R ( d i + 1 ) = R ( d i ) – α i Kp i
T
R ( di + 1 ) M ( R ( di + 1 ) – R ( di) )
Finite Element βi + 1 = -------------------------------
T
R ( d i ) MR ( d i )
Formulations
p i + 1 = MR ( d i + 1 ) + β i + 1 p i
Nonlinear
Finally, we are in a position to discuss a preconditioned matrix-free conjugate gradient
Equation structure for a fully nonlinear system. We might summarize such an algorithm as
< Go Back
ext int
R ( do) = F –F ( do)
SEACAS
p 0 = MR ( d o )
Library
α i obtained from nonlinear line search
d i + 1 = d i + αip i
ext int . (3.139)
Theory R ( di + 1 ) = F –F ( di + 1 )
Manuals
T
R ( di + 1 ) M ( R ( di + 1 ) – R ( di) )
βi + 1 = ----------------------------------
T
R ( d i ) MR ( d i )
Finite Element
Formulations
p i + 1 = MR ( d i + 1 ) + β i + 1 p i
Following the lead of the last section, one alternative for the line search would be the first
iteration for a Newton-Raphson strategy to find α i :
Nonlinear
Equation
T –1
< Go Back piM R ( di)
α i = -------------
. (3.140)
T
p i Kp i
To conclude, it should be remarked that the simplest choice for preconditioning, Jacobi
preconditioning, is accomplished when the matrix M is chosen as diagonal. One choice
might be to take M to be the inverse of the diagonal of the stiffness matrix.
Revisited
Basics of Element Design

Theory Introduction
Manuals
Convergence
Parameterization
Finite Element
Formulations
Shape Functions
Quadrature
< Go Back Local Arrays
i Blue text
indicates
a link to more
information.
Theory Manuals (2/19/99) Finite Element Formulations - Basics of Element Design 1

Introduction
SEACAS
Library In this section we explore the basic issues associated with the design of finite elements,
which are the building blocks of the methods we have discussed. In particular we will
discuss how definitions and manipulations are done at the local level to produce the
Theory e inte e
Manuals
elemental quantities like m , f , and k that are needed for assembly of the global
equations of motion. We concentrate in this section on one field problems (i.e., where only
the deformation mapping ϕ t is discretized). It will turn out that many nonlinear solid
Finite Element mechanics applications of interest, including nearly incompressible elasticity and metal
Formulations plasticity, require more sophisticated approximations in which other variables (like
pressure) must be explicitly included in the formulation. Discussion of such advanced
methods is deferred to Advanced Element Design Issues.
Basics of
Element Design
< Go Back
Theory Manuals (2/19/99) Finite Element Formulations - Basics of Element Design - Introduction 2
Convergence
SEACAS
Library Before introducing in detail the manipulations necessary at the element level, it is
worthwhile to discuss in general terms the general requirements usually placed upon
shape function definitions. It should be noted that these conditions are sufficient but not
Theory necessary, so that many formulations exist that violate one or more of them. However, it is
Manuals also fair to say that most finite elements in wide use satisfy the conditions we will now
place. The discussion we give now is brief and rather qualitative; the interested reader
should consult [Hughes, T.J.R., 1987] or [Strang, G. and Fix, G.J., 1973] for more
Finite Element technical discussions of these points (notably within the context of linear problems).
Formulations
We begin by defining m , which will denote the highest order of shape function spatial
derivative present in the expression for the stiffness matrix. For the class of problems we
have considered in this text, we recall from Newton Raphson Framework that the
Basics of element stiffness matrix takes the form
Element Design
< Go Back int

i ∂f p i
k pq  d n + 1  de  .
e e
= (3.141)
  e  n + 1
∂dq
The internal force vector required in (3.141) was given generically in Localization and
Assembly as
Theory Manuals (2/19/99) Finite Element Formulations - Basics of Element Design - Convergence 3
3
int e
∑
–1 h
f p = ∫ N a, j ( ϕ t ( x ) )T ij dv . (3.142)
h e
SEACAS ϕt( Ω ) j = 1
Library
Performing the differentiation indicated in (3.141) will produce no higher than first-order
derivatives of the shape functions; therefore, m = 1 .
Theory
Manuals The three general convergence requirements we wish to mention are as follows.
• The global shape functions N A should have global continuity of the order m – 1 . In
m–1 h
mathematical terms they should be C on Ω .
Finite Element
Formulations
• The restriction of the global shape functions to individual elements (i.e., the { N a } )
m
should be C on element interiors.
Basics of
• The elemental shape functions { N a } should be complete.
Element Design
m–1
< Go Back The first two of these requirements are fairly simple to understand. The first, the C
continuity requirement, simply means that all derivatives up to m – 1 of the shape
functions should not undergo jumps as element boundaries are crossed. In the current case
0
this means that all N A should be C . Since the approximation to the configuration
h
mapping ϕ t is a linear combination of these shape functions, we see that the physical
restriction placed by this condition amounts to no more than a requirement that the
displacement be single-valued throughout the domain (i.e., gaps and interpenetrations at
element boundaries may not occur).
SEACAS
Library
The second requirement on element interiors simply states that the shape functions should
be sufficiently smooth so that the element stiffness expression is integrable. Physically
speaking, the first derivatives of the configuration mapping produces strain measures, so
we simply require that the strains be well-behaved on element interiors by this restriction.
Theory
Manuals Note that global smoothness of the strains (and, therefore, stresses) is not required. This
point is of some importance in the reporting of results.
The third requirement, the completeness requirement, is somewhat more involved to
Finite Element explain and yet corresponds fairly directly to physical ideas. We say that a given element is
Formulations
complete when setting the element degrees of freedom according to a given low-order
h
polynomial forces the solution (in this case ϕ t ) to be interpolated according to the same
polynomial pointwise within the element. The degree of polynomials for which we place
Basics of
Element Design this requirement is referred to as the degree of completeness for the element.
< Go Back In the current case where we deal with solid continua, the usual degree of completeness
demanded is 1. This means that all polynomials, up to and including order 1, should be
exactly interpolated by the element. It is worthwhile to consider an example of this point.
Suppose we are in three dimensions, and set the element degrees of freedom via
e e e e
d a = c 0 + c 1 X ae x + c 2 Y ae y + c 3 Z ae z , (3.143)
e e e
where c 0 – c 3 are arbitrary constants and X a, Y a, Z a are the (reference) coordinates for
local node number a . The completeness condition requires that
SEACAS
Library
nen
∑
h e e e
ϕt ( X ) = N a ( X )d a
a=1 (3.144)
Theory e e e
Manuals = ( c0 + c1X ex + c2Y ey + c3Z ez)
e e
hold for all X ∈ Ω and for all values of the arbitrary constants.
Finite Element
Formulations As mentioned above, this requirement has a physical interpretation as well. In solid
mechanics we have already pointed out that the first spatial derivatives of the
displacements produce strains. Since we require that an element be able to reproduce
arbitrary linear polynomials, this also implies that any state where the first derivatives (i.e.,
Basics of strains) are constant should be exactly representable. Thus a complete element should be
Element Design
able to exactly represent any uniform strain state. A practical way to test for this condition
< Go Back is to impose a boundary value problem on an arbitrary patch of elements having a constant
strain (and thus stress) solution and then to demand exactness of the numerical solution.
Such a test is called a “patch test” and has become one of the standard benchmarks by
which any new proposed element formulation is tested and evaluated.
Parameterization
SEACAS
Library With the three criteria in hand for element definitions, we proceed to define a recipe
through which element definitions and manipulations can be systematically performed.
The most basic definition to be made toward this end is the concept of the local (or parent)
Theory parameterization of an element. In effect we seek to define a local coordinate system that
Manuals will be the same for every element in a problem, contributing in great part to the
modularity we will desire for element level operations.
We will denote a vector of these local variables by r , with r being a 2-vector in two
Finite Element
Formulations dimensions and a 3-vector in three dimensions. Specifically we define r as
 r (two dimensions)

Basics of  s
Element Design 
r =  . (3.145)
< Go Back  r
 s (three dimensions)

 t
The local variables r , s , and t are all assumed to range between -1 and 1, so that the
domain definition is likewise standardized among all elements in a given problem. The
domain of r is often referred to as the parent domain. In two dimensions it is a biunit
square, and in three dimensions it is a biunit cube (see Figure 3.6).
Theory Manuals (2/19/99) Finite Element Formulations - Basics of Element Design - Parameterization 7
s e
X (r) e
SEACAS
Library
Ω
r
e
Theory
s X (r) e
Manuals
Ω
r
Finite Element t
Formulations
Figure 3.6 Local parameterizations and coordinate mappings in two and three dimensions.
Of course, for this alternative element coordinate system to be of any use, its relationship
Basics of
with the global coordinate system must be defined. This is accomplished through a shape
Element Design function expansion via
< Go Back
nen
∑
e e
X (r) = Ñ a ( r )X a , (3.146)
a=1
e e
where X is the global (reference) coordinate mapping covering element e and where X a
are the element nodal (reference) coordinates, as before. Note also that in (3.146) the
shape functions have been written using the parent coordinates as the independent
variables. This is the reason for the superposed tilde on the shape function. One could
e
think of r as an material point label within the element, so that X and r are two
reference coordinate systems for the element that are related according to (3.146).
SEACAS
Library The most important generic class of finite elements is comprised of so-called
isoparametric elements. Such elements are defined by utilizing the same shape functions
h e e
for definition of ϕ t ( X ) (see Eq. (3.144)) as for the element coordinates X (as in
Theory
Manuals (3.146)). One can show (and the reader should) that providing all element shape functions
sum to one at any point in the element, an isoparametric element automatically satisfies
the completeness condition outlined in Convergence. Thus provided we choose shape
Finite Element
functions that sum to one, are suitably smooth on element interiors, and match
Formulations neighboring element descriptions on element boundaries, the convergence criteria are
automatically satisfied. We will concentrate on the mechanics of shape function definition
in the next section.
Basics of In the meantime let us consider the implications of the isoparametric approach for the
Element Design Lagrangian description of large deformation solid mechanics we now consider. So that we
< Go Back may distinguish carefully between mappings taking r as an argument and between those
taking X , we will use superposed tildes for the former (as in Eq. (3.146)). If an element is
isoparametric, then by definition the configuration mapping over an element is given by
nen
∑
h e
ϕ̃ t ( r ) = Ñ a ( r )d a , (3.147)
a=1
where the shape functions Ñ a ( r ) are exactly the same as in (3.146). However, it should
h
SEACAS
also be the case that the function ϕ̃ t ( r ) should be attainable from the composition of
Library
h e e
ϕ t ( X ) (defined according to (3.144)) with X ( r ) (defined according to (3.146)). Thus
we can write
Theory nen nen
Manuals
∑ ∑
e h e e
Ñ a ( r )d a = ϕ̃ t ( r ) = N a ( X ( r ) )d a . (3.148)
a=1 a=1
Finite Element Comparing the leftmost and rightmost expressions of (3.148) and realizing that the
Formulations e
equality must hold for any given combination of the element degrees of freedom d a , we
e
are led to conclude that the alternative shape function expressions Ñ a ( r ) and N a ( X )
Basics of
must be related by mere composition via
Element Design
e
< Go Back Ñ a = N a ° X . (3.149)
Thus we have the option of defining the shape functions over whatever domain is
convenient, and since the parent domain is the one that is standardized, we typically begin
with an expression for Ñ a and then derive the implied expression for N a according to
–1
e
N a = Ñ a ° X . (3.150)
Equation (3.150) has important implications in practice for, in general, we have no
–1
e e
guarantee that the inverse mapping X of X is well behaved, which it must be for the
SEACAS
Library shape functions N a to make sense. Fortunately, according to the implicit function
theorem, the inverse function to (3.146) is smooth and one-to-one, provided the Jacobian
of the indicated transformation is nonzero. This essentially amounts to a geometric
Theory restriction on elements in the reference domain. In two dimensions, using a four-noded
Manuals
element, the implication is that all interior angles in each element must be less than 180°
(see Figure 3.7).
Finite Element
Formulations
Basics of all interior angles ≤ 180°

Element Design
< Go Back
Figure 3.7 Geometric restrictions on a four-noded element to retain well-posedness of the

coordinate and configuration mappings.
)
Finally, let us introduce the notation N a for shape functions that take the current
e h e
coordinates x = ϕ t ( X ) as arguments. Such an expression is needed – for example, in
∂
)
Eq. (3.142) (note the abuse in notation) – where the spatial derivatives N a , i = eN a
∂xi
SEACAS
Library must be computed. Following similar reasoning to the above, one can conclude that the
)
functions N a must obey
–1
)
h
a ° ϕ̃ t
Theory
Manuals
Na = N . (3.151)
–1
h
Again for the needed function ϕ̃ t to be well-behaved, the Jacobian of the transformation
Finite Element (3.147) must be nonzero. This restriction amounts to:
Formulations
h h
∂ϕ̃ t ∂ϕ̃ t ∂X
e
det = det det ≠ 0. (3.152)
∂r ∂X
e ∂ r
Basics of
Element Design
Provided the original element definitions are not overly distorted, the second term on the
< Go Back right-hand side of (3.152) will be nonzero. Thus the well-posedness of the spatial shape
h
∂ϕ̃ t
)
functions N a requires that det e

be nonzero. The reader will recognize this as the
∂X
approximated determinant of the deformation gradient J , as defined in Measures of
Deformation. According to Eq. (2.10) J must be positive pointwise for the concept of
volume change to have any physical meaning. Thus provided the approximated
deformation mapping remains kinematically admissible (i.e., J > 0 ), the spatially defined
shape functions are guaranteed to be well-behaved.
SEACAS
Library With these arguments as background, we now turn our attention to definition of the N a ,
according to the parent domain. To keep notational complexity to a minimum, we will
)
drop the explicit distinction between N a , Ñ a , and N a , referring to all these objects as
Theory
Manuals simply N a in the sequel.
Finite Element
Formulations
Basics of
Element Design
< Go Back
Shape Functions
SEACAS
Library Most continuum-based finite elements rely on Lagrange polynomials for their shape
function definitions. Beginning with expressions appropriate for one-dimensional
domains, let us suppose that we have a one-dimensional element with nen nodes, which
Theory are equally spaced over the [ – 1, 1 ] domain. Use of nen – 1 order Lagrange polynomials
Manuals
nen – 1
La ( r ) for definition of each of the element shape functions N a leads to
nen
∏ ( r – rb)
Finite Element
Formulations
nen – 1 b = 1, b ≠ a
N a( r ) = L a ( r ) = -------------------
nen
, (3.153)
∏ ( ra – rb)
b = 1, b ≠ a
Basics of
Element Design
where the r b refer to the local (parent) coordinates of the individual element nodes. The
< Go Back
reader may care to verify that these shape functions have two useful properties. First, that
N a ( r b ) = δ ab (3.154)
and second, that
nen
∑ N a ( r ) = 1 for all r ∈ [ – 1, 1 ] . (3.155)
a=1
Theory Manuals (2/19/99) Finite Element Formulations - Basics of Element Design - Shape Functions 14
Equation (3.155) is noteworthy in that it provides an important ingredient of the
completeness argument (see Parameterization), whereas Eq. (3.154) ensures that the
SEACAS nodal degrees of freedom have the actual interpretation that
Library
e h e e
d a = ϕt( X a) = x a (3.156)
Theory
(see Eq. (3.144)).
Manuals
Before proceeding to the more interesting multidimensional case, some examples may be
useful. For nen = 2 the element coordinates are r 1 = – 1 and r 2 = 1 . The
1
Finite Element
corresponding element shape functions are computed from (3.153) as N 1 ( r ) = --(-1 – r )
Formulations 2
1
and N 2 ( r ) = --(-1 + r ) , thereby providing the basis for the one-dimensional linear finite
2
element. For nen = 3 the local nodal coordinates are r 1 = – 1 , r 2 = 0 , and r 2 = 1 ,
Basics of
Element Design
1 2
< Go Back with the shape functions turning out to be N 1 ( r ) = --r-( r – 1 ) , N 2 ( r ) = 1 – r , and
2
1
N 3 ( r ) = --r-( r + 1 ) , thereby defining the one-dimensional quadratic element. These
2
shape functions are plotted in Figure 3.8.
1 N1 1 N3
N2
N1 N2
SEACAS
Library
r1 = –1 r2 = 1
r1 = – 1 r2 = 0 r3 = 1
Theory
Manuals Figure 3.8 Low-order, one-dimensional Lagrange shape functions.
Next we turn to the generalization of these concepts to two and three dimensions. This can
be accomplished by building up “products” of one-dimensional shape functions, as
Finite Element
Formulations indicated schematically for the four-noded quadrilateral element depicted in Figure 3.9. In
general, let us suppose that local Node a in a two-dimensional element has local
coordinates ( r c, s d ) , where indices c and d refer to the node number in the r and s
directions, respectively. The two-dimensional shape function is given by
Basics of
Element Design
nenr – 1 nens – 1
< Go Back N a ( r, s ) = Lc (r) × Ld (s) , (3.157)
where nen r and nen s are the numbers of nodes in the r and s directions, respectively.
Taking the four-noded quadrilateral depicted in Figure 3.9 as an example, the shape
1 1
functions N a are found to be N 1 = --(-1 – r ) ( 1 – s ) , N 2 = --(-1 + r ) ( 1 – s ) ,
4 4
1 1
N 3 = --(-1 + r ) ( 1 + s ) , and N 4 = --(-1 – r ) ( 1 + s ) .
4 4
=
s 4 s 3
2
SEACAS ×
Library 1 2 r r
1 1 2
Theory Figure 3.9 Definition of element shape functions for two-dimensional, four-noded quadrilateral.
Manuals
The three-dimensional case can be treated analogously, and in doing so we use the trilinear
brick element depicted in Figure 3.10 as a template. Here we consider that local Node a
Finite Element has local coordinates ( r c, s d, t e ) and write the three-dimensional shape functions as:
Formulations
nenr – 1 nens – 1 nent – 1
N a ( r, s ) = Lc (r) × Ld (s) × Le (t) . (3.158)
Basics of
The appropriate shape functions for the trilinear brick turn out to be:
Element Design 1 1 1
N 1 = --(-1 – r ) ( 1 – s ) ( 1 – t ) , N 2 = --(-1 + r ) ( 1 – s ) ( 1 – t ) , N 3 = --(-1 + r ) ( 1 + s ) ( 1 – t ) ,
4 4 4
< Go Back
1 1 1
N 4 = --(-1 – r ) ( 1 + s ) ( 1 – t ) , N 5 = --(-1 – r ) ( 1 – s ) ( 1 + t ) , N 6 = --(-1 + r ) ( 1 – s ) ( 1 + t ) ,
4 4 4
1 1
N 7 = --(-1 + r ) ( 1 + s ) ( 1 + t ) , and N 8 = --(-1 – r ) ( 1 + s ) ( 1 + t ) .
4 4
s 4 3
s
2
1 8 7
SEACAS
Library
× × =
1 2 r 1 r2
1 2 t
t 5 6
Theory Figure 3.10 Definition of element shape functions for three-dimensional, eight-noded brick.
Manuals
Finite Element
Formulations
Basics of
Element Design
< Go Back
Quadrature
SEACAS
Library With the element parameterizations and shape functions now defined, we are in a position
e
to discuss how element-level calculations are performed to evaluate such quantities as m ,
inte e
Theory f , and k . Notably all these calculations involve integrals over the element domain,
Manuals
inte
as evidenced by the expression for f given in (3.142). One can evaluate these
integrals analytically only for a few highly specialized cases, meaning that numerical
integration (i.e., quadrature) is required for any type of generality to be present in the
Finite Element
Formulations element formulation. Accordingly, we are led in this section to consider the generic
problem of integrating a function, f , over the element domain via
e
Basics of
∫ f ( x ) dv , (3.159)
h e
Element Design ϕt( Ω )
< Go Back
where f could, in principle, be scalar, vector, or tensor-valued.
The first step in evaluation of (3.159) is generally to perform a change of variables,
h e
converting the integral in the current element physical domain ϕ t ( Ω ) to one over the
parent domain, which we shall denote by . This is accomplished using the standard
change-of-variables formula from multivariate calculus,
Theory Manuals (2/19/99) Finite Element Formulations - Basics of Element Design - Quadrature 19
e e
∫ f ( x ) dv = ∫ f ( x ( r ) )j ( r ) d , (3.160)
h e
ϕt( Ω )
SEACAS
Library
where j ( r ) is the Jacobian of the transformation from parent coordinates r to spatial
e
coordinates x :
Theory
e
Manuals ∂x
j ( r ) = det , (3.161)
∂r
e
Finite Element where x is as given in (3.147).
Formulations
The advantage of (3.160) over (3.159) is that the integration takes place over a
standardized domain, for which quadrature rules are readily tabulated. One typically
approximates the integral in (3.160) by applying quadrature via
Basics of
Element Design nint
∑
e e
< Go Back ∫ f ( x ( r ) )j ( r ) d ≈ f ( x ( r l ) )j ( r l )W l , (3.162)
l=1
where nint is the number of integration (quadrature) points in the element, r l is the
parent coordinate of quadrature Point l , and W l is the weight associated with quadrature
Point l . The choice of these quadrature point coordinates and weights effectively defines
the numerical integration scheme and the accuracy associated with it.
The most prevalent quadrature schemes in finite elements are based on Gaussian
quadrature rules, which may be derived in terms of Legendre polynomials. While this
SEACAS derivation is unnecessary for our present purposes, its result is that in one dimension the
Library Gaussian quadrature rules are optimally accurate in that no greater accuracy can be
achieved for lesser cost. By cost we mean the number of integration points used, whereas
by accuracy we mean the lowest order polynomial not integrated exactly by a given
Theory quadrature rule. The Gaussian rules in one dimension have the property that given nint
Manuals
integration points, 2 × nint order accuracy is achieved, meaning that a 2 × nint – 1
order polynomial in r will be integrated exactly. Below are listed the first few Gauss
integration rules over the domain [ – 1, 1 ] .
Finite Element
Formulations
• nint = 1 : r 1 = 0 , W 1 = 2 (second order accurate).
1 1
• nint = 2 : r 1 = – ---, r 2 = ---, W 1 = W 2 = 1 (fourth order accurate).
3 3
Basics of
Element Design
3 3 5 8
< Go Back • nint = 3 : r 1 = – --
5 , r2 = 0 , r3 = --,
5 1 W = W 3 = , W = (sixth order
9 2
-- --
9
accurate).
Now returning to the problem of interest, multidimensional quadrature, we can use very
similar reasoning to that used to define multidimensional shape functions in the last
section. Since integration in a multidimensional domain involves integrating with respect
to each variable separately while holding the others constant, we might expect that
numerical integration in successive directions is done in just the same way. The result of
this fact is that we can define multidimensional quadrature rules as products of one-
dimensional ones, just as was done for shape function definitions. It turns out the
optimality property present in the one-dimensional formulation is lost but that a highly
SEACAS systematic and effective integration procedure results.
Library
Beginning in two dimensions we refer to Figure 3.11, which depicts a four-point
quadrature scheme in two dimensions that we will use as a template. Let us consider
quadrature Point l , having local coordinates ( r l , s l ) , where indices l r and l s refer
Theory r s
Manuals
to the appropriate quadrature point number in the r and s directions, respectively. The
two-dimensional weight is given by the product of the appropriate weights from the one-
dimensional rules, i.e.:
Finite Element
Formulations
Wl = Wl × Wl . (3.163)
r s
s s
Basics of × × ×
Element Design
r×l × × =
< Go Back r × × r
r × sl
s
( rl , sl )
r s
Figure 3.11 Quadrature rule definition in two dimensions: four-point Gaussian quadrature.
Taking the two-dimensional, four-point quadrature rule depicted in Figure 3.11 as an

example, the appropriate parameters are found to be: W 1 = W 2 = W 3 = W 4 = 1 ,
( r 1, s 1 ) =  – ------
1 1  1 , – 1  , ( r , s ) =  1 , 1  , and
, – ----, ( r , s ) =
 3 -3 2 2  -------
3 ----
3
3 3  ------3-
3 ----
SEACAS
( r 4, s 4 ) =  – ------
1 1
Library , - .
 3 ---- 3
Three-dimensional quadrature rules are similarly conceived; the reader should consult
Theory
Manuals
Figure 3.12 for a template of the procedure. We consider quadrature Point l , having local
coordinates ( r l , s l , t l ) , where indices l r , l s , and l t refer to the appropriate
r s t
quadrature point number in the r , s , and t directions, respectively. The three-

Finite Element dimensional weight is given by the product of the appropriate weights from the one-
Formulations
dimensional rules, i.e.:
Wl = Wl × Wl × Wl . (3.164)
r s t
Basics of
Element Design s s
× ×
< Go Back ×
× × × × ××
× × = ×r
rl r
sl × ×t
× t×
r
s lt
( rl , sl , tl )
r s t
Figure 3.12 Quadrature rule definition in three dimensions: eight-point Gaussian quadrature.
Considering the case in Figure 3.12 as a specific example, we find the following
( r1, s1, t1 ) =  – ------
1
parameters: W 1 – W 8 = 1 , 1 1
, – - , – ----
 3 ------ 
, ( r2, s2, t2 ) =  ------
1 1
, – ------
1
, – ---- ,
3
- 3 - -
SEACAS
3 3 3
Library
( r3, s3, t3 ) =  ------
1
1 1
, -, – ----
 3 ------ 
, ( r4, s4, t4 ) =  – ------
1 1
, ------
1
, – ----

, ( r5, s5, t5 ) =  – ------
1 1 1
, – ------
, ----,
3 3-
3 - - - -
3 3 3 3 3
( r6, s6, t6 ) =  ------

1 1
1
, – - , ----
 3 ------ 
, ( r7, s7, t7 ) =  ------
1 1 1
, -, ----, and ( r8, s8, t8 ) =  – ------
1 1 1
, -, ----.
3 3 - 3 3-
3 ------  3 ------
3 3-
Theory
Manuals
Finite Element
Formulations
Basics of
Element Design
< Go Back
Local Arrays
SEACAS
Library The final task in the section is to give brief prescriptions for how element-level
e inte e
calculations are done to find m , f , and k for a given element, e . These quantities
int
Theory are needed by the global assembly algorithm to form M and F (as discussed in
Manuals
Localization and Assembly) and K (as discussed in Newton Raphson Framework).
e
Beginning first with m we recall the general expression for the element mass matrix
Finite Element e –1 –1
Formulations m pq = ∫ ρN a ( ϕ t ( x ) )δ ij N b ( ϕ t ( x ) ) dv . (3.165)
h e
ϕt( Ω )
One can apply a change of variables to the reference configuration to find that
Basics of
Element Design e e e
< Go Back
m pq = ∫ ρN a ( X )δ ij N b ( X )J dV
e
Ω
, (3.166)
e e
= ∫ ρ 0 N a ( X )δ ij N b ( X ) dV
e
Ω
where the second line of (3.166) holds because ρ 0 = Jρ (conservation of mass, see Eq.
(2.74)). Looking at the second line of (3.166), we see that the element mass matrix is
independent of the deformation. It is straightforwardly calculated using quadrature via
Theory Manuals (2/19/99) Finite Element Formulations - Basics of Element Design - Local Arrays 25
nint
∑
e
m pq ≈ ρ 0 ( r l )N a ( r l )δ ij N b ( r l )j 0 ( r l )W l , (3.167)
SEACAS l=1
Library
where
e
∂X
j 0 ( r l ) = det . (3.168)
Theory
Manuals
∂r
According to the discussion of the last section, the ordinary strategy in applying
quadrature would be to use a sufficiently accurate rule so that (3.167) is evaluated exactly
Finite Element (at least if the reference density ρ 0 is constant). This would lead one to employ a four-
Formulations
point Gauss quadrature rule for a four-noded quadrilateral in two dimensions and an eight-
point Gauss rule for an eight noded brick in three dimensions. Following this procedure
e
produces a “consistent” mass matrix m .
Basics of
Element Design The difficulty with a consistent mass matrix, however, is twofold. First, it is, in general,
< Go Back banded but not diagonal (as would be preferable for an explicit dynamics application, for
example), and second, experience shows that better accuracy is often exhibited in
dynamics problems with “lumped mass”, where the rows of (3.168) are actually summed
and the result placed on the diagonal of the mass matrix. Use of this row sum technique
produces the following alternative, more widely used expression for element mass:
nint
∑
e
m pq ≈ δ ij δ ab ρ 0 ( r l )N a ( r l )j 0 ( r l )W l . (3.169)
lumped
l=1
inte
Turning attention now to f , we begin by applying the change of variables to (3.142):
SEACAS nint
int e
∑
h
p≈ N a, j ( r l )T ij j 0 ( r l )W l .
Library
f (3.170)
l=1
h
Theory Two requirements of (3.170) are notable: the determination of the stress T ij (dependent,
Manuals
as indicated on the current element deformation field through the constitutive law), and the
need for the spatial derivatives N a, j of the shape functions. In fact, derivatives of this type
are also needed for the stress calculation, which will ordinarily involve the approximated
Finite Element
Formulations deformation gradient:
h
∂ϕ t nen ∂N a
∑
h e
F = e
= d .
e a
(3.171)
Basics of ∂X a = 1 ∂X
Element Design
< Go Back ∂N a ∂N a inte

Thus calculation of both e
and e
is typically necessary to obtain f .These
∂X ∂x
derivatives are usually produced by a shape function subroutine, called by the element
∂N a
subroutine for each quadrature point. Taking the spatial derivative e
as an example, the
∂x
chain rule can be invoked to obtain the appropriate expressions via:
      e –1
 ∂N a   ∂N a  ∂r k  ∂N a  ∂x j
  =   =   ∂r , (3.172)
 ∂rk  ∂rk
e
SEACAS
 ∂ x ej   ∂xj  k
Library
     
where ( r k ) = ( r, s ) in two dimensions, and ( r k ) = ( r, s, t ) in three dimensions. The
∂N a
Theory
Manuals
reader will recognize that can be computed through simple differentiation of the local
∂rk
e
∂x j
Finite Element shape functions, whereas is found by differentiation of (3.147). Calculation of the
Formulations ∂rk
required inverse is rather simple and is readily done in closed form, since the matrix
involved is 2x2 in two dimensions and 3x3 in three dimensions.
Basics of e
Element Design
Completing our discussion of element-level calculations, the element stiffness matrix k
is given generically by
< Go Back
int
i ∂f p i
k pq  d  =  de  ,
e e
(3.173)
  e  
∂dq
where the internal force vector is as given in (3.142). Some manipulation of (3.142) is
useful at this point:
int e e –1 h
f p = ∫ [ N a, J ( X )F Jj T ij ]J dV
e
Ω
SEACAS , (3.174)
e h
∫ [ N a, J ( X )P iJ ] dV
Library
=
e
Ω
h
Theory where in the second line of (3.174), the Piola stress P iJ has been introduced in
Manuals
accordance with (2.54). There are many ways in which the derivative indicated in (3.173)
can be expressed; here we do it by computing the derivative of the Piola stress with respect
to the deformation gradient (Eq. (3.171)) and invoking the chain rule:
Finite Element
Formulations
i ∂
k pq  d  =
e e e h
  ∂dq
e ∫ [ N a, J ( X )P iJ ] dV
e
Ω
h
. (3.175)
Basics of
e ∂P iJ ∂F kL
Element Design = ∫ N a, J ( X ) h e
dV
< Go Back Ω
e ∂ F kL ∂ d q
Simplification of the second line of (3.175) results in the following expression:
i
k pq  d  =
e e e h e
  ∫ N a, J ( X )C iJjL ( ϕ t )N b, L ( X ) dV , (3.176)
e
Ω
where the material moduli C iJjL are defined as
h
∂P iJ
C iJjL = . (3.177)
h
SEACAS
∂ F jL
Library
Application of the quadrature rule to (3.176) gives
i nint
k pq  d ≈
∑
e e
Theory N a, J ( r l )C iJjL ( r l )N b, L ( r l )j 0 ( r l )W l . (3.178)
Manuals  
l=1
The required reference coordinate shape function derivatives can be calculated as

Finite Element
discussed above for each quadrature point, as can the Jacobian j 0 . The material moduli
Formulations
C iJjL are, typically, the most difficult to compute, as they require linearization of the
tensor-valued constitutive relation with respect to a tensor-valued strain measure (in this
case the deformation gradient). It should be noted that (3.178) is given for illustrative
Basics of purposes only; the stress and strain measures conveniently utilized in the linearization vary
Element Design
widely, depending on the constitutive relation used. It should be noted, however, that
< Go Back provided the moduli are symmetric (in the major sense), then the element stiffness matrix
e
k will be as well.
Revisited
Advanced Element Design Issues

Theory Introduction
Manuals
Constrained Media and Locking
Selective/Fully Reduced Integration
Finite Element
Formulations
Hourglass Control
< Go Back
Blue text
i indicates
a link to more
information.
Theory Manuals (2/19/99) Finite Element Formulations - Advanced Element Design Issues 1
Introduction
SEACAS
Library In this section we discuss some advanced element design issues having particular
relevance to large deformation problems featuring inelastic response. We begin the
discussion with a specific example of how the standard element formulations discussed in
Theory Basics of Element Design can have difficulty in problems featuring near or complete
Manuals incompressibility, as is common, for example, in computational plasticity. Some basic
remedies for this situation are then discussed.
Finite Element
Formulations
Advanced
Element Design
< Go Back
Theory Manuals (2/19/99) Finite Element Formulations - Advanced Element Design Issues - Introduction 2
Constrained Media and Locking
SEACAS
Library The incompressibility dilemma can be motivated fairly simply by considering linear
elastic, isotropic behavior. Returning once more to the discussion in Linear Elastic IBVP,
we consider a stress strain relation of the form
Theory
Manuals T ij = C ijkl E kl , (3.179)
where the material moduli C ijkl are of the form
C ijkl = λδ ij δ kl + µ [ δ ik δ jl + δ il δ jk ] .
Finite Element
Formulations (3.180)
Plugging (3.180) into (3.179), one obtains
T ij = 2µu ( i, j ) + λδ ij u k, k . (3.181)
Advanced
Element Design
We are most interested in the volumetric response of this material; accordingly, let us
< Go Back
define the hydrostatic pressure p as
1
– p = --- T k, k (3.182)
3
and the dilitation (volume change) Θ as
Θ = u k, k . (3.183)
Theory Manuals (2/19/99) Finite Element Formulations - Advanced Element Design Issues - Constrained Media and Locking 3
Computing according to (3.181), we find that for an isotropic material,
1
SEACAS – p = --- ( 2µ + 3λ )Θ . (3.184)
Library 3
The coefficient relating – p and Θ is ordinarily called the bulk modulus K :
Theory K = 2µ + 3λ (3.185)
Manuals
and corresponds physically to the volumetric stiffness of the material. Recalling that µ is
the shear modulus, representing the resistance of the material to shearing motions, it
Finite Element proves useful to examine the ratio of K to µ as an indicator of the degree of
Formulations incompressibility of the material. Using the relationships between the Lamé parameters
and the more familiar elastic modulus and Poisson’s ratio (see (1.64) and (1.65)), we find
this ratio can be written solely in terms of the Poisson’s ratio:
Advanced K 2(1 + ν)
Element Design --- = ----------------------- . (3.186)
µ 3 ( 1 – 2ν )
< Go Back
1
Recalling that the thermodynamically admissible values for ν range between – 1 and --- ,
2
1
we see that the case where ν approaches --- from below causes (3.186) to grow without
2
1
bound, so that the bulk modulus becomes infinitely large when ν = --- . In this case the
2
volume change Θ is constrained to be zero pointwise in the medium, and the material is
said to be incompressible.
SEACAS
Library Let us now consider the behavior of a finite element discretization of the linear boundary
1
value problem described in Linear Elastic IBVP, where ν = --- . We consider the mesh
2
Theory
shown in Figure 3.13, comprised of linear triangles. This is an element not discussed to
Manuals this point but can be obtained formally by consideration of the four-noded quadrilateral
discussed in Basics of Element Design, with two of the nodal coordinates set to be the
same. Thus the displacement field is linearly interpolated with the result, in the case of
Finite Element triangles, that the strains are constant throughout the element.
Formulations
Advanced
Element Design II A
< Go Back I
Figure 3.13 Sample mesh illustrating mesh locking for the incompressible case.
Since the strains are constant in the element, the requirement that
h
u k, k = 0 (3.187)
pointwise causes the total volume change in the element to be zero also:
dA h
SEACAS
Library
------- =
A ∫ u k, k dA = 0 , (3.188)
A
meaning simply that each element in the mesh may not change area due to the
incompressibility constraint.
Theory
Manuals
Examining now the behavior of Element I in Figure 3.13, the constant area constraint
implies that Node A can only move in the horizontal direction, since the two lower nodes
of Element I are fixed. However, Element II places the restriction that Node A can only
Finite Element move vertically. Taken together, the isochoric constraint in each element prevents Node A
Formulations
from moving at all, in any direction. This argument can be repeated throughout the mesh,
h
to conclude that no node can move at all, so that u = 0 is the solution to the discrete
problem. However, it is clear that in the physical situation, the fact that the material is
Advanced incompressible does not preclude all deformation. Thus the finite element solution
Element Design
produces a solution that is nearly nonphysical because of the fact that the numerical
< Go Back approximation of the incompressibility condition overconstrains the numerical
representation of the physical system.
The phenomenon described for this admittedly specialized system is referred to, generally
h
speaking, as “mesh locking”. It will not, of course, always be the case that u = 0 for
every boundary value problem, but it does turn out that fully integrated elements of the
type discussed in Basics of Element Design will, in general, produce excessively stiff
solutions when the material is either nearly or completely incompressible. We are thus led
to consider techniques where the amount of constraint placed by the approximation of the
volumetric material response can be relaxed when appropriate.
SEACAS
Library
Theory
Manuals
Finite Element
Formulations
Advanced
Element Design
< Go Back
Selective/Fully Reduced Integration
SEACAS
Library One of the simplest techniques used to eliminate element locking is to deliberately
inte e
underintegrate the internal force vector f (and the element stiffness k in the case of
a quasistatic or implicit dynamic calculation). Selective reduced integration means that
Theory
Manuals only the troublesome volumetric terms are underintegrated, whereas fully reduced
integration means that all terms are underintegrated. The latter option is particularly
attractive in explicit dynamic and matrix-free quasistatic calculations, since the element
level calculations comprise a large proportion of total solution costs in these cases. Since
Finite Element
Formulations the cost of element calculations is directly proportional to the number of quadrature
points, fully reduced integration becomes very attractive when speed is of special concern.
Let us consider the case of the eight-noded hexahedron in three dimensions as an example,
Advanced inte
Element Design
and apply fully reduced integration in the calculation of f . The ordinary quadrature
rule for this element would be eight-point Gauss (two points in each direction), but it turns
< Go Back
out that this element locks. The reduced quadrature rule would then be one-point Gauss,
inte
which leads to the following expression for f :
int e h
f p≈ 8N a, j ( 0 )T ij ( 0 )j 0 ( 0 ) , (3.189)
where, as indicated, all quantities are evaluated at the origin 0 of the parent coordinates.
Theory Manuals (2/19/99) Finite Element Formulations - Advanced Element Design Issues - Selective/Fully Reduced Integration 8
Schemes such as this have the advantages of being cheap and of eliminating locking but
int e
come at a cost. They do not accurately integrate those parts of the integrand of f p that
SEACAS
Library come from deformation varying spatially in a superlinear fashion.
Theory
Manuals
Finite Element
Formulations
Advanced
Element Design
< Go Back
Theory Manuals (2/19/99) Finite Element Formulations - Advanced Element Design Issues - Selective/Fully Reduced Integration 9
Hourglass Control
SEACAS
Library Arguably the most important work done in this area was published by [Flanagan, D.P.
and Belytschko, T., 1981]. The development in this section closely follows their original
presentation. To understand more clearly the possible spurious behaviors enabled by
Theory reduced integration, we first note that the shape functions N a can be written in terms of
Manuals
some standardized element deformation modes as follows:
1 1 1 1
N a = --- Σ a + --- rΛ 1a + --- sΛ 2a + --- tΛ 3a
Finite Element 8 4 4 4
Formulations . (3.190)
1 1 1 1
+ --- stΓ 1a + --- rtΓ 2a + --- rsΓ 3a + --- rstΓ 4a
2 2 2 2
These modes are depicted schematically in Figure 3.14. As can be seen therein, Σ a
Advanced
Element Design represents a rigid body translation, the Λ ia represent constant strain deformation modes,
< Go Back and the Γ αa are referred to as hourglass modes.
If we consider the velocity field representable by this element, we find that it can be
written via
nen
∑
e e e
V = V i e i, Vi = V ia N a . (3.191)
a=1
Theory Manuals (2/19/99) Finite Element Formulations - Advanced Element Design Issues - Hourglass Control 10
SEACAS
Library
Theory Σa Λ 1a Λ 2a
Manuals
Finite Element
Formulations
Λ 3a Γ 1a Γ 2a
Advanced
Element Design
< Go Back
Γ 3a Γ 4a
Figure 3.14 Mode shapes for the eight-noded hexahedron element.
The fully linear portion of the velocity field is made up of the Σ a and Λ ia modes, so that
the hourglass portion of the velocity field can be written as
hg lin 1
V ia = V ia – V ia = ------- q ˙iα Γ αa , (3.192)
8
SEACAS
Library 1
where ------- is a normalizing factor, and q ˙iα are the hourglass normal velocities. It turns out
8
that the hourglass velocities are orthogonal to the element’s other modes in that
Theory
hg
Manuals
V ia Σ a = 0 (3.193)
and
Finite Element
hg
Formulations
V ia Λ ia = 0 . (3.194)
Basically the one-point integration scheme discussed in the last section fully controls the
linear modes of the system but provides no resistance at all to the hourglass velocities
Advanced hg
Element Design V ia . The objective of hourglass control therefore is to restore such control even in the
< Go Back context of one-point integration.
Flanagan and Belytschko wrote the hourglass nodal velocities in terms of the hourglass
shape vector γ αa via
1
q ˙iα = ------- V ia γ αa , (3.195)
8
with the shape vector γ αa found to be
1 ∂V e
γ αa = Γ αa – --- e X ib Γ αb . (3.196)
V ∂X
ia
SEACAS
Library
hg
Hourglass forces f ia are applied in these directions, so as to be orthogonal to the physical
modes of the system. One choice is
Theory
hg 1
f ia = --- Q iα γ αa ,
Manuals
(3.197)
2
where the generalized forces Q iα are given via

Finite Element
Formulations
µ̂ ∂V ∂V ˙ ,
Q iα = ---------- e q (3.198)
50V ∂ X ∂ X e iα
ja ja
Advanced with µ̂ being an effective shear modulus. The interested reader should refer to Hourglass
Element Design
Control Algorithm in Eight-Node Uniform Strain Element for the details of a specific
< Go Back
implementation of hourglass control.
Revisited
Element Examples
Theory Eight-Node Uniform Strain Element
Manuals
Finite Element
Formulations
< Go Back
Blue text
i indicates
a link to more
information.
Theory Manuals (2/19/99) Finite Element Formulations - Element Examples 14

Eight-Node Four-Node
SEACAS Uniform Corotational
Library Strain Element Shell
Eight-Node Uniform Strain Element

Theory Introduction
Manuals
Element Force Vector
Lumped Mass Matrix
Finite Element
Formulations
Finite Rotation Algorithm
Determination of Effective Moduli
Determination of the Stable Time Increment
Element Hourglass Control Algorithm
Examples
Artificial Bulk Viscosity
< Go Back
Theory Manuals (2/19/99) Finite Element Formulations - Element Examples: Eight-Node Uniform Strain Element 1
Introduction
SEACAS
Library The eight-node, three-dimensional isoparametric element is widely used in computational
mechanics. The determination of optimal integration schemes for this element, however,
presents a difficult dilemma. A one-point integration of the element underintegrates the
Theory element, resulting in a rank deficiency that manifests itself in spurious zero energy modes,
Manuals commonly referred to as hourglass modes (see Hourglass Control). A two-by-two-by-
two integration of the element, by contrast, overintegrates the element and can lead to
serious problems of element locking in fully plastic and incompressible problems (see
Finite Element Constrained Media and Locking). The eight-point integration also carries a tremendous
Formulations
computational penalty compared to the one-point rule. Particularly in explicit dynamic
applications (see Explicit Finite Element Methods), this added expense is extremely
undesirable.
Element In this section we present an element that is widely used in explicit analyses, wherein one-
Examples
point integration is utilized in combination with an hourglass control scheme that controls
the spurious modes. The implementation presented below follows directly from
[Flanagan, D.P. and Belytschko, T., 1981]. In particular the aspects of this element
Eight-Node
Uniform Strain pertaining to hourglass control have already been discussed in Hourglass Control.
Element
e
< Go Back The hexahedral element relates the spatial element coordinates x i to the nodal
coordinates d ia through the isoparametric shape functions N a as follows
e
x i = N a ( ξ, η, ζ )d ia , (3.199)
Theory Manuals (2/19/99) Finite Element Formulations - Eight-Node Uniform Strain Element - Introduction 2
where for convenience the summation convention on nodal indices a has been adopted.
Subscripts i have a range of three, corresponding to the spatial coordinate directions, and
SEACAS subscripts a have a range of eight.
Library
Our discussion here will focus primarily on explicit dynamics applications, where the
argument for the use of this element is most compelling. As such it is necessary to have
Theory expressions for the elemental velocity field
Manuals
v i = N a v ia , (3.200)
and for the acceleration field

Finite Element
Formulations
a i = N a a ia . (3.201)
In (3.200) and (3.201) the nodal velocities v ia and nodal accelerations a ia are the
Element localized global velocities and accelerations, which for central differences are produced by
Examples
the updates (3.75) given in Explicit Finite Element Methods.
The velocity gradient tensor L has been discussed in Rates of Deformation and was
Eight-Node specified in Eq. (2.30). It can be written in the domain of an element as
Uniform Strain
Element
L ij = v i, j = v ia N a, j . (3.202)
< Go Back
By convention a comma preceding a lowercase subscript denotes differentiation with
∂v i
respect to the spatial coordinates (e.g., v i, j denotes --------- ).
∂x j
The three-dimensional isoparametric shape functions map the unit cube in r -space ( r is
written explicitly as ( r, s, t ) to a general hexahedron in x -space, as depicted for the
SEACAS element reference configuration in Figure 3.6. The trilinear shape functions defined over
Library
this domain, as summarized in Shape Functions, can be conveniently expanded in terms
of an orthogonal set of base vectors, as was mentioned previously in Hourglass Control:
1 1 1 1
Theory
Manuals
N a = --- Σ a + --- rΛ 1a + --- sΛ 2a + --- tΛ 3a
8 4 4 4
. (3.203)
1 1 1 1
+ --- stΓ 1a + --- rtΓ 2a + --- rsΓ 3a + --- rstΓ 4a
2 2 2 2
Finite Element
Formulations The basis vectors represent the displacement modes of a unit cube, as was also discussed
in Hourglass Control. The first vector, Σ a , accounts for rigid body translation. The
vectors Λ ia may be readily combined to define three uniform normal strains and three
Element rigid body rotation modes for the unit cube. We refer to Λ ia as the volumetric base
Examples
vectors since, as we will illustrate below, they are the only base vectors which appear in
the element volume expression. The last four vectors, Γ αa , where Greek subscripts have a
Eight-Node range of four, give rise to linear strain modes which are neglected by uniform strain
Uniform Strain
Element integration (i.e., the one-point quadrature rule summarized in (3.189)). These vectors
< Go Back define the hourglass patterns for a unit cube, so that we refer to Γ αa as the hourglass base
vectors. The displacement modes represented by all these vectors are shown in Figure
3.14.
Element Force Vector
SEACAS
Library Recalling the development of Basics of Element Design, the generic expression for the
element internal force vector (see (3.142)) is
int e
Theory
Manuals
f ia = ∫ N a, j T ij dv . (3.204)
h e
ϕt( Ω )
In the element we consider, the one-point integration scheme neglects the nonlinear
Finite Element portion of the element displacement field, thereby considering a state of uniform strain and
Formulations
stress. The preceding expression is approximated by
int e
f ia ≈ T ij ∫ N a, j dv , (3.205)
h e
Element ϕt( Ω )
Examples
where T ij , the mean stress tensor, represents the assumed uniform stress field. By
neglecting the nonlinear displacements, we have assumed that the mean stresses depend
Eight-Node
Uniform Strain only on the mean strains. Mean kinematic quantities are defined by integrating over the
Element element as follows:
< Go Back
1
v i, j = ---
v ∫ v i, j dv . (3.206)
h e
ϕt( Ω )
We now define the discrete gradient operator as

Theory Manuals (2/19/99) Finite Element Formulations - Eight-Node Uniform Strain Element - Element Force Vector 5
B ia = ∫ N a, i dv . (3.207)
h e
ϕt( Ω )
SEACAS
Library
The mean velocity gradient, applying Eq. (3.202), is given by
1
v i, j = ---v ia B ja . (3.208)
Theory v
Manuals
Combining Equations (3.205) and (3.207), we may express the nodal forces by
int e
f ia = Tij Bja . (3.209)
Finite Element
Formulations
Computing nodal forces with this integration scheme requires evaluation of the gradient
operator and the element volume. These two tasks are linked since
x i, j = δ ij , (3.210)
Element
Examples
where δ ij is the Kronecker delta. Equations (3.199), (3.207), and (3.210) yield
Eight-Node
Uniform Strain
x ia B ja = ∫ ( x ia N a ) , j dv = vδ ij . (3.211)
e
Element v
< Go Back Consequently, the gradient operator may be expressed by
∂v
B ia = ----------------- . (3.212)
∂x ia
To integrate the element area in closed form, we use the Jacobian of the isoparametric
transformation to transform to an integral over the biunit cube:
SEACAS
Library +1 +1 +1
v = ∫ dv = ∫ ∫ ∫ j dr ds dt . (3.213)
v -1 -1 -1
Theory The Jacobian is given in terms of the permutation symbol e ijk as

Manuals
∂x ∂y ∂z
J = e ijk --------- --------- --------- . (3.214)
∂r i ∂r j ∂r k
Finite Element
Formulations
Therefore, Equation (3.213) can be written as
v = x a y b z c C abc , (3.215)
Element where
Examples
+1 ∂N
a ∂N b ∂N c
+1 +1
Cabc = eijk ∫ ∫ ∫ ---------- ----------- ------------- dr1 dr2 dr3 .
∂ri ∂rj ∂rk
(3.216)
-1 -1 -1
Eight-Node
Uniform Strain
Element Observe that the coefficient array C abc is identical for all hexahedrons. Furthermore, it
< Go Back possesses the following alternator properties:
C abc = C bca = C cab = – C acb = – C bac – C cba . (3.217)
Therefore, applying Equations (3.212) and (3.217) to (3.215) yields the following form for
evaluating the B-matrix:
SEACAS
Library ybzc
B ia = z b x c C abc . (3.218)
xbyc
Theory
Manuals
In light of Equation (3.203), it is evident that evaluating each component of C abc involves
integrating a polynomial that is at most biquadratic. However, since we are integrating
Finite Element over a symmetric region, any term with a linear dependence will vanish. The only terms
Formulations which survive the integration will be the constant, square, double square, and triple square
terms. Furthermore, the alternator properties cause half of these remaining terms to drop
out. The resulting expression for C abc is
Element
1
C abc = --------- e ijk ( 3 Λ ia Λ jb Λ kc + Λ ia Γ kb Γ jc
Examples
192 . (3.219)
+ Γ ka Λ jb Γ ic + Γ ja Γ ib Λ kc )
Eight-Node
Uniform Strain
Element Since C abc has the alternator properties given in Equation (3.217), only 56 (the
< Go Back combination of eight nodes taken three at a time) distinct nonzero terms are possible.
However, the volume must be independent of the selection of Node 1, which implies that
C abc must be invariant under all node numberings, preserving the relative orientation of
the element. Consequently, only 21 (the combination of seven nodes taken two at a time)
terms may be independent. Furthermore, once Node 1 is selected, three orientations of the
node numbering system are possible, so that only seven terms of C abc need be evaluated.
SEACAS
Library The seven sets of triples ( a, b, c ) giving rise to independent terms of C abc are: ( 1, 2, 3 ) ;
( 1, 2, 5 ) ; ( 1, 2, 6 ) ; ( 1, 2, 7 ) ; ( 1, 2, 8 ) ; ( 1, 3, 5 ) ; 1, 3, 6 . Of these, only the first three
terms do not vanish. All other nonzero terms of C abc are found by permutations and use
Theory
Manuals of the alternator properties summarized by Eq. (3.217).
With the C abc in hand, the first term of B ia is expressed using (3.218) as
Finite Element 1
Formulations B11 = --- [ y 2 ( ( z6 – z3 ) – ( z4 – z5 ) ) + y3 ( z2 – z4 )
2
+ y 4 ( ( z3 – z8 ) – ( z5 – z2 ) ) + y 5 ( ( z8 – z6 ) – ( z2 – z4 ) ) . (3.220)
+ y 6 ( z5 – z2 ) + y 8 ( z4 – z5 ) ]
Element
Examples
After permuting the nodes according to the transformations described above, other terms
of B ia are also evaluated using (3.218). The element volume is most easily computed by
contracting the B-matrix and nodal coordinates as per Eq. (3.211).
Eight-Node
Uniform Strain
Element
< Go Back
Lumped Mass Matrix
SEACAS
Library In order to reap the benefits of an explicit architecture, we must diagonalize the mass
matrix (see Explicit Finite Element Methods). We do this by integrating the inertial
terms as
Theory
e e e
Manuals ( m a ) ia = a ib m ab , (3.221)
where
Finite Element
e
Formulations m ab = ρvδ ab , (3.222)
and δ ab is the Kronecker delta. Clearly the assembly process for the global mass matrix
from the individual element matrices results in a global mass matrix that is diagonal and
Element
Examples can be expressed as a vector, M P , if desired.
Eight-Node
Uniform Strain
Element
< Go Back
Theory Manuals (2/19/99) Finite Element Formulations - Eight-Node Uniform Strain Element - Lumped Mass Matrix 10
Finite Rotation Algorithm
SEACAS
Library As discussed in more detail in Frame Indifference, an important factor in proper
formulation of large deformation problems is the assurance of material objectivity. In the
element we now consider, this is achieved by formulation of the constitutive updates in the
Theory rotated configuration depicted in Figure 2.2 and introduced more thoroughly in Rate of
Manuals Deformation Tensors. Of particular interest in constitutive modeling are quantities like
the rotated rate of deformation tensor D (see Eq. (2.39)), the rotated Cauchy stress T in
Eq. (2.57), and the Green-Naghdi rate of Cauchy stress T̃ defined in (2.122).
Finite Element
Formulations
Notably all of the above objects require the determination of R , the rotation tensor defined
by the polar decomposition summarized by Eq. (2.13). Here we describe an incremental
algorithm for determination of this tensor with emphasis on computational efficiency and
numerical accuracy. We begin by considering (2.42) as a first-order differential equation in
Element
Examples R
Ṙ = L R . (3.223)
Eight-Node
Uniform Strain The crux of integrating Eq. (3.223) for R is to maintain the orthogonality of R.
Element Unfortunately, if one merely applies a forward difference scheme, the orthogonality of R
< Go Back degenerates rapidly no matter how fine the time increments. Instead the algorithm of
Hughes and Winget ([Hughes, T.J.R. and Winget, J., 1980]) for integrating incremental
rotations can be adopted as follows.
A rigid body rotation over a time increment ∆t may be represented by

Theory Manuals (2/19/99) Finite Element Formulations - Eight-Node Uniform Strain Element - Finite Rotation Algorithm 11
x t + ∆t = Q ∆t x t , (3.224)
SEACAS where Q ∆t is a proper orthogonal tensor with the same rate of rotation as R, as given by
Library
Eq. (3.223). The total rotation R is updated via
R t + ∆t = Q ∆t R t . (3.225)
Theory
Manuals For a constant rate of rotation, the midpoint velocity and the midpoint coordinates are
related by
1 1
------- ( x t + ∆t – x t ) = --- L ( x t + ∆t + x t ) . (3.226)
Finite Element
Formulations
∆t 2
Combining Eqs. (3.224) and (3.226) yields
∆t
( Q ∆t – I )x t = ------- L ( Q ∆t + I )x t . (3.227)
Element 2
Examples
Since x t is arbitrary in Equation (3.227), it may be eliminated. We then solve for Q ∆t ,

which gives
Eight-Node
Uniform Strain
Element
 ∆t  –1  ∆t 
Q ∆t = I – ------- L I + ------ L . (3.228)
< Go Back  2   2 
The accuracy of this integration scheme is dependent upon the accuracy of the midpoint
relationship of Equation (3.226). The rate of rotation must not vary significantly over the
time increment. Furthermore, [Hughes, T.J.R. and Winget, J., 1980] showed that the
conditioning of Equation (3.228) degenerates as ∆t L grows.
SEACAS
Library Our complete numerical algorithm for a single time step can be summarized as below.
• Calculate the rate of deformation tensor D and the spin tensor W (see Eqs. (2.31) and
(2.32))
Theory
Manuals
• Determine L from W and V (the left stretch, see Eq. (2.13)), using the following
algorithm due to [Dienes, J.K., 1979]. Compute
z i = e ijk V jm D mk (3.229)
Finite Element
Formulations
–1
l = w – 2 ( V – Itr(V ) ) z (3.230)
1
L ij = --- e ijk l k , (3.231)
Element
2
Examples
where
w i = e ijk W jk . (3.232)
Eight-Node
Uniform Strain • Solve
Element
< Go Back
 I – ∆t   ∆t 
------- L R t + ∆t = I + ------ L R t . (3.233)
 2   2 
• Calculate
V̇ = ( D + W )V – V L . (3.234)
• Update
SEACAS
Library
V t + ∆t = V t + ∆tV̇ . (3.235)
• Compute the rotated rate of deformation (see (2.39))
Theory T
Manuals D = R DR . (3.236)
• Integrate the constitutive equations in the rotated frame of reference
Finite Element Ṫ = f ( D , T ) . (3.237)

Formulations
• Compute the Cauchy stress in the spatial configuration
T
T = RT R . (3.238)
Element Note that this algorithm requires that the tensors V and R be stored in memory for each
Examples
element.
Eight-Node
Uniform Strain
Element
< Go Back
Determination of Effective Moduli
SEACAS
Library Algorithms for calculating the stable time increment require effective moduli for each
element (see (3.73) in Explicit Finite Element Methods). Such calculations of
dilatational and shear moduli are also necessary for hourglass control, bulk viscosity, and
Theory nonreflecting boundaries. Here we present a procedure for adaptively determining the
Manuals effective dilatational and shear moduli of the material.
In an explicit integration algorithm, the constitutive response over a time step can be recast
a posteriori as a hypoelastic relationship. We approximate this relationship as isotropic.
Finite Element
Formulations This defines effective Lamé parameters, λ̂ and µ̂ , in terms of the hypoelastic stress
increment and strain increment (in the rotated frame of reference) as follows:
∆T ij = ∆t ( λ̂D kk δ ij + 2µ̂D ij ) . (3.239)

Element
Examples
Equation (3.239) can be rewritten in terms of volumetric and deviatoric parts as
∆T kk = ∆t ( 3λ̂ + 2µ̂ )D kk , (3.240)

Eight-Node
Uniform Strain
Element and
< Go Back
s ij = ∆t2µ̂e ij , (3.241)
where
Theory Manuals (2/19/99) Finite Element Formulations - Eight-Node Uniform Strain Element - Determination of Effective Moduli 15
1
s ij = ∆T ij – --- ∆T kk δ ij (3.242)
3
SEACAS
Library and
1
e ij = D ij – --- D kk δ ij . (3.243)
3
Theory
Manuals The effective bulk modulus K̂ follows directly from Equation (3.240) as
∆T kk
3K̂ = 3λ̂ + 2µ̂ = -----------------
-. (3.244)
Finite Element ∆tD kk
Formulations
Taking the inner product of Equation (3.241) with the deviatoric strain rate and solving for
the effective shear modulus 2µ̂ , gives
Element s ij e ij
Examples 2µ̂ = --------------------------
-. (3.245)
∆te mn e mn
Using the result of Equation (3.244) with Equation (3.245), we can calculate the effective
Eight-Node
Uniform Strain dilatational modulus λ̂ + 2µ̂ :
Element
< Go Back 1
λ̂ + 2µ̂ = --- ( 3K̂ + 2 ⋅ ( 2µ̂ ) ) . (3.246)
3
If the strain increments are insignificant, Equations (3.244) and (3.245) will not yield
numerically meaningful results. In this circumstance the dilatational modulus can be set to
an initial estimate, λ 0 + 2µ 0 . An initial estimate of the dilatational modulus is, therefore,
the only parameter which every constitutive model is required to provide to the time step
SEACAS control algorithm.
Library
In a case where the volumetric strain increment is significant but the deviatoric increment
is not, the effective shear modulus can be estimated by rearranging Equation (3.246) as
Theory follows:
Manuals
1
2µ̂ = --- ( 3 ( λ 0 + 2µ 0 ) – 3K̂ ) . (3.247)
2
Finite Element If neither strain increment is significant, the effective shear modulus can be set equal to the
Formulations
initial dilatational modulus.
Element
Examples
Eight-Node
Uniform Strain
Element
< Go Back
Determination of the Stable Time Increment
SEACAS
Library Flanagan and Belytschko [Flanagan, D.P. and Belytschko, T., 1984] provided eigenvalue
estimates for the uniform strain hexahedron described in this section. They showed that
the maximum eigenvalue is bounded by
Theory
Manuals λ + 2µ B ia B ia 2 8 λ + 2µ B ia B ia
- ≥ ω max ≥ --- ------------------ --------------------
8 ---------------- -------------------- -. (3.248)
ρ v
2 3 ρ v
2
Finite Element
Using the effective dilatational modulus from Determination of Effective Moduli with
Formulations the eigenvalue estimates of Equation (3.248) allows us to write the stability criteria of Eq.
(3.73) as
1
2 ( ρ0 V0 )
∆t̂ ≤ --- ----------------------------------------- . (3.249)
Element 2 ( λ + 2µ )B iI B iI
Examples
The stable time increment is determined from Equation (3.249) as the minimum over all
elements.
Eight-Node
Uniform Strain Equation (3.249) is numerically invalid if the effective dilatational modulus is less than or
Element
equal to zero. A negative modulus indicates a strain softening situation that renders the
< Go Back central difference operator unconditionally unstable. In practice, however, strain softening
is generally short lived, so that the calculations can continue in a stable manner once the
softening energy has been dissipated. To aid the user in controlling an unstable strain-
Theory Manuals (2/19/99) Finite Element Formulations - Eight-Node Uniform Strain Element - Determination of the Stable Time Increment 18
softening situation, the effective dilatational modulus can be adjusted with a strain-
softening scale factor, ssft, as follows
SEACAS
Library λ 0 + 2µ 0
If λ̂ + 2µ̂ < 0 then λ̂ + 2µ̂ = --------------------- . (3.250)
2
( ssft )
To avoid dividing by zero in Equation (3.249), one can enforce the following condition:
Theory
Manuals
–6
λ̂ + 2µ̂ ≥ ( λ 0 + 2µ 0 ) ⋅ 10 . (3.251)
The estimate of the critical time increment given in Equation (3.249) is for the case where
Finite Element
Formulations there is no damping in the system. If we define ε as the fraction of critical damping in the
highest element mode, the stability criteria of Eq. (3.249) becomes
2
∆t ≤ ∆t̂ ( 1 + ε – ε ) . (3.252)
Element
Examples
Conventional estimates of the critical time increment size have been based on the transit
time of a dilatational wave over the shortest dimension l of an element or zone. For the
Eight-Node
undamped case this gives
Uniform Strain
Element l
∆t ≤ --- , (3.253)
< Go Back c
where c is the dilatational wave speed.
There are two fundamental and important differences between the time increment limits
given by Equations (3.249) and (3.253). First, our time increment limit is dependent on a
characteristic element dimension, which is based on the finite element gradient operator
and does not require an ad hoc guess of this dimension. This characteristic element
SEACAS dimension, l , is defined by inspection of Equation (3.249) as
Library
1 v
l = --- ----------------------- . (3.254)
2 B B
iI iI
Theory
Manuals Second, the sound speed used in the estimate is based on the current response of the
material and not on the original elastic sound speed. For materials that experience a
reduction in stiffness due to plastic flow, this can result in significant increases in the
Finite Element critical time increment.
Formulations
It should be noted that the stability analysis performed at each time step predicts the
critical time increment for the next step. Our assumption is that the conservativeness of
this estimate compensates for any reduction in the stable time increment over a single time
Element step.
Examples
Eight-Node
Uniform Strain
Element
< Go Back
Hourglass Control Algorithm
SEACAS
Library The mean stress-strain formulation of the uniform strain element considers only a fully
linear velocity field. The remaining portion of the nodal velocity field is the so-called
hourglass field. Excitation of these modes may lead to severe, unresisted mesh distortion.
Theory The hourglass control algorithm described here is taken directly from [Flanagan, D.P.
Manuals and Belytschko, T., 1981]. The method isolates the hourglass modes so that they may be
treated independently of the rigid body and uniform strain modes.
A fully linear velocity field for the hexahedron can be described by
Finite Element
Formulations
LIN
vi = v i + v i, j ( x j – x j ) . (3.255)
The mean coordinates x i correspond to the center of the element and are defined as
Element
Examples
1
x i = --- x ia Σ a . (3.256)
8
Eight-Node The mean translational velocity is similarly defined by

Uniform Strain
Element
1
v i = --- v ia Σ a . (3.257)
< Go Back 8
The linear portion of the nodal velocity field may be expressed by specializing Eq. (3.255)
to the nodes as follows:
Theory Manuals (2/19/99) Finite Element Formulations - Eight-Node Uniform Strain Element - Hourglass Control Algorithm 21
LIN
v ia = v i Σ a + v i, j ( x ja – x j Σ a ) , (3.258)
SEACAS where Σ a is used to maintain consistent index notation and indicates that v i and x j are
Library
independent of position within the element. From Equations (3.211) and (3.258) and the
orthogonality of the base vectors, it follows that
LIN
v ia Σ a = v ia Σ a = 8v i ,
Theory
Manuals (3.259)
and
Finite Element LIN

Formulations
v ia B ja = v ia B ja = vv i, j . (3.260)
HG
The hourglass field v ia may now be defined by removing the linear portion of the nodal
velocity field:
Element
Examples
HG LIN
v ia = v ia – v ia . (3.261)
Eight-Node
Equations (3.259) through (3.261) prove that Σ a and B ja are orthogonal to the hourglass
Uniform Strain field:
Element
< Go Back HG
v ia Σ a = 0 (3.262)
HG
v ia B ja = 0 . (3.263)
Furthermore, it can be shown that the B matrix is a linear combination of the volumetric
base vectors Λ ia , so Eq. (3.263) can be written as
SEACAS
Library HG
v ia Λ ia = 0 . (3.264)
Equations (3.262) and (3.264) show that the hourglass field is orthogonal to all the base
Theory HG
Manuals vectors depicted in Figure 3.14 except the hourglass base vectors. Therefore, v ia may be
expanded as a linear combination of the hourglass base vectors as follows:
HG 1
Finite Element v ia = ------- q̇ iα Γ αa . (3.265)
Formulations 8
The hourglass nodal velocities are represented by q̇ iα above (the leading constant is
added to normalize Γ αa ). We now define the hourglass shape vector γ αa such that
Element
Examples
1
q̇ i = ------- u̇ ia γ αa . (3.266)
8
Eight-Node
Uniform Strain
By substituting Equations (3.258), (3.261), and (3.266) into (3.265), then multiplying by
Element Γ αa and using the orthogonality of the base vectors, we obtain the following:
< Go Back
u̇ iI Γ αI – u̇ i, j x jI Γ αI = u̇ iI γ αI . (3.267)
With the definition of the mean velocity gradient, Equation (3.208), we can eliminate the
nodal velocities above. As a result, we can compute γ αI from the following expression:
SEACAS
Library
1
γ αa = Γ αa – ---B ia x ib Γ αb . (3.268)
V
The difference between the hourglass base vectors Γ αa and the hourglass shape vectors
Theory
Manuals γ αa is very important. They are identical if and only if the hexahedron is a right-
parallelepiped. For a general shape Γ αa is orthogonal to B ja , whereas γ αa is orthogonal
LIN
Finite Element to the linear velocity field v ia . Γ αa defines the hourglass pattern, and γ αa is necessary
Formulations
to accurately detect hourglassing.
For the purpose of controlling the hourglass modes, we define generalized forces Q iα ,
Element which are conjugate to q̇ iα so that the rate of work is
Examples
HG 1
v ia f ia = --- Q iα q̇ iα (3.269)
2
Eight-Node
Uniform Strain for arbitrary u̇ iI . Using Equation (3.266) it follows that the contribution of the hourglass
Element
resistance to the nodal forces is given by
< Go Back
HG 1
f ia = --- Q iα γ αa . (3.270)
2
Two types of hourglass resistance are possible: artificial stiffness and artificial damping.
Considering the stiffness type as an example, we can define a tuneable hourglass stiffness,
SEACAS κ , and express the resistance by
Library
κ B jb B jb
Q iα = --- 2µ̂ --------------------- q̇ iα . (3.271)
2 v
Theory
Manuals Note that the stiffness expression must be integrated, which further requires that this
resistance be stored in a global array.
Observe that the nodal antihourglass forces of Equation (3.270) have the shape of γ αa
Finite Element
Formulations rather than Γ αa . This fact is essential since the antihourglass forces should be orthogonal
to the linear velocity field, so that no energy is transferred to or from the rigid body and
uniform strain modes by the antihourglassing scheme.
Element
Examples
Eight-Node
Uniform Strain
Element
< Go Back
Artificial Bulk Viscosity
SEACAS
Library Artificial viscosity may be desirable in numerical calculations for two reasons. First, high-
velocity gradients can collapse an element before it has a chance to respond if no viscosity
is employed. Second, viscosity is often useful in quieting truncation frequency “ringing”.
Theory
Manuals
Ideally one would like to add viscosity only to the highest mode of the element, but
isolating this mode is impractical. The standard technique is to simply add viscosity to the
volumetric or “bulk” response. This generates a viscous pressure in terms of the volume
strain rate as follows:
Finite Element
Formulations
V̇  V̇ 2
q = b 1 ρc --- – ρ b 2 l ----- , (3.272)
V  V
where b1 and b2 are coefficients for the linear and quadratic terms, respectively. The
Element
Examples quadratic term in Equation (3.272) is more important and is designed to “smear” a shock
front across several elements. This term yields a jump in energy as a smeared shock
passes, which simulates the shock heating. As a result, the smeared shock front can be
Eight-Node propagated as a steady wave.
Uniform Strain
Element
The linear term is intended to dissipate truncation frequency oscillations. The quadratic
< Go Back term is only applied to compressive strain rates, since an element cannot collapse in
expansion.
The preceding expression is simplified if we use the undamped stable time increment
defined by Equation (3.249) and write
Theory Manuals (2/19/99) Finite Element Formulations - Eight-Node Uniform Strain Element - Artificial Bulk Viscosity 26
2
l V ρ
∆t̂ = --- = ---------------------- ⋅ ---------------- , (3.273)
c 2B ia B ia λ + 2µ
SEACAS
Library
or
m V
∆t̂ = ---------------- ⋅ ---------------------- , (3.274)
Theory λ + 2µ 2B ia B ia
Manuals
where m is the element mass. We now define the factor ε such that the quadratic viscosity
term vanishes in expansion
Finite Element
2
Formulations
ε = b 1 – b 2 ∆t̂ min ( 0, D kk ) . (3.275)
This quantity is required for the damped stability criteria of Equation (3.252). Note that
the condition imposed by (3.251) prevents Equation (3.275) from yielding so large a value
Element
Examples
of ε that Equation (3.252) would numerically yield a zero value.
We will show below that ε can be used to estimate the fraction of critical damping in the
highest element mode. Using Equation (3.275) in Equation (3.274) allows us to write the
Eight-Node
Uniform Strain
viscous pressure as
Element
2
< Go Back q = ( b 1 – b 2 ∆t̂D kk ) ( λ + 2µ )∆t̂D kk . (3.276)
The bulk viscosity pressure is appended to the stresses during the internal force
calculations to yield the following forces:
f ia = qB ia . (3.277)
SEACAS
The above expression can be expanded using Equations (3.274) and (3.275) to yield
Library
l
f ia = ερc ---B jb B ia u̇ jb . (3.278)
V
Theory This form indicates that if B ia is an eigenvector, the modal damping is

Manuals
cV
ερ --------- . (3.279)
l
Finite Element
Formulations The critical damping estimate of the maximum element frequency is
ρV 2c cV
2mω = 2 --------- -------- = ρ --------- . (3.280)
8 l 2l
Element
Examples The two expressions above show that ε is half the fraction of critical damping in the
highest mode.
Eight-Node
Uniform Strain
Element
< Go Back
Eight-Node Four-Node
SEACAS Uniform Corotational
Library Strain Element Shell

Theory Introduction
Manuals
Shell Kinematics
Constitutive Assumptions
Finite Element
Formulations
Shell Element Coordinate Systems
Element Equations
Computation of Internal Force and Moment Resultants
Element Hourglass Control
Examples
Calculation of the Stable Time Increment
< Go Back Constitutive Models
Time-Stepping Algorithm
Theory Manuals (2/19/99) Finite Element Formulations - Element Examples: Four-Node Corotational Shell 1
Introduction
SEACAS
Library In this section we discuss in detail some of the implementational issues associated with a
frequently utilized structural element in nonlinear mechanics: the four-noded corotational
shell element depicted in Figure 3.15. In so doing, we will add some important detail to
Theory the very conceptual discussion of shells and other structural entities given in Structural
Manuals Components.
4
Finite Element
Formulations r 42
ê
3
ê r 31
2
1 3
Element ê
Examples
1
r 21
r 21
⊥
Four-Node
Corotational 2
Shell Figure 3.15 Four-noded, corotational shell element also showing the element coordinate system.
< Go Back Although much of the discussion is equally applicable to matrix-free quasistatic solution
strategies (see Conjugate Gradient Methods), We target our discussion here primarily to
explicit dynamic calculations (see Explicit Finite Element Methods). In such settings the
equations of motion of the deformable body thus become a system of uncoupled equations
Theory Manuals (2/19/99) Finite Element Formulations - Four-Node Corotational Shell - Introduction 2
governing the nodal motions in the discretized system. For continuum elements each
equation in the system is the equation for three-dimensional motion of a particle,
SEACAS
EXT INT
Library a iA = ( F iA – F iA ) ⁄ M A , (3.281)
EXT INT
where a iA is the acceleration, F iA and F iA are the external and internal forces, and
Theory
Manuals M A is the lumped mass, all associated with global Node A. For a continuum element with
displacement degrees of freedom d ia , a localized version of Eq. (3.281) completely
describes the motion of the nodes.
Finite Element
Formulations
By contrast, a shell element requires rotational degrees of freedom in addition to the

displacement degrees of freedom The additional equations governing these degrees of
Element freedom are Euler’s equations for the rotation of a rigid body about the principal axes,
Examples written here for an individual shell element
EXT INT
α 1b = [ ( m 1b – m 1b ) – ( I3b – I2b )ω 3b ω 2b ] ⁄ I1b
Four-Node EXT INT
Corotational α 2b = [ ( m 2b – m 2b ) – ( I1b – I3b )ω 1b ω 3b ] ⁄ I2b . (3.282)
Shell
EXT INT
α 3b = [ ( m 3b – m 3b ) – ( I2b – I1b )ω 2b ω 1b ] ⁄ I3b
< Go Back
In these equations α 1b is the angular acceleration for local Node b , ω ib is angular
velocity, and I 1b is the mass moment of inertia associated with local Node b in principal
Direction 1. Similarly numerical subscripts 2 and 3 designate the other principal
directions. A coordinate system is used at each node to track rotation of the principal axis
system.
SEACAS
Library
Theory
Manuals
Finite Element
Formulations
Element
Examples
Four-Node
Corotational
Shell
< Go Back
Shell Kinematics
SEACAS
Library In Mindlin shell theory (see [Mindlin, R.D., 1951]), the shell normals are assumed to
remain straight, although they are allowed to rotate. Rotation of the normals allows the
element to model transverse shear strains. Because displacements are assumed to vary
Theory linearly through the shell thickness, the velocity at any point can be expressed as
Manuals
v∗ = v + ẑê 3 × q̇ , (3.283)
Finite Element
in which v∗ is the velocity of a point in the shell body, v is the velocity of the point on the
Formulations midsurface lying on the same normal, and q̇ is the rotational velocity of the normal (see
Figure 3.16). ẑ is the coordinate through the shell thickness along the unit normal vector,
ê . Vector components in the corotational coordinate system are indicated by the ˆ
3
Element
Examples
symbol.
The components of velocity strain in the corotational coordinate system are
Four-Node
ˆ 1  ∂v̂ ∂v̂ 
Corotational d ij = ---  ---------i + ---------j  . (3.284)
2  ∂x̂ 
j ∂x̂ i
Shell
< Go Back
Strain-displacement relationships for the shell are obtained by substitution of Eq. (3.283)
into Eq. (3.284), giving
Theory Manuals (2/19/99) Finite Element Formulations - Four-Node Corotational Shell - Shell Kinematics 5
ˆd ∂v̂ 1 ∂q˙2
11 = --------- + ẑ ---------
SEACAS
∂ x̂ 1 ∂ x̂ 1
Library
ˆd ∂v̂ 2 ∂q˙1
22 = --------- – ẑ ---------
∂ x̂ 2 ∂ x̂ 2
Theory
Manuals ˆd 1 ∂v̂ 1 ∂v̂ 2  ∂q˙2 ∂q˙1 
12 = --- --------- + --------- + ẑ  --------- – ---------  , (3.285)
2 ∂ x̂
2 ∂ x̂ 1
 ∂ x̂ 2 ∂ x̂ 1 
 ∂v̂ 
Finite Element dˆ 23 = 1---  --------3- – q˙1
Formulations 2  ∂ x̂ 
2
ˆd 1  ∂v̂ 3 ˙ 
13 = ---  --------- + q 2
2  ∂ x̂ 
Element 1
Examples
where all quantities are expressed in the corotated coordinate system; the x̂ i being the
Four-Node
coordinates in this system, and the v̂ i being the spatial midsurface velocities expressed in
Corotational
Shell
this system.
< Go Back
ζm
t + ∆t
ẑ êt + ∆t ( P∗ )
SEACAS
3 t + ∆t
P
Library
u∗ = u + ζm t + ∆t
êt = êt + ∆t u∗
3
u
3
Theory
Manuals
u
t
Finite Element P
t ( P∗ )
Formulations
ẑ êt
t 3
Figure 3.16 Shell kinematics.
Element
Examples
Four-Node
Corotational
Shell
< Go Back
Constitutive Assumptions
SEACAS
Library If the velocity strain components in the corotational system are arranged in a column
matrix, d , as
d = [ dˆ 11, dˆ 22, 2dˆ 12, 2dˆ 23, 2dˆ 13 ] ,

T
Theory (3.286)
Manuals
then the conjugate Cauchy stress components in the corotated coordinate system can be
written as
Finite Element
T
Formulations
T = [ σ̂ 11, σ̂ 22, σ̂ 12, σ̂ 23, σ̂ 13 ] ,. (3.287)
Furthermore, the shell is assumed to be in a state of plane stress, so
Element σ̂ 33 = 0 . (3.288)
Examples
Note also the omission of dˆ 33 , the rate of through-the-thickness thinning, from (3.285)
Four-Node and (3.286).
Corotational
Shell The representations of deformation and stress given in (3.286) and (3.287) are conjugate,
< Go Back in the sense that the stress power P , discussed in general continuum mechanical terms in
Stress Power, can be expressed as
T
P = d T̂ . (3.289)
Theory Manuals (2/19/99) Finite Element Formulations - Four-Node Corotational Shell - Constitutive Assumptions 8
Shell Element Coordinate Systems
SEACAS
Library We consider now the element formulation of a four-node, quadrilateral shell (see Figure
3.15), where the midsurface velocities v̂ and rotation rates q˙ are interpolated using the
i i
bilinear shape functions described in Shape Functions. Thickness of the shell is handled
Theory
Manuals
as an element attribute. The element uses reduced integration with hourglass control and is
based on a corotational formulation (see [Belytschko, T., Lin, J.I. and Tsay, C.S., 1984]
and [Belytschko, T., Ong, J.S.-J., Liu, W.K. and Kennedy, J.M., 1984])
Finite Element Three coordinate systems are used in the shell element formulation. The translational
Formulations
equations of motion (Eq. (3.281)) are written in the global system with coordinates x i
and basis vectors e i . Strain-displacement relationships and constitutive equations are
enforced in a local element coordinate system that rotates and translates with each
Element
Examples element. This corotational coordinate system has orthonormal basis vectors denoted by ê i
and coordinates designated by x̂ i . The internal element force ( f̂ ) and moment ( m̂ )
Four-Node resultants are also computed in this element coordinate system. New element coordinate
Corotational
Shell
systems are computed at each time step using current nodal coordinates. The equations
governing the rotational motion are written in a local coordinate system at each node.
< Go Back
These nodal coordinate systems are assumed to rotate with the principal axes at each node
with the motion governed by Eqs. (3.282). The nodal systems have orthonormal basis
vectors, e i , with coordinates designated by x i . The nodal coordinate systems are
updated at each time step.
Theory Manuals (2/19/99) Finite Element Formulations - Four-Node Corotational Shell - Shell Element Coordinate Systems 9
The element (corotational) normal vector ê 3 is approximated by the normalized cross
product of the element diagonals (see Figure 3.15)
SEACAS
Library
r 31 × r 42
ê 3 = ---------------------------- (3.290)
r 31 × r 42
Theory In (3.290) r ab is the position vector of element Node a relative to element Node b in the
Manuals
current configuration. The direction of the normal vector is thus determined by the
element-nodal connectivity, and the positive side of the shell is the one for which the
nodes are numbered counterclockwise. Next the direction of ê 1 is taken as the portion of
Finite Element
Formulations the vector connecting nNodes 1 and 2 that is orthogonal to ê 3 .
r 21⊥ = r 21 – ( r 21 ê 3 )ê 3 (3.291)
Element
Examples
ê 1 = r 21⊥ ⁄ r 21⊥ (3.292)
The final basis vector for the element coordinate system is obtained from the cross product
of ê 3 and ê 1 ,
Four-Node
Corotational
Shell
ê 2 = ê 3 × ê 1 . (3.293)
< Go Back
Internal force resultants at the nodes are computed from the stress gradient in the local
element coordinate system (see (3.142) in Basics of Element Design); they are then
transformed to the global coordinate system via
 inte   inte 
 fx   f̂ 1 
  ê 1x ê 2x ê 3x  
SEACAS
 inte   inte 
Library
 fy  = ê 1y ê 2y ê 3y  f̂ 2 . (3.294)
   
 int 
e ê 1z ê 2z ê 3z  int e
 fz   f̂ 3 
Theory    
Manuals
Internal moment resultants are also computed in local element coordinates and
transformed to global coordinates using Equation (3.294). However, Euler’s equations for
the rotational accelerations (see (3.282))are written in nodal coordinates that rotate and
Finite Element
Formulations translate with the node in question (denoted here by superposed bars). The transformation
from global to nodal coordinates is accomplished by
 m1  e 1x e 1y e 1z  mx 
Element    
Examples  m 2  = e 2x e 2y e 2z  my . (3.295)
   
 m3  e 3x e 3y e 3z  mz 
Four-Node Therefore, the complete transformation for moments from the element to nodal
Corotational
Shell coordinates is
< Go Back T
{ m } = [ e ] [ ê ] { m̂ } . (3.296)
Element Equations
SEACAS
Library The shell element is based on a four-node quadrilateral with bilinear interpolations of
midsurface coordinates and of both translational and rotational velocities. Coordinates in
the midsurface of the shell are approximated as
Theory
Manuals x i = x ia N a ( r, s ) , (3.297)
where N a are the shape functions, and r and s are the parent domain coordinates.
Finite Element Repeated indices, a , indicate summation over the four nodes of the element. Similarly the
Formulations velocity of the midsurface and the angular velocity of the normal are interpolated as
v = v a N a ( r, s ) , (3.298)
Element and
Examples
q̇ = q˙a N a ( r, s ) , (3.299)
Four-Node using the same shape functions used for the midsurface coordinates.
Corotational
Shell
In a similar fashion as described for the three-dimensional, constant stress element (see
< Go Back Eight-Node Uniform Strain Element), the shape functions can be expanded in terms of
an orthogonal set of basis vectors as
1 1
φ a = Σ a + --- rΛ 1a + --- sΛ 2a + rsΓ a . (3.300)
2 2
Theory Manuals (2/19/99) Finite Element Formulations - Four-Node Corotational Shell - Element Equations 12
These basis vectors represent deformation modes of a unit square, analogous to their
three-dimensional counterparts shown in Figure 3.14. Rigid body motion is represented
SEACAS by the first vector Σ a . The volumetric basis vectors, Λ 1a and Λ 2a , can be combined to
Library
represent the normal and shear strain in an element. With the reduced integration
formulation, the element area involves only these two vectors. Since the last vector, Γ a , is
neglected in the uniform strain formulation, it represents spurious energy, or hourglass,
Theory
Manuals modes for the element. Substitution of Eqs. (3.298) and (3.299) into Eqs. (3.285) yields
the discretized strain-displacement relationships:
ˆ 1
Finite Element d 11 = --- [ B 1a v̂ 1a + ẑB 1a q̇ 2a ]
Formulations A
ˆ 1
d 22 = --- [ B 2a v̂ 2a – ẑB 2a q̇ 1a ]
A
ˆ 1
Element 2d 12 = --- [ B 2a v̂ 1a + B 1a v̂ 2a + ẑ ( B 2a q̇ 2a – B 1a q̇ 1a ) ] . (3.301)
Examples
A
ˆ 1
2d 23 = --- [ B 2a v̂ 3a – N a q̇ 1a ]
A
Four-Node ˆ 1
Corotational 2d 13 = --- [ B 1a v̂ 3a + N a q̇ 2a ]
Shell A
< Go Back The gradient operator, B αa , is defined as
∂φ a
B αa = ∫ ---------
- dA , (3.302)
A e ∂x̂ α
where Greek subscripts indicate the 1 and 2 coordinates in the plane of the element. For
the reduced integration element, the gradient operator is needed only at the point r = 0 ,
SEACAS s = 0 and can be expressed in closed form in terms of the corotational coordinates of the
Library
element nodes (see [Belytschko, T., Lin, J.I. and Tsay, C.S., 1984]) as
1 ŷ – ŷ 4 ŷ 3 – ŷ 1 ŷ 4 – ŷ 2 ŷ 1 – ŷ 3
Theory [ B αa ] = --- 2 . (3.303)
Manuals 2 x̂ – x̂ x̂ 1 – x̂ 3 x̂ 2 – x̂ 4 x̂ 3 – x̂ 1
4 2
Because the element is a quadrilateral, its area, A, can also be expressed in terms of nodal
Finite Element
coordinates:
Formulations
1
A = --- [ ( x̂ 3 – x̂ 1 ) ( ŷ 4 – ŷ 2 ) + ( x̂ 2 – x̂ 4 ) ( ŷ 3 – ŷ 1 ) ] . (3.304)
2
Element
Examples
Four-Node
Corotational
Shell
< Go Back
Computation of Internal Force and Moment
SEACAS
Library
Resultants
Representation of the rate of internal energy for a shell element in the global variational
principle requires derivation of the internal force and moment resultants. Let the velocity
Theory
Manuals vector for an element node be defined as
 v̂ 1a 
 
Finite Element  2a 
v̂
Formulations
 
v̂ a =  v̂ 3a  (3.305)
 
 1a 
q̇
 
Element  q̇ 2a 
Examples
in corotational coordinates of the element. Then the corresponding internal element force
and moment vector is
Four-Node
Corotational
Shell
< Go Back
Theory Manuals (2/19/99) Finite Element Formulations - Four-Node Corotational Shell - Computation of Internal Force and Moment Resultants 15
ˆ
 1a 
f
ˆ 
SEACAS  f 2a 
Library ˆ  
f a =  fˆ 3a  . (3.306)
 m̂ 
 1a 
 
Theory  m̂ 
2a
Manuals
The concept of stress power, discussed for a general continuum in Stress Power, can be
utilized here to define the vector f̂ ia via
Finite Element
Formulations
ˆd T Tˆ
Element stress power = ∫ T ij v i, j dV = ∫ dV . (3.307)
Ve Ve
Substitution of Eqs. (3.285) into (3.307) and using one-point quadrature leads to the
Element
Examples
following expressions for the element internal forces and moments:
Four-Node
Corotational
Shell
< Go Back
ˆ ˆ ˆ
f 1a = B 1a f 11 + B 2a f 12
ˆ ˆ ˆ
SEACAS f 2a = B 2a f 22 + B 1a f 12
Library
ˆ ˆ ˆ
f 3a = B 1a f 13 + B 2a f 23
A (3.308)
m̂ 1a = – B 2a m̂ 22 – B 1a m̂ 12 – --- fˆ yz
Theory 4
Manuals
A
m̂ 2a = B 1a m̂ 11 + B 2a m̂ 12 + --- fˆ xz
4
m̂ 3a = 0.
Finite Element
Formulations
with
ˆ ˆ
f αβ = ∫ T αβ dẑ
ˆ ˆ
Element f αz = κ ∫ T αz dẑ , (3.309)
Examples
m̂αβ = ẑ Tˆ dẑ
∫ αβ
Four-Node where κ = 5/6 is the shear factor from classical plate theory. The integrals of stress over the
Corotational shell thickness are computed analytically for linear elastic materials. For nonlinear
Shell
materials the internal force and moment resultants are computed by numerical integration
< Go Back through the element thickness (h) using the trapezoid rule. Currently the user may use
either three or five integration points; the first point is at ẑ = – h ⁄ 2 , the middle point is at
ẑ = 0 (the midsurface), the last point is at ẑ = +h ⁄ 2 .
Hourglass Control
SEACAS
Library Reduced integration of the stress divergence leads to spurious zero-energy modes
(hourglass modes) in the element, as discussed, in general, in Hourglass Control. The
hourglass control algorithm implemented for the shell element is one developed by
Theory [Flanagan, D.P. and Belytschko, T., 1981] and by [Belytschko, T., Ong, J.S.-J., Liu,
Manuals W.K. and Kennedy, J.M., 1984]. Removal of the linear portion of the nodal velocity
results in the definition of hourglass shape vectors γ a :
Finite Element 1
Formulations γ a = Γ a – ---B αa x αb Γ b , (3.310)
A
and the corresponding hourglass forces and moments
ˆ HG
Element
Examples
f αa = Q α γ a
ˆ HG
f 3a = Q 3 γ a (3.311)
HG
Four-Node m̂ αa = P α γ a
Corotational
Shell
with the generalized hourglass stresses given by
< Go Back
Theory Manuals (2/19/99) Finite Element Formulations - Four-Node Corotational Shell - Hourglass Control 18
Q α = C 1 γ b v̂ αb
Q 3 = C 2 γ b v̂ 3b , (3.312)
SEACAS
Library
P α = C 3 γ b q̇ αb
and the hourglass stiffnesses defined via

Theory
Manuals Eh
C 1 = r m --------- B αa B αa
8A
3
Gh
C 2 = r z -----------2-B αa B αa , (3.313)
Finite Element
Formulations 12A
3
Eh
C 3 = r θ ------------- B αa B αa
192A
Element where E and G are Young’s modulus and the shear modulus, respectively; and h is the
Examples
element thickness. In the preceding equations the fˆ αa represent the nodal hourglass
HG
membrane forces associated with in-plane velocities v̂ α ; fˆ 3a , the hourglass bending

HG
Four-Node
Corotational HG
Shell forces associated with out-of-plane velocity v̂ 3 ; and m̂ αa , the hourglass bending moments
< Go Back associated with rotational velocities q̇ α . The corresponding hourglass parameters r m , r z ,
and r θ are usually assigned values ranging from 0.01 to 0.05.
Theory Manuals (2/19/99) Finite Element Formulations - Four-Node Corotational Shell - Hourglass Control 19
Calculation of the Stable Time Increment
SEACAS
Library The central difference operator is conditionally stable with the stability limit for a system
with no damping given by
2
Theory ∆t ≤ ------------ , (3.314)
Manuals ω max
where ω max is the maximum frequency of the system (see also Explicit Finite Element
Finite Element
Methods). The maximum frequency in the system can be bounded by the maximum
Formulations element frequency, so the stability limit becomes
2
∆t ≤ ---------------------------------------------
e
-, (3.315)
MAXIMUM ( ω max )
Element
Examples
e
where MAXIMUM ( ω max ) is the maximum element frequency of all the elements in the
system.
Four-Node
Corotational Conservative estimates of the maximum frequency for quadrilateral shell elements were
Shell
developed by Belytschko and Lin [Belytschko, T. and Lin, J.I., 1985] and are given for
< Go Back membrane, bending, and shear deformation as
Theory Manuals (2/19/99) Finite Element Formulations - Four-Node Corotational Shell - Calculation of the Stable Time Increment 20
2 12D 2 2 2
ω max, m = -------------2- ( R 1 + R 1 – 16 ( 1 – ν )A )
MAh
SEACAS
Library
2
2 h 2
ω max, z = ---------- ω max, m (3.316)
12α
2 4  c s α c s A
ω max, θ = --- ---------- + ----------
Theory M A 4α 
Manuals
with
c s = κGh
Finite Element
Formulations
R 1 = 2B αa B αa
2 2
R2 = R 1 – 16A
, (3.317)
Element α1 = ( R1 + R2 ) ⁄ 4
Examples
3
Eh
D = -------------------------
2
12 ( 1 – ν )
Four-Node
Corotational where E is Young’s modulus, G is the shear modulus, ν is Poisson’s Ratio, h is the element
Shell
thickness, M is the element mass, and κ is the shear correction factor.
< Go Back
In the preceding expressions for maximum element frequency, α is a rotational inertia
scaling factor assumed to be
I
α = ------ , (3.318)
Ac
SEACAS
Library where I and A c are the mass moment of inertia and area, respectively, of the cross section.
α is approximated using the value for a flat, rectangular element:
2
Theory h +A
Manuals α = ---------------- . (3.319)
12
The maximum stable time step is computed using the maximum frequency over all shell
Finite Element
and brick elements.
Formulations
Element
Examples
Four-Node
Corotational
Shell
< Go Back
Constitutive Models
SEACAS
Library Two plane stress constitutive models that are widely implemented for the shell element
described here are described next: an elastic model and an elastic-plastic model with
combined linear hardening. In corotational coordinates the stress rate follows directly
Theory from the velocity strain. Since all stress and strain quantities are computed in corotational
Manuals
coordinates, the ˆ notation has been dropped.
For a plane stress, linear elastic model, the stress rate is computed from the velocity strain
Finite Element as
Formulations
Ṫ = λ' ( trd )I + 2µd , (3.320)
where µ is the shear modulus.
Element
E
Examples
µ = G = -------------------- , (3.321)
2(1 + ν)
λ' is the Lamé constant for plane stress

Four-Node
Corotational
Shell Eν
λ' = ------------------- (3.322)
2
< Go Back (1 – ν )
( trd ) = d ii , (3.323)
and I is the second-order identity tensor.

Theory Manuals (2/19/99) Finite Element Formulations - Four-Node Corotational Shell - Constitutive Models 23
Material properties required as input are Young’s modulus E and Poisson’s ratio ν, from
which the above Lamé parameters can be calculated. There are no internal state variables.
SEACAS
Library
In treatment of elastic-plastic materials, on the other hand, it should be first noted that
internal force and moment resultants for nonlinear materials are computed by numerical
integration through the element thickness, so the constitutive model must be evaluated at
each integration point. The general theory for elastic-plastic materials, with combined
Theory
Manuals isotropic and kinematic hardening, is discussed in detail in Constitutive Modeling.
Here we concentrate on the implementational details of plane stress radial return, as they
differ from the fully three-dimensional situation. For plane stress, d33 must be computed
Finite Element from the constraint on the constitutive model rather than directly from the finite element
Formulations
equations. This constraint requires iterations on the radial return algorithm. The secant
approach, presented by [Hallquist, J.O. and Benson, D.J., 1986] and assessed for
accuracy and efficiency by [Whirley, R.G., Hallquist, J.O. and Goudreau, G.L., 1988],
Element is described below.
Examples
Because the plane stress assumption only affects the volumetric strains, the trial shear
stresses can be computed outside the iteration loop. The trial shear stresses are
Four-Node
Corotational T
Shell
T 12 = T 12 + 2µ∆td 12 (3.324)
< Go Back T
T 23 = T 23 + 2µ∆td 23 (3.325)
T
T 23 = T 13 + 2µ∆td 13 , (3.326)
where dij is the rate of deformation at the integration point computed from Eq. (3.285).
The stress difference tensor ξij is defined as the difference between the deviatoric stress
SEACAS tensor Sij and the back stress tensor αij,
Library
ξ ij = S ij – α ij (3.327)
Theory S ij = T ij + pδ ij (3.328)
Manuals
1
p = – --- T ii (3.329)
3
Finite Element
Formulations The shear stress differences are
T T
ξ 12 = T 12 – α 12 (3.330)
T T
Element
Examples ξ 23 = T 23 – α 23 (3.331)
T T
ξ 13 = T 13 – α 13 . (3.332)
Four-Node
Corotational The iteration loop is entered for computation of the volumetric stress difference and the
Shell
equivalent plastic strain. For iteration i the volumetric velocity strain is
< Go Back
( i) ( i)
( trd ) = d 11 + d 22 + d 33 . (3.333)
The trial normal stresses and pressure then follow as:
T ( i)
T 11 = T 11 + λ∆t ( trd ) + 2µ∆td 11 (3.334)
SEACAS T ( i)
Library T 22 = T 22 + λ∆t ( trd ) + 2µ∆td 22 (3.335)
T ( i)
T 33 = λ∆t ( trd ) + 2µ∆td 33 (3.336)
Theory
Manuals T 1 T T T
p = – --- ( T 11 + T 22 + T 33 ) . (3.337)
3
The first two normal components of the stress difference computed from Eq. (3.327) are:
Finite Element
Formulations
T T T
ξ 11 = T 11 + p – α 11 (3.338)
T T T
ξ 22 = T 22 + p – α 22 . (3.339)
Element
Examples
Because the stress difference tensor is deviatoric, the third normal component is given by
T T T
ξ 33 = – ( ξ 11 + ξ 22 ) . (3.340)
Four-Node
Corotational
Shell The increment in equivalent plastic strain for the ith iteration is
< Go Back
( i) 1 T T
∆γ = ----------------------------- ( ξ ij ξ ij – R ) , (3.341)
2µ  1 + ------
H'
 3µ
where R is the radius of the yield surface from the previous time step,
SEACAS
R = ξ ij ξ ij . (3.342)
Library
The hardening slope H' depends on Young’s modulus and the plastic modulus Ep,
E pE
Theory H' = ---------------- . (3.343)
Manuals E – Ep
The normal stress difference in the thickness direction follows from the radial return
algorithm as
Finite Element
Formulations
( i) T T T T
T 33 = T 33 – 2µ∆γβξ 33 ⁄ ξ ij ξ ij , (3.344)
where β is a scalar parameter ranging from 0 to 1 that determines the relative amounts of
Element isotropic and kinematic hardening. For β = 1 all hardening is assumed to be isotropic. At
Examples
the other extreme β = 0 means only kinematic hardening is present. Finally, an estimate
(i + 1)
for d 33 is obtained from a secant update,
Four-Node
Corotational (i + 1) (i – 1) ( i) ( i) (i – 1) ( i) (i – 1)
Shell d 33 = d 33 + T 33 ( d 33 – d 33 ) ⁄ ( T 33 – T 33 ). (3.345)
< Go Back
The (i+1) iteration proceeds by substituting the new value of d33 into (3.333) and
repeating Eqs. (3.334) through (3.344) until σ33 has converged to zero.
Because the plane stress algorithm is based on a secant update, two starting values are
required for d33. The two starting values are
SEACAS
Library (0) ν
d 33 = – ----------------- ( d 11 + d 22 ) , (3.346)
(1 – ν)
which assumes a completely elastic step, and
Theory
Manuals (1)
d 33 = – ( d 11 + d 22 ) (3.347)
based on a completely plastic (incompressible) step.

Finite Element
Formulations Once the algorithm has converged, the yield surface radius, equivalent plastic strain, and
back stresses and are updated as
t + ∆t t 2
R = R + --- βH'∆γ (3.348)
Element 3
Examples
t + ∆t t
pl pl 2
d = d + --- ∆γ (3.349)
3
Four-Node
Corotational
t + ∆t t 2
Shell α ij = α ij + --- ( 1 – β )H'γ ξ ij . (3.350)
3
< Go Back
Finally, the stresses are updated using radial return
t + ∆t T T T T
α ij = T ij – 2µ∆γβξ ij ⁄ ξ kl ξ kl . (3.351)
Time-Stepping Algorithm
SEACAS
Library In a central difference algorithm to integrate the equations of motion, such as is used in
explicit dynamics, the translational variables are handled just as they would be for an
ordinary continuum. Once the nodal accelerations at time t are solved from Eq. (3.281),
Theory the velocities and displacements in global coordinates follow as
Manuals
∆t ∆t
t + ------- t – -------
2 2 t
v = v + ∆ta
. (3.352)
Finite Element ∆t
t + -------
t + ∆t t 2
Formulations
d = d + ∆tv
Angular accelerations at time t for nodes connected to shell elements are computed in
coordinate systems that are local to each node. Because these coordinate systems rotate
Element with the nodes, the angular accelerations cannot be integrated directly for updating the
Examples
configuration. Instead the procedure of Belytschko and coworkers (see [Belytschko, T.,
Lin, J.I. and Tsay, C.S., 1984]) is implemented.
Four-Node The nodal rotations are updated from the angular accelerations in the same manner as the
Corotational
Shell translational velocities,
< Go Back ∆t ∆t
t + ------- t – -------
2 2 t
ω = ω + ∆tα . (3.353)
The updated nodal basis vectors are then
Theory Manuals (2/19/99) Finite Element Formulations - Four-Node Corotational Shell - Time-Stepping Algorithm 29
t + ∆t t
ei = e i + ∆tė i
∆t . (3.354)
SEACAS
t  t + ------
2
-
t
Library = e i + ∆t  ω × e i
 
t + ∆t
For e 3 , the preceding equation becomes
Theory
Manuals
∆t ∆t
t + ∆t t  t + ------
2 t
- t + ------- 
2 t
e3 = e 3 + ∆t  ω y e1 – ωx e 2 . (3.355)
 
Finite Element
Formulations
t + ∆t t
The scalar product of e 3 and e 1 gives
∆t
t + -------
t + ∆t 2
Element e3 = ∆tω y . (3.356)
Examples x
t + ∆t t t + ∆t
Similarly the scalar product of e 3 and e 2 gives the y component of e 3 .
Four-Node
Corotational ∆t
Shell t + -------
t + ∆t 2
e3 = – ∆ tω x . (3.357)
< Go Back y
t + ∆t
The third component of e 3 is found by normality. The rotation over a single time step
is assumed to be smal,l so second-order terms are dropped; and the z component of the
SEACAS
Library updated vector is
t + ∆t t + ∆t 2 t + ∆t 2
e3 = 1 – ( e3 ) – ( e3 ) . (3.358)
z x y
Theory
Manuals
Next e 1 is updated. From Equation (3.354)
∆t ∆t
t + ∆t t  t + ------
2 t
- t + ------- 
2 t
Finite Element
Formulations
e1 = e 1 + ∆t  ω z e2 – ωy e 3 . (3.359)
 
t
The scalar product of Equation (3.359), with e 2 , gives the y component of e 1 as
Element
t + ∆t t + ∆t
Examples
e1 = ∆tω z . (3.360)
y
t + ∆t
Four-Node If it is assumed that e 1 is approximately one (i.e., small rotations over the time step)
x
Corotational
Shell t + ∆t t + ∆t
orthogonality of e 1 and e 3 yields
< Go Back
t + ∆t t + ∆t t + ∆t
( e3 + e1 e3 )
t + ∆t x y y
e1 = – --------------------------------------------------------------- , (3.361)
z t + ∆t
e3
z
t + ∆t
again neglecting the second-order terms. The third component of e 1 is then found by
imposing unit length on the vector, so
SEACAS
Library
t + ∆t t + ∆t 2 t + ∆t 2
e1 = 1– ( e1 ) – ( e1 ) . (3.362)
x y z
Theory t + ∆t t + ∆t t + ∆t
Manuals
The cross product of e 3 with e 1 gives e 2 to complete the update of the
basis vectors for the nodal coordinate system.
Finite Element
Formulations
Element
Examples
Four-Node
Corotational
Shell
< Go Back
Adley, M.D. and Moxley, R.E., 1996 PENCURV/ABAQUS: A Simply Coupled
< Go Back
Penetration Trajectory/Structural Dynamics Model for Deformable Projectiles
Impacting Complex Curvilinear Targets, Technical Report SL-96-6, U.S. Army
Engineer Waterways Experiment Station, Vicksburg, MS.
American National Standards Institute, 1978 American National Standards

Programming Language FORTRAN, ANSI X3.9-1978, American National
Standards Institute, New York.
American Society for Testing and Materials, 19xx Annual Book of ASTM Standards:
Section 3 - Metals Test Methods and Analytical Procedures, American Society for
Testing and Materials, Philadelphia, Pennsylvania.
Attaway, S.W., 1990 Update of PRONTO 2D and PRONTO 3D Transient Solid

Dynamics Program, SAND90-0102, Sandia National Laboratories, Albuquerque,
NM.
Balmer, H.A. and Witmer, E.A., 1964 Theoretical-Experimental Correlation of Large

Dynamic and Permanent Deformation of Impulsively Loaded Simple Structures,
FDP-TDR-64-108, Air Force Flight Dynamics Laboratory, Wright-Patterson AFB,
OH.
Bartels, R. and Daniel, J.W., 1973 “A conjugate gradient approach to nonlinear elliptic
boundary value problems in irregular regions,” in Lecture Notes in Mathematics, (A.
Dold and B. Eckmann, eds.), Springer-Verlag, New York, pp. 1-11.
SEACAS Library (2/19/99) Bibliography 1

Bathe, K.J. and Wilson, E.L., 1976 Numerical Methods in Finite Element Analysis,
< Go Back
Prentice-Hall, Inc., New Jersey.
Beisinger, Z.E., 1984 SEACO: Sandia Engineering Analysis Department Code

Output Data Base, SAND84-2004, Sandia National Laboratories, Albuquerque, NM.
Belytschko, T. and Bindelman, L.P., 1993 “Assumed strain stabilization of the eight node
hexahedral element,” Compter Methods in Applied Mechanics and Engineering,
Vol. 105, pp. 225-260.
Belytschko, T. and Lin, J.I., 1985 “Eigenvalues and Stable Time Step for the Bilinear
Mindlin Plate Element,” International Journal for Numerical Methods in
Engineering, Vol. 21, pp. 1729-1745.
Belytschko, T. and Lin, J.I., 1987 “A Three-Dimensional Impact-Penetration Algorithm

with Erosion,” Computers and Structures, Vol. 25, No. 1, pp. 95-104.
Belytschko, T., Lin, J.I. and Tsay, C.S., 1984 “Explicit Algorithms for the Nonlinear
Dynamics of Shells,” Computer Methods in Applied Mechanics and Engineering,
Vol. 42, pp. 225-251.
Belytschko, T. and Marchertas, A.H., 1974 “Nonlinear Finite Element Method for Plates
and its Application to the Dynamic Response of Reactor Fuel Subassemblies,” Journal
of Pressure Vessel Technology, ASME, pp. 251-157.

Belytschko, T. and Neal, M.O., 1991 “Contact-Impact by the Pinball Algorithm with
< Go Back
Penalty and Lagrangian Methods,” International Journal Numerical Methods of
Belytschko, T., Ong, J.S.-J., Liu, W.K. and Kennedy, J.M., 1984 “Hourglass Control in
Linear and Nonlinear Problems,” Computer Methods in Applied Mechanics and
Belytschko, T., Schwer, L. and Klein, M.J., 1977 “Large Displacement Transient
Analysis of Space Frames,” International Journal for Numerical Methods in
Benson, D.J. and Hallquist, J.O., 1990 “A Single Surface Contact Algorithm for the Post-
Buckling Analysis of Structures,” Computer Methods in Applied Mechanics and
Benz, W., 1988 “Applications of Smooth Particle Hydrodynamics (SPH) to Astrophysical

Problems,” Computational Physics Commuinications, Vol. 48, pp. 97-105.
Benz, W., 1990 “Smooth Particle Hydrodynamics: A Review,” in The Numerical

Modeling of Stellar Pulsation, (J.R. Buchler, ed.), Kluwer, Dordrecht, Netherlands, p.
269.
Bergmann, V.L., 1991 Transient Dynamic Analysis of Plates and Shells with
PRONTO 3D, SAND91-1182, Sandia National Laboratories, Albuquerque, NM.

Bertsekas, D.P., 1982 Constrained Optimization and Lagrange Multiplier Methods,
< Go Back
Academic Press, New York.
Biffle, J.H., 1981 JAC--A Two-Dimensional Finite Element Computer Program for
the Nonlinear Response of Solids with the Conjugate Gradient Method, SAND81-
0998, Sandia National Laboratories, Albuquerque, NM.
Biffle, J.H., 1984 JAC - A Two-Dimensional Finite Element Computer Program for
the Nonlinear Quasi-Static Response of Solids with the Conjugate Gradient
Method, SAND81-0998, Sandia National Laboratories, Albuquerque, NM.
Biffle, J.H., 1993 JAC3D - A Three-Dimensional Finite Element Computer Program

for the Nonlinear Quasi-Static Response of Solids with the Conjugate Gradient
Method, SAND87-1305, Sandia National Laboratories, Albuquerque, NM.
Biffle, J.H. and Gubbels, M.H., 1976 WULFF - A Set of Computer Programs for the
Large Displacement Dynamic Response of Three Dimensional Solids, SAND76-
0096, Sandia National Laboratories, Albuquerque, NM.
Biffle, J.H. and Stone, C.M., 1989 “Personal communications,” Applied Mechanics
Division, Sandia National Laboratories, Albuquerque, NM.
Bishop, R.F., Hill, R. and Mott, N.F., 1945 “The theory of indentation and hardness”,
Proceedings of the Royal Society, Vol. 57 (3), pp. 147-159.

Blacker, T.D., 1988 FASTQ Users Manual Version 1.2, SAND88-1326, Sandia National
< Go Back
Laboratories, Albuquerque, NM.
Boaz, M.L., 1966 Mathematical Methods in the Physical Sciences, John Wiley and
Sons, New York, NY.
Budiansky, B. and O’Connell, R.J., 1976 “Elastic moduli of a cracked solid,” Computer
Methods in Applied Mechanics and Engineering, Vol. 12, pp. 81-97.
CRAY Research, Inc., 1989 Volume 3: UNICOS Math and Scientific Library
Reference Manual, SR-2081 5.0.
Campbell, P.M., 1988 Some New Algorithms for Boundary Values Problems in
Smooth Particle Hydrodynamics, DNA-TR-88-286.
Carpenter, N.J., Taylor, R.L. and Katona, M.G., 1991 “Lagrange Constraints for Transient
Finite Element Surface Contact,” International Journal for Numerical Methods in
Engineering, Vol. 32, 103-128.
Chait, R., 1972 “Factors influencing the strength differential of high strength steels”,
Metalurgical Transactions A, Vol. 3, pp. 365-371.
Chakarabarty, J., 1987 Theory of Plasticity, McGraw-Hill Book Company, New York,
pp. 306-315, 342-350.

Chaudhary, A.B. and Bathe, K.J., 1986 “A Solution Method for Static and Dynamic
< Go Back
Analysis of Three-Dimensional Contact Problems with Friction,” Computers and
Structures, Vol. 24, No. 6, pp. 855-873.
Cheung, C.Y. and Cebon, D., 1997 “Experimental Study of Pure Bitumens in Tension,
Compression, and Shear,” Journal of Rheology, Vol. 41, No. 1, pp. 45-73.
Chou, T.S., 1989 The Dynamic Response at High-Explosive Inert-Solid Interface as

Predicted by Finite Difference/Element Computer Programs, Tech. Rep. MLM-
3562, EG&G Mound Applied Technologies, Miamisburg, OH.
Cloutman, L.D., 1990(a) Basics of Smoothed Particle Hydrodynamics, Lawrence

Livermore National Laboratory, Livermore, CA, report UCRL-ID-103698.
Cloutman, L.D., 1990(b) “An Evaluation of Smoothed Particle Hydrodynamics,”

Proceedings of The NEXT Free-Lagrange Conference, Jackson Lake Lodge,
Moran, WY, June 3-7.
Cohen, M. and Jennings, P.C., 1983 “Silent Boundary Methods,” in Computational

Methods for Transient Analysis, (Belytschko, T. and Hughes, T.J.R., eds.), North-
Holland.
Cole, R.H., 1965 UnderWater Explosions, Dover Publications.

Concus, P., Golub, G.H. and O’Leary, D.P., 1976 “A generalized conjugate gradient
< Go Back
method for the numerical solution of elliptic partial differential equations,” in Sparse
Matrix Computations, (J.R. Bunch and D.J. Rose, eds.), Academic Press, New York,
pp. 309-332.
Courant, R., Friedrichs, K.O. and Lewy, H., 1928 Mathematical Annotations, Vol. 100,
p. 32.
Crismann, J.M. and Zapas, L.J., 1979 “Creep Failure and Fracture of Polyethylene in
Uniaxial Extension,” Polymer Engineering Science, Vol. 19, pp. 99-103.
Curnier, A. and Alart, P., 1988 “A Generalized Newton Method for Contact Problems
with Friction,” Journal de Mecanique Theorique et Appliquee, Vol. 7, 67-82.
Daniel, J.W., 1967 “The conjugate gradient method for linear and nonlinear operator
equations,” SIAM Journal of Numerical Analysis, Vol. 4, no. 1, pp. 10-26.
Dennis, J.E. and Schnabel, R.B., 1996 Numerical Methods for Unconstrained
Optimization and Nonlinear Equations, Society for Industrial and Applied
Mathematics, Philadelphia.
Dienes, J.K., 1979 “On the analysis of rotation and stress rate in deforming bodies,” Acta
Mechanica, Vol. 32, pp. 217-232.

Dobratz, B.M., 1981 LLNL Explosives Handbook, Properties of Chemical Explosives
< Go Back
and Explosive Simulants, UCRL-52997, Lawrence Livermore National Laboratory,
Livermore, CA.
Duffey, T.A. and Macek, R.W., 1997 “Non-normal impact of earth penetrators”,
Proceedings of the International Symposium on Penetration and Impact Problems
(ICES'97), San Jose, Costa Rica.
Engleman, R. and Jaeger, Z., 1987 Theoretical Aids for Improvement of Blasting
Efficiencies in Oil Shale and Rocks, Tech. Rep AP-TR-12/87, Soreq Nuclear
Research Center, Yavne, Israel.
Ferry, 1950 ???
Flanagan, D.P. and Belytschko, T., 1981 “A Uniform Strain Hexahedron and
Quadrilateral with Orthogonal Hourglass Control,” International Journal for
Numerical Methods in Engineering, Vol. 17, pp. 679-607.
Flanagan, D.P. and Belytschko, T., 1984 “Eigenvalues and Stable Time Steps for the
Uniform Hexahedron and Quadrilateral,” Journal of Applied Mechanics, Vol. 51, pp.
35-40.
Flanagan, D.P., Mills-Curran, W.C. and Taylor, L.M., 1986 SUPES - A Software
Utilities Package for the Engineering Sciences, SAND86-0911, Sandia National

Flanagan, D.P. and Taylor, L.M., 1987 “An accurate numerical algorithm for stress
< Go Back
integration with finite rotations,” Computer Methods in Applied Mechanics and
Flanagan, D.P. and Taylor, L.M., 1987 “On the Analysis of Rotation and Stress Rate in
Deforming Bodies,” Computer Methods in Applied Mechanics and Engineering,
Vol. 62, pp. 305-320.
Fletcher, R. and Reeves, C.M., 1964 “Function minimization by conjugate gradients,”

The Computer Journal, Vol. 7, pp. 149-154.
Forrestal, M.J., Altman, B.S., Cargile, J.D. and Hanchak, S.J., 1994 “An empirical
equation for penetration depth of ogive-nose projectiles into concrete targets”,
International Journal of Impact Engineering, Vol. 15, pp. 395-405.
Forrestal, M.J., Brar, N.S. and Luk, V.K., 1991 “Penetration of strain-hardening targets
with rigid spherical-nose rods”, ASME Journal of Applied Mechanics, Vol. 58, pp.
7-10.
Forrestal, M.J. and Luk, V.K., 1991 “Penetration into soil targets”, International Journal
of Impact Engineering, Vol. 12, pp. 427-444.
Forrestal, M.J., Okajima, K. and Luk, V.K., 1988 “Penetration of 6061-T651 aluminum
targets with rigid long rods”, ASME Journal of Applied Mechanics, Vol. 55, pp. 755-
760.

Forrestal, M.J. and Tzou, D.Y., 1996 “A spherical cavity-expansion penetration model for
< Go Back
concrete targets”, International Journal of Solids Structures (accepted).
Forrestal, M.J., Tzou, D.Y., Askari, E. and Longcope, D.B., 1995 “Penetration into ductile
metal targets with rigid spherical-nose rods”, International Journal of Impact
Freudenthal, A.M., 1963 ???
Frew, D.J., Hanchak, S.J., Green, M.L. and Forrestal, M.J., 1997 “Penetration of concrete
targets with ogive-nose steel rods”, International Journal of Impact Engineering,
(submitted).
Frost and Ashby, 19?? Deformation Mechanisms Maps, ???.
Fung, Y.C., 1965 Foundations of Solid Mechanics, Prentice-Hall, Englewood Cliffs,

New Jersey.
Fung, Y.C., 1977 A First Course in Continuum Mechanics, 2nd edition, Prentice-Hall,
Englewood Cliffs, New Jersey.
Gallego, F.J. and Anza, J.J., 1989 “A Mixed Finite Element Model for the Elastic Contact
Problem,” International Journal for Numerical Methods in Engineering, Vol. 28,
1249-1264.

Giannakopoulos, A.E., 1989 “The Return Mapping Method for the Integration of Friction
< Go Back
Constitutive Relations,” Computers and Structures, Vol. 32, 157-167.
Gilkey, A.P., 1988 ALGEBRA - A Program that Algebraically Manipulates the

Output of a Finite Element Analysis (EXODUS Version), SAND88-1431, Sandia
National Laboratories, Albuquerque, NM.
Gilkey, A.P. and Glick, J.H., 1989 BLOT - A Mesh and Curve Plot Program for the
Output of a Finite Element Analysis, SAND88-1432, Sandia National Laboratories,
Albuquerque, NM.
Gilkey, A.P. and Sjaardema, G.D., 1989 GEN3D: A GENESIS Database 2D to 3D

Transformation Program, SAND89-0485, Sandia National Laboratories,
Albuquerque, NM.
Gingold, R.A. and Monaghan, J.J., 1982 “Kernel Estimates as a Basis for General
Particle Methods in Hydrodynamics,” Journal Computational Physics, Vol. 46, pp.
429-453.
Glowinski, R. and LeTallec, P., 1989 Augmented Lagrangian and Operator Splitting
Methods in Nonlinear Mechanics, SIAM, Philadelphia.
Goodier, J.N. 1965 “On the mechanics of indentation and cratering in the solid targets of
strain-hardening metal by impact of hard and soft spheres”, Proceedings of the 7th
Symposium on Hypervelocity Impact III, pp. 215-259.

Grady, D., 1983 “The Mechanics of Fracture Under High-rate Stress Loading,” in
< Go Back
William Prager Syposium on Mechanics of Geomaterials: Rocks, Concretes and
Soils, (Bazant, Z.P., ed.).
Guenther, C., Hicks, D.L. and Swegle, J.W., 1994 Conservative Smoothing versus
Artificial Viscosity, SAND94-1853, Sandia National Laboratories, Albuquerque, NM.
Gurtin, M.E., 1981 An Introduction to Continuum Mechanics, Academic Press, Inc.
Hallquist, J.O., 1981 User’s Manual for DYNA3D and DYNAP, Lawrence Livermore
National Laboratory, Livermore, CA.
Hallquist, J.O., 1982 User’s Manual for DYNA2D - An Explicit Two-Dimensional

Hydrodynamic Finite Element Code with Interactive Rezoning, Lawrence
Livermore National Laboratory, Livermore, CA.
Hallquist, J.O., 1984 NIKE3D: An Implicit, Finite Deformation, Finite Element Code
for Analyzing the Static and Dynamic Response of Three Dimensional Solids,
UCID-18822, Lawrence Livermore National Laboratory, Livermore, CA.
Hallquist, J.O., 1984 User’s Manual for DYNA2D: An Explicit Two-Dimensional

Hydrodynamic Finite Element Code with Interactive Rezoning, Rev. 2, UCID-
18756, Lawrence Livermore National Laboratory, Livermore, CA.

Hallquist, J.O., 1986 NIKE2D: A Vectorized Implicit, Finite Deformation Finite
< Go Back
Element Code for Analyzing the Static and Dynamic Response of 2-D Solids with
Interactive Rezoning and Graphics, UCID-19677, Lawrence Livermore National
Laboratory, Livermore, CA.
Hallquist, J.O., Goudreau, G.L. and Benson, D.J., 1985 “Sliding Interfaces with Contact-
Impact in Large-Scale Lagrangian Computations,” Computer Methods in Applied
Mechanics and Engineering, Vol. 51, 107-137.
Hallquist, J.O. and Benson, D.J., 1986 “A comparison of an Implicit and Explicit
Implementation of the Hughes-Liu Shell,” in Finite Element Methods for Plate and
Shell Structures, Volume 1: Element Technology, (Hughes, T.J.R. and Hinton, E.,
eds.), Pineridge Press International, Swansea, U.K., pp. 394-430.
Hallquist, J.O. and Benson, D.J., 1987 User’s Manual for DYNA3D: Nonlinear
Dynamic Analysis of Structures, Rev. 3, UCID-19592, Lawrence Livermore National
Harding, D.C., Attaway, S.W., Neilsen, J., Blacker, T.D. and Pierce, J., 1992 Evaluation
of Four Multiple Package Crush Environment to the Common Package, Model 1,
Plutonium Air Transport Container, SAND92-0278, Sandia National Laboratories,
Albuquerque, NM.
Harlow, F.H., 1988 “PIC and its Progeny,” Computational Physics Communications,
Vol. 48, pp. 1-10.

Hearmon, R. F. S., 1961 An Introduction to Applied Anisotropic Elasticity, Oxford
< Go Back
University Press, Oxford, England.
Herrman and Peterson, 1968 ???
Heinstein, M.W., Attaway, S.W., Mellow, F.J. and Swegle, J.W., 1993 A General-
Purpose Contact Detection Algorithm for Nonlinear Structural Analysis Codes,
SAND92-2141, Sandia National Laboratories, Albuquerque, NM.
Hestenes, M.R. and Stiefel, E., 1952 “Methods of conjugate gradients for solving linear
systems,” Journal of Research of the National Bureau of Standards, Vol. 49, pp.
409-436.
Hibbitt, Karlsson and Sorensen, Inc., 1989 ABAQUS Users Manual, Version 4-8,
Hibbitt, Karlsson and Sorensen, Inc., Providence, Rhode Island.
Hibbitt, Karlsson and Sorensen, Inc., 1992 Contact Calculations with ABAQUS -
ABAQUS Explicit Users Manual, Hibbitt, Karlsson and Sorensen, Inc., Providence,
Rhode Island.
Hicks, D.L., 1978 “Stability Analysis of WONDY (A Hydrocode Based on the Artificial
Viscosity Method of von Neumann and Richtmyer) for a Special Case of Maxwell’s
Law,” Mathematics of Computation, Vol. 32, pp. 1123-1130.

Hilber, Hughes and Taylor, 1977 “Improved Numerical Dissipation for Time Integration
< Go Back
Algorithms in Structural Dynamics,” Earthquake Engineering and Structural
Dynamics, Vol. 5, pp.283-292.
Hill, R., 1948 A Theory of Earth Movement Near a Deep Underground Explosion,
Memo No. 21-48, Armament Research Establishment, Fort Halstead, Kent, UK.
Hill, R., 1950 The Mathematical Theory of Plasticity, Oxford University Press,
London.
Holcomb, D.J., 1985 “Results of the 55 Day Consolidation Test,” Internal Memorandum
to J. Stormont, 6332, Sandia National Labs, Albuquerque, NM.
Holcomb, D.L and Hannum, D.W., 1982 Consolidation of Crushed Salt Backfill Under
Conditions Applicable to the WIPP Facility, SAND82-0630, Sandia National Labs,
Albuquerque, NM.
Holcomb, D.J. and Shields, M.F., 1987 Hydrostatic Consolidation of Crushed Salt
with Added Water, SAND87-1990, Sandia National Labs, Albuquerque, NM.
Holden, J.T., 1972 “On the finite deflections of thin beams,” International Journal of
Solids and Structures, Vol. 8, pp. 1051-1055.
Hopkins, H.G., 1960 “Dynamic expansion of spherical cavities in metals”, Progress in

Solid Mechanics, Vol. 1, (Editors I. Sneddon and R. Hill), North Holland, New York,
pp. 85-164.
Huang, J., 1969 “Transient Interaction of Plane Acoustic Waves with a Spherical Elastic
< Go Back
Shell,” Journal Acoustical Society of America, Vol. 45, pp. 661-670.
Hughes, T.J.R., 1987 The Finite Element Method: Linear Static and Dynamic Finite
Element Analysis, Prentice-Hall, Inc., Englewood Cliffs, New Jersey.
Hughes, T.J.R. and Winget, J., 1980 “Finite rotation effects in numerical integration of
rate constitutive equations arising in large-deformation analysis,” International
Journal for Numerical Methods in Engineering, Vol. 15, No. 12, pp. 1862-1867.
Ide,Y. and White, J.L., 1977 “Investigation of Failure During Elongational Flow of
Polymer Melts,” Journal of Non-Newtion Fluid Mechanics, Vol. 2, pp. 281-298.
Ide,Y. and White, J.L., 1979 “Experimental Study of Elongational Flow and Failure of
Polymer Melts,” Journal of Applied Polymer Science,Vol. 22, pp. 1061-1079.
International Nickel Company, Inc., The, 1964 18% Nickel Maraging Steels, New York,
NY.
Irons, B. and Elsawaf, A., 1977 “The conjugate Newton algorithm for solving finite
element equations,” Formulations and Computation Algorithms in Finite Element
Analysis, (K.J. Bathe and W. Wunderlich, eds.), MIT Press, pp. 655-672.
Johnson, G.C. and Bammann, D.J., 1984 “A discussion of stress rates in finite
deformation bodies,” International Journal of Solids and Structures, Vol. 20, no. 8,
pp. 725-737.
Johnson, G.C. and Bammann, D.J., 1984 “On the Analysis of Rotation and Stress Rate in
< Go Back
Deforming Bodies,” International Journal of Solids and Structures, Vol. 20, No. 8,
pp. 725-737.
Johnson, G.R., 1983 Development of Strength and Fracture Models for

Computations Involving Severe Dynamic Loading, Vol. I: Strength and Fracture
Models, Tech. Rep. AFATL-TR-83-05, Air Force Armament Laboratory.
Johnson, G.R., 1988 “Implementation of simplified constitutive models in large computer

codes,” in Dynamic Constitutive/Failure Models, (A. M. Rajendran and T. Nichols,
eds.), pp. 409-426, Dayton, OH.
Johnson, G.R. and Holmquist, T.J., 1989 Test Data and Computational Strength and
Fracture Model Constants for 23 Materials Subjected to Large Strains, High
Strain Rates, and High Temperatures, Tech. Rep. LA-11463-MS, Los Alamos
National Laboratory, Los Alamos, NM.
Johnson, G.R., Stryk, R.A., Holmquist, T.J. and Beissel, S.R., 1996 User Instructions
for the 1996 version of the EPIC Code, Alliant Techsystems Inc., Hopkins, MN.
Johnson, G.R., Stryk, R.A., Holmquist, T.J. and Beissel, S.R., 1997 EPIC 97 (software to
be released), Alliant Techsystems Inc., Hopkins, MN.
Johnson, J.B., 1989 “Personal communication,” USACRREL, Building 4070, Ft.

Wainwrite, Alaska.

Ju, J.-W. and Taylor, R.L., 1988 “A Perturbed Lagrangian Formulation for the Finite
< Go Back
Element Solution of Nonlinear Frictional Contact Problems,” Journal of Theoretical
and Applied Mechanics, Vol. 7, 1-14.
Kamei, E. and Onogi, 1975 “Extensional and Fractural Properties of Monodisperse

Polystyrenes at Elevated Temperatures,” Applied Polymer Symposium, Vol 27, pp.
19-46.
Kaplan, W., 1959 Advanced Calculus, Addison-Wesley Publishing, Reading MA.
Kelley, C.T., 1995 Iterative Methods for Linear and Nonlinear Equations, Society for
Industrial and Applied Mathematics, Philadelphia.
Kerighan, B.W. and Ritchie, D.M., 1978 The C Programming Language, Prentice-Hall,
Inc., New Jersey.
Key, S.W., Beisinger, Z.E. and Kreig, R.D., 1978 HONDO II -A Finite Element
Computer Program for the Large Deformation Dynamics Response of
Axisymmetric Solids, SAND78-0422, Sandia National Laboratories, Albuquerque,
NM.
Kikuchi, N. and Song, Y.J., 1981 “Penalty/Finite Element Approximations of a Class of

Unilateral Problems in Linear Elasticity,” Quarterly of Applied Mathematics, Vol.
39, 1-22.

Kikuchi, N. and Oden, J.T., 1988 Contact Problems in Elasticity: A Study of
< Go Back
Variational Inequalities and Finite Element Methods, SIAM, Philadelphia.
Kipp, M.E. and Grady, D.E., 1978 Numerical Studies of Rock Fragmentations,
Kipp, M.E., Grady, D.E. and Chen, E.P., 1980 “Strain-rate Dependent Fracture Initiation,”
International Journal of Fracture, Vol. 16, pp. 471-478.
Kipp, M.E. and Lawrence, R.J., 1982 WONDY V - A One-Dimensional Finite-

Difference Wave Propagation Code, SAND81-0930, Sandia National Laboratories,
Albuquerque, NM.
Krieg, R.D., 1978 A Simple Constitutive Description for Soils and Crushable Foams,
SC-DR-72-0883, Sandia National Laboratories, Albuquerque, NM.
Krieg, R.D. and Key, S.W., 1976 “Implementation of a Time Dependent Plasticity Theory
into Structural Computer Programs,” Constitutive Equations in Viscoplasticity:
Computational and Engineering Aspects, Vol. 20, ASME, New York, p. 125.
Kreig, R.D. and Krieg, D.B., 1977 “Accuracies of numerical solution methods for the
elastic-perfectly plastic model,” ASME Journal of Pressure Vessel Technology, Vol.
99, pp. 510-515.

Kuszmaul, J.S., 1987 “A New Constitutive Model for Fragmentation of Rock Under
< Go Back
Dynamic Loading,” in Proceedings of the Second International Symposium on
Fragmentation by Blasting, pp. 412-423, Keystone, CO.
Kuszmaul, J.S., 1987 “A Technique for Predicting Fragmentation and Fragment Sizes
Resulting from Rock Blasting,” in Proceedings of the 28th U.S., Symposium on
Rock Mechanics, Tucson, AZ.
Lai, W.M, Rubin, D. and Krempl, E., 1993 Introduction to Continuum Mechanics, 3rd
edition, Pergamon Press.
Laursen, T.A., 1994 “The Convected Description in Large Deformation Frictional

Contact Problems,” International Journal of Solids and Structures, Vol. 31, pp. 669-
681.
Laursen, T.A. and Simo, J.C., 1993 “A Continuum-Based Finite Element Formulation for
the Implicit Solution of Multibody, Large Deformation Frictional Contact Problems,”
International Journal for Numerical Methods in Engineering, Vol. 36, 3451-3485.
Laursen, T.A. and Simo, J.C., 1993 “Algorithmic Symmetrization of Coulomb Frictional
Problems Using Augmented Lagrangians,” Computer Methods in Applied
Mechanics and Engineering, Vol. 108, 133-146.
Leaderman, 1943 ???

Lee, E., Finger, M. and Collins, W., 1973 JWL Equation of State Coefficients for High
< Go Back
Explosives, UCID-16189, Lawrence Livermore National Laboratory, Livermore, CA.
Lee, Radok and Woodward, 1959 ???
Lenard, M.L., 1976 “Convergence conditions for restarted conjugate gradient methods
with inaccurate line searches,” Mathematical Programming, Vol. 10, pp. 32-51.
Libersky, L. and Petschek, A.G., 1990 “Smooth Particle Hydrodynamics with Strength of
Materials,” Proceedings of The NEXT Free-Lagrange Conference, Jackson Lake
Lodge, Moran, WY, June 3-7.
Longcope, D.B., 1991 Coupled Bending/Lateral Load Modeling of Earth

Penetrators, SAND90-0789, Sandia National Laboratories, Albuquerque, NM.
Longcope, D.B., 1996 Oblique Penetration Modeling and Correlation with Field
Tests into a Soil Target, SAND96-2239, Sandia National Laboratories, Albuquerque,
NM.
Longcope, D.B. and Tabbara, M.R., 1998 Modeling of Oblique Earth Penetration
Using Cavity Expansion Loading with Surface Effects, SAND98-xxxx, Sandia
National Laboratories, Albuquerque, NM, in preparation.
Lovejoy, S.C. and Whirley, R.G., 1990 DYNA3D Example Problem Manual, UCRL-
MA-105259, Lawrence Livermore National Laboratory, Livermore, CA.

Luenberger, D.G., 1984 Linear and Nonlinear Programming, 2nd ed., Addison-
< Go Back
Wesley, Reading, MA.
Lucy, L.B., 1977 “A Numerical Approach to the Testing of the Fission Hypothesis,”
Astrophysics Journal, Vol. 82, pp. 1013-1024.
Luk, V.K. and Piekutowski, A.J., 1991 “An analytical model on penetration of eroding
long rods into metallic targets,” International Journal of Impact Engineering, Vol.
11, pp. 323-340.
Lysmer, J. and Kuhlemeyer, R.L., 1979 “Finite Dynamic Model for Infinite Media,”
Journal of the Engineering Mechanics Division of ASCE, pp. 859-877.
Malkin, A. Y. and Petrie, C.J.S., 1997 “Some Conditions for Rupture of Polymer Liquids
in Extension,” Journal of Rheology, Vol. 41 (1), pp. 1-25.
Malvern, L.E., 1969 Introduction to the Mechanics of a Continuous Medium,

Prentice-Hall, Inc., New Jersey, pp. 226-228.
Marsden, J.E. and Hughes, T.J.R., 1983 Mathematical Foundations of Elasticity,

Prentice-Hall, Englewood Cliffs, New Jersey.
Matthies, H. and Strang, G., 1979 “The Solution of Nonlinear Finite Element Equations,”
International Journal for Numerical Methods in Engineering, Vol. 14, 1613-1626.
McClelland, 19?? ???

McClure, C., 1993 Preliminary Report on Explosive Field Tests in Support of the
< Go Back
Hull Deformation/Rupture Study, NSWC Report NSWCDD/TN-93/94.
Mendelson, A., 1968 Plasticity: Theory and Application, The Macmillan Company,
New York, pp. 138-156.
Metzinger, K., Attaway, S. and Mello, F., 1991 “Bobbin Stresses Generated by Wire
Winding,” First International Conference of Web Handling, Oklahoma State
University, Stillwater, OK.
Michalowski, R. and Mroz, Z., 1978 “Associated and Non-Associated Sliding Rules in
Contact Friction Problems,” Archives of Mechanics, Vol. 30, 259-276.
Mills-Curran, W.C., 1988 EXODUS: A Finite Element File Format for Pre- and Post-
Processing, SAND87-2997, Sandia National Laboratories, Albuquerque, NM.
Mindlin, R.D., 1951 “Influence of Rotatory Inertia and Shear on Flexural Motions of
Isotropic, Elastic Plates,” Journal of Applied Mechanics, Vol. 18, pp. 31-38.
Monaghan, J.J., 1982 “Why Particle Methods Work,” SIAM Journal Scientific
Statistical Computing, Vol. 3, pp. 422-433.
Monaghan, J.J., 1985 “Particle Methods for Hydrodynamics,” Computational Physics

Reports, Vol. 3, pp. 71-124.

Monaghan, J.J., 1988 “An Introduction to SPH,” Computational Physics
< Go Back
Communications, Vol. 48, pp. 89-96.
Monaghan, J.J. and Gingold, R.A., 1983 “Shock Simulation by the Particle Method
SPH,” Journal of Computational Physics, Vol. 52, pp. 374-389.
Morland and Lee, 1960 ???
Mungiza, A., Owen, D.R.J. and Bicanic, N., 1995 “A Combined Finite-Discrete Element
Method in Transient Dynamics of Fracturing Solids,” Engineering Computations,
Vol. 12, pp. 145-174.
Muki and Sternberg, 1961 ???
Neilsen, M.K., Morgan, H.S. and Krieg, R.D., 1986 A Phenomenological Constitutive
Model for Low Density Polyurethane Foams, SAND86-2927, Sandia National
Nemat-Nasser, S., Chung, D. and Taylor, L.M., 1989 “Phenomenological modeling of

rate-dependent plasticity for high strain rate problems,” Mechanics of Materials, Vol.
7, No. 4, pp. 319-344.
Newmark, N.M., 1959 “A Method of Computation for Structural Dynamics,” Journal of

the Engineering Mechanics Division, ASCE, 67-94.

Ogden, R.W., 1987 “Recent advances in the phenomenological theory of rubber
< Go Back
elasticity,” Rubber Chemistry and Technology, Vol. 59, pp. 361-383.
Parisch, H., 1989 “A Consistent Tangent Stiffness Matrix for Three Dimensional Non-
Linear Contact Analysis,” International Journal for Numerical Methods in
Engineering, Vol. 28, 1803-1812.
Patel, N.R. and Finnie, I., 1969 ???, Report UCRL-13420, Lawrence Livermore
Pearson, G.H. and Connelly, R.W., 1982 “The Use of Extensional Rheometry to Establish
Operating Parameters for Stretching Processes,” Journal of Applied Polymer Science,
Vol. 27, pp. 969-981.
Perzyna, P., 1966 “Fundamental Problems in Viscoplasticity,” in Recent Advances in

Applied Mechanics, Academic Press, New York, pp. 243-377.
Pfeiffle, T.W. and Senseny, P.E., 1985 Permeability and Consolidation of Crushed Salt
from the WIPP Site, Topical Report RSI-0278, RE/SPEC Inc., Rapid City, SD.
Pilkey, W.D. and Wunderlich, W., 1994 Mechanics of Structures: Variational and
Computational Methods, CRC Press.
Plimpton, S.J., 1990 “Molecular Dynamics Simulations of Short-Range Force Systems on

1024-Node Hypercubes,” Proceedings of the Fifth Distributed Memory Computing
Conference, Charleston, SC.
Powell, M.J.D., 1977 “Restart procedures for the conjugate gradient method,”
< Go Back
Mathematical Programming, Vol. 12, pp. 242-254.
Ralston, A. and Wilf, H.S., 1960 Mathematical Methods of Digital Computers, John
Wiley and Sons Inc., New York, pp. 62-72.
Ratigan and Wagner, 19?? ???
Rebelo, N., Nagtegaal, J.C. and Hibbitt, H.D., 1990 “Finite Element Analysis of Sheet
Forming Processes,” International Journal for Numerical Methods in Engineering,
Vol. 30, 1739-1758.
Reedy, J.E.D., 1990 “Memo: Code corrections to pronto’s soils and crushable foams
model,” Applied Mechanics Division, Sandia National Laboratories, Albuquerque,
NM.
Richtmyer, R.D. and Morton, K.W., 1967 Difference Methods for Initial Value
Problems, Interscience, New York.
Rivlin, R.S., 1948 “???”, Philosophical Transactions of the Royal Society of London,
A, pp. 459-490.
Schoof, L.A. and Yarberry, V.R., 1994 EXODUS II: A Finite Element Data Model,

Schreyer, H.L., Kulak, R.F. and Kramer, J.M., 1979 “Accurate numerical solutions for
< Go Back
elastic-plastic models,” ASME Journal of Pressure Vessel Technology, Vol. 101, pp.
226-234.
Schwarzl and Staverman, 1952 ???
Silling, S. A., 1991 CTH Reference Manual: Viscoplastic Models, SAND91-0292,

Sandia National Laboratories, Albuquerque, NM.
Simo, J.C., 1992 “Algorithms for Static and Dynamic Multiplicative Plasticity that
Preserve the Classical Return Mapping Schemes of the Infinitesimal Theory,”
Computer Methods in Applied Mechanics and Engineering, Vol. 99, 61-112.
Simo, J.C., Marsden, J.E. and Krishnaprasad, P.S., 1988 “The Hamiltonian Structure of
Nonlinear Elasticity: The Material and Convective Representations of Solids, Rods and
Plates,” Archive for Rational Mechanics, Vol. 104, 125-183.
Simo, J.C. and Miehe, C., 1992 “Associative Coupled Thermoplasticity at Finite Strains:
Formulation, Numerical Analysis and Implementation,” Computer Methods in
Applied Mechanics and Engineering, Vol. 98, 41-104.
Simo, J.C. and Taylor, R.L., 1985 “Consistent Tangent Operators for Rate-Independent
Elastoplasticity,” Computer Methods in Applied Mechanics and Engineering, Vol.
48, 101-118.

Simo, J.C., Taylor, R.L. and Wriggers, P., 1991 “A Note on Finite Element
< Go Back
Implementation of Pressure Boundary Loading,” Communications in Applied
Numerical Methods, Vol. 50, 163-180.
Sjaardema, G.D., 1989 NUMBERS: A Collection of Utilities for Pre- and Post-
Processing Two-and Three-Dimensional EXODUS Finite Element Models,
Sjaardema, G.D., 1992 GJOIN: A Program for Merging Two or More GENESIS
Databases Version 1.4, SAND92-2290, Sandia National Laboratories, Albuquerque,
NM.
Sjaardema, G.D., 1993 GENSHELL: A GENESIS Database Shell Transformation

Program, Sandia National Laboratories, Albuquerque, NM.
Sjaardema, G.D., 1993 Overview of the Sandia National Laboratories Engineering

Analysis Code Access System, SAND92-2292, Sandia National Laboratories,
Albuquerque, NM.
Sjaardema, G.D. and Krieg, R.D., 1987 A Constitutive Model for the Consolidation of
WIPP Crushed Salt and Its Use in Analyses of Backfilled Shaft and Drift
Configurations, SAND87-1977*UC-70, Sandia National Laboratories, Albuquerque,
NM.
Stone, C.M., 1989 “Personal communication,” Applied Mechanics Division, Sandia

Stone, C.M., 1995 SANTOS: A Two-Dimensional Finite Element Program faor the
< Go Back
Quasistatic Large Deformation, Inelastic Response of Solids, SAND90-0543,
Stone, C.M., Krieg, R.D. and Beisinger, Z.E., 1988 SANCHO - A Finite Element
Computer Program for the Quasistatic, Large Deformation, Inelastic Response of
Two-Dimensional Solids, SAND84-2618, Sandia National Laboratories,
Albuquerque, NM.
Stone, C. M. and Wellman, G. W., 1993 “Implementation of Ductile Failure in

PRONTO2D and PRONTO3D”, Memo to Distribution, Sandia National Laboratories,
Albuquerque, NM.
Stone, C.M., Wellman, G.W. and Krieg, R.D., 1990 A Vectorized Elastic/Plastic Power
Law Hardening Material Model Including Luders Strain, SAND90-0153, Sandia
Strang, G. and Fix, G.J., 1973 An Analysis of the Finite Element Method, Prentice-
Hall, Inc., Englewood Cliffs, New Jersey.
Swegle, J.W., 1978 TOODY IV - A Computer Program for Two-Dimensional Wave

Propagation, SAND78-0552, Sandia National Laboratories, Albuquerque, NM.
Swegle, J.W., 1992 “Search Algorithm,” Memo to Distribution, Sandia National


Swegle, J.W., Attaway, S.W., Heinstein, M.W., Mello, F.J. and Hicks, D.L., 1993 An
< Go Back
Analysis of Smoothed Particle Hydrodynamics, SAND93-2513, Sandia National
Laboratoaries, Albuquerque, NM.
Swenson, D.V. and Taylor, L.M., 1983 “A Finite Element Model for the Analysis of
Tailored Pulse Stimulation of Boreholes,” International Journal for Numerical and
Analytical Methods in Geomechanics, Vol. 7, pp. 469-484.
Takaki, T. and Bogue, D.C., 1975 “The Extensional and Failure Properties of Polymer
Melts,” Vol. 19, pp. 419-433.
Taylor, L.M. and Becker, E.B., 1983 “Some Computational Aspects of Large
Deformation, Rate-Dependent Plasticity Problems,” Computer Methods in Applied
Mechanics and Engineering, Vol. 41, No. 3, pp. 251-277.
Taylor, L.M., Chen, E.P. and Kuszmaul, J.S., 1986 “Microcrack-Induced Damage
Accumulation in Brittle Rock Under Dynamic Loading,” Computer Methods in
Applied Mechanics and Engineering, Vol. 55, No. 3, pp. 301-320.
Taylor, L.M., Flanagan, D.P. and Mills-Curran, W.C., 1986 The GENESIS Finite
Element Mesh File Format, SAND86-0910, Sandia National Laboratories,
Albuquerque, NM.
Taylor, L.M. and Flanagan, D.P., 1987 PRONTO 2D: A Two Dimensional Transient
Solid Dynamics Program, SAND86-0594, Sandia National Laboratories,
Albuquerque, NM.
Taylor, L.M. and Flanagan, D.P., 1989 PRONTO 3D: A Three Dimensional Transient
< Go Back
Solid Dynamics Program, SAND 87-1912, Sandia National Laboratories,
Albuquerque, NM.
Thompson, S.L., 1977 CSQII - An Eulerian Finite Difference Program for Two-
Dimensional Material Response - Part 1, Material Selections, SAND77-1339,
Thorne, B.J., 1990 A Damage Model for Rock Fragmentation and Comparison of
Calculations with Blasting Experiments in Granite, SAND90-1389, Sandia
Thorne, B.J. and Preece, D.S., 1989 “Personal communications,” Geoenergy Technology
Department, Sandia National Laboratories, Albuquerque, NM.
Thrun, R., Goertner, J.F. and Harris, G.S., 1993 Underwater Explosion Bubble
Collapse Against a Flat Plate, Seneca Lake Test Series Data Report, NSWC Report,
NSWCDD/TR-92/482.
Timoshenko, S. and MacCollough, 1958 ???
Timoshenko, S. and Goodier, J.N., 1970 Theory of Elasticity, 3rd ed., McGraw-Hill,
New York.
Treloar, 1994 ???

Truesdell, C., 1966 The Elements of Continnum Mechanics, Springer Verlag, New
< Go Back
York.
Truesdell, C., 1977 A First Course in Rational Continuum Mechanics, Vol. 1,

General Concepts, Academic Press, Inc., New York, p. 162.
Truesdell, C. and Noll, W., 1965 “Non-Linear Field Theories”, in the Handbook of
Physics by Flugge, Springer-Verlag, Berlin.
Valanis, K. and Lnadel, R.F., 1967 “The strain energy function of a hyperelastic material
in terms of the extension ratios,” Journal of Applied Physics, Vol. 38.
Von Neumann, J. and Richtyer, R.D., 1950 “A Method for the Numerical Calculation of
Hydrodynamic Shocks,” Journal of Applied Physics, Vol. 21, p. 232.
Warren, T.L. and Forrestal, M.J., 1997 “Effects of strain hardening and strain-rate
sensitivity on the penetration of aluminum targets with spherical-nosed rods”,
International Journal of Solids Structures (submitted).
Warren, T.L. and Tabbara, M.R., 1997 Spherical Cavity-Expansion Forcing Function
in PRONTO3D for Application to Penetration Problems, SAND97-1174, Sandia
Weaver, Jr., W. and Johnson, P. R., 1984 Finite Elements for Structural Analysis,
Prentice-Hall, Englewood Cliffs, New Jersey.

Wellman, G. W., 1993 “Investigation of Mesh Dependencies in Ductile Failure for
< Go Back
Transient Dynamics”, Memo to Distribution, Sandia National Laboratories,
Albuquerque, NM.
Wen, Y., Hicks, D. L. and Swegle, J. W., 1994 Stabilizing S.P.H. with Conservative
Smoothing, SAND94-1932, Sandia National Laboratories, Albuquerque, NM.
Whirley, R. G., Engelmann, B. E. and Hallquist, J. O., 1991 DYNA3D Users Manual,
Lawrence Livermore Laboratory, Livermore, CA.
Whirley, R.G., Hallquist, J.O. and Goudreau, G.L., 1988 An Assessment of Numerical
Algorithms for Plane Stress and Shell Elastoplasticity on Supercomputers, UCRL-
99690, Lawrence Livermore National Laboratory, Livermore, CA.
Wilkins, M.L. and Guinam, M.W., 1973 “Impact of CYlinders on a Rigid Boundary”,
Journal of Applied Physics, Vol. 44.
Williams, Landell and Ferry, 1955 ???
Wriggers, P. and Simo, J.C., 1985 “A Note on Tangent Stiffness for Fully Nonlinear
Contact Problems,” Communications in Applied Numerical Methods, Vol. 1, 199-
203.
Wriggers, P., Vu Van, T. and Stein, E., 1990 “Finite Element Formulation of Large
Deformation Impact-Contact Problems with Friction,” Computers and Structures,
Vol. 37, 319-331.
Zeuch, Holcomb and Lauson, 19?? ???
< Go Back
Zhang, P. and Geers, T. L., 1993 “Excitation of a fluid-filled, submerged spherical shell
by a transient acoustic wave,” Journal Acoustical Society of America, Vol. 93, pp.
696-705.
Zhong, Z.H. and Nilsson, L., 1990 “A Contact Searching Algorithm for General 3D
Contact-Impact Problems,” Computers and Structures, Vol. 34, No. 2, pp. 327-335.

SEACAS Theory PRONTO 3D JAS 3D ACCESS Author’s
Library Manuals Manuals Manuals Manuals Manual
< Go Back
Theory Manuals
Blue text
i indicates Formulation of Nonlinear Problems
a link to more
information. Nonlinear Continuum Mechanics
Constitutive Modeling
Finite Element Formulations
SPH Formulation
Contact Surfaces
Theory Manuals (2/19/99) Main Menu 1

< Go Back SEACAS Library
Blue text
Theory Manuals
i indicates
PRONTO3D Manuals
a link to more
information.
JAS3D Manuals
ACCESS Manuals
Author’s Manual
SEACAS Library (2/19/99) Main Menu 1

Fin Elem

Uploaded by

Copyright:

Available Formats

Fin Elem

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Fin Elem

Uploaded by

Copyright:

Available Formats

SEACAS Formulation of Nonlinear Constitutive Finite Element SPH Contact

Library Nonlinear Continuum Models Formulations Formulation Surfaces

Finite Element Formulations

Theory Manuals (2/19/99) Finite Element Formulations - Main Menu 1

Theory Manuals (2/19/99) Finite Element Formulations - Introduction 1

Theory Manuals (2/19/99) Finite Element Formulations - Introduction - Introduction 3

Weak Form Revisited

Theory Manuals (2/19/99) Finite Element Formulations - Weak Form Revisited 1

The problem is also subject to initial conditions of the form

< Go Back ϕ∗ = 0 on Γ u (3.6)

SEACAS for any x ∈ ϕ t ( X ) . This causes the condition

Given the boundary conditions t on ϕ t ( Γ σ ) , ϕ t on ϕ t ( Γ u ) , the initial conditions ϕ 0

for all admissible w , where S t is as defined as

In addition, the solution ϕ is subject to the following conditions at t = 0 :

< Go Back and

both of which must hold for all ϕ∗ ∈ W .

Theory Manuals (2/19/99) Finite Element Formulations - Discretization 1

As referred to in Figure 3.1, the reference domain Ω is subdivided into a number of

Theory Manuals (2/19/99) Finite Element Formulations - Discretization - Introduction 2

Theory Manuals (2/19/99) Finite Element Formulations - Discretization - Introduction 3

In other words, we require members of the discrete solution space to (approximately)

Finite Element nnp

related to the material variations ϕ∗ ∈ W via

< Go Back The second term of (3.29) can be simplified via

Finally, the last two terms of (3.29) can be treated as

The results of Eqs. (3.30)-(3.35) can now be summarized as follows:

< Go Back a=1 a=2

Step 2: For each element e , e = 1, …, n el :

and Eq. (3.42).

int int int e

Theory Manuals (2/19/99) Finite Element Formulations - Quasistatics 1

Theory Manuals (2/19/99) Finite Element Formulations - Quasistatics - Introduction 2

Equation (3.53) as written implies no dependence of either t i or f i upon ϕ t ( x ) (and

Finite Element N–1

Given the solution d n corresponding to time level t n , find d n + 1 corresponding to

Solution of (3.56), therefore, amounts to finding the root of the equation

Theory Manuals (2/19/99) Finite Element Formulations - Dynamics 1

to be solved for t ∈ [ 0, T ] , subject to the initial conditions

Theory Manuals (2/19/99) Finite Element Formulations - Dynamics - Introduction 2

Theory Manuals (2/19/99) Finite Element Formulations - Dynamics - Time-Stepping Procedures 4

Theory Manuals (2/19/99) Finite Element Formulations - Dynamics - Time-Stepping Procedures 5

where the first equation has been written as solved for a n + 1 .

Subtracting Eq. (3.66) from (3.65) gives

However, evaluation of (3.64) during the time step [ t n – 1, t n ] reveals that

Nonlinear Equation Solving

Theory Manuals (2/19/99) Finite Element Formulations - Nonlinear Equation Solving 13

followed by the update

so that Eq. (3.79) can be rewritten as

This causes (3.81) to take the form

< Go Back least controlling growth in) the residual.

Find s such that

a typical algorithm to find s could be outlined as:

Iterate for s ∈ ( 0 , 1 ] such that G ( s ) < ( TOL ) G ( 0 ) (3.97)

In generalizing this concept to multiple directions, we seek to approximate so-called

< Go Back or in terms of the inverse:

The next search direction ∆d i is then computed via

where α i is given by (cf. (3.93)):

Subsequent search directions are defined via

where the coefficients β ik must be found to ensure K-orthogonality of p i with previous