NONLINEAR GAP AND MINDLIN SHELL ELEMENTS
FOR THE ANALYSIS OF CONCRETE STRUCTURES
By
KOOKJOON AHN
A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
U::i: L "At; 7 OF FLORIN LIfRPMR17
1990
ACKNOWLEDGEMENTS
I would like to express my deep gratitude to professor
Marc I. Hoit for his invaluable guidance and support. I also
thank professors Clifford 0. Hays, Jr., Fernando E. Fagundo,
John M. Lybas, and Paul W. Chun for being on my committee. I
also express my gratitude to professor Duane S. Ellifritt
for his help as my academic advisor at the start of my Ph.D.
program.
Thanks are also due to all the other professors not
mentioned above and my fellow graduate students, Alfredo,
Jose, Mohan, Prasan, Prasit, Shiv, Tom, Vinax, and Yuhyi.
Finally, I am thankful to every member of my family,
especially my wife and son, for their patience and support
in one way or another.
The work presented in this dissertation was partially
sponsored by the Florida Department of Transportation.
TABLE OF CONTENTS
page
ACKNOWLEDGEMENTS ...................................... ii
ABSTRACT ...................................... ........ v
CHAPTERS
1 INTRODUCTION ....................... ............. 1
1.1 General Remarks ............................ 1
1.2 Link Element ............................... 2
1.3 Shell Element .............................. 5
1.4 Literature Review .......................... 5
2 GENERAL THEORIES OF NONLINEAR ANALYSIS .......... 13
2.1 Introduction ............................... 13
2.2 Motion of a Continuum ...................... 14
2.3 Principle of Virtual Work .................. 16
2.4 Updated Lagrangian Formulation ............. 18
2.5 Total Lagrangian Formulation ............... 22
2.6 Linearization of Equilibrium Equation ...... 26
2.7 StrainDisplacement relationship
Using von Karman Assumptions ............... 28
3 THREEDIMENSIONAL LINK ELEMENT .................. 34
3.1 Element Description ...................... 34
3.2 Formation of Element Stiffness ............. 43
3.3 Solution Strategy .......................... 51
3.4 Element Verification ....................... 52
4 LINEAR SHELL ELEMENT ............................ 59
4.1 Introduction ............................... 59
4.2 Formulation of Shape Functions ............. 59
4.3 The Inverse of Jacobian Matrix ............. 64
4.4 Membrane Element ........................... 66
4.5 Plate Bending Element .................... 73
5 NONLINEAR SHELL ELEMENT ........................ 92
5.1 Introduction ............. ................. 92
5.2 Element Formulation ........................ 93
iii
5.3 Finite Element Discretization .............. 100
5.4 Derivation of Element Stiffness ............ 113
5.5 Calculation of Element Stiffness Matrix .... 115
5.6 Element Stress Recovery .................... 119
5.7 Internal Resisting Force Recovery ......... 122
6 NONLINEAR SHELL ELEMENT PERFORMANCE ............ 126
6.1 Introduction ............................... 126
6.2 Large Rotation of a Cantilever ............ 126
6.3 Square Plates .............................. 133
7 CONCLUSIONS AND RECOMMENDATIONS.................. 143
APPENDICES
A IMPLEMENTATION OF LINK ELEMENT .................. 146
B IMPLEMENTATION OF LINEAR SHELL ELEMENT ......... 170
C IMPLEMENTATION OF NONLINEAR SHELL ELEMENT ....... 219
REFERENCES ............................................ 230
SUPPLEMENTAL BIBLIOGRAPHY ............................ 238
BIOGRAPHICAL SKETCH ................................... 240
Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy
NONLINEAR GAP AND MINDLIN SHELL ELEMENTS
FOR THE ANALYSIS OF CONCRETE STRUCTURES
By
KOOKJOON AHN
August, 1990
Chairman: Marc I. Hoit
Major Department: Civil Engineering
Segmental posttensioned concrete box girders with
shear keys have been used for medium to long span bridge
structures due to ease of fabrication and shorter duration
construction.
Current design methods are predominantly based on
linear elastic analysis with empirical constitutive laws
which do not properly quantify the nonlinear effects, and
are likely to provide a distorted view of the factor of
safety.
Two finite elements have been developed that render a
rational analysis of a structural system. The link element
is a twodimensional friction gap element. It allows opening
and closing between the faces of the element, controlled by
the normal forces. The Mindlin flat shell element is a
combination of membrane element and Mindlin plate element.
This element considers the shear responses along the element
thickness direction. The shell element is used to model the
segment itself.
The link element is used to model dry joints and has
shown realistic element behavior. It opens under tension and
closes under compression. The link element has shown some
convergence problems and exhibited a cyclic behavior.
The linear Mindlin shell element to model the concrete
section of the hollow girder showed an excellent response
within its small displacement assumption.
The nonlinear Mindlin flat shell element has been
developed from the linear element to predict large
displacement and initial stress (geometric) nonlinearities.
The total Lagrangian formulation was used for the
description of motion. The incrementaliterative solution
strategy was used. It showed satisfactory results within the
limitation of moderate rotation.
Three areas of further studies are recommended. The
first is the special treatment of finite rotation which is
not a tensorial quantity. The second is the displacement
dependent loadings commonly used for shell elements. The
third is the material nonlinearity of concrete which is
essential to provide realistic structural response for safe
and cost effective designs.
CHAPTER 1
INTRODUCTION
1.1 General Remarks
In the past few decades segmental posttensioned
concrete box girders have been used for medium to long span
bridge structures. Highway aesthetics through long spans,
economy due to ease of fabrication, shorter construction
duration are some of the many advantages of precast segment
bridge construction.
The segments are hollow box sections, match cast with
shear keys in a casting yard, then assembled in place,
leaving the joints entirely dry. The shear keys are meant to
transfer service level shears and to help in alignment
during erection.
Current design methods are heavily based on linear
elastic analysis with empirically derived constitutive laws
assuming homogeneous, isotropic materials. The behavior
under load of the bridge system is very complex. Analyses
which do not properly quantify the nonlinear effects
including the opening of joints in flexure, are likely to
provide a distorted view of the factor of safety existing in
a structural system between service loads and failure. The
potential sliding and separation at the joints due to shear,
and by deformations generated by temperature gradients over
the depth and width of the box further complicate the
problem [1].
Two finite elements have been developed that render a
rational analysis of the system. The link element is a two
dimensional friction gap element. It allows sliding between
the faces of the element, controlled by a friction
coefficient and the normal forces. It also accounts for zero
stiffness in tension and a very high stiffness under
compression. This link element was borrowed from rock
mechanics and newly applied to this problem to model the dry
joint between the segments. The Mindlin flat shell element
is a combination of membrane element and Mindlin plate
element. This element considers the shear responses along
the element thickness direction. The shell element was used
to model the segment itself. This element can handle large
displacement and geometric nonlinearities.
1.2 Link Element
A link element is a nonlinear friction gap element used
to model discontinuous behavior in solid mechanics. Some
examples are interfaces between dissimilar materials and
joints, fractures in the material, and planes of weakness.
These have been modeled using constraint equations, discrete
springs and a quasicontinuum of small thickness [2]. The
following characteristics of prototype joints were
considered.
1. Joints can be represented as flat planes.
2. They offer high resistance to compression in the
normal direction but may deform somewhat modeling
compressible filling material or crushable
irregularities.
3. They have essentially no resistance to a net tension
force in the normal direction.
4. The shear strength of joints is frictional. Small
shear displacements probably occur as shear stress
builds up below the yield shear stress.
A model for the mechanics of jointed rocks was
developed by Goodman [3]. The finite element approximation
was done as a decomposition of the total potential energy of
a body into the sum of potential energies of all component
bodies. Therefore, element stiffness is derived in terms of
energy.
The Goodman element was tested for several modeled
cases.
1. Sliding of a joint with a tooth.
2. Intersection of joints.
3. Tunnel in a system of staggered blocks.
A problem with the Goodman's two dimensional model is
that adjacent elements can penetrate into each other.
Zienkiewicz et al. [4] advocate the use of continuous
isoparametric elements with a simple nonlinear material
property for shear and normal stresses, assuming uniform
strain in the thickness direction. Numerical difficulties
may arise from ill conditioning of the stiffness matrix due
to very large offdiagonal terms or very small diagonal
terms which are generated by these elements in certain
cases. A discrete finite element for joints was introduced
which avoids such theoretical difficulties and yet is able
to represent a wide range of joint properties, including
positive and negative dilatency (expansion and compaction
accompanying shear) [3].
The element uses relative displacements as the
independent degrees of freedom. The displacement degrees of
freedom of one side of the slip surface are transformed into
the relative displacements between the two sides of the slip
surface. This element has been incorporated into a general
finite element computer program [5]. The use of relative
displacement as an independent degree of freedom to avoid
numerical sensitivity is discussed in detail [6]. An
isoparametric formulation is given by Beer [2]. A fournode,
twodimensional link element and a eightnode plate bending
element were used to model the dry jointed concrete box
girder bridge with shear keys [7].
1.3 Shell Element
The shell element is formulated through the combination
of two different elements, the membrane element and the
Mindlin plate bending element.
The Mindlin plate element is different from the
Kirchhoff plate element in that the former allows transverse
shear deformation while the latter does not.
The nonlinearities included in the formulation of the
flat shell element is for large displacement and geometric
nonlinearity due to initial stress effects. The large
displacement effects are caused by finite transverse
displacements. These effects are taken into account by
coupling transverse displacement and membrane displacements.
The initial stress effects are caused by the actual stresses
at the start of each iteration. These stresses change the
element stiffness for the subsequent iteration. These
effects are evaluated directly from the stresses at the
start of each iteration and included in the element
stiffness.
1.4 Literature Review
The purpose of nonlinear analysis is to develop the
capability for determining the nonlinear loaddeflection
behavior of the structures up to failure so that a proper
evaluation of structural safety can be assured. There are
two general approaches for nonlinear analysis. The first
approach is a linearized incremental formulation by reducing
the analysis to a sequence of linear solutions. The second
approach is mathematical iterative techniques applied to the
governing nonlinear equations [8].
The advantage of the incremental approach results from
the simplicity and generality of the incremental equations
written in matrix form. Such equations are readily
programmed in general form for computer solutions [9].
A generalized incremental equilibrium equation for
nonlinear analysis can be found in [10, 11, 12]. The
formulation is valid for both geometrical and material
nonlinearities, large displacements and rotations,
conservative and displacement dependent (nonconservative)
loads.
There are two frames for the description of motion. The
difference lies in the coordinate systems in which the
motion is described. These are the total Lagrangian
formulation which refers to the initial configuration [10,
11] and the updated Lagrangian formulation which refers to
the deformed configuration [12]. There have evolved two
types of notations in the description of motion. A
correlation is given these two notations, the Bnotations
and the Nnotations, currently used in the Lagrangian
formulation of geometrically nonlinear analysis [13]. A
short history of early theoretical development of nonlinear
analysis can be found in [9, 14].
One form of updated Lagrangian formulation is the
corotational stretch theory [15].
Shell elements are often derived from governing
equations based on a classical shell theory. Starting from
the field equations of the threedimensional theory, various
assumptions lead to a shell theory. This reduction from
three to two dimensions is combined with an analytical
integration over the thickness and is in many cases
performed on arbitrary geometry. Static and kinematic
resultants are used. These are referred to as classical
shell elements. Alternatively, one can obtain shell elements
by modifying a continuum element to comply with shell
assumptions without resorting to a shell theory. These are
known as degenerated shell elements. This approach was
originally introduced by Ahmad, Irons, and Zienciewicz [16,
17]. Other applications can be found in [8, 1825].
In large rotation analysis, the major problems arise
from the verification of the kinematic assumptions. The
displacement representation contains the unknown rotations
of the normal in the arguments of trigonometric functions.
Thus additional nonlinearity occurs. Further difficulties
enter through the incremental procedure. Rotations are not
tensorial variables, therefore, they cannot be summed up in
an arbitrary manner [17]. One of the special treatment of
finite rotation is that the rotation of the coordinate
system is assumed to be accomplished by two successive
rotations, an outofplane rotation followed by an inplane
rotation using updated Lagrangian formulation [26, 27].
Usually the loadings are assumed to be conservative,
i.e., they are assumed not to change as the structure
deforms. One of the well known exceptions is pressure
loading which can be classified as conservative loading or a
nonconservative loading [28]. Another is the concentrated
loading that follows the deformed structure. For example, a
tip loading on a cantilever beam will change its direction
as the deformation gets larger. As loading is a vector
quantity, the change in direction means that the loading is
not conservative. Sometimes this is called a follower
loading.
The governing equation for large strain analysis can be
used for small increments of strain and large increments of
rotations [29]. This can be regarded as a generalization of
nonlinearity of small strain with large displacement. If
large strain nonlinearity is employed, an important question
is which constitutive equation should be used [9].
The degree of continuity of finite element refers to
the order of partial differential of displacements with
respect to its coordinate system. Order zero means
displacement itself must be continuous over the connected
elements. Order one means that the first order differential
of displacement must be continuous. Thus the higher order
the continuity requirement, the higher the order of assumed
displacement (shape, interpolation) function.
MindlinReissner elements require only Co continuity,
so that much lower order shape functions can be used,
whereas in KirchhoffLove type elements, high order shape
functions must be used to satisfy the C1 continuity.
Furthermore, since MindlinReissner elements account for
transverse shear, these elements can be used for a much
larger range of shell thickness. The relaxed continuity
requirements which permit the use of isoparametric mapping
techniques gives good computational efficiency if formulated
in the form of resultant stresses [30]. Unlike compressible
continuum elements, which are quite insensitive to the order
of the quadrature rule, curved Co shell elements require
very precisely designed integration scheme. Too many
integration points result in locking phenomena, while using
an insufficient number of quadrature points results in rank
deficiency or spurious modes [30]. While Gauss point stress
results are very accurate for shallow and deep, regular and
distorted meshes, the nodal stresses of the quadratic
isoparametric Mindlin shell element are in great error
because of the reduced integration scheme which is necessary
to avoid locking [31].
The degenerate solid shell element based on the
conventional assumed displacement method suffers from the
locking effect as shell thickness becomes small due to the
condition of zero inplane strain and zero transverse shear
strain. Element free of locking for linear shell analysis
using the formulation based on the HellingerReissner
principle with independent strain as variables in addition
to displacement is presented in [32].
Shear locking is the locking phenomenon associated with
the development of spurious transverse shear strain.
Membrane locking is the locking phenomenon associated with
the development of nonzero membrane strain under a state of
constant curvature. Machine locking is the locking
phenomenon associated with the different order of dependence
of the flexural and real transverse shear strain energies on
the element thickness ratio, and it is therefore strictly
related to the machine finite word length [33]. Some of the
solutions are as follows:
1. Assumed strain stabilization procedure using the Hu
Washizu or HellingerReissner variational principles
[33].
2. The assumed strain or mixed interpolation approach [34,
35].
3. Suppressing shear with assumed stress/strain field in a
hybrid/mixed formulation [30]. Suppression of zero
energy deformation mode using assumed stress finite
element [36].
4. Coupled use of reduced integration and nonconforming
modes in quadratic Mindlin plate element [37].
5. Higher order shallow shell element, with 17 to 25 nodes
[38, 39].
6. Global spurious mode filtering [40].
7. Artificial stiffening of element to eliminating zero
energy mode, special stabilizing element [41].
In the faceted elements, due to the faceted
approximation of the shell surface, coupling between the
membrane and the flexural actions is excluded within each
individual element, the coupling is, however, achieved in
the global model through the local to global coordinate
transformation for the elements [39].
In geometrically nonlinear analysis with flat plate
elements, it is common to use the von Karman assumptions
when evaluating the straindisplacement relations. The
assumption invoked is that the derivatives of the inplane
displacements can be considered to be small and hence their
quadratic variations neglected. However, this simplification
of the nonlinear straindisplacement relationship of the
plate, when used in conjunction with the total Lagrangian
approach, implies that the resulting formulation is valid
only when the rotation of the element from its initial
configuration is moderate. Thus for the total Lagrangian
approach to handle large rotations, simplifications of the
kinematic relationship using the von Karman assumptions is
not permitted [39].
Some of the special solution strategies to pass the
limit point are given in references [25, 4248]. A limit
point is characterized by the magnitude of tangential
stiffness. It is zero or infinite at a limit point. Thus
conventional solution strategies fail at the limit point.
Arc length method was introduced in reference [42], and
applied in the case of cracking of concrete [43]. This was
improved with line search and accelerations in references
[44, 45]. Line search means the calculation of an optimum
scalar step length parameter which scales the standard
iterative vector. This can be applied to load and
displacement control and arc length methods [44].
The traditional solution strategies are iterative
solutions, for example, NewtonRaphson, constant stiffness,
initial stiffness, constant displacement iteration, load
increment [46] along with Cholesky algorithm with shifts for
the eigensolution of symmetric matrices [47] for element
testing for spurious displacement mode.
The vector iteration method without forming tangent
stiffness for the postbuckling analysis of spatial
structures is also noted [48].
The linearized incremental formulation in total
Lagrangian description has been used for this study of large
displacement nonlinearity including initial stress effects.
The special treatment of finite rotation is not included in
the current study. Material nonlinearity is also excluded.
CHAPTER 2
GENERAL THEORIES OF NONLINEAR ANALYSIS
2.1 Introduction
The incremental formulations of motion in this chapter
closely follow the paper by Bathe, Ramm, and Wilson [11].
Other references are also available [9, 10, 12, 14, 15, 49,
50, 51].
Using the principle of virtual work, the incremental
finite element formulations for nonlinear analysis can be
derived. Time steps are used as load steps for static
nonlinear analysis. The general formulations include large
displacements, large strains and material nonlinearities.
Basically, two different approaches have been pursued
in incremental nonlinear finite element analysis. In the
first, Updated Lagrangian Formulation, static and kinematic
variables, i.e., forces, stresses, displacements, and
strains, are referred to an updated deformed configuration
in each load step. In the second, Total Lagrangian
Formulation, static and kinematic variables are referred to
the initial undeformed configuration.
It is noted that using either of two formulations
should give the same results because they are based on the
same continuum mechanics principles including all nonlinear
effects. Therefore, the question of which formulation should
be used merely depends on the relative numerical
effectiveness of the methods.
2.2 Motion of a Continuum
Consider the motion of a body in a Cartesian coordinate
system as shown in Fig. 21. The body assumes the
equilibrium positions at the discrete time points 0, dt,
2dt, ..., where dt is an increment in time. Assume that the
solution for the static and kinematic variables for all time
steps from time 0 to time t, inclusive, have been solved,
and that the solution for time t+dt is required next.
The superscript on left hand side of a variable shows
the time at which the variable is measured, while the
subscript on left hand side of a variable indicates the
reference configuration to which the variable is measured.
Thus the coordinates describing the configuration of the
body using index notation are
At time 0 = x
At time t = txi
At time t+dt = t+dtx
P t+dt
P( x )
t
P( Xi)
0
P( Xi)
Fig. 21 Motion of a Body
The total displacements of the body are
At time 0 = ui
At time t = tui
At time t+dt = t+dtu
The configurations are denoted as
At time 0 = C
At time t = tc
At time t+dt = t+dtc
Thus, the updated coordinates at time t and time t+dt are
txi = Oxi + tui
t+dtx = Oxi + t+dtu
The unknown incremental displacements from time t to
time t+dt are denoted as (Note that there is no superscript
at left hand side.)
u = t+dtui tui (2.1)
2.3 Principle of Virtual Work
Since the solution for the configuration at time t+dt
is required, the principle of virtual work is applied to the
equilibrium configuration at time t+dt. This means all the
variables are those at time t+dt and are measured in the
configuration at time t+dt and all the integration are
performed over the area or volume in the configuration at
time t+dt. Then the internal virtual work (IVW) by the
corresponding virtual strain due to virtual displacement in
t+dtC is
S t+dt t+dt t+dt
t+dt ij t+dt eij (+dt dV) (2.2)
where,
t+dt rj = Stresses at time t+dt measured in the
t+dt j configuration at time t+dt.
= Cauchy stresses.
= True stresses.
t+dt
t+dt eij = Cauchy's infinitesimal(linear) strain tensor
referred to the configuration at time t+dt.
= Virtual strain tensor.
6 = Delta operator for variation.
and the external virtual work (EVW) by surface tractions and
body forces is
EVW = t+dt t ] 6 [ t+dt u (t+dt dA)
t+dt k t+dt k (t+dt
St+dt t+dt b 6 t+dt u (t+dt dV)
St+dt t+dt k t+dt uk
(2.3)
where,
t+dt tk = Surface traction at time t+dt measured in
t+dt the configuration at time t+dt.
t+dt Uk = Total displacement at time t+dt measured in
t+dt the configuration at time t+dt.
6 t+dt Uk = Variation in total displacement at time
t+dt t+dt measured in configuration at time t+dt
= Virtual displacement.
t+dt p = Mass density per unit volume.
t+dt
t+dt bk = Body force per unit mass.
t+dt k
and all the integration is performed over the area and the
volume at time t+dt.
2.4 Updated Lagrangian Formulation
In this formulation all the variables in Eqs. (2.2) and
(2.3) are referred to the updated configuration of the body,
i.e, the configuration at time t. The equilibrium position
at time t+dt is sought for the unknown incremental
displacements from time t to t+dt.
The internal virtual work, the volume integral in Eq.
(2.2) measured in the configuration at time t+dt can be
transformed to the volume integral measured in the
configuration at time t in a similar manner that is given in
reference [52]
IVW t+dt t+dt e (t+dt dV)
t+dt ij t+dt i (t dV)
t Sij 6 Et ij (t dV) = EVW (2.4)
I ~t t~d
where,
t+dt Sj = Second PiolaKirchhoff (PKII) stresses
t measured in the configuration at time t.
6 t+dt ej = Variations in GreenLagrange (GL) strain
t tensor measured in the configuration at
time t.
The PKII stress tensor at time t+dt, measured in the
configuration at time t can be decomposed as
t+dt t t
Sij =t Sij + t Sij j + t Sij (2.5)
because the second PKII stress at time t measured in the
configuration at time t is the Cauchy stress.
From Eq. (2.1), the total displacements at time t+dt
measured in the configuration at time t is
t+dt t
t Ui = t Ui + Ui = t ui (2.6)
This is true because the displacement at time t measured in
the configuration at time t is zero. In other words, the
displacement at time t+dt with respect to the configuration
at time t is the incremental displacement itself.
And the GL strain is defined in terms of displacement as
Eij = 1 (Ui,j + Uji + Uk,iUk,j) (2.7)
E and U are used in the places of e and u to avoid confusion
between general strain and incremental strain, and between
general displacement and incremental displacement used in
this formulation. It is noted that these finite strain
components involve only linear and quadratic terms in the
components of the displacement gradient. This is the
complete finite strain tensor and not a second order
approximation to it. Thus this is completely general for any
threedimensional continuum [52].
Then the GL strain tensor at time t+dt measured at time
t can be calculated as
+dt = i [(ui+ i)+ u) + ( j + tuj),i
+{(uk + tuk),i)( tuk + uk),j)]
= tuij + tuj,i+ tuki tuk,j]
= teij + tij
= tij (2.8)
where,
ttij = tei + t ij
= Incremental GL strain in tC.
tej = (ui,j + tuji )
= Linear portion of incremental GL strain in t.
= This is linear in terms of unknown incremental
displacement.
= Linearized incremental GL strain in tc.
ttij = t (uk,i tuk,j)
= Nonlinear portion of incremental GL strain.
The variations in GreenLagrange strain tensor at time
t+dt measured in the configuration at time t can be shown as
using Eq. (2.8).
t+dt t
6 t ij = 6 ( tij + Ei ) = 6 t (2.9)
6teij = 0 because teij is known. There is no variation in
known quantity.
Then using the Eqs. (2.5), (2.8) and (2.9), the
integrand of Eq. (2.4) becomes
t+dt t+dt t
t Sij 6 qt =ij ( t i + tSij ) 6 tCij
= trij + tSij)(6 teij + 6tij)
=tij(6eij + 6t 1ij) + t~ij 6 te + rij 6 tij
t t
=tij 6 ti ij + ti 6eJ 6 (2.10)
tj t t ttij 6 t tii j` t
The constitutive relation between incremental PKII
stresses and GL strains are
tSij = tCijkj t'kl
(2.11)
Finally the equilibrium Eq. (2.4) from the principle of
virtual work using Eqs. (2.10) and (2.11) is
I tCijkl tkl 6 tcij tdv
+ tri 6 t Vij tdV
= EVW J ij 6 eij tdV
t 1t
(2.12)
where, the external virtual work must be transformed from
t+dtc to tC. This is not applicable to conservative loading,
i.e., loading that is not changed during deformation.
EVW= t+dt tk [ t+dt uk (tdA)
+ t+dt t+dt b t+dt td) (2 )
and this is the general nonlinear incremental equilibrium
equation of updated Lagrangian formulation.
2.5 Total Lacrangian Formulation
Total Lagrangian formulation is almost identical with
the updated Lagrangian formulation. All the static and
kinematic variables in Eqs. (2.2) and (2.3) are referred to
the initial undeformed configuration of the body, i.e, the
configuration at time 0. The terms in the linearized strain
are also slightly different from those of updated Lagrangian
formulation.
The volume integral in Eq. (2.2) measured in the
configuration at time t+dt can be transformed to the volume
integral measured in the configuration at time 0 as [52]
t+dt 'i 6 t+dt eij t+dt
t+dt ]i t+dt
o t Sij t+dt j (dV) (2.14)
where,
t+dtS = Second PiolaKirchhoff stress tensor
0 measured in the configuration at time 0.
6 t+dtij = Variations in GreenLagrange (GL) strain
0 tensor measured in the configuration at
time 0.
The PKII stress tensor at time t+dt, measured in the
configuration at time 0 can be decomposed as
t+dtsij= tsij+ *ij (2.15)
o1 ) o) 1 o01j (2.15)
From Eq. (2.1), the total displacements at time t+dt
measured in the configuration at time 0 is
t+dtu = u + Ui (2.16)
o 0i o0 o
Then the GL strain tensor at time t+dt measured at time
0 can be calculated as
t+dt = [(tui + oui)j + (tuj + ou i
=tt+
+(( e*k+ oUk),i)(( Uk+ Uk) j)]
= Eij + e + i
= ij + Eij (2.17)
where,
tt tutu tu tu)
oij = I (t u + + )
0 o o rJ o k,i k,j
= GL strain at time t in oC.
ij = oeij + oij
= Incremental GL strain in oC.
oeij = (oi, + oUj1 + uk,i ouk,kj+ k,j ouk,i)
= Linear portion of incremental GL strain in oC.
=This is linear in terms of unknown incremental
displacement.
= Linearized incremental GL strain in oC.
oij = i (oUk, ouk,j)
= Nonlinear portion of incremental GL strain.
The variations in GreenLagrange strain tensor at time
t+dt measured in the configuration at time t can be shown as
using Eq. (2.17).
Stdtij = 6 ( ij+ ij ) = 6 (2.18)
o i= oec olu olt
6 te, = 0 because ej is known. There is no variation in
o ij oJ
known quantity.
Then using the Eqs. (2.15), (2.17) and (2.18), the
integrand of Eq. (2.14) becomes
t+dt t+dt ts
t+dt*i 6 t+dti ( tSi + S ) 6 e
o o o o0J o'J
= (ij oSij) oeij + 6 oiij)
= Sij(6 eij + 6 ij) + oSij 6 e 6 +
= Si 6 eij + tS i 6 ei + ij 6 o7ij (2.19)
The constitutive relation between incremental PKII
stresses and GL strains are
Sij = oCijkj ockl (2.20)
Finally the equilibrium Eq. (2.14) from the principle
of virtual work using Eqs. (2.19) and (2.20) is
I oCijkl oekl 6 j dV
+ tsi 6 oij OdV
= EVW tij 6 e dV
o 1 oe ij
(2.21)
where, the external virtual work must be transformed from
t+dtc to oC. This is not applicable to conservative loading,
that is, loading that is not changed during deformation.
EVW = t+dt t t+dt uk (dA)
+ t+dt p t+dt b t+dt uk (odV) (2.22)
and this is the general nonlinear incremental equilibrium
equation of total Lagrangian formulation.
2.6 Linearization of Equilibrium Equation
The incremental strain from time t to t+dt is assumed
to be linear, i.e., ekl = ekl in Eqs. (2.11), (2.12),
(2.20), and (2,21).
For the updated Lagrangian formulation,
tSij = Cijkj tekl (2.23)
and,
ItCijkl tekl teij d
tI t
+ tij 6 tij tdV
= EVW rj 6 e tdV
(2.24)
For the total Lagrangian formulation,
oSij = oCijkj oekl (2.25)
and,
IoCijkl oekl 6 oeij odV
ddV
+ [sij ij dV
= EVW Sj e odV
o 0J o ij
(2.26)
It should be noted that the surface tractions and the
body forces in the calculation of external virtual work may
be treated configuration dependent when the structure
undergoes large displacements or large strains. If this is
the case, the external forces must be transformed to the
current configuration at each iteration [10, 11, 12].
2.7 StrainDisplacement Relationship
Using the von Karman Assumptions
The nonlinear strain terms can be simplified for the
plate or shell type structures using von Karman assumption
of large rotation.
In the mechanics of continuum the measure of
deformation is represented by the strain tensor Eij [52] and
is given by using index notation.
2Eij = ( ui,j + uj,i + k,iuk,j ) (2.27)
where,
ui = Displacement in idirection.
uij = aui / axj
xi = Rectangular Cartesian coordinate axes, i=1,2,3.
uk,iuk,j = ul,iul,j + u2,iu2,j + u3,iu3,j
The von Karman theory of plate is a nonlinear theory
that allows for comparatively large rotations of line
elements originally normal to the middle surface of plate.
This plate theory assumes that the strains and rotations are
both small compared to unity, so that we can ignore the
changes in geometry in the definition of stress components
and in the limits of integration needed for work and energy
considerations [53]. It is also assumed that the order of
the strains is much less than the order of rotations.
If the linear strain eij and the linear rotation rij
are defined as
2eij = uij + uji (2.28)
2rij = u ui (2.29)
Then the sum of Eqs. (2.28) and (2.29) gives
2(eij + rij) = 2uij (2.30)
and the subtraction of Eq. (2.29) from Eq. (2.28) gives
2(eij rij) = 2uj (2.31)
From Eqs. (2.30) and (2.31), it is concluded that
uk,j = ekj + rkj (2.32)
uk,i = eik rik (2.33)
Eq. (2.33) can be rewritten as
uk,i = eki + rki (2.34)
since eik = eki from the symmetry of linear strain terms and
rik = rki from the skew symmetry of the linear rotation
terms.
The straindisplacement Eq. (2.27) now becomes
2Eij = 2eij + (eki + rki)(ekj + rkj)
(2.35)
by substituting Eqs. (2.30) through (2.34) into Eq. (2.27).
Thus the nonlinear strain terms have been decomposed into
linear strain terms and linear rotation terms.
From the assumption on the order of strains and
rotations
eki << rki and ekj << rkj (2.36)
Thus Eq. (2.35) can be simplified as by ignoring eki and
ekj.
2Eij = 2eij + rkirkj (2.37)
The straight line remains normal to the middle surface
and unextended in the Kirchhoff assumption, but it is not
necessarily normal to the middle surface for the Mindlin
assumption. For both assumptions the generic displacements
u,v,w can be expressed by the displacements at middle
surface.
For the Kirchhoff plate [20],
u(x,y,z) = uo(x,y) z[Wo(x,y),x]
v(x,y,z) = vo(x,y) z[wo(x,y),y] (2.38)
w(x,y,z) = Wo(x,y)
where,
uo, Vo, wo = Displacements of the middle surface
in the direction of x, y, z.
u, v, w = Displacements of an arbitrary point
in the direction of x, y, z.
Now the linear strain components eij and the linear rotation
components rij can be calculated using Eqs. (2.28) and
(2.29).
ell = (ul,1 + Ul,1) = ul,1 = ,x
el2 = (l,2 + u2,1) = 1(uly + Vx)
el3 = (u1,3 + U3,1) = (Wox + Wx)
22 = i(u2,2 + u2,2) = u2,2 = Vy (2.39)
e23 = (u2,3 + u3,2) = (woy + Wy)
e33 = (u3,3 + u3,3) = u3,3 =
The rotation terms r12, r13, r23 are the rotation quantities
about the axes 3(z), 2(y) and l(x), respectively. For the
plate located in the xy plane, the rotation about z axis rl2
is much smaller than rotation about x axis r23 and y axis
r13 and therefore rl2 is assumed to be zero here. And it is
noted further that wo(x,y) is the same as w(x,y) and is a
function of only x and y so that w,3 = w,z = 0.
lrl21 << lr231 or 1r13 (2.40)
rl = (ul,1 ul,1) = 0
r12 = 1(ul,2 u2,1) = i(u, Vx) = 0
r13 = I(ul,3 u3,1) = (Wox Wx) = Wx
r22 = 1(u2,2 u2,2) = 0 (2.41)
23 = 1(u2,3 u3,2) = (Wo'y W'y) = Wy
r33 = i(u3,3 3,3) = 0
The linear strain component eij is symmetric and the linear
rotation component rij is antisymmetric.
eij = eji
rij = rji (2.42)
The strain components from Eq. (2.37) can be rewritten using
Eqs. (2.39) and (2.41).
Ex = ell + (r112 +
yy = e22 + (r122 +
Ezz = e3 + 1(r132 +
Exy = el2 + 2(rllrl2
Exz = e3 + (rllrl3
Eyz = e23 + U(rl2r13
r21 + r31 ) = ell + r31
r22 + r322) = e22 + r32
r232 + r332) = 1(r132 + r232) = 0
+ r21r22 + r31r32) = el2 + r31r32
+ r21r23 + r31r33)
+ r22r23 + r32r33) (2.43)
Egz term is assumed to be zero because it does not have the
linear term. Exz and Eyz terms are transverse shear terms
which can be ignored for thin plate. Then Eq. (2.43) can be
rearranged as follows using Eqs. (2.41) if all the zero
terms are removed.
Exx = ell + r312 = ell + (W,x)2
Eyy = e22 + r322 = 22 + )2
Ey = el2 + Ir31r32 = e2 + (W'x) (,y)
Exz = e13
Eyz = e23
(2.44)
33
Thus the decomposition of exact strain components has been
done using the Kirchhoff plate assumptions (2.38) and the
von Karman assumption (2.40) on the magnitude of rotation.
It is noted that all the inplane displacement gradients in
nonlinear strain terms are ignored through von Karman
assumptions [20]. This fact will be applied in chapter 5.
CHAPTER 3
THREE DIMENSIONAL LINK ELEMENT
3.1 Element Description
The link element used here is based on the two
dimensional element developed by Cleary [54]. The link
element is based on the following assumptions.
Any normal compressive force is transferred to the
other side of the link without any loss. To facilitate this,
a very limited amount of loss through displacement should be
allowed. Currently, this limited displacement is defaulted
to .001 units, while it is a input parameter. The link
separates in response to any net tension, losing its normal
stiffness.
To discuss the shear force transfer, some definitions
for friction are needed. The force to start one body sliding
along the other body is called the static friction force.
The force to keep it moving is the kinetic friction force.
There are two corresponding coefficients of friction, static
friction coefficient and dynamic friction coefficient, where
the static friction coefficient will generally be greater
than the dynamic friction coefficient.
Two laws of friction were used in the link element. The
first law is that the frictional force is proportional to
the normal force, with the constant of proportionality being
the friction coefficient. The second law is that friction
does not depend on the apparent area of the connecting
solids, i.e., it is independent of the size of the bodies.
The shear force is transferred through friction. The
uncertainty in friction is the factor which limits the
overall accuracy of the calculation. Therefore, it is
assumed that the static friction coefficient is proportional
to the dynamic friction coefficient. For nonmetallic
materials, the ratio of dynamic coefficient to static
coefficient is about 0.75.
The link element is composed of two surfaces. If the
shear force is less than or equal to the static friction
force, i.e., coefficient of friction times the normal force,
the shear force is balanced by the friction force and the
total force is transferred. This is shown in Fig. 31. But
if the shear force is greater than the static friction
force, one surface of the link element will move along the
other surface. In this case there will be a dynamic friction
force which is less than the shear force. This dynamic
friction force can only resist a portion of the shear and
the system is not in static equilibrium. Therefore, if the
shear force is greater than the static friction force, the
link element will lose its shear stiffness. This can also be
modeled with a body on roller and spring as shown in Fig. 3
2. The spring model of the link element is shown in Fig. 3
3.
The link element here has four nodes and each node has
three translational degrees of freedom in local u, n, and
wdirections. The total number of element degrees of freedom
is 12. The element degrees of freedom are shown in Fig. 34.
The equivalent "strain" for the link element is defined
as the average deformation at the center of the element. The
average deformation corresponding to the translational
degrees of freedom, i.e., uo, vo, and wo, can be directly
calculated from the joint displacements by averaging the
difference in nodal displacements at the ends of element in
local u, n, and wdirection in turn. The relative rotation
at the center of the element, ro, can be found using nodal
displacements in local ndirection and the element length as
shown in Fig. 35. This angle is not an "average" value but
the "relative" change in angle of the center line due to
rotation.
The two joint parameters must be introduced. These are
kn, the unit stiffness normal to the joint, and ks, the unit
stiffness along the joint.
The offdiagonal term kns to account for dilatation
during shearing is ignored because this joint element will
model the dry joint between concrete box girder bridge
segments. No significant dilation is expected in this case.
Some values of kn and ks were reported in geotechnique
area [3]. As the values are those for natural joints, they
do not directly apply to this case.
From the test results [7], it can be seen that the
shear stiffness of dry joint ranges from 70,000 to 286,000
psi per inch at ultimate. In case of frictional strength,
this can be interpreted as linear behavior between the
origin and the ultimate point.
It seems reasonable that the normal stiffness of the
element, kn, is assumed to be stiffer than the connected
material by the order of 103 to transfer the normal force
without any significant loss. The forces are either totally
transferred in compression or totally lost in tension. The
latter case has no problem related to the value of kn.
The shear stiffness parameter is more difficult to
define. The data available is so limited that even a
statistical treatment cannot be done. But in the analysis of
structural behavior up to the ultimate, these properties do
not have great influence because the forces are transferred
through friction.
The shear stiffness becomes zero upon sliding. But
there is some 'residual' shear force. This 'residual' force
is equal to friction force. Therefore, if shear displacement
is more than the displacement just before the sliding the
shear stiffness is set to be zero.
Ff
P, N = External forces.
F = Friction force.
f
m = Friction coefficient
1) P < or = mN then P = Friction Force. In Equilibrium.
2) P > mN then the body moves but the frictional force
mN is acting against the other body.
Fig. 31 Friction Force
FRICTIONAL SPRING WITH SHEAR STIFFNESS
Force in spring = mN
N
%AAN
BEFORE
SLIDING
SFRICTIONAL SPRING WITH ZERO STIFFNESS
N
A F
mN
mN < F
AFTER SLIDING
Fig. 32 Spring Model of Friction Force
__ __
Kn = Zero under tension.
= Very large under
compression.
SPRING MODEL FOR NORMAL FORCES
SPRING MODEL FOR SHEAR FORCES
Fig. 33 Spring Model of Link Element
11
K
Fig. 34 Element Degrees of Freedom of Link Element
L
12
I
34
CENTER OF ELEMENT
K
UK
NoLNo
uo = [(UK + UL) (UI + UJ)] / 2
ro= p a
S[(VKVL) (VJVI)] / L
VK
O
VI
o o
Fig. 35 Element "Strain"
uo_ZL
3.2 Formation of Element Stiffness
There are four nodes per element. Each node has three
degrees of freedom corresponding the translational
displacements in u, n, and wdirection resulting in 12
element degrees of freedom as shown in Fig. 34. The element
stiffness is derived directly from the physical behavior of
the element described in section 3.1.
The mathematical symbol {} is used for a column vector
and [] for a matrix. The nodal displacement column vector
{q)(12) is composed of 12 translational nodal displacements
corresponding to the 12 element degrees of freedom.
(q} = { ui vi wi Uj Vj Wj uk vk wk ul vl wl )T
The "strain" is defined as the average deformation at
the center of the element as shown in Fig. 35. The "strain"
column vector {e)(4) is
{e} = { uo vo wo ro )T
where,
uo = ( uk + ul ) / 2 ( u + j ) / 2
o = vk + 1 ) / 2 ( + Vj ) / 2
w = (wk + 1 ) / 2 ( i + j ) / 2
ro = ( vk v ) / L ( vj vi ) / L
where,
L = The length of the element.
uo, vo, wo = Average nodal displacements in local
u, n, wdirections.
ro = The relative angle change about local z axis.
Therefore the relationship between "strain" and nodal
displacements is
(e}(4) = [B](4,12) (q)(12)
The [B](4,12) matrix which gives strains due to unit
values of nodal displacements is
0.5 0.0 0.0 0.5 0.0 0.0 0.5 0.0 0.0 0.5 0.0 0.0
0.0 0.5 0.0 0.0 0.5 0.0 0.0 0.5 0.0 0.0 0.5 0.0
0.0 0.0 0.5 0.0 0.0 0.5 0.0 0.0 0.5 0.0 0.0 0.5
0.0 1/L 0.0 0.0 1/L 0.0 0.0 1/L 0.0 0.0 1/L 0.0
The "stress" is defined as the normal and shear stress
per unit of area. {s} is the average stress on the surface
due to the two nodal forces exerted in the plane of the
surface. This stress is in equilibrium with the stress on
the other surface of the element as shown in Fig. 36. m is
the moment of the nodal forces on one surface in local n
direction about the center of the element. This moment is
also balanced by the moment of the nodal forces on the other
surface of the element. This moment is used to define the
distribution of the normal stress of the element as shown in
Fig. 37.
Pnl
fl1
P...'
uk
~ ~
~ ~ ~
~ ~ ~
~ ~
uj nj
Pj J ~pnj
W Stresses
 Nodal Forces
Local Coordinate System
Fig. 36 Nodal Forces and Stresses of Link Element
/
ni
W. . w
0o 00 p
FORCE TRANSFER THROUGH ONE EDGE
OF THE LINK ELEMENT
P
L K
Vo ro
2 Vo
CENTER OF ELEMENT
p J
Fig. 37 Element "Strain", m
The "stress" column vector (s}(4) is
(s) = { Sx, sn, Sz, m )
The "stressstrain" relationship is
(s)(4) = [E](4,4) {e)(4)
where,
kx 0 0 0
[E] = 0 kn 0 0
0 0 kz 0
0 0 0 km
where km can be related to kn using the definition of
the moment m, i.e.,
m = (Sn)(L)(t)(0.5)(L)
= kmro = (km)(vo/(0.5)(L))
Thus, km = (knVo)(L)(t)(0.5)(L) / [Vo/(0.5)(L)]
= (0.25)(t)(kn)(L3)
where, t = Element thickness.
This assumes that there is no coupling between the
shear stress and normal stress.
The element nodal force column matrix (P)(12) is
composed of the 12 nodal forces shown in Fig. 36.
(P) Pui Pni Pwi Puj Pnj Pwj Puk Pnk Pwk
Pul nl Pwl T
Stress can then be related to nodal forces using the
definition of stress and force equilibrium between the two
surfaces of the element.
By the definition of stress,
sn = (1/Lt)( Pnk + Pnl ) (3.1)
sx = (1/Lt)( Puk + Pul ) (3.2)
sz = (1/Lt)( Pwk + Pwl ) (3.3)
m = Pnk(0.5)(L) Pnl(0.5)(L) (3.4)
where, Lt = (L)(t)
By force equilibrium of the two surfaces,
Pi = P and Pj = Pk (3.5)
To express the element nodal forces in terms of the
stress, we use Eqs. (3.1) through (3.5) to find the force
recovery matrix [FR]. [FR] gives the nodal forces in
equilibrium with the element stresses.
From (Eq. (3.1) + Eq. (3.4)),
2Pnk = (L)(t) (sn) + 2(m)/L
Pnk = 0.5(L) (t) (sn) + (l/L)(m)
From Eq. (3.5),
Pj = Pk
nj = Pnk
= 0.5(L)(t)(sn) (1/L)(m)
From Eq. (3.1),
Pn1 = (L)(t)(sn) Pnk
= (L)(t)(sn) ( 0.5(L)(t)(sn) + (1/L)(m))
= 0.5(L) (t) (n) (1/L) (m)
From Eq. (3.5),
Pni = Pn1
= 0.5(L)(t)(sn) + (1/L)(m)
From the assumption that Puk = Pul and Eq. (3.2),
Puk = (L)(t)(sx)/2
Pul = (L)(t)(sx)/2
From Eq. (3.5),
Pui = ul = (L)(t)(sx)/2
Puj = uk = (L)(t)(sx)/2
From the assumption that Pwk = Pwl and Eq. (3.3),
Pwk = (L)(t)(sz)/2
Pwl = (L)(t)(sz)/2
From eqn 5,
Pwi = wl = (L)(t)(sz)/2
Pwj = Pwk = (L)(t) (s)/2
Therefore, the forcestress relationship is
{P)(12) = [FR](12,4) (s)(4)
where the force recovery matrix [FR](12,4) is
[FR] =
Lt/2
0
0
Lt/2
0
0
Lt/2
0
0
Lt/2
0
0
0
Lt/2
0
0
Lt/2
0
0
Lt/2
0
0
Lt/2
0
0
0
Lt/2
0
0
Lt/2
0
0
Lt/2
0
0
Lt/2
0
1/L
0
0
1/L
0
0
1/L
0
0
1/L
0
And this relationship is further expanded using the
stressstrain relationship and the strainnodal displacement
relationship as follows.
(P}(12) = [FR](12,4)
= [FR](12,4)
= [Bt](12,4)
[E](4,4) {e)(4)
[E](4,4) [B](4,12)
[E](4,4) [B](4,12)
Then finally this can be symbolized as equilibrium
equation.
(q)(12)
{q)(12)
{P)(12)= [Ke](12,12) (q}(12)
where [Ke] = [Bt][E][B]
Here it is noted that [FR] = [Bt] and [Ke] = [Bt][E][B] just
as in the case of common finite element method.
The final element stiffness matrix [Ke] is
(L)t
4
kx
0
0
kx
0
0
kx
0
0
kx
0
0
0
2kn
0
0
0
0
0
0
0
0
2kn
0
0
0
kz
0
0
kz
0
0
kz
0
0
kz
kx
0
0
kx
0
0
kx
0
o
kx
0
0
0
0
0
0
2kn
0
0
2kn
0
0
0
0
0
0
kz
0
0
kz
0
0
kz
0
0
kz
kx
0
0
kx
0
0
kx
0
0
kx
0
0
0
0
0
0
2kn
0
0
2kn
0
0
0
0
0
0
kz
0
0
kz
0
0
kz
0
0
kz
kx
0
0
kx
0
0
kx
0
0
kx
0
0
0
2kn
0
0
0
0
0
0
0
0
2kn
0
0
0
kz
0
0
kz
0
0
kz
0
0
kz
This matrix can be
standard rotation.
rotated to any direction using the
3.3 Solution Strategy
The structural stiffness changes because of the slip
and debonding of the link. Therefore, the process of the
resistance of the total structure physically becomes
nonlinear. Correspondingly, special solution techniques for
nonlinear behavior are needed.
This can be done using the iterative solution technique
with initial stiffness or tangent stiffness. The latter can
be formed by assembling the structural stiffness at the
beginning of each iteration and this converges faster than
the initial stiffness.
A third solution strategy for this case is eventto
event technique which is usually employed for the linear
stiffnesses between any two "events," which are defined as
the intersection point between two linear segments. This
also provides means of controlling the equilibrium error.
Any significant event occurring within any element
determines a substep. The tangent stiffness is modified in
each substep, and hence, the solution closely follows the
exact response.
3.4 Element Verification
3.4.1 SIMPAL
The finite element analysis program SIMPAL [55], is
used to implement and verify the element formulation. SIMPAL
was chosen for the initial implementation because that was
the original implementation done by Cleary [54]. This way,
the 3D aspects could be implemented and verified using
Cleary's original program. A table of the element
verification is shown in Fig. 38 and Fig. 39.
LOADING
10
10
Y
1 2
NODE THEORY SIMPAL ERROR
DISP 2 .1333 .1333 .000
DISP 4 .1333 .1333 .000
STRESS N/A 80 80 .000
NODE THEORY SIMPAL ERROR
DISP 2 .1333 .1333 .000
DISP 4 .1333 .1333 .000
STRESS N/A 80 80 .000
NODE THEORY SIMPAL ERROR
DISP 2 .1333 .1337 .003
DISP 4 .1333 .1337 .003
STRESS N/A 80 80 .000
* NODE 2 Y DISP =0.101704 Z DISP = 0.868405
NODE 4 Z DISP =0.868405 Y DISP =0.101704
2 2
SQRT((.1017) + (.08684) )=0.1337
Fig. 38 Link Element Test Using SIMPAL
3
3
THICKNESS = .25
Ks = 3E6
Kn = 6E6
RESULTS
LOADING
4 6
Aft Am
.P ......
3 5
1 _~~
RESULTS
8
Aft
Ii
MP, 
Adh 
10
1 I
1 3 5 7 Y
1 3 5 7 Y
Z
2
4 6
3 5
Fig. 39 Combined Test Model for SIMPAL
Z
2
~araooa a;a8oo~ uru 
iam i io n

HtS~ ^~ k
10
8
10
3.4.2 ANSR
The test examples used are the same as those used in
the initial element verification using SIMPAL. The results
from ANSR [56] are exactly the same as those from SIMPAL.
The link element was tested further using a modeled membrane
element composed of 22 truss elements as a membrane element
was not available at the time of element verification in
ANSR. The results are shown in Table 31 and the structures
used are shown in Fig. 310 and Fig. 311.
Table 31 Displacements of
for ANSR
Truss Model
Node Truss Truss Diff.
No. only w/ LINK (%)
10x .1027e4 .1049e4 2.2
10y .1990e5 .2010e5 0.9
llx .9017e5 .9211e5 2.2
11y .4906e6 .4973e6 1.4
12x .9915e5 .1049e4 2.2
12y +.9742e6 +.9703e6 0.4
LINK ELEMENT
20
.2
20
.2
20
Fig. 310 Combined Test Model for ANSR
Fig. 311 Truss Model for ANSR
CHAPTER 4
LINEAR SHELL ELEMENT
4.1 Element Description
The shell element is formulated through a combination
of two different elements, the membrane element and the
Mindlin plate bending element [57].
The Mindlin plate element is different from the
Kirchhoff plate element in that the former allows transverse
shear deformation while the latter does not.
The common portions of the formulation of two elements
are
1. Formation of the shape functions.
2. Formation of the inverse of Jacobian matrix.
These processes can be done at the same time. The four
to ninenode shape functions and their derivatives in rs
space can be formed and then transformed into xyspace
through the inverse of Jacobian matrix.
4.2 Formulation of Shape Functions
The formulation of shape functions starts with three
basic sets of shape functions shown in Fig. 41.
1. The bilinear shape functions for fournode element.
2. The linearquadratic shape functions for nodes
five to eight of the eightnode element.
3. The bubble shape function for node nine of
ninenode element.
These shape functions can be formulated directly from
the local coordinates of the element nodes through the
multiplication of the equations of the lines which have zero
values in the assumed displacement shapes and the scale
factor to force the shape function value to one at the node
for which the shape function is formed. The derivative of
each shape function with respect to r and s is then
evaluated from the shape function expressed in terms of r
and s.
If node nine exists, the value at node nine of shape
functions one to eight must be set to zero. The value of the
bilinear shape functions for a fournode element at node
nine is one fourth and the value of the linearquadratic
shape functions for the five to eightnode element at the
node nine is one half. This can be forced to zero using the
bubble shape function of the ninenode element because this
shape function has the value of one at node nine and zero at
all other nodes. Therefore the modification is the
subtraction of one fourth of the value the bubble shape
function has at node nine from the each shape function for
the corner nodes and the subtraction of one half of the
value of the bubble shape function of node nine for the
nodes five to eight, whichever exists.
If any of the center nodes on the edge of the element
(any one of nodes five to eight) exists, the bilinear shape
functions of fournode element must be modified further
because the value at center of the edge is one half in those
bilinear shape functions. This can be done by subtracting
one half of the linearquadratic shape function for the
newly defined center node on the edge of the element from
the bilinear shape functions of the two adjacent corner
nodes. The value of any five node shape functions at the
corner node is zero. Therefore, no further consideration is
needed except for the shifting of the shape functions in the
computer implementation. These processes are shown in Fig.
42.
If any of the linearquadratic shape functions of nodes
five to eight is missing, all the linearquadratic shape
functions thereafter and the bubble shape function must be
shifted to the proper shape function number. For example, if
linearquadratic shape function five is missing, then the
shape functions six to eight must be shifted to five through
seven and the bubble shape function must be shifted to the
node eight because all of the linearquadratic shape
functions have been defined and numbered as shape functions
for the nodes five through eight and the bubble shape
function for the node nine.
Four Node Element
Shape Function for Corner Node
Five Node Element
Shape Function for Edge Center Node
Nine Node Element
Shape Function for Element Center Node
Fig. 41 Three Basic Shape Functions
SF 1
SF 2
SF 3
SF 3
SF 4 = (SF 1) (1/4) (SF 3)
SF 5 = (SF 2) (1/2) (SF 3)
Fig. 42 Formation of Shape Functions
4.3 The Inverse of Jacobian Matrix
While the generic displacements are expressed in terms
of rscoordinate, the partial differential with respect to
the xycoordinate is needed for the calculation of strain
components. Thus the inverse of the Jacobian matrix must be
calculated. This can directly be found from the chain rule
using the notation (a,b) defined as the partial differential
of function a with respect to the variable b for simplicity.
f,x = (f,r)(r,x) + (f,s)(s,x)
f,y = (f,r)(r,y) + (f,s)(s,y)
In matrix form,
f,x r,x s,x f,r Jll1 J121 f,r
f,y r,y s,y f,s J211 J221 f,s
The inverse of Jacobian matrix
But the terms in the inverse of the Jacobian matrix are not
readily available because the rscoordinate cannot be solved
explicitly in terms of xycoordinate. On the other hand, for
the isoparametric formulation, the geometry is interpolated
using the nodal coordinate values(constants) and the
displacement shape functions in terms of r and s. Thus the
generic coordinate x and y can be expressed in r and s
easily and explicit partial differentials of x and y with
respect to r and s can be performed. Therefore the Jacobian
matrix is computed and then inverted.
The Jacobian matrix is derived by the chain rule.
f,r = (f,x)( x,r) + (f,y)( y,r)
f,s = (f,x)( x,s) + (f,y)( y,s)
In matrix form,
f,r x,r y,r f,x Jll J12 f,x
f,s x,s y,s f,y J21 J22 f,
Jacobian matrix
nn
Let s be .
i=1
where nn = number of nodes (4 to 9).
From geometric interpolation equations,
x = Z fi*xi
y = Z fi*Yi
The terms in the Jacobian matrix are
J11 = x,r = (z fi*xi),r = Z ((fi,r) xi)
J12 = y,r = (Z fi*yi),r = Z ((fi,r) yi)
J21 = x,s = (Z fi*xi),s = Z ((fi,s) xi)
J22 = y,s = (Z fi*yi),s = Z ((fi,s) yi)
xi, yi are coordinate values of the element and are
constants and therefore can be taken out of the partial
differentiation.
The inverse of twobytwo Jacobian matrix can be found
i
Jl1 = r,x = J22 / det(J)
i
J12 = s,x = J12 / det(J)
I
J21 = r,y = J21 / det(J)
i
J22 = s,y = J11 / det(J)
where det(J) = J11J22 12J21
4.4 Membrane Element
The formulation of the membrane element used for the
implementation follows the procedure shown on pages 115
through 118 in reference [57]. The ( ) symbol will be used
for the column vectors.
Nodal displacements are the nodal values of two in
plane translations and denoted as {ui vi}T. The generic
displacements are defined as two translational displacements
at a point and denoted as { u v )T. By the word generic it
is meant that the displacement is measured at an arbitrary
point within an element. The generic displacements u and v
can be calculated using shape functions. The shape function
is a continuous, smooth function defined over the closed
element domain and is differentiable over the open domain of
the element. The shape function is also the contribution of
displacement of a node for which the shape function has been
defined to the generic displacement. Thus the generic
displacement at an arbitrary point can be found by summing
up all the contributions of all the nodes of the element.
The displacement interpolation equations are
u = Z fi Ui
v = Z fi Vi
In the isoparametric formulation the geometry is
interpolated using the same shape function assumed for the
displacement interpolation.
Therefore, the geometry interpolation is
x = Z fi xi
y = Z fi Yi
where,
fi = Shape function for node i.
xi, Yi = Coordinates of node i.
ui, vi = Displacements at the node i.
u, v = Displacements at an arbitrary point within
an element.
The three inplane strain components for a membrane
element are
{ E ) = ( ex Cy 7xy )
These strain components can be found through the
partial differentials of the generic displacements with
respect to xycoordinates.
ex = ux
Ey = v,y
7xy = uy + VX
Using the inverse of the Jacobian matrix, the strain
components can be evaluated.
eX
= u,x
= (u,r)(r,x) + (u,s)(s,x)
= (u,r)(J111) + (u,s)(J121)
= ((Efiui),r)(J111) + ((EfiUi),s)(J121)
= Z[(fi,r)(r,x) + (fi,s)(s,x)] ui
y
= v,y
= (v,r)(r,y) + (v,s)(s,y)
= (v,r)(J211) + (v,s)(J221
= ((.fivi),r)(J211) + ((Zfivi),s)(J221)
= [(fi,r)(r,y) + (fi,s)(s,y)] vi
7xy
= u,y + v,x
= [(u,r)(r,y) + (u,s)(s,y)] + [(v,r)(r,x) + (v,s)(s,x)]
= [((Zfiui),r)(J211) + ((Zfiui),s)(J221)]
+ [((Zfivi),r)(J11) + ((Zfivi),s)(J121)]
= Z[(fi,r)(r,y) + (fi,s)(s,y)] ui
+ Z[(fi,r)(r,x) + (fi,s)(s,x)] vi
New notations are introduced here to simplify the
equations. These are ai and bi and defined as follows:
ai = (r,x)(fi,r) + (s,x)(fi,s) = fi,x
bi = (r,y)(fi,r) + (s,y)(fi,s) = fiY
Then the strain terms above become
ex = Zaiui = Zf,x ui
Ey = Zbivi = Zfiy vi
7xy = Zbiui + Zaivi = Zfi,y ui + Zfi,x vi
In matrix form,
ex ai 0 ui
yE = 0 bi vi
7xy bi ai
In symbolic form,
[E] = E[Bi]Cqi]
where,
ai 0 fix 0
[Bi] = 0 bi 0 fi',
bi ai fi'y fix
and,
ui
[q] = Vi
vi
Therefore the strain at an arbitrary point within an
element is
[e] = [Bl][q1] + [B2][q2] + ... + [B9][q99
= [ B1 B2 B3 B4 B5 B6 B7 B8 B9 ]
ql
q2
q3
q4
95
q6
q7
q8
q9
The size of the vectors and matrix are
[e(3,1)] = [B(3,18)][q(18,1)]
In the actual calculation, this can be done by summing
up the [Bij[qi] over all the nodes for the given coordinates
of the point under consideration, i.e., the coordinates of
one of the integration points.
The stresses corresponding to the strains are
{ a } = { ax ay rxy )T
The stressstrain relationship of an isotropic material
is
Ell E12 0
[E] = E21 E22 0
0 0 E33
where,
E11 = E22 = E / ( 1 2 )
E12 = E21 = pE / ( 1 p2)
E33 = G
where,
E = Young's modulus
p = Poisson's ratio
G = shear mQdulus = E / ( 2*(1+p))
The integration over volume must be introduced for the
calculation of the element stiffness and equivalent nodal
loads of distributed loads or temperature effects. As the
thickness of the element is constant, the integration over
volume can be changed to the integration over area.
The element stiffness related to the degrees of freedom
of the node i can be calculated through the volume
integration of Bi(2,3)E(3,3)Bi(3,2).
[Ki] = BT E Bi dV
As [B] and [E] are constant about z, the integration
through the thickness from t/2 to t/2 can be performed on z
only and yields
T
[Ki] = [ Bi(t)EBi] dA
A B
T 
= [ Bi E Bi] dA
A
where,
E = tE
The size of membrane element stiffness is 18 by 18.
[K] =
[ E ] [ B1 B2 B3 ... B8 B9 ] dV
(3,3)
(3,18)
B8
B9
(18,3)
Equivalent nodal loads due to body forces on the
membrane element are calculated as
Pb Jv fTbdV =11 11
fTbIJI dr ds
in which {b} = ( 0 0 bz T or { 0 b 0 }T or { bx 0 0 } in
accordance with the direction of gravity in the coordinate
system used. The nonzero quantities bx, by, or bz represent
the body force per unit area in the direction of
application.
Equivalent loads caused by initial(temperature) strains
are
PO = BTEEO dV
JV
1 1
BTEcIJI dr ds
S1 Ji
where,
(C0)= exxO yyO 0 0 0 )T
= {aAT aAT 0 0 0 )
4.5 Plate Bending Element
The formulation of the plate bending element used for
the implementation has followed the procedures shown on
pages 217 through 221 in reference [57]. The { } symbol will
be used for the column vectors.
Many plate bending elements have been proposed. The
most commonly used are Kirchhoff plate elements and Mindlin
plate elements.
Kirchhoff theory is applicable to thin plates, in which
transverse shear deformation is neglected. The assumptions
made on the displacement field are
1. All the points on the midplane(z = 0) deform only
in the thickness direction as the plate deforms in
bending. Thus there is no stretching of midplane.
2. A material line that is straight and normal to the
midplane before loading is to remain straight and
normal to the midplane after loading. Thus there is
no transverse shear deformation (change in angle
from the normal angle).
3. All the points not on the midplane have displacement
components u and v only in the x and y direction,
respectively. Thus there is no thickness change
through the deformation.
Strain energy in the Kirchhoff plate is determined entirely
by inplane strains ex, Cy, and 7xy which can be determined
by the displacement field w(x,y) in the thickness direction.
The interelement continuity of boundarynormal slopes is not
preserved through any form of constraint.
Mindlin theory considers bending deformation and
transverse shear deformation. Therefore, this theory can be
used to analyze thick plates as well as thin plates. When
this theory is used for thin plates, however, they may be
less accurate than Kirchhoff theory because of transverse
shear deformation. The assumptions made on the displacement
field are
1. A material line that is straight and normal to the
midplane before loading is to remain straight but
not necessarily normal to the midplane after
loading. Thus transverse shear deformation (change
in angle from normal angle) is allowed.
2. The motion of a point on the midplane is not
governed by the slopes (w,x) and (w,y) as in
Kirchhoff theory. Rather its motion depends on
rotations Ox and 0 of the lines that were normal to
the midplane of the undeformed plate. Thus 0, and 0
are independent of the lateral displacement w, i.e.,
they are not equal to (w,x) or (w,y).
It is noted that if the thin plate limit is approached, xz
= 7yz = 0 because there is no transverse shear deformation.
In this case the angles 0x and 0y can be equated to the
(w,x) and (w,y) numerically but the second assumption still
holds.
The stiffness matrix of a Mindlin plate element is
composed of a bending stiffness [kb] and a transverse shear
stiffness [ks]. [kb] is associated with inplane strains ex,
cy, and 7xy. [ks] is associated with transverse shear
strains 7xz and 7yz. As these two groups of strains are
uncoupled, i.e., one group of the strains do not produce the
other group of strains, the element stiffness can be shown
as [82]
[k] = ( BbEBb ) dA + (BsEBs) dA
because BbEBs = BsEBb = 0 from uncoupling (corresponding E =
0). Each integration point used for the calculation of [kg]
places two constraints to a Mindlin plate element,
associated with two transverse shear strains lyz and yzx. If
too many integration points are used, there will be too many
constraints in transverse shear terms, resulting in locking.
Therefore, a reduced or selective integration can prevent
shear locking. Or, the transverse shear deformation can be
redefined to avoid such locking.
For example, a bilinear Mindlin plate element responds
properly to pure bending with either reduced or selective
integration. But full twobytwo integration is used for
pure bending, shear strains appear at the Gauss points as
shown in Fig 43. As the element becomes thin, its stiffness
is due almost entirely to transverse shear. Thus, if fully
integrated, a bilinear Mindlin plate element exhibits almost
no bending deformations, i.e., the mesh "locks" against
bending deformations.
Nodal displacements for the plate bending consist of
one outofplane translation and two outofplane rotations
and are denoted as { wi 0xi 8yi )T. The rotations are chosen
independently of the transverse displacement and are not
related to it by differentiation. Thus the transverse shear
strains jxz and 7yz are considered in the formulation
resulting in five strain components. The generic
displacements are defined as three translational
displacements and denoted as { u v w }T. By the word generic
it is meant that the displacement is measured at an
arbitrary point within an element. These generic
displacements are different quantities from the nodal
displacements and therefore must be related to the nodal
displacements.
The generic displacements u and v can be calculated as
functions of the generic outofplane rotations using the
small strain(rotation) assumption. The relationship between
generic displacements and rotation is shown in Fig 44.
u = zBy
v = ZBx
The generic displacements Ox and By can be found using the
assumed displacement shape functions and the corresponding
nodal displacements Oxi and 0yi.
The generic displacement w does not need any conversion
because it corresponds to the nodal displacement wi.
In the isoparametric formulation the geometry is
interpolated using the same shape function assumed for the
displacement interpolation.
The displacement interpolation is
8x = Z fi 0xi
6y = Z fi 0yi
w = Z fi Wi
Similarly, the geometric interpolation is
x = Z fi xi
y = Z fi Yi
Zero Shear
Strain
One Point Gauss Integration
I Nonzero
Shear Strain
Two Point Gauss Integration
Fig. 43 Shear Strains at Gauss Point(s)
i~Amlh
Z
+ u
y
X
Positive small rotational angle about yaxis gives
positive generic displacement in xdirection ( u ).
Shown is xzplane.
Positive small rotational angle about xaxis gives
negative generic displacement in ydirection ( v
Shown is yzplane.
Fig. 44 Displacements due to Rotations
where,
fi = shape function for node i
xi, Yi = coordinates of node i
Therefore,
u = Zoy = z z fi 6yi
v = zOx = z Z fi Oxi
W = E fi Wi
The five strain components for plate bending element
are { ex ey 7xy 7xz Iyz )T. These strain components can be
found through the partial differentials of the generic
displacements with respect to xycoordinates.
e = uI,
Ey = v,y
7xy = uy + V,X
7XZ = u,z + W,X
7yz = vz + w,y
Using the inverse of the Jacobian matrix found, the
strain components can be evaluated.
Ex
= u,X = (ZBy),X
= (u,r)(r,x) + (u,s)(s,x)
= (u,r)(J111) + (u,s)(J121)
= ((ZzfiOyi),r)(J~l1) + ((zzfiyi) ,) (J121)
= z Z[(fi,r)(r,x) + (fi,s)(s,x)] Oyi
= v,y = (zOx),y
= (v,r)(r,y) + (v,s)(s,y)
= (v,r)(J211) + (v,s)(J221)
= ((zzfioxi),r)(J211) + ((zZfi0xi),s)(J221)
= z Z[(fi,r)(r,y) + (fi,s)(s,y)] Oxi
7xy
= u,y + v,x = (zoy),y + (zox),X
= [(u,r)(r,y) + (u,s)(s,y)] + [(v,r)(r,x) + (v,s)(s,x)]
= [((zzfioyi),r)(J211) + ((zZfiyi),s)(J221)]
+ [((zZfisxi),r)(Jll1) + ((zzfixi),s) (J121)]
= z Z[(fi,r)(r,y) + (fi,s)(s,y)] Oyi
z Z[(fi,r)(r,x) + (fi,s)(s,x)] exi
7xz = (u,z) + (w,x)
= (ZOy),Z + w,x = 9y + w,x
= Zfisyi + [(w,r)(r,x) + (w,s)(s,x)]
= Zfioyi + [((Zfiwi),r)(J111) + ((Zfiwi),s) (J121)
= Zfieyi + [((Zfi,r)wi) (J111) + ((Zfi,s)wi)(J121)]
= ZfiOyi + z[(fi,r)(r,x) + (fi,s)(s,x)] wi
lyz = (v,z) + (w,y)
= (z0x),z + w,y = (OX) + w,y
= (Zfioxi) + [(w,r)(r,y) + (w,s)(s,y)]
= (Zfioxi) + [((Zfiwi),r)(J211) +((Zfiwi),s)(J221)]
= (Zfioxi) + [((Zfi,r)wi)(J211) +((Zfi,s)wi)(J221)]
= (Zfi0xi) + Z[(fi,r)(r,y) +(fi,s)(s,y)] Wi
New notations are introduced here to simplify the
equations. These are ai and bi and defined as follows:
ai = (r,x)(fi,r) + (s,x)(fi,s) =
bi = (r,y)(fi,r) + (s,y)(fi,s) =
Then the strain terms above become
f.l
i1,x
ex = z aioyi
Cy = z Zbixi
7xy = z Zbiyi z Zaixi
7xz = fi0yi + Zaiwi
Tyz = fioxi + Zbiwi
In matrix form,
eX
fy
7xy =
7XZ
Lyz
In symbolic form,
[e] = E[Bi][qi]
0
zbi
zai
0
fi
i
zai
0
zbi
fi
0
wi
exi
Oyi
where,
[Bi] =
or,
[Bi] =
0
 zbi
 zai
0
fi
i
0
0
0
fi,x
fi'y
zai
0
zbi
fi
0
0
 zfi,y
 zfi,x
0
fi
zfi,x
0
zfi,Y
fi
0
and,
[qi]
wi
ixi
Oyi
Therefore the strain at an arbitrary point within an
element is
[e] = [B][ql] ] + ... + [B2[ + + 9][q9]
= [ B1 B2 B3 B4 B5 B6 B7 Bg B9 ]
q91
92
93
94
95
96
97
98
99
The size of the vectors and matrix are
[E(5,1)] = [B(5,27)][q(27,1)]
In the actual calculation, this can be done by summing
up the [Bi][qi] over all the nodes for the given coordinates
of the point under consideration, i.e., the coordinates of
one integration point.
The stresses corresponding to the strains are
{ o ) = ( ax ay rxy rxz yz )TT
The stressstrain relationship of an isotropic material
is
E11 E12 0 0 0
E21 E22 0 0 0
[E] = 0 E33 0 0
0 0 0 E44 0
0 0 0 0 E55
where,
E11 = E22 = E / ( 1 p2 )
E12 = E21 = pE / ( 1 p2)
E33 = G
E44 = E55 = G / 1.2
where,
p = Poisson's ratio
G = shear modulus = E / ( 2*(1+p))
The form factor 1.2 for the E44 and E55 terms is
provided to account for the parabolic distribution of the
transverse shear stress rzx over a rectangular section. This
form factor 1.2 can be shown from the difference in
deflections of a cantilever beam at its free end [58].
Let a beam have a rectangular cross section of
dimensions b by t with a length of L. If P is the transverse
shear force, then the parabolic distribution of the
transverse shear stress rz, is
r'X = (3P/2bt3)(t2 4z2)
where z = 0 at the neutral axis.
Then the transverse shear strain energy from the parabolic
distribution can be calculated by
Us = (1/2)V (zx 2/G) dV
= (1/2) [((3P/2bt3)(t2 4z2))2 / G] dV
= [(1/2)(3P/2bt3)2]/G (t2 4z22 dAdz
= (area)[(1(/)(3P/2bt3)2]/G (t2 4z2)2 dz
= bL[(1/2)(3P/2bt3)2]/G (t2 4z22 dz
= 1.2(P2L/btG)/2
While the transverse shear strain energy from the constant
distribution is
Us = (1/2) I(rx2/G) dV
= (1/2) ((P/bt)2 /G) dV
= [(1/2)(P/bt)2/G] (btL)
= (p2L/btG)/2
This result suggests the view that a uniform stress
P/bt acts over a modified area bt/1.2, so that the same Us
results. Therefore the deflection at the free end of a
cantilever beam with parabolically distributed transverse
shear stress will be 1.2 times that with constantly
distributed transverse shear stress.
The Mindlin plate is generalized from the Mindlin beam.
Thus the higher transverse shear stiffness from the
assumption of constant transverse shear stress has been
reduced by dividing the corresponding elastic constants by
the factor 1.2 for flat plate element. The reduced stiffness
will produce the more flexible response in shear that is
expected from the actual parabolic distribution.
The integration over volume must be introduced for the
calculation of the element stiffness and equivalent nodal
loads of distributed loads or temperature effects. As the
thickness of the element is constant, the integration over
volume can be changed to the integration over area. This can
be accomplished through the partition of the strainnodal
displacement matrix [Bi] as follows:
0 0 zfj,x
0 zfi,y 0
[Bi] = 0 zfi,x zfi,y
fi,x 0 fi
fi,Y fi 0
The submatrices are named as follows:
BiA ZBiA
[Bi] =
BiB BiB
The element stiffness can be calculated through the
T
volume integration of Bi (3,5)E(5,5)Bi(5,3). Thus the [E]
matrix is to be partitioned as follows:
Ell E12 0 0 0
E21 E22 0 0 0
[E] = 0 0 E33 0 0
0 0 0 E44 0
0 0 0 0 E55
The submatrices are named as follows:
E EA 0
E =
0 EB
Then the stiffness of the element is
FT B T T [EA 0 ziA 1
[Ki] =I B E Bi dV = zBiA [ A B dV
0 EB BiB
2 T T
= [ BA EA BA] + [BB EB BB]
As [B] and [E] are constant about z, the integration
through the thickness from t/2 to t/2 can be performed on z
only and yields
= T 3
[Ki] = [ BA(t3/12)EABA + BB(t)EBBB] dA
I T 
= [ BAEABA + BBEBBB] dA
where,
EA = (t3/12)EA and EB = tEB
Then [Ki] can be rewritten as matrix equation as
follows:
T T EA 0 BiA
[Ki] = [ BiA BiB A dA
0 E B B iB d
T 
= BiE Bi dA
The size of plate element stiffness will be 27 by 27.
[K] =
[ E ]
(5,5)
[ B1 B2 B3 ... Bg Bg ] dV
(5,27)
B8
B9J
(27,5)
The strainnodal displacement matrix from which
the constant thickness is taken out is defined as [Bi].
[Bi] =
0
0
0
fi,x
fi,Y
0
 fi,y
 fi,x
0
fi
fi,x
0
fiy
fi
0
BiA]
BiB
Equivalent nodal loads due to body forces on the plate
element are calculated as
Pb = fTb dV = 1
V 1 i
fTblJI dr ds
in which {b} = { 0 0 bz )T or ( 0 by 0 }T or { bx 0 0
T in accordance with the direction of the gravity in the
coordinate system used. The nonzero quantities bx, by, or bz
represent the body force per unit area in the direction of
application.
Equivalent loads caused by initial strains are
PO = BTEeo dV
JV
1 1
= BTEO#IJI dr ds
where,
{40)= ( xx0 Oyy0 xy 0 0 )T
= { aAT/2 aAT/2 0 0 0 T
The stresses can be calculated from the equation
[a] = [E][E]
The corresponding generalized stresses, if desired, may
be computed from
M = ( Mxy My y QM Q y T = E ( B q 40 )
It is noted that the generalized stresses are actually
moment and shear forces applied per unit length of the edge
of the plate element. Therefore these can also be turned
into common stresses using the formulation for the bending
stress calculation. The moment of inertia for the unit
length of the plate is t3 / 12. Then the inplane stress at
a point along the thickness can be calculated as
a = Mz / I = M(t/2) / (t3/12) = 6M / t2
91
The transverse shear stresses can be found as
r =Q/ t
But this may be multiplied by a factor of 1.5 to get the
maximum shear stress at a point on a neutral surface because
the transverse shear stresses show parabolic distribution
while the calculated stresses are average stresses coming
from the assumption of a constant transverse strain along
the element z axis.
CHAPTER 5
NONLINEAR SHELL ELEMENT
5.1 Introduction
The nonlinearities included in the formulation of the
Mindlin flat shell element are those due to large
displacements and those due to initial stress
effects(geometric nonlinearity). The large displacement
effects are caused by finite transverse displacements. These
effects are taken into account by coupling transverse
displacement and membrane displacements. The initial stress
effects are caused by the stresses at the start of each
iteration. These stresses change the element stiffness for
the current iteration. These effects are evaluated directly
from the stresses at the start of each iteration and are
included in the element stiffness formation.
The total Lagrangian formulation is used. If the
updated Lagrangian formulation is used, the element
coordinate system cannot be easily formed for the next
iteration because the deformed shell is not usually planar
[26]. The symbol {} is used for a column matrix (a vector)
and the symbol [] is used for a matrix of multiple columns
and rows throughout the chapter.
5.2 Element Formulation
The generic displacements of Mindlin type shell element
are translational displacements {u v w)T and denoted as {U}.
The displacements and rotations at a point on the midplane
are (uo vo wo 9x 0y)T and denoted as (Uo}. The generic
displacements can be expressed in terms of the midplane
displacements and z as
u = Uo(x,y) + zOy(x,y)
v = Vo(x,y) zOx(x,y) (5.1)
w = Wo(x,y)
The linearized incremental strain from Eq. (2.17) is
e = ( ui,j + uj,i + uk,i uk,j + k,j k,i)
(5.2)
This equation can be written out for the strain terms to be
used for shell element using the generic displacements {u,
v, w)T
exx =
eyy=
exy =
u,
U,y
i (
2'
ez = i (
eyz = (
+ tu u, + tv, v, + tw,x ,
+ t,y U,y + t,y V,y + tw,y Wy
U,y + V,x + tU, U,y + t,x V,y + tx w,y
+ u, tuy + V,x tvy + W,x ty )
,z + w,x + tU u,z + tVx ,z + tW,x w,Z
+ u,x tu,z + V,x tv, + ,x tw,z )
V,z + w,y + tu,y u,z + t,y v,z + tWy w,z
+ u,y tU, + V,y tV, + w,y tW, )
(5.3)
The derivatives of inplane displacements u and v with
respect to x, y, and z are assumed to be small and thus the
second order terms of these quantities can be ignored
through von Karman assumption from Eqs. (2.44) [20, 21].
Furthermore the transverse displacement w is independent of
z for the shell element which means that w,z is zero. Then
Eqs. (5.3) can be reduced to
xx = u, + tw,x w,x
eyy = Uy + ty Wy
exy = ( Uy + Vx + tw,x Wy + ,x wy )
exz = ( u,z + wx )
eyz = 1 ( v,z + W,y )
(5.4)
The incremental Green's strains, sometimes called
engineering strains, can then be shown as
x = exx = Ux + tw,x w,x
4e = eyy = Uy + tw,y Wy
Ixy = 2exy = Uy + v,x + tw,x Wy + Wx ty
7xz = 2exz = u,z + wx
yz e = e = + Wy (5.5)
It is noted that the linearized nonlinear strains are left
only for inplane strain terms.
By substituting Eqs. (5.1) into Eqs. (5.5), the Green's
strain can be expressed in terms of midplane displacements.
