FINITE ELEMENT ANALYSIS DEPARTMENT OF

MECHANICAL
ENGINEERING

SUBJECT: FEM
TOPIC : PLATE AND
SHELL ANAYSIS

DR. PINANK A PATEL FACULTY:

DR.PINANK A PATEL
INTRODUCTION
When a flat plate is subjected to both in-plane and transverse or normal loads as shown in
Figure any point inside the plate can have displacement components u, v, and w parallel to
x, y, and z axes, respectively. In the small deflection (or linear) theory of thin plates, the
transverse deflection w is uncoupled from the in-plane deflections u and v. Consequently,
the stiffness matrices for the in-plane and transverse deflections are also uncoupled and
they can be calculated independently. Thus, if a plate is subjected to in-plane loads only, it
will undergo deformation in its plane only. In this case, the plate is said to be under the
action of “membrane” forces. Similarly, if the plate is subjected to transverse loads (and/or
bending moments), any point inside the plate experiences essentially a lateral
displacement w (in-plane displacements u and v are also experienced because of the
rotation of the plate element). In this case, the plate is said to be under the action of
bending forces. If the plate elements are used for the analysis of three-dimensional
structures, such as folded plate structures, both in-plane and bending actions have to be
considered in the development of element properties. This
TRIANGULAR MEMBRANE ELEMENT

The triangular membrane element is considered to lie in the xy plane of a local xy

coordinate system as shown in Figure (b). By assuming a linear displacement variation
inside the element, the displacement model can be expressed as
u(x, y) = α1 + α2x+ α3y
v(x, y) = α4 + α5x+ α6y (1)
By considering the displacements ui and vi as the local degrees of freedom of node i(i =
1, 2, 3), the constants α1,…, α6 can be evaluated. Thus, by using the conditions
u(x, y) = u1 = q1 and v(x, y) = v1 = q2 at (x1, y1)
u(x, y) = u2 = q3 and v(x, y) = v2 = q4 at (x2, y2)
u(x, y) = u3 = q5 and v(x, y) = v3 = q6 at (x3, y3)
we can express the constants α1,…, α6 in terms of the nodal degrees of
freedom
This leads to the displacement model:
 (15)

(16)

where V(e) denotes the volume of the element. If the plate thickness is taken as
a constant (t), the evaluation of the integral in Eq. (16) presents no difficulty
since the elements of the matrices [B] and [D] are all constants (not functions of
x and y). Hence, Eq. (10.16) can be rewritten as

(17)
Although the matrix products involved in Eq. (17) can be performed
conveniently on a computer, the explicit form of the stiffness matrix is given
below for convenience:
(18)
 TRANSFORMATION MATRIX

To generate the transformation matrix [λ], the direction cosines of lines 0x and 0y
with respect to the global X, Y, and Z axes are required. Since the direction cosines of
the line 0y are the same as those of line ij, we obtain

where the distance between the points i and j(dij) is given

by

Then the direction cosines of the line 0x will be the same as those of the line pk :
where dpk is the distance between the points p and k. The coordinates (Xp, Yp, Zp) of
the point p in the global coordinate system can be computed as

where dip is the distance between the points i and p. To find the distance dip, we use
the condition that the lines ij and pk are perpendicular to each other:
lijlpk +mijmpk + nijnpk = 0

The transformation matrix [λ] can now be constructed by using the direction cosines of
lines ij and pk as
Finally, the stiffness matrix of the element in the global XYZ coordinate system can be
computed as

CONSISTENT LOAD VECTOR

The total load vector in the local coordinate system is thus given by
This load vector, when referred to the global system, becomes

CHARACTERISTICS OF THE ELEMENT

1. The CST element was the first finite element developed for the analysis of plane stress
problems [1.7]. Because the displacement model is linear (Eq.1), the element is called a
linear triangular element. From Eqs. (10.11) and (10.12), we find that the [B] matrix is
independent of the position within the element and hence the strains are constant throughout
the element. This is the reason why this element is often referred to as a CST element
(constant strain triangular element).
2. The displacement model chosen (Eq. 10.1) guarantees continuity of displacements with
adjacent elements because the displacements vary linearly along any side of the triangle.
3. From Eq. (10.13), we can notice that the stresses are also constant inside an element.
Hence, the element is also called a CST (constant stress triangular) element. Since the
stresses are independent of x and y, the equilibrium equations are identically satisfied inside
the element since there are no body forces.
4. If the complete plate structure being analyzed lies in a single (e.g., XY) plane, the vector
Q will also contain six components. In such a case, the matrices [λ] and [K(e)] will each be
of the order 6 × 6.
5. In problems involving bending, this element overestimates the bending stiffness (normal
stresses). More accurate normal stresses can be obtained by using smaller size elements.
However, the convergence to the correct solution will be very slow. To illustrate the
numerical performance of the element, a uniform beam with rectangular cross section
subjected to nodal forces at the free end.
The nodal forces indicated will produce bending in the beam. If the beam is modeled
using CST elements as indicated in Figure, each element gives a constant value of σx
in the depth or y direction instead of a linear variation predicted by the exact solution.
In fact, the exact solution along the x axis would be σx = 0 while the CST model
would predict a constant value of σx with alternate signs as we move along the x axis
from one element to the next as shown in Figure. In addition, the CST element
predicts a spurious shear stress. The element predicts a constant, nonzero, value of the
shear stress σxy in each element. This means that the shear strain εxy is constant in the
beam while it should be zero. Also, in some cases, the CST element can exhibit a
phenomenon known as locking. The term locking denotes excessive stiffness in one or
more deformation modes. The numerical results obtained with the CST element are
given in Table 1. The results indicate that even 512 elements could not predict the
bending behavior of the beam (deflection and stress) very accurately.
 EXAMPLE 01
A triangular membrane element of thickness t = 0.1 cm, with the (x, y) coordinates of
nodes indicated besides the node numbers, is shown in Figure. If the material of the
element is steel with Young’s modulus E = 207 GPa and Poisson ratio v = 0.3, determine
the following:
1. Shape functions of the element, Ni(x, y), Nj(x, y), and Nk(x, y).
2. Matrix [B] that relates the strains to the nodal displacements.
3. Elasticity matrix [D].
4. Element stiffness matrix, [K(e)].
EXAMPLE 2
If the element described in Example 1 undergoes a temperature increase of 80°C and the
coefficient of thermal expansion of the material is 10.8 × 10−6 m/m-°C, find the nodal load
vector due to thermal loading.
Solution
The load vector due to initial strain (thermal loading) as
EXAMPLE 3
If distributed body forces of magnitude ϕ0x = 20 N/m2 and ϕ0y = 30 N/m2 act along x
and y directions throughout the element described in Example1, determine the
corresponding nodal load vector of the element.
 EXAMPLE 4
For the triangular membrane element considered in Example1, distributed surface tractions
of magnitude Φx = p0x = 50 N/m2 and Φy = p0y = –30 N/m2 act on the edge (face)
connecting the nodes 1 and 3. Find the resulting nodal load vector of the element.
Solution
Using an equation similar to Eq. (42), we obtain the load vector of the element due to the
specified distributed surface tractions as

where S13 is the surface area of the edge (face) 13 = t d13 with d13 denoting the distance between the nodes 1
and 3:
EXAMPLE 5
A concentrated load, with components P0x = 1000 N and P0y = –500 N, acts at the point
(x0, y0) =(3, 5) cm of the plate described in Example 1. Determine the corresponding
nodal load vector of the element.
The equivalent nodal load vector of the element corresponding to the point load,
can be expressed as

At the given point (x0 = 3 cm, y0 = 5 cm), the shape functions, given by Eqs. (E.1) to (E.3) ,
assume the values
NUMERICAL RESULTS WITH MEMBRANE ELEMENT
A Plate under Tension
The uniform plate under tension, shown in Figure (a), is analyzed by using the CST
elements. Due to symmetry of geometry and loading, only a quadrant is considered for
analysis. The finite element modeling is done with eight triangular elements as shown in
Figure (b). The total number of nodes is nine and the displacement unknowns are 18.
However, the x components of displacement of nodes 3, 4, and 5 (namely Q5, Q7, and Q9)
and the y components of displacement of nodes 5, 6, and 7 (namely Q10, Q12, and Q14) are
set equal to zero for maintaining symmetry conditions. After solving the equilibrium
equations, the global displacement components can be obtained as
 COMPUTATION OF STRESSES
For finding the stresses inside any element “e,” shown in Figure (c), the following
procedure can be adopted:
Step 1: Convert the global displacements of the nodes of element e into local
displacements as

where

and [λ] is the transformation matrix of the element given

by
 Here (lpk, mpk) and (lij, mij) denote the direction cosines of lines pk (x axis) and ij (y axis)
with respect to the global (X, Y) system.

Step 2: Using the local displacement vector q of element e, find the stresses inside the
element in the local system by using Eqs. (13) and (11) as

where [D] and [B] are given by Eqs. 15 and 12 respectively.

Step 3: Convert the local stresses σxx, σyy, and σxy of the element into global stresses
σXX, σYY, and σXY by using the stress transformation relations
 A Plate with a Circular Hole
The performance of membrane elements for problems of stress concentration due to
geometry is studied by considering a tension plate with a circular hole (Figure). Due to
the symmetry of geometry and loading, only a quadrant was analyzed using four
different finite element idealizations as shown in Figure 8. The results are shown in
Table 3. The results indicate that the stress concentration is predicted to be smaller
than the exact value consistently.
A Cantilevered Box Beam
The cantilevered box beam shown in Figure 10.9 is analyzed by using CST elements.
The finite element idealization consists of 24 nodes, 72 degrees of freedom (in global
XYZ system), and 40 elements as shown in Figure 10. The displacement results
obtained for two different load conditions are compared with those given by simple
beam theory in Table 4. It can be seen that the finite element results compare well with
those of simple beam theory.
QUADRATIC TRIANGLE ELEMENT

A triangle element with a quadratic displacement model has six nodes—three at the vertices
and three at the mid points of the sides—as shown in Figure 11. The displacement
components of a point in the element parallel to the x and y axes are assumed as

where the constants or generalized degrees of freedom αi, i = 1, 2, …, 12 can be

expressed in terms of the vector of nodal displacements of the element
The strains in the element are given by

FIGURE 10.7
Plate with a Circular Hole under
Uniaxial Tension
FIGURE 10.8
Finite Element Idealization of the Plate with a Circular Hole
RECTANGULAR PLATE ELEMENT (IN-PLANE FORCES)
Consider a rectangular plate undergoing in plane displacements due to in plane
forces as shown in Figure 12. The variations of displacements inside the element
are assumed as

It can be seen that although the displacement distribution is represented by a second

degree surface, the displacement u(x, y), for example, varies linearly along the x (or y)
direction for any constant value of y (or x) as shown in Figure 10.13. By using the nodal
values, u(x1 = −a, y1 = −b) = u1, u(x2 = a, y2 = −b) = u2,…, v(x4 = −a, y4 = b) = v4, Eqs.
(10.54) can be exprese as
The global stiffness matrix of the element ½K(e) in three-dimensional space can be generated as
and [0] is a zero matrix of size 2 × 3. In above Eq., lpq, mpq, npq and lps, mps, nps
denote, respectively, the direction cosines of the lines pq (x-axis) and ps (y-axis).
1. The rectangular element described in this section is called the Q4 element. It can
be observed from the assumed displacement field that the strain εxx is constant in
the x direction and varies linearly in the y direction. Similarly, the strain εyy is
constant in the y direction, but varies linearly in the x direction. On the other hand,
the shear strain εxy varies linearly in both x and y directions.
2. If the element undergoes a constant temperature change, T, the stresses in the
element are given by

Each stress component varies linearly with x and y that, in general, violates the stress
equilibrium equations within the element.
3. This element, as in the case of CST element, cannot exhibit pure bending. In
problems involving bending, this Q4 element displays not only the expected bending
(normal) strain, but also spurious shear strain. Thus, the element exhibits shear locking
behavior under bending deformation.
 The term "nonconforming elements" is used when the trial functions do not form.
a subspace of the corresponding energy space. Nonconforming elements are.
reported to have been successfully used in some cases in the displacement analysis.
of plate bending

When the finite element space is a subspace of the solution space, the method is
called CONFORMING. ... It means that the shape function in this conforming
finite element space is continuous together with its m−1 order derivatives.

The basic assumptions employed in degenerated shell elements are similar to

these used for Mindlin plate elements. Ahmad degenerated a three dimensional
brick element to a general curved shell element which has nodes only at the mid-
surface.
The most popular procedures are uniform or selectively reduced integration. ...
This may result in undesired internal mechanism which produces hour-glass modes.
Even if these zero energy modes do not exist at the beginning, they may show up in a
later deformed stage.
Shear locking, caused by the so-called “parasitic shear” in the bending of a thin
shell, is a typical locking phenomenon. Due to the incompetence of the shell
elements in capturing the bending deformation appropriately, the strain energy is
absorbed by the shear mode erroneously.
Shear locking is an error that occurs in finite element analysis due to the linear nature
of quadrilateral elements. The linear elements do not accurately model the curvature
present in the actual material under bending, and a shear stress is introduced. The
additional shear stress in the element (which does not occur in the actual beam) causes
the element to reach equilibrium with smaller displacements, i.e., it makes the element
appear to be stiffer than it actually is and gives bending displacements smaller than
they should be.
In areas where linear elements are loaded by in plane bending, shear locking is
prevented by using preferably 3 elements over the height. This is illustrated in the
example below. In the area of interest ensure that the elements are as rectangular as
possible (preferably square), to give the most accurate results.

Hour glass phenomenon : Consider a 4 noded Rectangular shell element called as S4

element. They may have any number of Integration points (depending upon full or
reduced integration). Stresses and strains at this element are captured by integrating
the values got at these integration points. It is to be remembered that values are more
accurate at the exact location of the Integration point and not so accurate as we move
away from it. For sake of lesser computational time, we can use S4R element which
means Reduced Integration where in the one Integration point is used and it is ideally
at the centre of the element.
In the case of these Reduced Integration elements, while post processing results, we
can see a false or fake deformation mode (zig zag lines). Such a thing is the Hour
Glass phenomenon
The above is an example of a Single Integration point shell element. The values
captured at the single integration point at the centre remains same as the dotted lines
have not changed in magnitude or in angle. However, the element as a whole has shown
a change in configuration. This is because elements with lesser integration points are
less stiffer and as there is no stiffness in this mode, these elements have shown this kind
of change in Configuration.
In actual physical scenario or in Experiments, this does not happen. But in FEA this
happens. This is just a non physical zero energy mode which leads to excessive element
distortions. The more structured and less distorted your elements are, the better the results.
To solve this problem,
1. Mesh Refining is a solution. Hourglassing is more prone in coarse mesh.
2. Fully Integrated elements are better though computational time increases vastly and fully
integrated elements are prone to shear locking.
3. Providing a small arbitrary stiffness in the mode in which hourglass deforms. This is
only a FEA point of view and in actual physical system, there is no such stiffness. This is
what is called as Built in hourglass controls.
4. Viscous damping.
h-method: This method increases the number of elements and hence
decreases the element size while keeping the polynomial order of the constant.

p-method: This method increases the polynomial order of the interpolation

function while keeping the number of elements constant.

Aspect ratio is the measure of a mesh element's deviation from having all
sides of equal length. A high aspect ratio occurs with long, thin elements.
Entering an overly large value for the Minimum Element Size mesh control
may cause the mesh generator to create solid elements with high aspect
ratios. It is the ratio of length to thickness of the element.