A STUDY OF THE APPLICATION OF MINDLIN
PLATE ELEMENTS TO THIN PLATES
by
Piyoros Jirawattana
B.Eng., Chiang Mai University, 1994
A thesis submitted to the
University of Colorado at Denver
in partial fulfillment
of the requirements for the degree of
Master of Science
Mechanical Engineering
1997
{777
This thesis for the Master of Science
degree by
Piyoros Jirawattana
has been approved
/
Samuel W
. Welch
Date
Jirawattana, Piyoros (M.S., Mechanical Engineering)
A Study of the Application of Mindlin Plate Elements to Thin Plates
Thesis directed by Assistant Professor Samuel W. J. Welch
ABSTRACT
This thesis is a study of Mindlin plate theory in finite element analysis with
particular emphasis on locking and instability behaviors of different element
formulations. A computer model is used to simulate the bending action of thin plates
and to study locking and instabilities associated with the different element
formulations.
We start with the bilinear Mindlin plate element. Using techniques of numerical
integration we demonstrate that shear locking can be eliminated and a very
competitive convergence results. This is done by using a technique known as selective
reduced intergation. However, the selectivereduced integration method can lead to
an unstable mode in the finite element solution.
Finally, we consider the development of a Mindlin plate element that does not
lock and does not contain unstable modes. This is done by approximating the
transverse shear strain using shape functions called mixed interpolation functions.
The detection of locking and instability will be presented by imposing spatial
Fourier modes and comparing the strain energy from finite elements to the strain
energy of the continuous Mindlin plate theory. This idea can also be used in the
analysis of convergence of the solution of finite element analysis.
This abstract accurately represents the content of the candidates thesis
I recommend its publication. /
Signed
[/ Samuel W. V W
in
ACKNOWLEDGMENTS
I would like to thank Samuel W. J. Welch for the thesis topic and his guidance. I
also wish to thank John A. Trapp for the fundamental plate theory from advanced
strength class.
IV
CONTENTS
Acknowledgments........................................................iv
CHAPTER
1. INTRODUCTION AND PLATE THEORETICAL PRELIMINARIES......................1
Scope................................................................1
Kirchhoff and Mindlin Plate Bending Theories.........................2
Introduction.......................................................2
StressStrain Relationship for Flat Plate............................2
Kirchhoff Theory...................................................3
MomentCurvature Relation.......................................6
Mindlin Theory.....................................................7
MomentCurvature Relation......................................10
2. FINITE ELEMENTS FOR PLATES...........................................11
Kirchhoff Plate Elements............................................11
Mindlin Plate Elements..............................................13
Stiffness Matrix.................................................14
Discussion of Kirchhoff and Mindlin Plate Elements.................16
3. INSTABILITY AND LOCKING..............................................18
Introduction........................................................18
Locking.............................................................19
Locking in Mindlin Beam..........................................19
Numerical Examples Using Mindlin Plate Elements...................23
Instability.........................................................27
Zero Strain Energy in Mindlin Beam...............................27
Eigenvalue Analysis of the Zero Strain Energy Modes..............30
Analysis of Instability and Locking of Bilinear Mindlin Plate Elements
Using Spatial Fourier Modes......................................31
Typical Instability Mode in Mindlin Plate Element................38
4. THE DEVELOPMENT OF A MINDLIN PLATE ELEMENT...........................41
Approach...........................................................41
Numerical Examples...............................................46
The Accuracy of the Mixed Interpolation Element..................48
5. CONCLUSION........................................................50
APPENDIX
A. Stress Resultants in a Flat Plate ...............................51
Force and Moment Balance of Flat Plates..........................55
Force Balance for Plates........................................55
Moment Balance for Plates.......................................56
B. Derivation of FourNode Quadrilateral Mindlin Plate Element......58
Numerical Integration............................................65
Two Dimensional Integrals......................................65
Stiffness Integration..........................................66
REFERENCES...........................................................68
VI
CHAPTER 1
INTRODUCTION AND PLATE THEORETICAL
PRELIMINARIES
Scope
Plate bending has been studied in the field of classical strength of materials for a
long time. The problem of plate bending can be considered as an extension of either
elementary beam theory or Mindlin beam theory depending on which theory is
appropriate. Both plates and beams support transverse loads by bending action. Two
theories of classical plate bending will be discussed, KirchhofF theory and Mindlin
theory. The difference between these plate bending theories can be classified by the
magnitude of the thickness compared to the magnitude of the other dimensions and
also the magnitude of the transverse displacement compared to the thickness. The
plate in this context will be considered as a relatively thin plate with small
deformation.
The classical KirchhofF plate bending theory makes an assumption that the
transverse shear deformation of thin plates in bending can be neglected. Information
about the resulting deformation can be assumed to be contained in the function
w(xi,x2) and in its derivatives. Therefore, finite elements based on this theory must
use C'continuous elements. Even though this classical theory is applicable to thin
plates, its use in finite elements is problematic as it requires higherorder nodal
interpolation function to enforce interelement compatibility. These higher order
interpolations can be quite complicated and impractical for general use.
In Mindlin plate bending theory, the transverse shear deformation is taken into
account so that this theory is applicable to thick plate bending. The deformation of
Mindlin plate theory consists of one transverse deformation and two small angle
rotations. The resulting finite element implementation requires only Ccontinuous
elements so that lower order interpolation functions may be used. In addition,
because the deformation in Mindlin plate theory is more physical, it is easier to apply
boundary conditions. However, the application of Mindlin plate elements to thin
plates may not be as accurate as KirchhofF plate elements because the transverse
shear strains from Mindlin theory, if approximated inaccurately, may make thin plates
1
too stiff. We will see that by the use of special techniques, Mindlin plate elements can
be made accurate for thin plate problems.
Thus, the context of this study is the application of Mindlin plate elements to the
bending of thin plates. We will begin with general background theory of Kirchhoff
plate bending and Mindlin plate bending in chapter 1. Next, the application of both
theories in finite elements will be discussed in chapter 2. Chapter 3 introduces the
concept of shear locking and how methods to prevent this problem result in the
side effect known as instability. Finally the development of a Mindlin plate element
for thin plates that does not lock and does not exhibit instability will be outlined in
chapter 4. The numerical results presented in this study were obtained using software
written by the author. The software is not included as part of this thesis.
Kirchhoff and Mindlin Plate Bending Theories
This section introduces, discusses and compares the Kirchhoff and Mindlin plate
bending theories. A more comprehensive discussion of these theories is contained in
reference [15] and in the references therein.
Introduction
A flat plate supports transverse loads by a bending action similar to a straight
beam but the deformation modes are two dimensional and thus are more
complicated. Plate theories assume that the plate thickness is relatively small
compared to the plate dimensions in the plane of the middle surface. Thus the
bending behavior of a plates in both Kirchhoff and Mindlin theory depends on the
plate thickness. The former theory makes the assumption that straight fibers normal
to the plane of the middle surface remain normal to the plane of the middle surface
after bending. The latter theory does not make this assumption. Before discussing the
details of these theories, we will start by discussing the general stress strain
relationships for flat plates.
StressStrain Relationship for Flat Plate
The stressstrain relationship may be obtained from Hookes law for isotropic
elasticity.
2
ai =777*1 *'
Ev
'^ij(Â£l\ + Â£22 + f3j)
\ + v'IJ (l + vXl2v) ~JJ' (11)
In both Kirchhoff and Mindlin theories a33 is considered negligible. From equation
(1.1) with 033 = 0, we obtain
v
Â£=~
\v
(Â£u +Â£22)
(1.2)
substituting this into Hookes law we obtain
vE
Gap =
1 + V
&afl(Â£ 11 +Â£22 ) + 2.GÂ£ap
(1.3)
equation (1.3) is the general stressstrain relation for a flat plate and from equation
(1.1) we get the equation for transverse shear stress
(1.4)
where a and /? range from 1 to 2 Sap= 1 if a = and Sap = 0 if a /?
Kirchhoff Theory
Figure 1.1 Kirchhoff plate after loading
3
In Kirchhoff theory, it is assumed that the straightline normals to the
undeformed middle plane of the plate remains straight and inextensional under the
deformation of the plate. In addition, it is assumed that the normals in the deformed
plate remain normal to the deformed middle plane which is a reasonable assumption
for thin plates. Under this assumption, it follows that Â£& = 0 which implies that =
0, hence the shear forces Q\ and Q2 are implied to be zero. Note that is neglected
in the stress strain laws but Qa is kept in the equations of motion. Despite this
contradiction, the theory is more accurate than the membrane theory of plates [6],
If we expand transverse displacement u3(xi,xzx3) about the mid plane using
power series in x3 we get [14]
u3(xi,x2,Xi) = w(xi,x2) + x3a(xux2) + x3 ......
if the plate is thin
!]
and we can approximate
u3(x i,x2,x3) = w(x,,x2)
(1.5)
from the assumption that = 0 we get
e
a 3
0 2
V<%3
\
y
and if the plate is thin we can introduce equation (1.5) into the equation above, then
0zd*a Mx{,x2)
<%3 <%a
integrating with respect to x3 we get
M*i.*2.*s) = *3^X^ + tf(*i*2) (16>
Ua (*1, x2) is constant of the x3 integration which has no effect on the moment Map
and thus may be neglected. The resulting kinematics assumptions for Kirchhoff
theory (neglecting
4
(1.7)
wi =
x
u2 =
d*
3 ~aT
dxx
3i
x
3 A
OX,
u3 = w
from f =
op j
&am+AA
\&fi &a
we get
f22 X3
OC,
*12 = "*3
(1.8)
(1.9)
(1.10)
(1.11)
(1.12)
Equations (1.7) and (1.8) are the kinematics assumptions of KirchhofF theory
mentioned earlier. Considering equations (1.7), (1.8) and (1.9), we can see that the
deformation in X\ and x2 direction are related to w by first order differential
equations. This leads to the secondorder differential equations for strain in equations
(1.10), (1 11) and (1.12).
To get the expression for bending moments in term of the transverse
displacement, we substitute equations (1.10), (1.11) and (1.12) into equation (1.3),
and then use the definition of bending moments (see the details of this definition in
appendix A)
h
2
MaP = J X&apfol
_h
2
to obtain
5
(1.13)
Mn=D
^ d2 w dlw'^
v
a, ac
2 /
m22=
+ V
\ 3c2 3c l j
(1.14)
Mn=D(\v)
d1
w
3cx3c2
(1.15)
where D is the flexural rigidity (define in the next section) which is approximately
equivalent to E times the moment of inertia of a unit width plate.
MomentCurvature Relation. In Kirchhoff theory we neglect Ojj and also
transverse shear strain, so that we can consider the in plane stressstrain relation. By
assumption the material is isotropic. Let xu x2 be the principle axis of material then
we can write
1v2
1 v
v 1
0 0
0
(1v)
22
(1.16)
We obtain the momentcurvature relation by substitution of equations (1.10),
(1.11) and (1.12) into equation (1.16) and the result into the definition of bending
moments, yields [12]
\M] D vD 0
m22 m]2 * vD 0 D 0 0 (1 v)D 2
(1.17)
6
where D =
Eh3
12(1 v2 )
E is the elastic modulus and
v is Poissons ratio.
Mindlin Theory
In KirchhofFtheory normals to the midsurface remain normal during deformation
and hence transverse shear deformation is neglected. Mindlin theory makes an
attempt to keep the transverse shear deformation which by assuming that fibers
originally normal to the midsurface can rotate independently of the normal direction.
Figure 1.2 Mindlin plate after loading
We define the transverse shear strain, y =6 and note that in KirchhofF
3c
theory this transverse shear strain is zero. In Mindlin plate theory we assume that the
deformation in the Xi and x2 directions is a function of small angle rotations so that
we can write,
7
W1 = *30,,(*,*2) (1.18)
2 =*3^(*.*2) (1.19)
ui = w(xx,x2) (1.20)
Equations (1.18) and (1.19) have included some rotation of fibers originally normal
to the midsurface to orientations not normal to the midsurface direction. To obtain
the expression for strains in Mindlin theory, introduction of equations (1.18), (1.19)
and (1.20) into the definition of =
1 f dij
di,
K&j
dc,
yields
i y
= "*3
(1.21)
and
' a 3
1 d4,  + _1 (
2 dc 3 J 2l
dv
dc
0.
(1.22)
from equations (1.21) and (1.22) we can write the expression for strains as
Â£u X30X) ,Xj
= ~XA,
, dda
where 6n,g = 
a p ee
p
Yu =^(^x,,xj+<5,x2.x1)
rn=w,Xl0Xl
Yu =^.x,^x,
(1.23)
8
Considering equation (1.22), cr^ is independent of x3 direction, so it is constant
over the plate thickness. In reality has a parabolic profile over plate thickness.
Instead of developing a very complicated theory to approximated in parabolic
manner with x3 (zero at h'lhQ. and maximum at x3 = 0). A k factor is added in the
transverse shear force constitutive equation. We can evaluate the correct k shear
correction factor by forming the strain energy for plate theory corresponding to a 3
dimensional energy with a parabolic shear stress. [14]
To get the expression of bending moment and the transverse shear force in
Mindlin theory in terms of the displacement function, substitute strain equations
(1.23) into the stressstrain relations (1.3) and (1.4). Then using the relationships
h
2
h
2
h
2
and
h
2
we obtain
(1.24)
(1.25)
It should be noted as transverse shear strain has been included, the transverse shear
stress can be written as
MomentCurvature Relation. Let *i and x2 be the principle axis of the material.
The MomentCurvature relation for Mindlin theory includes the shear stressstrain
relations on = Gy23 and cr3i = Gy3i The MomentCurvature relation is abbreviated
as {M)=[DKj\{K}. Again, we assume that the material is isotropic, from equations
(1.24) and (1.25) we can explicitly write
Mu 0 0 0 ' K ,jci
m22 [Dk] 0 0 0 ,Z1
Mn . =  0 0 0
02 0 0 0 Gh 0
.G . 0 0 0 0 Gh H 1 h
(1.27)
Where [Z\] is the same as equation (1.17). Compare equation (1.27) with equation
(1.17), as we mentioned ,Kirchhoff theory neglects transverse shear deformations so
that in Kirchhoff theory the transverse shear forces Q\= Q2 = 0.
10
CHAPTER 2
FINITE ELEMENTS FOR PLATES
The preceding chapter is devoted to the analytical derivation of Kirchhoff and
Mindlin plate bending theories. We might observe that numerical approximations of
Mindlin theory should be easier because the governing equations are of lower order.
We will discuss finite elements using these two theories.
Kirchhoff Plate Elements
As we have discussed in Chapter 1, Kirchhoff theory neglects transverse shear
deformation and is applicable to thin plates. Strain energy is determined by inplane
strains E\\, Â£22 and yn From equation (1.3), strains are determined by the lateral
displacement field w = h'(xiix2). Nodal d.o.f. of Kirchhoff elements consist of lateral
deflections w, and small rotations w,_ and w,_ of midsurface normals. We now
Jl ,xl
express the displacement field in terms of independent shape function interpolations as
[12]
* II M (2.1)
.x, = HN (2.2)
,2 = LA>,,x2 (2.3)
In order to formulate an element stiffness matrix [Ar]e, we now form the strain energy.
U = l{e}T[EUeidV (2.4)
These strains of Kirchhoff theory in equations (1.10), (1.11) and (1.12), may be
written in matrix notation as
ll
= [
*3w>Vl
2x3w,IiJ2
Figure 2.1 Plate 12 d.o.f. at node 3 of KirchhofF plate element
We consider that the thickness constant over the plate, so that we have dV = hdA and
finally dA = dxdy. Now using equation (1.16), equation (2.4) can be rewritten as
= (2.5)
A ^
where f/r}7 = In',,, 2tf, is known as the curvature for KirchhofF
K 9 [ XJXJ X2X2 9 xtx2 J
elements and \DK\ is known as the flexural rigidity
From equations (2.1), (2.2) and (2.3), we can write the nodal displacement matrix
of a KirchhofF element as
12
(d) = [w, w ,
(2.6)
.
WN W^N Y
Equation (2.6) can be written via the nodal interpolation shape functions as
w = [AG{J} (2.7)
equation (2.7) allows us to write the curvature using nodal shape functions as
{K) = [B]{d)
Where [B] is the differentiation matrix that contains information about the curvature.
Substituting equation (2.7) into equation (2.5), we obtain the element stiffness
matrix, [Â£]
where
u=^d)rikm
(2.8)
Mindlin Plate Elements
Mindlin plate theory includes transverse shear deformation, therefore it is a theory
that is reasonable for thick plates.
Nodal d.o.f. of Mindlin plate elements consist of one lateral deflection wi and two
rotations, 9xi and 0 of midsurface normals. If all interpolations use the same
polynomial, then for an element of N nodes, we can write
w V 0 0
8* II M 0 N, 0 **
KJ i=i 0 0 N._
(2.9)
where {
0x,n 0x,Jr
(2.10)
13
Consider equations (1.23) and (2.9), the field for w is coupled to the fields for 6
and 0Xj only via the transverse shear strains yn and yu [12],
A
X3
*2
Figure 2.2 12 d.o.f. at nodes 2 ofMindlin plate element.
Stiffness Matrix
In a like manner as we did for the KirchhofF plate element. We can write the
formula for the strain energy U [ 12]
h
U = ~f /{Â£}[Â£]{Â£}&,dA (2.11)
2
14
where A is the area of the plate midsurface and
{Â£} = [ex,x, V, rVi y xtxl ] Expressions for each strain are stated in
equation (1.23). Integration through the thickness yields
U = (2.12)
Â£ A
where [DM ] and {*} are stated in equation (1.27). We can write
{^)=[L]{W)
(2.13)
where
[L] =
0 d 0
0 0 0 0C2
n 0 0
u 0c2 0XX
d 0 1
d 1 0
Equations (2.9) and (2.13) yield
{K) = [B]{d)
where
15
[5] = [L][N\
0 AVx, 0
0 0
0 A\2 Ar..x,
0 A\
0
(2.14)
And finally from equations (2.12) and (2.14)
U = \{d)T{k]{d)
where [k} = \[B]T[DU][BYA
Discussion of Kirchhoff and Mindlin Plate Elements
Sometimes Mindlin theory is called thick plate theory and Kirchhoff theory is
called thin plate theory. Consider the d.o.f. in the displacement matrix, {
Kirchhoff and Mindlin plate elements. In Kirchhoff theory, the d.o.f. at node / are
one lateral displacement w and two rotations, w,z and w, .. In Mindlin theory, the
d.o.f. at node / are one lateral displacement wf and two rotations 9Zi, and 0XjI We
can see that the displacement matrix of Mindlin has lower order partial differentials
that make it more physical and easier to apply numerical approximations. In addition,
the higher order partial differentials in Kirchhoff theory requires at least biquadratic
shape function that lead to complicated numerical evaluation. Mindlin theory needs
only bilinear shape functions.
Considering equations (1.7), (1.8) and (1.9), Kirchhoff plate elements must display
interelement continuity of w, w,Zt and w (compatibility). Consider the connection
between elements in the x2 direction, continuity of w forces the continuity of the
tangential slope, w but not continuity of the boundarynormal slope w, Now for
Mindlin theory, consider equations (1.18), (1.19) and (1.20). Compatibility requires
the continuity of the fields, w, 9Z and 9Xj at the interelement boundaries. Mindlin
plate elements are easier to enforce compatibility because the lateral displacement
fields are not related to the rotation fields by differentiation.
16
In the next chapter, we will show that by using selectivereduced integration,
Mindlin plate elements can capture behavior of thin or Kirchhoff plates. These are
the reasons that Kirchhoff plate elements are not widely used.
17
CHAPTER 3
INSTABILITY AND LOCKING
Introduction
So far, we have discussed the derivation of Mindlin plate theory and its
implementation in finite element theory. In this chapter we compare results obtained
after selectivereduced integration and full integration are used for the integration of
the element shear stiffness matrix. This leads to the study of two important concepts
in finite element analysis that can make practical results unrealizable, instability and
locking. Numerous researchers have made attempts to prevent these problems.
P
Figure 3.1 Cantilever beam
18
Locking
Mindlin plate theory make plate bending theory more physical by including
transverse shear deformations into account. However in finite element analysis
approximation of these transverse shear strains can lead to inaccurate results.
Physically the deformation of thin plates should be dominated by the bending
deformation but if it is dominated by transverse shear deformation we call this shear
locking.
Locking in Mindlin Beam
To explain how this shear locking occurs we consider the cantilever beam of
length L with load P in Figure 3.1 and we use two different integration methods, full
integration and selectivereduced integration.
We employ Mindlin plate theory with the nodal d.o.f. wi, ft w2, ft We separate
the stiffness matrix [Â£] into a component with bending stiffness [Â£fc] and a
component with transverse shear stiffness [Â£J. We do so in order to integrate the
bending and the shear terms separately.
We obtain [ft] that is associated with inplane strain
0 ML 0
0 0 0
ML
0
and [ft] that is associated with transverse shear strain
0 0 0 0
ML x/Ll ML x/L
(3.1)
(3.2)
case 1) Full integration
Introduction of equations (3.1) and (3.2) into the expression for the stiffness
matrix and using exact integration for both the bending stiffness [&6] and the
transverse shear stiffness [ks ] we obtain the stiffness matrices
19
0 0 0 0
0 101
0 0 0 0
0101
(3.3)
[*,] =
1 L 1 L
2 2
L L r l2\
kGA 2 l 3 J 2 l 6 ,
L 1 L 1 L
2 2
L L iL2)
2 l 6 J 2 { 3 J
(3.4)
With k =
, we solve for the tip displacement
*2
for a thick beam, h L, the above solution becomes
6 PL
w. =
2 5GA
(3.5)
(3.6)
Equation (3.6) is a reasonable result for a thick beam with deformation caused by
shear. Compare equation (3.6) with the deformation caused by shear from
PL
Castiglianos theory, w2 = 1.33 [10], Thus, the Mindlin beam provides a
AG
reasonably accurate solution for thick beams.
20
Now for a thin beam, h L, from equation (3.5) we obtain
24 PL
5GA
(3.7)
In the latter case the deformation should be dominated by bending (compare to the
pi}
analytical solution from elementary beam theory w2 but it is not. It is
dominated by shear. This is an example of shear locking. Now we consider what
happens if selectivereduced integration is applied.
Case 2) Selective integration
In this case, we integrate bending stiffness with two Gauss points and use one
Gauss point for the shear stiffness. This yields the same bending stiffness as before but
we get a different shear stiffness.
1 L 1 L
2 2
L (Ln L iL2)
kGA 2 l 4 J 2 l 4 J
L 1 L 1 L
2 2
L iL2) L
2 l 4 J 2 L 4 J
(3.8)
Solving for the deflection at the free end, we get
w2
P_ _L_ 11)
G khb 2bh3
(3.9)
The first term and the second are deformation from shear and bending, respectively.
When the beam becomes thin, the second term contributes and we can neglected the
first term. The solution from elementary beam theory ( using Poissons ratio of 0.3 )
21
is w, =. Selectivereduced integration with a Mindlin beam finite element
2 1.3 Gbh3
produces a result of 97.5 % accuracy. It should be noted that for a thick beam
solution, either full integration or selectivereduced integration can capture the
behavior of shear deformation because the bending deformation can be neglected.
It is counter to our intuition that an exact integration produces an inaccurate
result but an underintegration might introduce an acceptable answer. In order to
explain, we consider the cantilever beam in Figure 3.1, obtained by simply reducing a
Mindlin plate element into one dimension. We can write the equation for axial strain
eu and transverse shear strain y 31 as
Â£ii=x30,x, (3.10)
Yi\ =w>x,^x, (311)
where the notation is the same as that of the Mindlin plate element.
We separate the element strain into the bending strain and the transverse shear
strain. We can write
{[*,]+[*.]}<
Our discussion will be concerned with the transverse shear approximation in equation
(3.12).
= w2w] x(6>20,) e
L L 1
(3.13)
If we consider equation (3.13) we see that the terms associated with the transverse
deformation are constant. The terms associated with the rotational degree of freedom
are linear in x. If we use selectivereduced integration we will evaluate (3.13) at x =
L/2, we are in effect replacing (3.13) with
Yi\ =
_ *2
w, 6, + 0,
(3.14)
22
In the literature it is common to see the above process referred to as eliminating the
excess monomial. This excess monomial results from using inappropriate order
polynomials for the displacement and rotation fields.
When the beam becomes thin, the transverse shear terms should become
infinitesimal. The shear resulting from evaluating equation (3.13) at two Gauss points
cannot, in general, be infinitesimal at both points and the resulting error becomes
magnified in thin beams as the stiffness matrix divides this error by the thickness
squared. Thus, in the case of thin beams, if full integration is used with the shear
stiffness matrix, the deformation will be dominated by [Â£J not by [^]. This is called
Locking. In the literature the use of full integration on the shear stiffness terms is
sometimes referred to as over constraining the element to have zero shear at two
points thus causing locking.
Numerical Examples Using Mindlin Plate Elements
We now demonstrate the concepts we have been discussing with numerical
simulations using various Mindlin plate elements. We will numerically demonstrate
that full integration leads to locking that can be suppressed by selectivereduced
integration. This section will present numerical examples of Mindlin plate element
used on a square plate with various boundary conditions and loads. We model the
square plate with 10x10 elements as shown in Figure 3.2
23
Figure 3.2 Square plate
The data for the plate consists of the following
E = 3 x 107 psi
v= 0.3
a = b = 1 inch.
Considerations similar to those in Mindlin beam elements can be applied to the
Mindlin plate element. If the plate is thin, the transverse shear strain in both directions
(equations (3.15)) should be small.
7 2 3 =0
r*\ =v%0x, = 0
(3.15)
Using bilinear shape function for Mindlin plate elements, twopoint Guassian
quadrature is required to exactly integrate the shear stiffness. This full integration
brings in the higher order terms of the rotational degrees of freedom that the terms
24
corresponding to the transverse deformation do not have. Thus, full integration can
lead to shear locking unless both transverse shear strains are evaluated at
x, = x2 = 0 (selectivereduced integration). The point being that the excess
monomials in the terms corresponding to the rotational degrees of freedom are
eliminated. Using selective integration with Mindlin plate elements can prevent
locking and gives acceptable results for thin plates as shown in Figure 3.3. We can see
that when the plate become thin, exact integration gives unacceptable results due to
locking.
a) Clampededges with distributed load.
25
fiumcrtcaUanalytlcBl numtricalfanalytlcal
1.2
b) Simple supports with distributed load
c) Simple supports with point load at the center of plate
Figure 3.3 Numerical/Analytical Ratio of center deflection of 10x10 element, square
plate. Thin plates correspond to large ab/t.
26
Instability
Generally, finite elements do not give exact results. We have demonstrated that
the use of a more accurate stiffness integration might make an element too stiff
causing locking. We saw that using selectivereduced integration eliminates locking.
However, using selectivereduced integration can also cause something worse than
locking called instability.
Zero Strain Energy in Mindlin beam
Instability is sometimes referred to as a zero strain energy mode. This is because
the global displacement matrix results in zero strain at the Gauss points so that the
mode is easily excited. For this reason, instability is often associated with rank
deficiencies of the element stiffness matrix which is caused by lowerorder integration.
To explain the term zeroenergy mode, we form an expression for the element
strain energy.
u.=UdYmw=\ iisViDmdA
2 2 A.
where \k]= \[B]T[D][ByA
A.
and {e) = [B]{d}
When we evaluate [it] numerically, it contains only the information that can be
sensed at the sampling points of the quadrature rule. If the strains vanish at all Gauss
points for a certain mode {d}, then Ut must be equal to zero, and [Â£] must be a
zerostiffness matrix. We know that the strain energy must vanish for only rigid body
motion. If it vanishes for other motions, an instability has happened.
To support the above discussion, we consider the inplane twisting mode of the
two Mindlin beam elements shown in Figure 3.4. We will compare the strain energy
using stiffness matrices derived with both full integration and selectivereduced
integration.
27
^=0, 0i=l
^3=0, 01=1
o e o
^2=0, 02= 1
Figure 3.4 Instability in the Mindlin beam element, dashedline shows inplane
twisting instability mode.
The global displacement matrix of the 3 nodes is
{
w,
e2
e
3
The displacement matrix of the inplane unstable twisting mode is
o'
1
0
1
0
1
(3.16)
(3.17)
28
The global stiffness matrix, assembled from the element shear stiffness in equation
(3.8) (obtained by selective reduced integration), is
, kGA
g ~ L
1
Z
2
1
Z
2
0
Z
2
Z2
4
Z
2
j}_
4
0
1
Z
2
0
1
0
1
0 
Z
2
Z2
4
1
_Z
2
L
2
0
L
2
Z?
4
L
2
4
(3.18)
where k*g is global shear stiffness matrix obtained by selectivereduced integration.
Now by applying equations (3.17) and (3.18) into the expression of strain energy, we
obtain
= 0
where UT denotes the total strain energy.
This result corresponds to the definition of zerostrain energy mode. In a like manner
we use full integration and obtain the exactly integrated shear stiffness matrix, kÂ£
29
1
1
0
0
f kGA
g L
1
L L2 L L2 0 n
2 3 2 6
L L
1 2 0 1
2 2
L L2 n 2L2 L L2
2 6 3 2 6
1 L L
0 0 1
2 2
n fi L L2 L L2
u 2 6 2 3
(3.19)
Now calculating the total strain energy, we obtain
UT=\{d)TlkI,\{d)
= 21}
3
We can see that the zerostrain energy mode is not present with exact integration.
However we should be concerned about exact integration because of shear locking.
Eigenvalue Analysis of the Zero Strain Energy Modes
There is another approach to detect instability in the finite element method. This
approach performs an eigenvalue analysis of the element stiffness matrix. Consider the
element equilibrium equation [12]
md) = {/)
Now if the load matrix {/} is proportional to element nodal displacements {d}
through a factor A, we have
md) = {f)=m
30
or
mmw = {0}
(3.20)
This is an eigen problem with the eigen values of {Â£} being the X. Each X, corresponds
an eigenvector {
(3.20), we obtain
{<#[*]{,=*,
or 2Ut = X, (3.21)
where Â£/, is the strain energy in the element when its nodal d.o.f. are normalized
displacements {d},.
[Â£] is complete element stiffness matrix of the element in an infinite media.
Equation (3.21) can be interpreted in the same way as was the expression of strain
energy, that is [it] should yield Xt = 0 for the rigid body motion case only. If it is not
a rigid body motion and A, = 0 then instability is present. Generally a solution with
an instability is easy to detect as it usually contains short wavelength, high amplitude
unphysical deformations.
Analysis of Instability and Locking of Bilinear Mindlin Plate Elements
Using Spatial Fourier Modes
Corresponding to the above discussion we might say that instability happens when
strain energy vanishes for deformation modes not related to rigid body motion. Next
we will examine the instability and locking by applying Fourier modes on an infinite
media. We begin by writing the expression for strain energy for a linearly elastic
isotropic Mindlin plate
31
u(w,ex ,ex)= Eh3. ff
x' Xl 24(1 v2):
ee, \
dx
+ 2v
i
d0x ddx
X1______^2_
dxx dx2
+
h 1 (aex eex > xi xj 1 fN
l ) 2 h
dx
Eh
4(1 + v)
A, '
+
3>v
\ck2
e,
dx^dx2
(3.22)
We will use the definition above to calculate a continuous strain energy to compare
with the strain energy calculated from the finite element methods. We write the nodal
deflections w 0 and 0Xjl as general harmonic waves in 2 dimensions
w, = cos(*Xi x,) cos(*X2 x2) (3.23)
For thin plates we can neglect the transverse shear deformations, so that
0*,i = "*x, sin(*Xix,)cos(*Xjx2) (3.24)
9Xj = kH cos(^xl)sin(^x2) (3.25)
We note that in Mindlin theory, other modes may be present but we expect these to
be the dominant modes. These three equations above can be written via nodal shape
functions as
wt = cos(kxdx^n)cos(kxdx2n) (3.26)
0*,, = sin(Arx,dxxri)cos(Â£Xjdx2m) (3.27)
= ~kXi cos(kx dx]n)s\n(kx dx2m) (3.28)
where w is the node element deflection
6Xi is the small angle rotation in x\ direction
0Xj is the small angle rotation in x2 direction
32
kT kT is the wave number in each direction
dx\ is the horizontal length of element
dx2 is the vertical length of element
n is the element index in horizontal direction
m is the element index in vertical direction
Equations (3.26), (3.27) and (3.28) allow us to write the displacements as a function
of the wave number and then calculate the strain energy. We will compare the
continuous and finite element strain energies calculated using these modes.
case 1) kx = k. =0.
/ x, x2
This case corresponds to rigid body motion. It is certain that strain energy must
vanish. This is true for any method.
case 2) maximum number of wave, kXi x, = n, k x2 = n
The maximum admissible wave number that can be represented by nodes of finite
elements for the transverse deformation is presented in Figure 3.5. This is the
characteristic instability called whourglass mode.
wi = l, #i=0
h3 = 1, 0,=O
Figure 3.5 The maximum admissible wave number.
33
The wave number corresponding to the wave length is
k
2 n
~T
(3.29)
where k is the length of one wave. At the minimum admissible wave length k, which
is equal to the length of two elements, 2dx, we have the maximum wave number
k
max
n
dx
(3.30)
The deformation at the maximum admissible wave number is not a rigid body motion.
If we apply the maximum wave number into equations (3.26), (3.27) and (3.28) and
calculate the strain energy we will obtain, in some cases, the zero strain energy mode.
We will show that the zero strain energy mode only appears, in some cases, at the
maximum admissible wave number. We will do this calculation numerically with a
single one inch square element. In the case of Mindlin plate elements, only the w
hourglass mode is possible. We will discuss this in more detail later. We can write the
displacement
w, * cos(Â£Xi dxji) cos( kxdx2m)
{d}= > = k s\n(k dxxn)cos{k dx2rri)
K\ k cos(kxdx{n)siT\(kXidx2Tn)
Using equation (3.31) in the expression for strain energy, we can calculate the strain
energy as a function of the wave number using stiffness matrices obtained using full
integration and selectivereduced integration. We compare these solutions with the
continuous solution using the same Fourier modes. Figure 3.5 presents the strain
energy of a single one inch square Mindlin plate element. We set the shortest wave
length equal to 2 inches so that (Â£x *,), = (kX2X2)mm n. Consider Figure 3.5 a
and c, The locking behavior makes the strain energy obtained from full integration too
big, even at a very small wave number. In Figure 3.5 a and c we plot the logarithm of
strain energy as the strain energy obtained from continuous case is very small
34
compared to that obtained using full integration. Compare these results with the
results from selectivereduced integration in Figure 3.5 b and d, we can see that the
smaller the wave number, the more accurate these elements are. Practical use of these
elements does not generally consider deformations near the maximum admissible
wave number and the accuracy of these elements would be acceptable were it not for
the behavior at the maximum admissible wave number. At the maximum admissible
wave number the strain energy using selectivereduced integration vanishes indicating
the presence of an unstable mode.
35
b) Strain energy from selectivereduced integration and analytical from kAx = 0
to n
36
d) Strain energy from selectivereduced integration and analytical from kAx = 0
to ^lO
Figure 3.5 The strain energy of a single one inch square Mindlin plate element in
an infinite media, thin plate limit.
We note that Figure 3.5 indicates the expected result that for a given wavelength,
2n,
the numerical results become more accurate as Ax > 0.
37
Typical Instability Mode in Mindlin Plate Element
We split [A:] into the bending stiffness [Â£*], and the transverse shear stiffness [Â£,].
We integrate the bending stiffness fully and we under intergrate the shear stiffness
which makes only two modes of instability possible, whourglass and inplane twist
mode. The latter can not happen with two or more element because it can not be
communicated between adjacent elements, so that only the whour glass mode
appears. In general the global stiffness matrix is singular by itself and zeroenergy
modes can be prevented after boundary conditions are applied. From the assumption
that the boundary conditions preclude the rigid body mode from appearing in one
element, it also precludes the mode in the remainder of the mesh [8], A critical test is
to clamp one node in the middle of the square plate, and the instability appears as
shown in Figure 3.6. Now we fix one node that is adjacent to the center node. The
instability is gone. To retain symmetries as shown in Figure 3.7, we fix 3 adjacent
clamped nodes in the middle of the plate.
38
x 10
10
Figure 3.6 whour glass mode with clamped center node of 10x10 Mindlin
plate elements with selective integration and distributed load.
39
15
Figure 3.7 Deflection of 3 clamped center nodes of 10x10 Mindlin plate
elements with selective integration and distributed load.
40
CHAPTER 4
THE DEVELOPMENT OF A MINDLIN PLATE ELEMENT
The preceding chapter was devoted to the study of instability and locking in
Mindlin plate elements with bilinear shape functions. We showed that using exact
integration of the element stiffness matrix resulted in shear locking. Shear locking can
be suppressed by adopting selectivereduced integration, however for some specific
situations selectivereduced integration introduces instability. In this chapter we will
present an reliable element that does not lock and contains no unstable modes[4].
Approach
The development of Mindlin plate element obtained from the method of selective
reduced integration will be tried with a different point of view. Let us consider the
curvature of Mindlin plate element.
M =
*1
0*. .*,+0,,.
0r
(4.1)
The first three terms in the curvature matrix are bending stressstrain relations, cth,
022 and 012, respectively. The last two terms on the bottom are the transverse shear
deformation. For the isoparametric element we used the standard nodal shape
interpolations as follow,
N,= UlWi7)
4
(4.2a)
41
^2= la+MJv) 4 (4.2b)
M= ) 4 (4.2c)
N, I(/K/+7) 4 (4.2d)
Where each node is located in Figure 4.2
0,1)
(1,1)
Figure 4.1 The quadrilateral element in Tj coordinate (the master element)
In the previous chapter we showed that selectivereduced integration solves the shear
locking problem and gives good convergence behavior, but leads to unstable modes.
To motivate the development of a better Mindlin plate element we consider the
positive characteristics of selectivereduced integration. The one Gauss point located
at the center of an element results in constant transverse shear strains over the
element in each direction. The element presented in this chapter has similar properties
in that one of the shear strains is constant along element edges. The shear strain is
obtained by approximating the slope of the transverse deformation with out strict
adherence to the element shape functions. The resulting expression allows us to use
two Gauss points and integrate the resulting polynomials exactly.
42
From the definition of transverse shear strain
Y 23=^^ (43)
Y*\=w,xl0xl (44)
We can define these transverse shear strains with different interpolation functions
from the in plane stressstrain terms [13]
r23 =^a + ^)(rf3)+^aa(r?3) (45)
rn =^(1+'7)(rfi)+o'7)(r3i) (46)
where yn and yi\ are defined in Figure 4.2, and y\^ y^, y\x and yc%x are shear
strains at points B, D, A, C respectively
4 L / *1 3
i 3 1 n c
i \ i 4 t D 4 B
4
1 2
PI 1 *
Â£
T
h
jk_
Figure 4.2 Convention for the transverse shear deformation
43
Approximating the transverse shear strain from equations (4.3) and (4.4) and using
equations (4.5) and (4.6) we obtain
(
w2w
Unlike the shear strains resulting from strict observance to the element shape
functions these shear strains may have infinitesimally small values at two points and
thus will not lock when full integration is used. We are now ready to modify the
stiffness matrix by using these difference interpolation functions with linear shape
functions, the curvature can be modified to be
Where ^ and yu are defined in equations (4.7) and (4.8). The [5,] matrix can be
straight forwardly constructed as
(4.9)
Yl 3
Yi i
44
(10 07)
A A
0 07) 2
00 0
2
0 + 0 (17)
A A
0 07) 2
00 0
2
0 + 0 (1 + 7)
A A
0 0 + 7) 2
0 + 0 0
2
00 0 + 7)
A A
0 1? + 1
00 0
2
(4.10)
Then from the shear stiffness matrix is calculated
nDu\{B,w
A
Using this new mixed interpolation function, Mindlin plate elements do not lock
and because every integration is done exactly, it is certain that there are no unstable
modes. It should be noted that these mixed interpolation functions may be under
integrated as well as with the bilinear interpolation function but we must beware of
the possible unstable modes. For example, if this new mixed interpolation function is
under integrated and applied to the clampedcenter node problem shown at the end if
chapter 3 also leads to unstable modes. However, we can use full integration to solve
45
this problem without locking or instability and this will be demonstrated by using
Fourier modes to calculate the strain energy.
Numerical Examples
We used the same material and dimensions presented before but with the mixed
interpolation with full integration discussed above. Note that the locking behavior
exhibited with full integration in the previous chapter does not appear here.
Simple supportdistributed load
a) simple supportdistribute load
46
Simple Supportpoint load
b) simple supportpoint load
Clamped edges distributed load
c) clampdistribute load
Figure 4.3 Numerical/Analytical Ratio of center deflection of 10x10 element,
square plate. Thin plates correspond to large ab/t.
47
The Accuracy of the Mixed Interpolation Element
The preceding chapter, we used strain energy calculations to analyze the accuracy
of selectivereduced integration. We will use the same to analyze the mixed
interpolation element discussed in this chapter. The purpose of the element introduced
in this chapter is to prevent both instability and locking by adopting a different nodal
interpolation function with full integration of both the bending stiffness matrix and
the shear stiffness matrix. Considering Figures 4.4 a and b we can see that for small
wave numbers the element produces acceptable results. There is no shear locking and
no zerostrain energy mode.
a) Strain energy from mixed interpolation and analytical from kAx = 0
to 7t
48
b) Strain energy from mixed interpolation and analytical from kAx = 0 to itf\0
Figure 4.4 The strain energy of a single one inch square Mindlin plate element
in an infinite media, thin plate limit.
49
CHAPTER 5
CONCLUSION
A computer model was developed to simulate the bending action of Mindlin plate
elements. This program was also dedicated to the study of shear locking, instability
and strain energy in thin plate deformations. The selectivereduced integration with
standard bilinear shape functions is introduced to eliminate shear locking successfully.
This technique requires simple programming and produces a very good solution for
thin plates (more than 99% accurate compared to classical analytical solutions).
Unfortunately, this selectivereduced integration leads to unstable modes in the finite
element solution. These unstable modes can be precluded via an appropriate boundary
condition.
Finally, by considering the plate as an extension of twonode beam elements to
twodimensions, the different interpolation functions for transverse shear deformation
can be constructed. The mixed interpolation produced even better convergence
results than the selectivereduced integration method did. This mixed interpolation
function does not lock and also does not contain any unstable modes. It is obtained by
a slight modification of the shear stiffness matrix for an element in cartesian
coordinates. However, this element is quite complicated in general geometries
element compared to the previous techniques. This is because it involves tensorial
shear strains in arbitrary directions.
We presented two different viable element formulations, selectivereduced
integration and mixed interpolation. Each has practical advantages and disadvantages.
The instability from selectivereduced integration results from use of boundary
conditions not usually seen in practical problems. The mixed interpolation formulation
for arbitrary quadrilateral shapes can be cumbersome. We let the reader decide which
technique is appropriate to their problems with optimization of computer work and
the workers time.
50
APPENDIX A
Stress Resultants in a Flat Plate
Consider Figure A.l by assuming that the material is isotropic. Let Mu be the
bending moment per unit length of x2 axis. We obtain [6]
h
2
Mn = jx3audx 3
_h
2
and let M\2 be the twisting moment per unit length of x2 axis, we get
h
2
M\2 = *30'l2^1C3
h
~ 2
Similar development in x\ axis, we obtain
h
2
M22 = {*3022^3
h
~ 2
h
2
and M2l = J x3 ondx2 = Mn
_h
2
note: positive direction is indicated by the righthand rule for moments.
If we let Qi and Q2 be the transverse shears per unit length in x2 axis and Xi
direction. We can find equations for transverse shear forces per unit length.
(A. 3)
(A.4)
(A.2)
(A.l)
51
h
and
2
Ql= S CT13^3
_h
2
h
2
Q2 = f a23dx3
_h
2
(A.5)
(A.6)
Resultant moments and transverse shear forces are shown in Figure A.2
Next, is the derivation of the equation for the traction force on a cross sectional
surface of the element. Let Mi be the traction force of the element per unit length in
the *2 direction on the middle surface. From Figure A. 1 we get
A
Figure A.l Stresses that act on cross sections of a differential element of a
homogeneous linearly elastic plate. q(xi,x2) is the distributed force per unit area on
lateral surface
52
h
2
h
2
h
2
N
12
h
2
Traction forces are shown in Figure A.3
A
<<
> *2
Q2
O>
M\2
f m22
Xi
Q10
A Ml
Figure A.2 Resultant moments and shears
(A. 7)
(A. 8)
(A. 9)
53
Nl2
Nu
N22
Xi
f
X\
*
Nn
>
12
N2i
N22
x2
Figure A.3 Resultant tractions
Normal stresses on and 022 vary linearly with x3, are maximum at h/2 and zero
at x3 = 0, and are associated with bending moments Mu and M22. The shear stress
also varies linearly with x3 and is associated with twisting moment Mn. Transverse
shear stresses oi3 and 023 vary quadratically with x3. Integration of equations (A.l)
(A.6) results in
a
11
hi3 /12
(A. 10)
a22
^22X3
h3 /12
(A.l 1)
54
(A. 12)
al2 =
Mnx3
h3 /12
The maximum magnitude of 023 and <7,3 is at x3 = 0
15 Q2
<7 23 
CT]3 =
t
150,
(A. 13)
(A. 14)
Force and Moment Balance of Flat Plates
We will start the derivation of the force and moment balance equations for a flat
plate by considering the equations of motion
da 11 da 1 da1
" 12+= pUx
dx, dx2 dx
da2x , da22 da.
dx\ dx2 dx.
da 3\ ( da 32 da H =
dx, dx2 dx
= pu2
= pu3
(A. 15)
(A. 16)
(A. 17)
Force Balance for Plates
by using equation (A. 15) divided by thickness h and integrating with respect to x3
[14]
h
_1_
h
h
1 f {3(7 w  da 12
h \ \ dxx dx2
~2
da
dx
= f\
dx.
we obtain
55
(A. 18)
A parallel development using equation (A. 16) results in
(A. 19)
Equations (A. 18) and (A. 19) are force balance equations for the plate element per
unit thickness as shown in Figure A. 3 X\ and X2 are the force per unit area of
loading on top and bottom of the plate in the X\ and x2 direction, respectively.
Moment Balance for Plates
We now consider
Equation (A.21) is the vertical force balance for a plate element, as shown in Figure
h
2
integrating equation (A. 20) yields
(A.21)
A. 2
h
2
We now form Jxi (A.l5)dx3 and obtain
h
2
(A. 22)
56
h
Similarly we
2
form ^x3{A\6)dx3
to obtain
dc\
cM22
ac2
 + C2Q2 =phx3u3
(A.23)
Equations (A.22) and (A.23) are the balance of moments about the xi and x2 axes,
respectively. The physical picture is shown in figure A.2.
57
APPENDIX B
Derivation of the FourNode Quadrilateral Mindlin Plate Element
,1)
0,1)
Figure B.l The quadrilateral element in Â£ 77 coordinates (the master element)
We first develop the standard, bilinear shape function on the master element
(Figure B.l) which is defined in 77, % coordinates. The shape functions are defined
such that Ni is equal to unity at node / and equal to zero at other nodes. The four
shape functions can be written as
i(l3(l/7) (B.la)
4
Ni 7O+0OIJ) (B ib)
4
7(1+00+17) (B.lc)
4
58
(Bid)
JV4= 7(10(1+17)
4
or M=i(l,)(l+77,) (B.le)
4
where and 77, are the coordinates of node /
We now express the displacement field within the element in terms of the nodal
values.
4
w= ^Ntwt 1=1 (B.2a)
1=1 (B.2b)
(B.2c)
1=1
which can be written in matrix from as
{"HAW (B.3)
where
59
M
W1
L*VJ
(B4)
and
0 0 0 0 0 0 n4 0 0
N = 0 0 0 N2 0 0 *3 0 0 N4 0
0 0 0 0 Ni 0 0 N, 0 0 x4
0.5)
Next, we express the derivatives of a function of x, in Â£ q coordinates, using the
chain rule,
4L=Â£.*l+4L*2.
dE, dcx dE, dc2 dE,
Â§L = Â£.*l + Â£.*i
dq dcx dq ck2 dq
(B.6)
or
rr
rJ >=j< * f
Â§_
A.
(B.7)
60
where J is the Jacobian matrix
From equations (B. 1) and (B.2) ,we obtain
j_ J\\ Jn
_^21 *^22 _
where
Jn = 7 [(l7)(^i)i + (17)(^i)2 + (1+7)(*i)3 (l+7X*i)4]
4
j\2 = 7 [(l7)(^2)l + (1 n)(x2)2 + (\ + Tj)(x2)i (1+ 7X^)4]
4
J21 = \ KlfX*,), (l+fX*,)2 + (l + ^)(X,)3 + (lfl(*l)4]
Ja \[OflWi (1+5X*j)2 + (1+Sfe), + (lafe).]
4
from equation (1.27), the curvature is defined as
M
^x, >x,
9** ,x:
e*x +0*
0I2
0. w.
by applying equation (B.7) we obtain, after inversion
1 ^22 Jn
detJ ~^2\ J\\ . .^7.
(B.8)
fay
ac, 1 J22 J\2
detJ 1 1 h t
.&2 . . #n .
(B.9)
<*> 1 J22 ~Jn
. ^2 . dety L^2, Jn J ay .
(BIO)
62
from equations (B.8), (B.9) and (B.10), we get {/c} in the form
M =
1
dety
0 0 Jn ~Jn 0 0 0 0
0 0 0 0 Ax Ax 0 0
0 0 Ax Ax J* ~Jn 0 0
Ax ~Ax 0 0 0 0 0 dety
1 > K> Jn 0 0 0 0 dety 0
d#
an
an
%
an
(B. 11)
Now from the interpolation equation, we have
3>v
av
an
aexx
%
aex,
an
aa>c2
%
aox2
an
flr,
(k2
= [G]{d)
(B. 12)
63
where
(17) (15) 0 0 0 0 0 0
0 0 (17) d5) 0 0 d5)d7) 0
0 0 0 0 d7) d5) 0 (15X17)
(17) (1+5) 0 0 0 0 0 0
0 0 d7) d+5) 0 0 /"v C* 1 W + 0
0 0 0 0 (17) d+5) 0 (1 + 5X17)
0 + 7) 0+5) 0 0 0 0 0 0
0 0 d + 7) d+5) 0 0 d + 5)d+ 7) 0
0 0 0 0 (1+7) d+5) 0 (1+5X1+7)
(1 + 7) d5) 0 0 0 0 0 0
0 0 d + 7) d5) 0 0 (15X1+ 7) 0
0 0 0 0 d + 7) d5) 0 (15X1+ 7)
equations (B.l 1) and (B.12) yields
[B] = [A][G]
where
[A] =
1
detJ
0 0 Jr. ~Jn 0 0 0 0
0 0 0 0 ~Jr J\\ 0 0
0 0 Jr Ju Jj2 ~Jn 0 0
Jr 0 0 0 0 0 dety
Jn Jn 0 0 0 0 detJ 0
Finally we can calculate the element stiffness ffom
A
64
Numerical Integration
Consider the problem of numerically evaluating a onedimensional integral of the
form [11]
l
The Gaussian quadrature approach for evaluating / with /7point approximation is
given below
1
/= ) + w2/(f2)+ +"/(
1
where w2,..........,wn are the weights
and 4\, g2,........,Â£n are Gauss points
From equation (B.13), //point Guassian quadrature provides an exact answer for
polynomial of order (2nl) or less.
Two Dimensional Integrals
We can write the equation to evaluate the twodimensional integrals of the form
i i
f=[f/({.vWv (Bi4)
11
Now, from equation (B.13)
1U=1
i=l )=l
(B 15)
65
Stiffness Integration
Introduction of [5] into the expression for the stiffness matrix, we may now
numerically evaluate the element stiffness matrix [11],
1
A
0.57735
0.57735
4
Â£, = 0.57735 Â£2 = 0.57735
Figure B.2 twodimensional Gaussian quadrature with 2x2 rule.
[DMteUd&i! (B.16)
11
where [5] and det/ are the functions of Â£ and tj.
Let ftt,rj) = [B]T[D][B]detJ (B.17)
Then, If we use a 2point Gaussian rule in each direction, we get
66
[k .] = ,TJi)+^^2
where W] = w2 = 1.0
i,=*~ Yji
b = m= Xf}
and the Gauss points are shown above diagrammatically
67
REFERENCES
1. T. J. L. Hughes, R. L. Taylor, and Kanoknukulchai, A Simple and Efficient
Element for Plate Bending , International Journal for Numerical Methods in
Engineering, 11, no. 10(1977), 15291543.
2. R. H. Macneal, A Simple Quadrilateral Shell Element, Computers and
Structures, 8(1978), 175183.
3. T. J. R. Hughes, Generalization of Selective Integration Procedures to
Anisotropic and Nonlinear Media, International Journal of Numerical Methods in
Engineering, Vol. 15, 1980, 14131418
4. T. J. R. Hughes and T. E. Tezduyar, Finite Elements Based Upon Mindlin Plate
Theory With Particular Reference to the Fournode Bilinear Isoparametric
Element, Journal of Applied Mechanics, (September 1981), 587596.
5. K. C. Park and D. L. Flaggs, A Fourier Analysis of Spurious Mechanisms and
Locking in The Finite Element Method, Computer Methods in Applied
Mechanics and Engineering, 46 (1984), 6581.
6. A. B. Boresi, O. M. Sidebottom, F. B. Seeley and J. Smith, Advanced Mechanics
of Materials, 3rd ed., New York: John Wiley & Son Inc., 1985.
7. A. C. Ugural and S. K. Fenster, Advanced Strength and Applied Elasticity, New
Jersey: PrenticeHall Inc., 1987.
8. T. J. R. Hughes, The Finite Element Method: Linear Static and Dynamic Finite
element analysis, New Jersey: PrenticeHall, 1987.
9. R. D. Cook, D. S. Malkus and M. E. Plesha, Concepts and Applications of Finite
Element Analysis, 3rd., New York: John Wiley, 1989.
10. J. E. Shigley and C. R. Mischke, Mechanical Engineering Design, 5th ed., New
York: McGrawHill, 1989.
68
11. T. R. Chandrupatla and A. D. Belegundu, Introduction to Finite Elements in
Engineering, New Jersey: PrenticeHall Inc., 1991
12. R. D. Cook, Finite Element Modeling for stress Analysis, New York: John
Wiley & Son Inc., 1995.
13. K. J. Bathe, Finite Element Procedures, New Jersey: PrenticeHall Inc., 1996
14. J. A. Trapp, Course notes, ME 5153., 1996
15. I. H. Shames and C. L. Dym, Energy and Finite Element Methods in Structural
Mechanics, USA: Hemisphere Publishing Corporation., 1985.
69

PAGE 1
A STUDY OF THE APPLICATION OF MINDLIN PLATE ELEMENTS TO THIN PLATES by Piyoros Jirawattana B.Eng., Chiang Mai University, 1994 A thesis submitted to the University of Colorado at Denver in partial fulfillment of the requirements for the degree of Master of Science Mechanical Engineering 1997 r
PAGE 2
This thesis for the Master of Science degree by Piyoros Jirawattana has been approved Date
PAGE 3
Jirawattana, Piyoros (M.S., Mechanical Engineering) A Study of the Application ofMindlin Plate Elements to Thin Plates Thesis directed by Assistant Professor Samuel W. J. Welch ABSTRACT This thesis is a study of Mindlin plate theory in finite element analysis with particular emphasis on locking and instability behaviors of different element formulations. A computer model is used to simulate the bending action of thin plates and to study locking and instabilities associated with the different element formulations. We start with the bilinear Mindlin plate element. Using techniques of numerical integration we demonstrate that "shear locking" can be eliminated and a very competitive convergence results. This is done by using a technique known as selective reduced intergation. However, the selectivereduced integration method can lead to an unstable mode in the finite element solution. Finally, we consider the development of a Mindlin plate element that does not lock and does not contain unstable modes. This is done by approximating the transverse shear strain using shape functions called mixed interpolation functions. The detection of locking and instability will be presented by imposing spatial Fourier modes and comparing the strain energy from finite elements to the strain energy of the continuous Mindlin plate theory. This idea can also be used in the analysis of convergence of the solution of finite element analysis. This abstract accurately represents the content of the candidate's _thesi,' I recommend its publication. Signed ,. / Samuel W. W lch iii
PAGE 4
ACKNOWLEDGMENTS I would like to thank Samuel W. J. Welch for the thesis topic and his guidance. I also wish to thank John A. Trapp for the fundamental plate theory from advanced strength class. iv
PAGE 5
CONTENTS Acknowledgments ............................ ............. ....... ..... ........................................ .... iv CHAPTER 1. INTRODUCTION AND PLATE THEORETICAL PRELIMINARIES ................ 1 Scope ................................................................................................................. 1 Kirchhoff and Mindlin Plate Bending Theories .................................................. 2 Introduction .................................................................................. ............... 2 StressStrain Relationship for Flat Plate ............................................................ 2 Kirchhoff Theory .......................................................................................... 3 MomentCurvature Relation ................... ......................... ...... ... ... ........ 6 Mindlin Theory ............................................................................................. 7 MomentCurvature Relation .................................................................... 1 0 2. FINITE ELEMENTS FOR PLATES .................................................................. 11 KirchhoffPlate Elements ................................................................................. 11 Mindlin Plate Elements .................................................................................... 13 Stiffness Matrix .......................................................................................... 14 Discussion of Kirchhoff and Mindlin Plate Elements ....................................... 16 3. INSTABILITY AND LOCKING ...................................................................... 18 Introduction ......................................................................................... .......... 18 Locking ............................................................................. ............................ 19 Locking in Mindlin Beam ........................................................................... 19 Numerical Examples Using Mindlin Plate Elements ...................................... 23 Instability ........................................................................................................ 27 Zero Strain Energy in Mindlin Beam .......................................... ............... 27 Eigenvalue Analysis of the Zero Strain Energy Modes ................................ 30 Analysis oflnstability and Locking ofBilinear Mindlin Plate Elements Using Spatial Fourier Modes ...................................................................... 31 Typical Instability Mode in Mindlin Plate Element ....................................... 38 4. THE DEVELOPMENT OF A MINDLIN PLATE ELEMENT ......................... .41 Approach ....................................................................................................... 41 v
PAGE 6
Numerical Examples .................................................................................... 46 The Accuracy ofthe Mixed Interpolation Element ...................................... .48 5. CONCLUSION ., ........................................................................................... 50 APPENDIX A. Stress Resultants in a Flat Plate ................................................................... 51 Force and Moment Balance ofFlat Plates .................................................... 55 Force Balance for Plates ........................................................................... 55 Moment Balance for Plates ....................................................................... 56 B. Derivation ofFourNode Quadrilateral Mindlin Plate Element ....................... 58 Numerical Integration .................................................................................. 65 Two Dimensional Integrals ...................................................................... 65 Stiffness Integration ................................................................................ 66 REFERENCES ................................................................................................. 68 vi
PAGE 7
CHAPTER 1 INTRODUCTION AND PLATE THEORETICAL PRELIMINARIES Scope Plate bending has been studied in the field of classical strength of materials for a long time. The problem of plate bending can be considered as an extension of either elementary beam theory or Mindlin beam theory depending on which theory is appropriate. Both plates and beams support transverse loads by bending action. Two theories of classical plate bending will be discussed, Kirchhoff theory and Mindlin theory. The difference between these plate bending theories can be classified by the magnitude of the thickness compared to the magnitude of the other dimensions and also the magnitude of the transverse displacement compared to the thickness. The plate in this context will be considered as a relatively thin plate with small deformation. The classical Kirchhoff plate bending theory makes an assumption that the transverse shear deformation of thin plates in bending can be neglected. Information about the resulting deformation can be assumed to be contained in the function and in it's derivatives. Therefore, finite elements based on this theory must use C1continuous elements. Even though this classical theory is applicable to thin plates, it's use in finite elements is problematic as it requires higherorder nodal interpolation function to enforce interelement compatibility. These higher order interpolations can be quite complicated and impractical for general use. In Mindlin plate bending theory, the transverse shear deformation is taken into account so that this theory is applicable to thick plate bending. The deformation of Mindlin plate theory consists of one transverse deformation and two small angle rotations. The resulting finite element implementation requires only CJ continuous elements so that lower order interpolation functions may be used. In addition, because the deformation in Mindlin plate theory is more physical, it is easier to apply boundary conditions. However, the application of Mindlin plate elements to thin plates may not be as accurate as Kirchhoff plate elements because the transverse shear strains from Mindlin theory, if approximated inaccurately, may make thin plates
PAGE 8
too stiff. We will see that by the use of special techniques, Mindlin plate elements can be made accurate for thin plate problems. Thus, the context of this study is the application of Mindlin plate elements to the bending of thin plates. We will begin with general background theory of Kirchhoff plate bending and Mindlin plate bending in chapter 1. Next, the application of both theories in finite elements will be discussed in chapter 2. Chapter 3 introduces the concept of "shear locking" and how methods to prevent this problem result in the side effect known as "instability". Finally the development of a Mindlin plate element for thin plates that does not lock and does not exhibit instability will be outlined in chapter 4. The numerical results presented in this study were obtained using software written by the author. The software is not included as part of this thesis. Kirchhoff and Mindlin Plate Bending Theories This section introduces, discusses and compares the Kirchhoff and Mindlin plate bending theories. A more comprehensive discussion of these theories is contained in reference [ 15] and in the references therein. Introduction A flat plate supports transverse loads by a bending action similar to a straight beam but the deformation modes are two dimensional and thus are more complicated. Plate theories assume that the plate thickness is relatively small compared to the plate dimensions in the plane of the middle surface. Thus the bending behavior of a plates in both Kirchhoff and Mindlin theory depends on the plate thickness. The former theory makes the assumption that straight fibers normal to the plane of the middle surface remain normal to the plane of the middle surface after bending. The latter theory does not make this assumption. Before discussing the details of these theories, we will start by discussing the general stress strain relationships for flat plates. StressStrain Relationship for Flat Plate The stressstrain relationship may be obtained from Hooke's law for isotropic elasticity. 2
PAGE 9
(I. I) In both Kirchhoff and Mindlin theories CT33 is considered negligible. From equation (I. I) with CT33 = 0, we obtain (I.2) substituting this into Hooke's law we obtain (I.3) equation (1.3) is the general stressstrain relation for a flat plate and from equation (I. I) we get the equation for transverse shear stress CTaJ = 2GEa3 (1.4) where a and Prange from I to 2 8ap= I if a= P and 8ap = 0 if a :t: p Kirchhoff Theory iW u1 =x3 ......... ... ...... Figure 1.1 Kirchhoff plate after loading 3
PAGE 10
In Kirchhoff theory, it is assumed that the straightline normals to the undeformed middle plane of the plate remains straight and inextensional under the deformation of the plate. In addition, it is assumed that the normals in the deformed plate remain normal to the deformed middle plane which is a reasonable assumption for thin plates. Under this assumption, it follows that EaJ = 0 which implies that aaJ = 0, hence the shear forces Q 1 and Q2 are implied to be zero. Note that aaJ is neglected in the stress strain laws but Qa is kept in the equations of motion. Despite this contradiction, the theory is more accurate than the membrane theory of plates [6]. If we expand transverse displacement u3(x1,x:z,x3 ) about the mid plane using power senes m x3 we get [ 14] X ifthe plate is thin 3 << 1 and we can approximate I from the assumption that Ea3 = 0 we get (1.5) and if the plate is thin we can introduce equation (1.5) into the equation above, then integrating with respect to x3 we get (1.6) Ua(x1,x2 ) is constant of the XJ integration which has no effect on the moment Map and thus may be neglected. The resulting kinematics assumptions for Kirchhoff theory (neglecting 0'13 and a23 in the stressstrain law) are 4
PAGE 11
(1.7) (1.8) (1.9) 1 (iii iii J from & all = _a + _fJ we get 2 tJxp tJxa (1.10) (1.11) (1.12) Equations ( 1. 7) and ( 1. 8) are the kinematics assumptions of Kirchhoff theory mentioned earlier. Considering equations (1.7), (1.8) and (1.9), we can see that the deformation in x1 and x2 direction are related to w by first order differential equations. This leads to the secondorder differential equations for strain in equations (1.10), (1.11) and (1.12). To get the expression for bending moments in term of the transverse displacement, we substitute equations (1.10), (1.11) and (1.12) into equation (1.3), and then use the definition of bending moments (see the details of this definition in appendix A) to obtain h 2 M afJ = J x3u afJdxJ h 2 5
PAGE 12
(1.13) (1.14) (1.15) where D is the flexural rigidity (define in the next section) which is approximately equivalent to E times the moment of inertia of a unit width plate. MomentCurvature Relation. In Kirchhoff theory we neglect u33 and also transverse shear strain, so that we can consider the in plane stressstrain relation. By assumption the material is isotropic. Let x1 x2 be the principle axis of material then we can write v 0 v 1 0 0 0 (1v) 2 (1.16) We obtain the momentcurvature relation by substitution of equations (1.1 0), ( 1.11) and ( 1. 12) into equation ( 1.16) and the result into the definition of bending moments, yields [ 12] vD D 0 6 (1.17)
PAGE 13
Eh3 where D = 2 E is the elastic modulus and v is Poisson's ratio. 12(1 v ) Mindlin Theory In Kirchhoff theory normals to the midsurface remain normal during deformation and hence transverse shear deformation is neglected. Mindlin theory makes an attempt to keep the transverse shear deformation which by assuming that fibers originally normal to the midsurface can rotate independently ofthe normal direction. Ut=x (} 3 .... ,' Figure 1.2 Mindlin plate after loading iW We define the transverse shear strain, y =(} and note that in Kirchhoff a theory this transverse shear strain is zero. In Mindlin plate theory we assume that the deformation in the x1 and x2 directions is a function of small angle rotations so that we can write, 7
PAGE 14
(1.18) ( 1.19) (1.20) Equations (1.18) and (1.19) have included some rotation of fibers originally normal to the midsurface to orientations not normal to the midsurface direction. To obtain the expression for strains in Mindlin theory, introduction of equations ( 1.18), ( 1.19) and (1.20) into the definition of E;J =_!_(iii; + iUJJ yields 2 Ox} Ox; ((}a.fJ+(}fJ.a) EaPx3 2 and E = _! ( CU 3 + CYI a ) = _!_ (} ) a3 2 Ox Ox 2 Ox a a 3 a from equations (1.21) and (1.22) we can write the expression for strains as where (}a fJ = Of} a 0(}/l 8 (1.21) (1.22) (1.23)
PAGE 15
Considering equation ( 1.22), C1a3 is independent of x3 direction, so it is constant over the plate thickness. In reality CTa3 has a parabolic profile over plate thickness. Instead of developing a very complicated theory to approximated CTa3 in parabolic manner with x 3 (zero at h'2,hl2. and maximum at x3 = 0). A k factor is added in the transverse shear force constitutive equation. We can evaluate the correct k shear correction factor by forming the strain energy for plate theory corresponding to a 3dimensional energy with a parabolic shear stress.[14] To get the expression of bending moment and the transverse shear force in Mindlin theory in terms of the displacement function, substitute strain equations (1.23) into the stressstrain relations (1.3) and (1.4). Then using the relationships and h 2 Map = I X3C1 apdx3 h 2 h we obtain M h 8 B +B + 2G+3 [ vE { ) ( 1 ( 00 a Of) p JJ] afJ 1v2 a/1 .... ..... rl 'rl 2 Oxp Oxa (1.24) (1.25) It should be noted as transverse shear strain has been included, the transverse shear stress can be written as {1.26) 9
PAGE 16
MomentCurvature Relation. Let x1 and x2 be the principle axis of the material. The MomentCurvature relation for Mindlin theory includes the shear stressstrain relations 0"23 = Gy23 and 0"31 = GyJJ. The MomentCurvature relation is abbreviated as {M}=[Du] { K}. Again, we assume that the material is isotropic, from equations (1.24) and (1.25) we can explicitly write Mil 0 0 0 (} x. '.r M22 [DK] 0 0 0 (} %1 'Zl Ml2 = 0 0 0 (} +B .1'1 'Zl %2 'Zt (1.27) Q2 0 0 0 Gh 0 e.r2 w,.r2 Ql 0 0 0 0 Gh e.r w,.r I I Where [DK] is the same as equation (1.17). Compare equation (1.27) with equation (1.17), as we mentioned ,Kirchhofftheory neglects transverse shear deformations so that in Kirchhofftheory the transverse shear forces Q1 = Q2 = 0. 10
PAGE 17
CHAPTER2 FINITE ELEMENTS FOR PLATES The preceding chapter is devoted to the analytical derivation of Kirchhoff and Mindlin plate bending theories. We might observe that numerical approximations of Mindlin theory should be easier because the governing equations are of lower order. We will discuss finite elements using these two theories. Kirchhoff Plate Elements As we have discussed in Chapter 1, Kirchhoff theory neglects transverse shear deformation and is applicable to thin plates. Strain energy is determined by inplane strains &11, &22 and YJ2 From equation (1.3), strains are determined by the lateral displacement field w = w(x1,x2). Nodal d.o.f. of Kirchhoff elements consist of lateral deflections w, and small rotations w,... and w,... of midsurface normals. We now I 2 express the displacement field in terms of independent shape function interpolations as [12] (2.1) (2.2) (2.3) In order to formulate an element stiffness matrix [k]e, we now form the strain energy. (2.4) These strains of Kirchhoff theory in equations (1.10), (1.11) and (1.12), may be written in matrix notation as II
PAGE 18
4 .x2 Figure 2.1 Plate 12 d.o.f. at node 3 of Kirchhoff plate element We consider that the thickness constant over the plate, so that we have dV = hdA and finally dA = dxdy. Now using equation (1.16), equation (2.4) can be rewritten as (2.5) where {K"f = [w,.r1.r1 w,.r2.r2 2w,.r,.r2 ] is known as the curvature for Kirchhoff elements and [DK] is known as the flexural rigidity From equations (2.1), (2.2) and (2.3), we can write the nodal displacement matrix of a Kirchhoff element as 12
PAGE 19
{d} = [w1 w,.r 1 w,.r 1 I 1 (2.6) Equation (2.6) can be written via the nodal interpolation shape functions as w = [N]{d} (2.7) equation (2. 7) allows us to write the curvature using nodal shape functions as {I(}= [B]{d} Where [B] is the differentiation matrix that contains information about the curvature. Substituting equation (2. 7) into equation (2. 5), we obtain the element stiffness matrix,[k] (2.8) where [k] = j[Bf[DK ][B]dA A Mindlin Plate Elements Mindlin plate theory includes transverse shear deformation, therefore it is a theory that is reasonable for thick plates. Nodal d.o.f. of Mindlin plate elements consist of one lateral deflection w, and two rotations, (} .r 1 and (} .r 1 of midsurface normals. If all interpolations use the same I 1 polynomial, then for an element of N nodes, we can write (2.9) where (2.1 0) 13
PAGE 20
Consider equations (1.23) and (2.9), the field for w is coupled to the fields for 8z I and 8z only via the transverse shear strains 'Y23 and [12]. 2 4 3 w 8 1 2 z, Figure 2.2 12 d.o.f at nodes 2 ofMindlin plate element. Stiffness Matrix In a like manner as we did for the Kirchhoff plate element. We can write the formula for the strain energy U [ 12] h 1 2 U =J J {e}[E]{e}ch3d4 2 A h 2 14 (2.11)
PAGE 21
where A is the area of the plate mid surface and { E} = [c X X EX X r X X r X X r Xr ] Expressions for each strain are stated in I I l l I l l l rl equation {1.23). Integration through the thickness yields (2.12) where [DM] and {K'} are stated in equation (1.27). We can write {K}= [L]{u} (2.13) 0 0 0 a1 0 0 0 a2 where 0 0 0 [L]= a2 a1 0 0 &2 0 1 0 a1 Equations (2.9) and (2.13) yield {K} = [B]{d} where 15
PAGE 22
0 N":r.. 0 0 0 Nl.r l [B] = [L][N] = 0 Nl.r Nl.r (2.14) l I NI.rl 0 Nl ........... Nl.r Nl 0 I And finally from equations (2.12) and (2.14) where Discussion of Kirchhoff and Mindlin Plate Elements Sometimes Mindlin theory is called thick plate theory and Kirchhoff theory is called thin plate theory. Consider the d.o.f in the displacement matrix, {d}, of Kirchhoff and Mindlin plate elements. In Kirchhoff theory, the d.o.f at node i are one lateral displacement w, and two rotations, w,.r, and w,.r:t. In Mindlin theory, the d.o.f at node i are one lateral displacement w,. and two rotations (}.r, and (}.r, We I l can see that the displacement matrix of Mindlin has lower order partial differentials that make it more physical and easier to apply numerical approximations. In addition, the higher order partial differentials in Kirchhoff theory requires at least biquadratic shape function that lead to complicated numerical evaluation. Mindlin theory needs only bilinear shape functions. Considering equations (1.7), (1.8) and (1.9), Kirchhoffplate elements must display interelement continuity of w, w .r and w .r (compatibility). Consider the connection I l between elements in the x2 direction, continuity of w forces the continuity of the tangential slope, w,.r but not continuity ofthe boundarynormal slope w,.r Now for l I Mindlin theory, consider equations (1.18), (1.19) and (1.20). Compatibility requires the continuity of the fields, w (} .r and (} .r at the interelement boundaries. Mindlin I l plate elements are easier to enforce compatibility because the lateral displacement fields are not related to the rotation fields by differentiation. 16
PAGE 23
In the next chapter, we will show that by using selectivereduced integration, Mindlin plate elements can capture behavior of thin or Kirchhoff plates. These are the reasons that Kirchhoff plate elements are not widely used. 17
PAGE 24
CHAPTERJ INSTABILITY AND LOCKING Introduction So far, we have discussed the derivation of Mindlin plate theory and it's implementation in finite element theory. In this chapter we compare results obtained after selectivereduced integration and full integration are used for the integration of the element shear stiffness matrix. This leads to the study of two important concepts in finite element analysis that can make practical results unrealizable, instability and locking. Numerous researchers have made attempts to prevent these problems. p Figure 3.1 Cantilever beam 18
PAGE 25
Locking Mindlin plate theory make plate bending theory more physical by including transverse shear deformations into account. However in finite element analysis approximation of these transverse shear strains can lead to inaccurate results. Physically the deformation of thin plates should be dominated by the bending deformation but if it is dominated by transverse shear deformation we call this "shear locking". Locking in Mindlin Beam To explain how this shear locking occurs we consider the cantilever beam of length L with load P in Figure 3. 1 and we use two different integration methods, full integration and selectivereduced integration. We employ Mindlin plate theory with the nodal d.o.f w2 Oz. We separate the stiffness matrix [ k] into a component with bending stiffness [ k b ] and a component with transverse shear stiffness [ k.] We do so in order to integrate the bending and the shear terms separately. We obtain [Bb] that is associated with inplane strain B [0 1 I L 0 1 I L] [ b]0 0 0 0 (3.1) and [B.] that is associated with transverse shear strain [ 0 0 0 0 ] [B] = 11L xiL1 IlL xiL (3.2) case 1) Full integration Introduction of equations (3 .1) and (3 .2) into the expression for the stiffness matrix and using exact integration for both the bending stiffness [kb] and the transverse shear stiffness [k.] we obtain the stiffness matrices 19
PAGE 26
I L [k.] 2 I L 2 0 0 0 I 0 0 0 1 0 L I L 2 2 ( L 2 L I L 2 2 (:) L ( 2 With k = we solve for the tip displacement 6 12(%r + 20 (6PL) (h)2 5GA I2 +5 L for a thick beam, h >> L the above solution becomes 6PL w ==25GA (3.3) (3.4) (3.5) (3.6) Equation (3.6) is a reasonable result for a thick beam with deformation caused by shear. Compare equation (3.6) with the deformation caused by shear from Castigliano's theory, w2 = 1.33 PL [IO]. Thus, the Mindlin beam provides a AG reasonably accurate solution for thick beams. 20
PAGE 27
Now for a thin beam, h << L, from equation (3.5) we obtain 24PL w =25GA (3.7) In the latter case the deformation should be dominated by bending (compare to the analytical solution from elementary beam theory w2 = PL3 ) but it is not. It is 3/ dominated by shear. This is an example of shear locking. Now we consider what happens if selectivereduced integration is applied. Case 2) Selective integration In this case, we integrate bending stiffitess with two Gauss points and use one Gauss point for the shear stiffitess. This yields the same bending stiffitess as before but we get a different shear stiffitess. 1 L 1 L 2 2 L L [k.] 2 2 (3.8) L L 1 1 2 2 L L 2 2 Solving for the deflection at the free end, we get (3.9) The first term and the second are deformation from shear and bending, respectively. When the beam becomes thin, the second term contributes and we can neglected the first term. The solution from elementary beam theory (using Poisson's ratio of 0.3 ) 21
PAGE 28
2PL3 IS w = Selectivereduced integration with a Mindlin beam finite element 2 1.3Gbh3 produces a result of 97.5 % accuracy. It should be noted that for a thick beam solution, either full integration or selectivereduced integration can capture the behavior of shear deformation because the bending deformation can be neglected. It is counter to our intuition that an exact integration produces an inaccurate result but an underintegration might introduce an acceptable answer. In order to explain, we consider the cantilever beam in Figure 3.1, obtained by simply reducing a Mindlin plate element into one dimension. We can write the equation for axial strain c 11 and transverse shear strain r 31 as (3 .1 0) (3 .11) where the notation is the same as that of the Mindlin plate element. We separate the element strain into the bending strain and the transverse shear strain. We can write (3.12) Our discussion will be concerned with the transverse shear approximation in equation (3.12). (3.13) If we consider equation (3.13) we see that the terms associated with the transverse deformation are constant. The terms associated with the rotational degree of freedom are linear in x. If we use selectivereduced integration we will evaluate (3 .13) at x = L/2, we are in effect replacing (3. 13) with (3.14) 22
PAGE 29
In the literature it is common to see the above process referred to as eliminating the excess monomial. This excess monomial results from using inappropriate order polynomials for the displacement and rotation fields. When the beam becomes thin, the transverse shear terms should become infinitesimal. The shear resulting from evaluating equation (3.13) at two Gauss points cannot, in general, be infinitesimal at both points and the resulting error becomes magnified in thin beams as the stiffuess matrix divides this error by the thickness squared. Thus, in the case of thin beams, if full integration is used with the shear stiffuess matrix, the deformation will be dominated by [ k sl not by [ k b ] This is called "Locking". In the literature the use of full integration on the shear stiffuess terms is sometimes referred to as over constraining the element to have zero shear at two points thus causing locking. Numerical Examples Using Mindlin Plate Elements We now demonstrate the concepts we have been discussing with numerical simulations using various Mindlin plate elements. We will numerically demonstrate that futl integration leads to locking that can be suppressed by selectivereduced integration. This section will present numerical examples of Mindlin plate element used on a square plate with various boundary conditions and loads. We model the square plate with lOxlO elements as shown in Figure 3.2 23
PAGE 30
a I b Figure 3.2 Square plate The data for the plate consists of the following E = 3 x 107 psi v= 0.3 a= b = 1 inch. r I X Considerations similar to those in Mindlin beam elements can be applied to the Mindlin plate element. If the plate is thin, the transverse shear strain in both directions (equations (3.15)) should be small. YzJ = w, .. l O .. l = 0 y31=w, .. O .. :::0 I I (3.15) Using bilinear shape function for Mindlin plate elements, twopoint Guassian quadrature is required to exactly integrate the shear stiffitess. This full integration brings in the higher order terms of the rotational degrees of freedom that the terms 24
PAGE 31
corresponding to the transverse defonnation do not have. Thus, full integration can lead to shear locking unless both transverse shear strains are evaluated at x1 = x2 = 0 (selectivereduced integration). The point being that the excess monomials in the tenns corresponding to the rotational degrees of freedom are eliminated. Using selective integration with Mindlin plate elements can prevent locking and gives acceptable results for thin plates as shown in Figure 3.3. We can see that when the plate become thin, exact integration gives unacceptable results due to locking. r 'I 0.6 L<<}\(\ .. c 0.4 0.2 a) Clampededges with distributed load. 25
PAGE 32
.. 0.8 u t I o.s u i E .. c 0.4 0.2 0 0 50 100 150 ... b) Simple supports with distributed load 0 50 100 200 150 c) Simple supports with point load at the center of plate 250 300 200 250 Figure 3.3 Numerical/ Analytical Ratio of center deflection of 1 Ox 10 element, square plate. Thin plates correspond to large ab/t. 26
PAGE 33
Instability Generally, finite elements do not give exact results. We have demonstrated that the use of a more accurate stiffuess integration might make an element too stiff causing locking. We saw that using selectivereduced integration eliminates locking. However, using selectivereduced integration can also cause something worse than locking called instability. Zero Strain Energy in Mindlin beam Instability is sometimes referred to as a zero strain energy mode. This is because the global displacement matrix results in zero strain at the Gauss points so that the mode is easily excited. For this reason, instability is often associated with rank deficiencies of the element stiffuess matrix which is caused by lowerorder integration. To explain the term "zeroenergy mode", we form an expression for the element strain energy. U. = {df[k]{d} = f {&f[D]{&}d4 A, where [k] = j[Bf[D][B]dA A, and {E} = [B]{d} When we evaluate [ k] numerically, it contains only the information that can be sensed at the sampling points of the quadrature rule. If the strains vanish at all Gauss points for a certain mode {d}, then U. must be equal to zero, and [k] must be a zerostiffuess matrix. We know that the strain energy must vanish for only rigid body motion. If it vanishes for other motions, an instability has happened. To support the above discussion, we consider the inplane twisting mode of the two Mindlin beam elements shown in Figure 3. 4. We will compare the strain energy using stiffuess matrices derived with both full integration and selectivereduced integration. 27
PAGE 34
. . . Figure 3.4 Instability in the Mindlin beam element, dashedline shows inplane twisting instability mode. The global displacement matrix of the 3 nodes is WI OJ {d} = w2 82 w3 83 The displacement matrix ofthe inplane unstable twisting mode is {d} = 0 1 0 1 0 28 (3.16) (3.17)
PAGE 35
The global stiffuess matrix, assembled from the element shear stiffuess in equation (3.8) (obtained by selective reduced integration), is L 1 L 0 0 2 2 L L2 L L2 0 0 2 4 2 4 1 L 2 0 1 L ks = kGA 2 2 (3.18) g L L L2 L2 L L2 0 2 4 2 2 4 0 0 1 L L 2 2 0 0 L L2 L L2 2 4 2 4 where k; is global shear stiffuess matrix obtained by selectivereduced integration. Now by applying equations (3 .17) and (3 .18) into the expression of strain energy, we obtain UT = {d}1[k;]{d} =0 where U 1 denotes the total strain energy. This result corresponds to the definition of zerostrain energy mode. In a like manner we use full integration and obtain the exactly integrated shear stiffuess matrix, k { 29
PAGE 36
L 1 L 0 2 2 L L2 L L2 0 2 3 2 6 1 L 2 0 1 kr = kGA 2 g L L L2 2L 2 L 0 2 6 3 2 0 0 1 L I 2 0 0 L L2 L 2 6 2 Now calculating the total strain energy, we obtain UT = {d}T[k{ ]{d} 2L2 = 3 0 0 L 2 (3 .19) L2 6 L 2 L2 3 We can see that the zerostrain energy mode is not present with exact integration. However we should be concerned about exact integration because of shear locking. Eigenvalue Analysis of the Zero Strain Energy Modes There is another approach to detect instability in the finite element method. This approach performs an eigenvalue analysis of the element stiffness matrix. Consider the element equilibrium equation [ 12] [k]{d} = {f} Now if the load matrix {f} is proportional to element nodal displacements { d} through a factor A., we have [k]{d} = {f} = A.{d} 30
PAGE 37
or ([k]A[/]){d} = {0} (3.20) This is an eigen problem with the eigen values of { k }" being the A. Each A; corresponds an eigenvector { d} ;. If each {d); is normalized so that {d}; {d}, = 1, from equation (3.20), we obtain or 2U; =A, (3.21) where U; is the strain energy in the element when its nodal d.o.f are normalized displacements { d} ;. [k] is complete element stiffness matrix of the element in an infinite media. Equation (3 .21) can be interpreted in the same way as was the expression of strain energy, that is [k] should yield A;= 0 for the rigid body motion case only. If it is not a rigid body motion and A; = 0 then instability is present. Generally a solution with an instability is easy to detect as it usually contains short wavelength, high amplitude unphysical deformations. Analysis of Instability and Locking of Bilinear Mindlin Plate Elements Using Spatial Fourier Modes Corresponding to the above discussion we might say that instability happens when strain energy vanishes for deformation modes not related to rigid body motion. Next we will examine the instability and locking by applying Fourier modes on an infinite media. We begin by writing the expression for strain energy for a linearly elastic isotropic Mindlin plate 31
PAGE 38
+ {} + {} dx Eh )2 J2f 4(1 + v) If acl I ac2 2 I 2 (3.22) We will use the definition above to calculate a continuous strain energy to compare with the strain energy calculated from the finite element methods. We write the nodal deflections w,, ()"'; and ()"' ; as general harmonic waves in 2 dimensions I 1 (3.23) For thin plates we can neglect the transverse shear deformations, so that (3.24) (3.25) We note that in Mindlin theory, other modes may be present but we expect these to be the dominant modes. These three equations above can be written via nodal shape functions as (3.26) (3.27) (3.28) where w is the node element deflection () "' is the small angle rotation in x1 direction ()"' is the small angle rotation in x2 direction 1 32
PAGE 39
k"' k"' is the wave number in each direction I 2 dx1 is the horizontal length of element dx2 is the vertical length of element n is the element index in horizontal direction m is the element index in vertical direction Equations (3.26), (3.27) and (3.28) allow us to write the displacements as a function of the wave number and then calculate the strain energy. We will compare the continuous and finite element strain energies calculated using these modes. case 1 ) k"' = k"' = 0. I 2 This case corresponds to rigid body motion. It is certain that strain energy must vanish. This is true for any method. case 2) maximum number of wave, k"' x, = 1r, k"' x2 = 1r I 2 The maximum admissible wave number that can be represented by nodes of finite elements for the transverse deformation is presented in Figure 3.5. This is the characteristic instability called whourglass mode. . . . Figure 3.5 The maximum admissible wave number. 33
PAGE 40
The wave number corresponding to the wave length is (3.29) where A. is the length of one wave. At the minimum admissible wave length A, which is equal to the length of two elements, 2dx, we have the maximum wave number 7i kmax = dx (3.30) The deformation at the maximum admissible wave number is not a rigid body motion. If we apply the maximum wave number into equations (3 .26), (3 .27) and (3 .28) and calculate the strain energy we will obtain, in some cases, the zero strain energy mode. We will show that the zero strain energy mode only appears, in some cases, at the maximum admissible wave number. We will do this calculation numerically with a single one inch square element. In the case of Mindlin plate elements, only the whourglass mode is possible. We will discuss this in more detail later. We can write the displacement (3.31) Using equation (3.31) in the expression for strain energy, we can calculate the strain energy as a function of the wave number using stiffness matrices obtained using full integration and selectivereduced integration. We compare these solutions with the continuous solution using the same Fourier modes. Figure 3.5 presents the strain energy of a single one inch square Mindlin plate element. We set the shortest wave length equal to 2 inches so that (k" x1)m"" = (k" x2)m"" = 7i. Consider Figure 3.5 a I 2 and c, The locking behavior makes the strain energy obtained from full integration too big, even at a very small wave number. In Figure 3.5 a and c we plot the logarithm of strain energy as the strain energy obtained from continuous case is very small 34
PAGE 41
compared to that obtained using full integration. Compare these results with the results from selectivereduced integration in Figure 3.5 b and d, we can see that the smaller the wave number, the more accurate these elements are. Practical use of these elements does not generally consider deformations near the maximum admissible wave number and the accuracy of these elements would be acceptable were it not for the behavior at the maximum admissible wave number. At the maximum admissible wave number the strain energy using selectivereduced integration vanishes indicating the presence of an unstable mode. Strain Energy Log (lb.inch) "'Numerical Solution 5 0 Analytical Solution 5 10 15 0 0.5 1 1.5 2 2.5 3 3.5 Wave Number (kAx) a) Strain energy from full integration and analytical from ktu = 0 to ;r 35
PAGE 42
Strain Energy (lb.inch) Numerical Solution 2000 1500 1000 500 Analytical Solution 0 0 0.5 1 1.5 2 2.5 3 3.5 Wave Number (k.1x) b) Strain energy from selectivereduced integration and analytical from kta = 0 to 1r Strain Energy Log (lb.inch) 1 5 10 15 20 Numerical Solution Analytical Solution 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 Wave Number (k.1x) c) Strain energy from full integration and analytical from kta = 0 to w'IO 36
PAGE 43
Strain Energy (lb.inch) 0.03 0.02 0.01 0.05 0.1 Numerical Solution \ 0.15 0.2 0.25 0.3 0.35 Wave Number (k!u) d) Strain energy from selectivereduced integration and analytical from kllx = 0 to wlO Figure 3.5 The strain energy of a single one inch square Mindlin plate element in an infinite media, thin plate limit. We note that Figure 3.5 indicates the expected result that for a given wavelength, the numerical results become more accurate as llx + 0 37
PAGE 44
Typical Instability Mode in Mindlin Plate Element We split [k] into the bending stiffness [kb] and the transverse shear stiffness [ks]. We integrate the bending stiffness fully and we under intergrate the shear stiffness which makes only two modes of instability possible, whourglass and inplane twist mode. The latter can not happen with two or more element because it can not be communicated between adjacent elements, so that only the whour glass mode appears. In general the global stiffness matrix is singular by itself and zeroenergy modes can be prevented after boundary conditions are applied. From the assumption that the boundary conditions preclude the rigid body mode from appearing in one element, it also precludes the mode in the remainder of the mesh [8]. A critical test is to clamp one node in the middle of the square plate, and the instability appears as shown in Figure 3.6. Now we fix one node that is adjacent to the center node. The instability is gone. To retain symmetries as shown in Figure 3. 7, we fix 3 adjacent clamped nodes in the middle of the plate. 38
PAGE 45
0 15 Figure 3.6 whour glass mode with clamped center node of IOxiO Mindlin plate elements with selective integration and distributed load 39 15
PAGE 46
0.01 0 0.01 0.02 0.03 0.04 15 0 0 Figure 3. 7 Deflection of 3 clamped center nodes of lOxlO Mindlin plate elements with selective integration and distributed load. 40 15
PAGE 47
CHAPTER4 THE DEVELOPMENT OF A MINDLIN PLATE ELEMENT The preceding chapter was devoted to the study of instability and locking in Mindlin plate elements with bilinear shape functions. We showed that using exact integration of the element stiffness matrix resulted in shear locking Shear locking can be suppressed by adopting selectivereduced integration, however for some specific situations selectivereduced integration introduces instability In this chapter we will present an reliable element that does not lock and contains no unstable modes[4]. Approach The development of Mindlin plate element obtained from the method of selective reduced integration will be tried with a different point of view Let us consider the curvature of Mindlin plate element. (4.1) (} w .1'1 'XI The first three terms in the curvature matrix are bending stressstrain relations, a11, 0"22 and 0"12. respectively. The last two terms on the bottom are the transverse shear deformation For the isoparametric element we used the standard nodal shape interpolations as follow, (4.2a) 41
PAGE 48
1 N3 = 4 Where each node is located in Figure 4.2 ( 1,1) ......_,(1, 1) (1,1) (1,1) Figure 4.1 The quadrilateral element 77 coordinate (the master element) (4.2b) (4.2c) (4.2d) In the previous chapter we showed that selectivereduced integration solves the shear locking problem and gives good convergence behavior, but leads to unstable modes. To motivate the development of a better Mindlin plate element we consider the positive characteristics of selectivereduced integration. The one Gauss point located at the center of an element results in constant transverse shear strains over the element in each direction. The element presented in this chapter has similar properties in that one of the shear strains is constant along element edges. The shear strain is obtained by approximating the slope of the transverse deformation with out strict adherence to the element shape functions. The resulting expression allows us to use two Gauss points and integrate the resulting polynomials exactly. 42
PAGE 49
From the definition oftransverse shear strain (4.3) (4.4) We can define these transverse shear strains with different interpolation functions from the in plane stressstrain terms [ 13] (4.5) (4.6) where Yl3 and are defined in Figure 4.2, and r:l and are shear strains at points B, D, A, C respectively /3 ; ,., 4 c 3 r 1 B D A L 1 /1i!J 2 Figure 4.2 Convention for the transverse shear deformation 43
PAGE 50
Approximating the transverse shear strain from equations (4.3) and (4.4) and using equations (4.5) and (4.6) we obtain (4.7) (4.8) Unlike the shear strains resulting from strict observance to the element shape functions these shear strains may have infinitesimally small values at two points and thus will not lock when full integration is used. We are now ready to modify the stiffness matrix by using these difference interpolation functions with linear shape functions, the curvature can be modified to be (4.9) Where Y2J and YJJ are defined in equations (4.7) and (4.8). The [Bs] matrix can be straight forwardly constructed as 44
PAGE 51
(1;) (1T/) /4 /I 0 (1T/) 2 (1;) 0 2 (1 +;) (1T/) /2 /I 0 (1T/) 2 (1;) 0 [BsY = 2 (I+;) (1 + T/) /2 /3 0 (1 + T/) 2 (1+;) 0 2 (1;) (1 + T/) /4 /3 0 (I+ T/) (4.10) 2 (1;) 0 2 Then from the shear stiffness matrix is calculated [ks] = j[Bsf[DM ][Bs]d4 A Using this new mixed interpolation function, Mindlin plate elements do not lock and because every integration is done exactly, it is certain that there are no unstable modes. It should be noted that these mixed interpolation functions may be under integrated as well as with the bilinear interpolation function but we must beware of the possible unstable modes. For example, if this new mixed interpolation function is under integrated and applied to the clampedcenter node problem shown at the end if chapter 3 also leads to unstable modes. However, we can use full integration to solve 45
PAGE 52
this problem without locking or instability and this will be demonstrated by using Fourier modes to calculate the strain energy. Numerical Examples We used the same material and dimensions presented before but with the mixed interpolation with full integration discussed above. Note that the locking behavior exhibited with full integration in the previous chapter does not appear here. "i i "i c: ::a u "i E :!1 c: Simple load 1.12 1.1 108 106 Bilinear: Selectivereduced integration 1.04 1.02 Mix interpolation:Full integration 0.98 0 50 100 150 8/t 200 250 300 a) simple supportdistribute load 46
PAGE 53
1.18 1.16 1.14 11.12 E 1.05 :II c: 1.04 1.02 0.98 Simple Supportpoint load ablt b) simple supportpoint load Clamped edgea diatributld load 1.1 108 1.05 I 1.04 l! :!: 1.02 I :II c: 0.98 0.96 0.94 ablt c) clampdistribute load Figure 4.3 NumericaV Analytical Ratio of center deflection of 1 Ox 10 element, square plate. Thin plates correspond to large ab/t. 47
PAGE 54
The Accuracy of the Mixed Interpolation Element The preceding chapter, we used strain energy calculations to analyze the accuracy of selectivereduced integration. We will use the same to analyze the mixed interpolation element discussed in this chapter. The purpose of the element introduced in this chapter is to prevent both instability and locking by adopting a different nodal interpolation function with full integration of both the bending stiffness matrix and the shear stiffness matrix. Considering Figures 4.4 a and b we can see that for small wave numbers the element produces acceptable results. There is no shear locking and no zerostrain energy mode. Strain 15 Energy Log(lb.inch) 1 0 5 0 Analytical Solution 5 10 15 0 0.5 1 1.5 2 2.5 3 Wave Number (kAx) a) Strain energy from mixed interpolation and analytical from k!lX = 0 to ;r 48 3.5
PAGE 55
Strain Energy (lb.inch) 0.05 0.04 0.03 Numerical Solution +1 0.02 0.01 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 Wave Number (kAx) b) Strain energy from mixed interpolation and analytical from kllx = 0 to w'l 0 Figure 4.4 The strain energy of a single one inch square Mindlin plate element in an infinite media, thin plate limit. 49
PAGE 56
CHAPTERS CONCLUSION A computer model was developed to simulate the bending action of Mindlin plate elements. This program was also dedicated to the study of shear locking, instability and strain energy in thin plate deformations. The selectivereduced integration with standard bilinear shape functions is introduced to eliminate shear locking successfully. This technique requires simple programming and produces a very good solution for thin plates (more than 99% accurate compared to classical analytical solutions). Unfortunately, this selectivereduced integration leads to unstable modes in the finite element solution. These unstable modes can be precluded via an appropriate boundary condition. Finally, by considering the plate as an extension of twonode beam elements to twodimensions, the different interpolation functions for transverse shear deformation can be constructed. The mixed interpolation produced even better convergence results than the selectivereduced integration method did. This mixed interpolation function does not lock and also does not contain any unstable modes. It is obtained by a slight modification of the shear stiffness matrix for an element in cartesian coordinates. However, this element is quite complicated in general geometries element compared to the previous techniques. This is because it involves tensorial shear strains in arbitrary directions. We presented two different viable element formulations, selectivereduced integration and mixed interpolation. Each has practical advantages and disadvantages. The instability from selectivereduced integration results from use of boundary conditions not usually seen in practical problems. The mixed interpolation formulation for arbitrary quadrilateral shapes can be cumbersome. We let the reader decide which technique is appropriate to their problems with optimization of computer work and the worker's time. 50
PAGE 57
APPENDIX A Stress Resultants in a Flat Plate Consider Figure A. 1 by assuming that the material is isotropic. Let M11 be the bending moment per unit length of x2 axis. We obtain [6] h 2 MII = J xJaiidxJ h 2 and let M12 be the twisting moment M12 per unit length of x2 axis, we get h 2 Ml2 = f XJU 12dx3 h 2 Similar development in x1 axis, we obtain and h 2 M22 = J X3U22dx3 h h 2 2 M21 = J X3a12dx3 = M12 h 2 note: positive direction is indicated by the righthand rule for moments. (A. I) (A.2) (A.J) (A.4) If we let Q1 and Q2 be the transverse shears per unit length in X2 axis and X1 direction. We can find equations for transverse shear forces per unit length. 51
PAGE 58
(A.5) 2 h and (A.6) 2 Resultant moments and transverse shear forces are shown in Figure A.2 Next, is the derivation of the equation for the traction force on a cross sectional surface of the element. Let N11 be the traction force of the element per unit length in the x2 direction on the middle surface. From Figure A. I we get q l I. Figure A.l Stresses that act on cross sections of a differential element of a homogeneous linearly elastic plate. q(x.,x 2 ) is the distributed force per unit area on lateral surface 52
PAGE 59
h 2 Nu = J h 2 h 2 h 2 Nl2 = J h 2 Traction forces are shown in Figure A.3 0 q Figure A.2 Resultant moments and shears 53 (A.7) (A.8) (A.9) + X2 Q2 + Mn f M22
PAGE 60
)I Figure A.3 Resultant tractions Nu )I Normal stresses a,, and Oi2 vary linearly with X3, are maximum at h/2 and zero at x3 = 0, and are associated with bending moments M11 and M22 The shear stress 0"12 also varies linearly with x3 and is associated with twisting moment Mt2 Transverse shear stresses o13 and Oi3 vary quadratically with X3. Integration of equations (A.l)(A.6) results in (A.IO) (A. II) 54
PAGE 61
The maximum magnitude of a23 and a13 is at x3 = 0 1.5Q2 (723 =t Force and Moment Balance of Flat Plates (A.12) (A.l3} (A.l4} We will start the derivation of the force and moment balance equations for a flat plate by considering the equations of motion (A.l5} (A.l6} (A.l7) Force Balance for Plates by using equation (A.15) divided by thickness h and integrating with respect to x3 [14] we obtain h _!_ f (A.l)dr3 h h 2 2 55
PAGE 62
(A. IS) A parallel development using equation (A.l6) results in (A.I9) Equations (A. IS) and (A.I9) are force balance equations for the plate element per unit thickness as shown in Figure A.3 X1 and X2 are the force per unit area of loading on top and bottom of the plate in the x1 and x2 direction, respectively. Moment Balance for Plates We now consider (A.20) 2 integrating equation (A.20) yields (A.21) Equation (A.21) is the vertical force balance for a plate element, as shown in Figure A.2 We now form h 2 J x 3 (A .15)dt3 and obtain h 2 (A.22) 56
PAGE 63
h 2 Similarly we form J x3(A.16)dx3 to obtain h 2 (A.23) Equations (A.22) and (A.23) are the balance of moments about the x1 and x 2 axes, respectively. The physical picture is shown in figure A.2. 57
PAGE 64
APPENDIX B Derivation of the FourNode Quadrilateral Mindlin Plate Element (1,1) (1,1) Figure B. I The quadrilateral element 71 coordinates (the master element) We first develop the standard, bilinear shape function on the master element (Figure B.1) which is defined in 71, coordinates. The shape functions are defined such that N; is equal to unity at node i and equal to zero at other nodes. The four shape functions can be written as (B. Ia) (B.1b) (B.1 c) 58
PAGE 65
or I M= (t;;,)(I+7777,) 4 where ;, and 71, are the coordinates of node i (B.ld) (B.Ie) We now express the displacement field within the element in terms of the nodal values. 4 w=I:N1w1 i=l 4 (}x, = LN;8xJ i=l 4 (}x, = r=l which can be written in matrix from as {u}=[N]{d} where 59 (B.2a) (B.2b) (B.2c) (B.3)
PAGE 66
WI (} :r I I (}:rl 1 {d} = (B.4) wn (} :r,n (} :r n 1 and [N' 0 0 N2 0 0 NJ 0 0 N4 0 j.J N= Nl 0 0 N2 0 0 NJ 0 0 N4 (B.S) 0 Nl 0 0 N2 0 0 NJ 0 0 Next, we express the derivatives of a function of x, in Tf coordinates, using the chain rule, (B.6) or (B.7) 60
PAGE 67
where J is the Jacobian matrix From equations (B. I) and (B.2) ,we obtain where from equation (1.27), the curvature is defined as 61
PAGE 68
by applying equation (B 7) we obtain, after inversion (B.S) (B.9) (B.IO) 62
PAGE 69
from equations (B.8), (B.9) and (8.10), we get {K} in the form 0 0 J22 J12 0 0 0 0 0 J21 1 {K}=0 0 J21 J11 J22 detJ J21 Jll 0 0 0 J22 J12 0 0 0 Now from the interpolation equation, we have cw cw a, O{;kl O{;kl =[G]{d} a, iJ(}x2 iJ(}x2 a, Brl 6r2 63 0 0 0 Jll 0 0 J12 0 0 0 0 detJ 0 detJ 0 cw cw a, iJ8 .... iJ(} .... (B. II) a, iJ(} ... 2 iJ(} ... 2 a, e .... () ... 2 (B.12)
PAGE 70
where (1'7) 0 0 0 0 0 0 0 0 (1'7) (1;) 0 0 0 0 0 0 0 (1'7) 0 (1'7) (1'7) 0 0 0 0 0 0 0 0 (1'7) (1 +;} 0 0 0 (Gf = .!_ 0 0 0 0 (1'7) 0 4 (1 + '7) 0 0 0 0 0 0 0 0 (1 + '7) 0 0 (1 + '7) 0 0 0 0 0 (1 + '7) 0 '7) (1 + '7) 0 0 0 0 0 0 0 0 (1 + '7) 0 0 (1'7) 0 0 0 0 0 (1 + '7) 0 (1'7) equations (B. II) and (B.12) yields [B] = [A][G] where 0 0 122 112 0 0 0 0 0 0 0 0 121 111 0 0 [A]1 0 0 121 111 122 112 0 0 det1 121 111 0 0 0 0 0 det1 122 112 0 0 0 0 det1 0 Finally we can calculate the element stiffness from [k] = j[Bf[DM ][B]dA A 64
PAGE 71
Numerical Integration Consider the problem of numerically evaluating a onedimensional integral of the form [11] I I= f f(f)d; I The Gaussian quadrature approach for evaluating I with npoint approximation is given below where and I I= J f(f)d; +w2J(;2)+ ...... +wnJ(;n) I w 1 w 2 ..... ,wn are the weights ;1, ;2, ..... ,;n are Gauss points (B.13) From equation (B.13), npoint Guassian quadrature provides an exact answer for polynomial of order (2n1) or less. Two Dimensional Integrals We can write the equation to evaluate the twodimensional integrals of the form I 1 I= f f J(;, q)dgiq (B.14) 11 Now, from equation (B.13) n n LLw;w1J(;;,T]1 ) (B.15) i=lj=l 65
PAGE 72
Stiffness Integration Introduction of [B] into the expression for the stiffness matrix, we may now numerically evaluate the element stiffness matrix [II ] 172 = 0.57735 0 0 __.. 1]1 = 0.57735 0 0 ;I= 0.57735 ;2 = 0.57735 Figure B.2 twodimensional Gaussian quadrature with 2x2 rule. I I [k]= J f[Bf[DM][B]detJd;dq (B.I6) 11 where [B] and detJ are the functions of ; and 7]. Let J(;,q) = [Bf[D][B]detJ (B.I7) Then, If we use a 2point Gaussian rule in each direction, we get 66
PAGE 73
where [k .] = WJ2 rp(;l, 771) +WI w2;(;1, 772) + Wz wlrp(;2, 771) + w; ;(;2, 772) (B.18) ;I = 711 = YJJ ;2 = T/2 = y../3 and the Gauss points are shown above diagrammatically 67
PAGE 74
REFERENCES 1. T. J. L. Hughes, R. L. Taylor, and Kanoknukulchai, "A Simple and Efficient Element for Plate Bending ," International Journal for Numerical Methods in Engineering, 11, no. 10(1977), 15291543. 2. R. H. Macneal, "A Simple Quadrilateral Shell Element," Computers and Structures, 8(1978), 175183. 3. T. J. R. Hughes, "Generalization of Selective Integration Procedures to Anisotropic and Nonlinear Media, International Journal of Numerical Methods in Engineering, Vol. 15, 1980, 14131418 4. T. J. R. Hughes and T. E. Tezduyar, "Finite Elements Based Upon Mindlin Plate Theory With Particular Reference to the Fournode Bilinear Isoparametric Element," Journal of Applied Mechanics, (September 1981 ), 587596. 5. K. C. Park and D. L. Flaggs, "A Fourier Analysis of Spurious Mechanisms and Locking in The Finite Element Method," Computer Methods in Applied Mechanics and Engineering, 46 ( 1984 ), 6581. 6. A. B. Boresi, 0. M. Sidebottom, F. B. Seeley and J. Smith, Advanced Mechanics of Materials, 3rd ed., New York: John Wiley & Son Inc., 1985. 7. A. C. Ugural and S. K. Fenster, Advanced Strength and Applied Elasticity, New Jersey: PrenticeHall Inc., 1987. 8. T. J. R. Hughes, The Finite Element Method: Linear Static and Dynamic Finite element analysis, New Jersey: PrenticeHall, 1987. 9. R. D. Cook, D. S. Malkus and M. E. Plesha, Concepts and Applications of Finite Element Analysis, 3rd., New York: John Wiley, 1989. 10. J. E. Shigley and C. R. Mischke, Mechanical Engineering Design, 5th ed., New York: McGrawHill, 1989. 68
PAGE 75
11. T. R. Chandrupatla and A. D. Belegundu, Introduction to Finite Elements in Engineering, New Jersey: PrenticeHall Inc., 1991 12. R. D. Cook, Finite Element Modeling for stress Analysis, New York: John Wiley & Son Inc., 1995. 13. K. J. Bathe, Finite Element Procedures, New Jersey: PrenticeHall Inc., 1996 14. J. A. Trapp, Course notes, ME 5153., 1996 15. I. H. Shames and C. L. Dym, Energy and Finite Element Methods in Structural Mechanics, USA: Hemisphere Publishing Corporation., 1985. 69
