Topics To Be Discussed: Quadrilateral Elements

Quadrilateral elements
Asst.Prof.Dr. A. Ahmed
(awais.ahmed@uetpeshawar.edu.pk)
Department of Civil Engineering, University of Engineering & Technology Peshawar, Pakistan
Asst.Prof.Dr. A. Ahmed Finite element method 1
Topics to be discussed
1 Local coordinates/Natural coordinates
2 Isoparametric elements
3 Four Node Quadrilateral elements
4 Derivates, strain, N and B matrices

Completeness of trial solution
1 Is the following a complete linear trial function ??
θ(x, y) = a1 + a2 x + a3 y
2 Is the following a complete quadratic trial function ??
θ(x, y) = a1 + a2 x + a3 y + a4 x2 y 2 + a5 xy + a6 y 3
⇒ The answer can be determined by examining Pascal’s triangle

⇒ Trial functions of higher degree can be systematically developed
with the help of Pascal’s triangle

Figure: Pascal’s triangle
Pascal’s triangle can help to find out a complete polynomial

expansion for trial solution


Figure: Pascal’s triangle
Each row of the triangle gives the monomials that must be

included in a finite element approximation to provide an
element with the required order of completeness
Higher order element shape fucntions

Exercise: Steps required for the development of shape functions
◮ Choose the interpolation function

X
n
θe (x) = aei pi = p(x)ae
i=1
◮ Evaluate interpolation functions at nodal points
θ̂e (x) = Me (x)ae
◮ Interpolation function is then expressed as
θe (x) = p(Me )−1 θ̂e

e
θ (x) = Ne θ̂e
in which Ne = p(Me )−1 is the shape function matrix

Coordinate systems
A coordinate system defined for the entire body or structure is
called the Global coordinate system
A coordinate system that is defined for a particular element and

not necessarily for the entire body is called the Local coordinate
system
A local coordinate system which permits identification /

referencing of a point within an element in terms of
dimensionless numbers, whose magnitude never exceeds unity is
called Natural coordinate system
Coordinate systems

Coordinate systems

It is possible to establish a mapping between Global coordinate

system and parent element natural coordinate system
1 1
x= (a + b) + ξ(b − a)
2 2
re-arranging the terms
1−ξ 1+ξ
x = a +b
2 2
= x̂1 N1 (ξ) + x̂2 N2 (ξ) = N x̂
Coordinate systems

Iso-parametric Finite elements
Iso-parametric formulation
Finite elements are usually developed with shape functions
expressed in terms of parent element coordinates, often called
as element coordinates or natural coordinates
Examples of element coordinates are triangular coordinates and

isoparamteric coordinates
We will next describe the use of shape functions expressed in

terms of element coordinates

What is Iso-parametric element
Elements can utilizes different parametrization for geometry
and displacements
Elements for which geometry and displacements are given by

the same parametric representation are called as iso-parametric
elements
◮ Development of iso-parametric elements: Key idea is to use the
same shape functions to represent both the geometry and the
problem unknown, which in structural mechanics are
displacements. Hence the name isoparametric element (“iso”
means equal)
◮ This property may be generalized to arbitrary elements. Thus,

for triangular element, replacing the term ”triangular
coordinates” by the more general one ”natural coordinates”
Elements which utilizes higher order shape function for

displacement field compared to geometry are called
sub-parametric
What is Iso-parametric element

There are also super-parametric elements whose geometry is
more refined than the displacement expansion

General Iso-parametric formulation
Element geometry is interpolated in terms of shape functions as
n
X
x = x̂I NIe (ξ)
I=1
n
X
y = ŷI NIe (ξ)
I=1
n
X
with NIe = 1
I
Displacements are interpolated as
3
X
ux = ûxI NI (ξ)
I=1
3
X
uy = ûyI NI (ξ)
i=1
Derivatives of functions
The spatial derivatives of displacement field is obtained by
implicit differentiation. By chain rule we have
∂ ∂ ∂ξ ∂ ∂η ∂ ∂ζ
= + +
∂x ∂ξ ∂x ∂η ∂x ∂ζ ∂x
∂ ∂ ∂ξ ∂ ∂η ∂ ∂ζ
= + +
∂y ∂ξ ∂y ∂η ∂y ∂ζ ∂y
∂ ∂ ∂ξ ∂ ∂η ∂ ∂ζ
= + +
∂z ∂ξ ∂z ∂η ∂z ∂ζ ∂z
In matrix form
 ∂   ∂ξ ∂η ∂ζ   ∂ 
 ∂x   ∂x ∂x ∂x   ∂ξ 
    
    
 ∂   ∂ξ ∂η ∂ζ   ∂  
 =  
 ∂y   ∂y ∂y ∂y   ∂η 
    
     
∂ ∂ξ ∂η ∂ζ  ∂ 
   
∂z ∂z ∂z ∂z ∂ζ
Requires evaluation of ∂ξi /∂xj
∂ξi /∂xj are difficult to evaluate by direct method. However
using the chain rule we write
∂ ∂x ∂y ∂z  ∂ 
   
 ∂ξ   ∂ξ ∂ξ ∂ξ  ∂x 
   
   


 ∂   ∂x ∂y ∂z   ∂ 
 
 =
   
 ∂η   ∂η ∂η ∂η   ∂y 

   
 
    
 ∂   ∂x ∂y ∂z  ∂

∂ζ ∂ζ ∂ζ ∂ζ ∂z

∂ ∂ ξ ∂xj
= Fξ in which Fji =
∂ξ ∂x ∂ξi
In indicial form
∂ ∂xj ∂
=
∂ξi ∂ξi ∂xj
Fξ is Jacobian matrix relating the natural coordinate derivatives
to the global coordinates derivatives
Jacobian matrix can be found by taking derivatives of

coordinate interpolation relations
To compute the derivatives in global coordinates we invert the

equation
∂ ξ −1 ∂
= F
∂x ∂ξ
Thus the computation of the derivatives with respect to x

requires inverse of the Jacobian between current and parent
element coordinates
Derivatives of shape functions are

Ni,x = [Fξ ]−1 Ni,ξ
Area/volume transformation
A differential area and volume in global coordinates is related
to the differential area and volume in natural coordinates by
following relations
dx dy = det(Fξ ) dξ dη = J dξ dη
dx dy dz = det(Fξ ) dξ dη dζ = J dξ dη dζ
in which J = det(Fξ )
Integration and nodal forces

Integrals on the current configuration are related to integrals
over the parent domain by
Z Z
g(x)dΩ = g(ξ)Jd∆
Ωe ∆
When internal forces are computed by integration over the

parent domain, it is computed as
∂NI ∂NI
Z Z
fint = σji dΩ = σji Jd∆
Ωe ∂x j ∆ ∂xj

Iso-parametric Quadrilateral element
Isoparametric Q4 element - Geometry

Consider the quadrilateral element as shown

Spatial coordinates are expressed as
x = x̂I NI
y = ŷI NI
 
N1
x x1 x2 x3 x4 
N2 

=
y y1 y2 y3 y4  N 3 
N4
Isoparametric Q4 element - Displacement
Displacement field is approximated as
ux = ûxI NI
uy = ûyI NI
 
  N1
ûx1 ûx2 ûx3 ûx4  
ux   N2 
=
uy  N3 
ûy1 ûy2 ûy3 ûy4
N4
or  
ûx1
ûy1 
 
ûx2 
 

ux N1 0 N2 0 N3 0 N4 0  ûy2 

=
uy 0 N1 0 N2 0 N3 0 N4  û
 x3 

ûy3 
 
ûx4 
ûy4
Isoparametric Q4 element - Shape functions

Assume a linear displacement approximation is
ux = a0 + a1 ξ + a2 η + a3 ξη
uy = a4 + a5 ξ + a6 η + a7 ξη
in matrix form
 
a0
a1 
 
a2 
 

ux 1 ξ η ξη 0 0 0 0 a3 

= or u = pa
uy 0 0 0 0 1 ξ η ξη a4 


a5 
 
a6 
a7

Imposition of nodal conditions yield
    
ûx1 1 −1 −1 1 0 0 0 0 a0
ûy1  0 0 0 0 1 −1 −1 1  a1 
 
  
ûx2  1 1 −1 −1 0 0 0 0 a2 
 
  
ûy2  0 0 0 0 1 1 −1 −1 a3 
 
 = or û = Ma
ûx3  1 1 1 1 0 0 0 0 a4 
 
  
ûy3  0 0 0 0 1 1 1 1 a5 
 
  
ûx4  1 −1 1 −1 0 0 0 0  a6 
ûy4 0 0 0 0 1 −1 1 −1 a7
This implies
a = M−1 û =⇒ u = pM−1 û =⇒ u = Nû

in matrix form
 
ûx1
ûy1 
 
ûx2 
 

ux N1 0 N2 0 N3 0 N4 0  ûy2 

=
uy 0 N1 0 N2 0 N3 0 N4 ûx3 


ûy3 
 
ûx4 
ûy4
in which
1 1
N1 = (1 − ξ)(1 − η), N2 = (1 + ξ)(1 − η)
4 4
1 1
N3 = (1 + ξ)(1 + η), N4 = (1 − ξ)(1 + η)
4 4

N1 N2
1 1
0.5 0.5
0 0
−1 1 −1 1
−0.5 0.5 −0.5 0.5
0 0 0 0
0.5 −0.5 0.5 −0.5
1 −1 1 −1
N3 N4
1 1
0.5 0.5
0 0
−1 1 −1 1
−0.5 0.5 −0.5 0.5
0 0 0 0
0.5 −0.5 0.5 −0.5
1 −1 1 −1

In a more compact form, shape functions are written as
1
NI = (1 + ξI ξ)(1 + ηI η)
4
Node I ξI ηI
1 -1 -1
2 1 -1
3 1 1
4 -1 1

Isoparametric Q4 element - Derivatives
It is impossible to write explicit expressions for the element
coordinates in x and y, and the derivatives of the shape
functions are evaluated by using implicit differentiation

T
ξ,x ξ,y
= NTI,ξ Fξ
−1
NI,x = [NI,x NI,y ] = NI,ξ NI,η
η,x η,y
Jacobian of the current configuration with respect to the

element coordinates is given by
ξ ∂xj
Fξ = Fji =
∂ξi

x,ξ x,η
=
y,ξ y,η
∂NI

= xiI , ∵ x = x I NI
∂ξj

xI x N xI NI,η
NI,ξ NI,η = I I,ξ

=
yI yI NI,ξ yI NI,η

For the four node quadrilateral element
4
ξ
X xI ξI (1 + ηI η) xI ηI (1 + ξI ξ)
F =
yI ξI (1 + ηI η) yI ηI (1 + ξI ξ)
I=1
or more explicitly

ξ1 (1 + η1 η)
 
η1 (1 + ξ1 ξ)

1 ξ (1 + η η) 1 η (1 + ξ ξ)
x,ξ = [x1 x2 x3 x4 ]  2
 2 
, x,η = [x1 x2 x3 x4 ]  2
 2 
4 ξ3 (1 + η3 η) 4 η3 (1 + ξ3 ξ)
ξ4 (1 + η4 η) η4 (1 + ξ4 ξ)

ξ1 (1 + η1 η)
 
η1 (1 + ξ1 ξ)

1 ξ2 (1 + η2 η) 1 η2 (1 + ξ2 ξ)
y,ξ = [y1 y2 y3 y4 ] 
ξ (1 + η η) , y,η = [y1 y2 y3 y4 ] 
η (1 + ξ ξ)
4 3 3 4 3 3
ξ4 (1 + η4 η) η4 (1 + ξ4 ξ)

The inverse of Fξ is given by

1 y,η −x,η
Fξ
−1
=
J −y,ξ x,ξ
in which J = x,ξ y,η − x,η y,ξ

The gradients of shape functions for the 4-node quadrilateral
with respect to the element coordinates are
 
∂N1 /∂ξ ∂N1 /∂η
∂N2 /∂ξ ∂N2 /∂η 
NT,ξ = [∂NI /∂ξi ] = 
∂N3 /∂ξ ∂N3 /∂η 

∂N4 /∂ξ ∂N4 /∂η

 
ξ1 (1 + η1 η) η1 (1 + ξ1 ξ)
1ξ2 (1 + η2 η) η2 (1 + ξ2 ξ)

=
4 3
 ξ (1 + η 3 η) η 3 (1 + ξ 3 ξ) 
ξ4 (1 + η4 η) η4 (1 + ξ4 ξ)

The gradients of shape functions for the 4-node quadrilateral
with respect to the spatial coordinates are
NTI,x = NTI,ξ Fξ
−1
BI =
 
ξ1 (1 + η1 η) η1 (1 + ξ1 ξ)
1ξ2 (1 + η2 η) η2 (1 + ξ2 ξ) 1 y,η −y,ξ
=
4 ξ3 (1 + η3 η) η3 (1 + ξ3 ξ) J −x,η x,ξ
ξ4 (1 + η4 η) η4 (1 + ξ4 ξ)
Displacement gradients are
u,x = uI BI = uI NTI,x
Higher order elements

Higher order elements - Quadrilateral element
Higher order isoparametric elements provide one of the most
attractive features of finite elements, the ability to model
curved boundaries

It is possible to use classical method to derive shape functions
for higher order elements
u = pa =⇒ a = M−1 û =⇒ u = pM−1 û =⇒ u = Nû
However, closed form symbolic inversion is very cumbersome
An alternate method is to construct shape functions by the

Tensor product method
The method is based on

◮ taking products of lower dimensional shape functions and
◮ exploiting the Kronecker delta property of shape functions
For example two-dimensional shape functions for a four node

quadrilateral element are obtained as the product of the
one-dimensional shape functions


N[I,J] (ξ, η) = NI (ξ)NJ (η)
1 1
N1 = N[1,1] = N1 ∗ N1 = (1 − ξ)(1 − η) N2 = N[2,1] = N2 ∗ N1 = (1 + ξ)(1 − η)
4 4
1 1
N3 = N[2,2] = N2 ∗ N2 = (1 + ξ)(1 + η) N4 = N[1,2] = N1 ∗ N2 = (1 − ξ)(1 + η)
4 4
Higher order elements - Q9

One-dimensional quadratic element shape functions are

1
N1 = − ξ(1 − ξ)
2
N2 = (1 − ξ 2 )
1
N3 = ξ(1 + ξ)
2

1
N1 = N1 (ξ) ∗ N1 (η) = ξη(1 − ξ)(1 − η)
4
1
N2 = N3 (ξ) ∗ N1 (η) = − ξη(1 + ξ)(1 − η)
4
1
N3 = N3 (ξ) ∗ N3 (η) = ξη(1 + ξ)(1 + η)
4
1
N4 = N1 (ξ) ∗ N3 (η) = − ξη(1 − ξ)(1 + η)
4
1
N5 = N2 (ξ) ∗ N1 (η) = − η(1 − ξ 2 )(1 − η)
2
1
N6 = N3 (ξ) ∗ N2 (η) = ξ(1 + ξ)(1 − η 2 )
2
1
N7 = N2 (ξ) ∗ N3 (η) = η(1 − ξ 2 )(1 + η)
2
1
N8 = N1 (ξ) ∗ N2 (η) = − ξ(1 − ξ)(1 − η 2 )
2
N9 = N2 (ξ) ∗ N2 (η) = (1 − ξ 2 )(1 − η 2 )

N1 N5
1 1
0.5 0.5
0 0
−0.5 −0.5
−1 −1
1 1
0 0.5 0 0.5
0 0
−0.5 −0.5
1 −1 1 −1
N9
0.8
0.6
0.4
0.2
0
−1
1
0 0.5
0
−0.5
1 −1


Exercise: Find out the shape functions of the Q12 element
Higher order elements

Isoparametric elements in two (or three) dimensions
constructed by a tensor product of one-dimensional element
shape functions are called Lagrange elements.
Some Lagrange elements possess internal nodes that do not
contribute to the interelement compatibility
These nodes can be condensed out at the element level to
reduce size of global matrices
Commercial softwares usually employs the formulation of higher
order element without internal nodes; these are called
serendipity elements
Shape functions for serendipity elements cannot be constructed
by a tensor product of one-dimensional shape functions
Serendipity element shape functions are obtained by a tensor
product of carefully selected functions to satisfy the Kronecker
delta property of the shape functions

Higher order elements - Q8 Serendipity element

Develop shape functions for the eight-noded quadrilateral

element (Q8)
Higher order elements - Q8 Serendipity element

Step 1

Lets take Lagrangian interpolation of a quadratic x linear for

the mid-side node

Higher order elements- Q8 Serendipity element
Step 2

For corner node, start with bi-linear Lagrangian family

Successive subtract of fraction of mid-side node shape function
to ensure zero at these nodes
Three-dimensional elements
The parent domain of the eight-node hexahedral element is a
biunit cube with element coordinates ξ, η and ζ

Three-dimensional elements
Geometry is approximated as
x = N8H (ξ, η, ζ)x̂

y = N8H (ξ, η, ζ)ŷ
z = N8H (ξ, η, ζ)ẑ
The shape functions can be constructed by tensor-product of

one-dimensional linear shape functions
NL8H (ξ, η, ζ) = NI2L (ξ) NJ2L (η) NK

2L
(ζ)
For isoparametric element
ux = N8H (ξ, η, ζ)ûx

uy = N8H (ξ, η, ζ)ûy
uz = N8H (ξ, η, ζ)ûz

Topics To Be Discussed: Quadrilateral Elements

Uploaded by

Copyright:

Available Formats

Topics To Be Discussed: Quadrilateral Elements

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Topics To Be Discussed: Quadrilateral Elements

Uploaded by

Copyright:

Available Formats

Quadrilateral elements

Department of Civil Engineering, University of Engineering & Technology Peshawar, Pakistan

Asst.Prof.Dr. A. Ahmed Finite element method 1

Asst.Prof.Dr. A. Ahmed Finite element method 2

2 Is the following a complete quadratic trial function ??

⇒ The answer can be determined by examining Pascal’s triangle

Asst.Prof.Dr. A. Ahmed Finite element method 3

Completeness of trial solution

Figure: Pascal’s triangle

Pascal’s triangle can help to find out a complete polynomial

Asst.Prof.Dr. A. Ahmed Finite element method 4

Figure: Pascal’s triangle

Each row of the triangle gives the monomials that must be

Higher order element shape fucntions

◮ Choose the interpolation function

◮ Evaluate interpolation functions at nodal points

θ̂e (x) = Me (x)ae

◮ Interpolation function is then expressed as

θe (x) = p(Me )−1 θ̂e

in which Ne = p(Me )−1 is the shape function matrix

Asst.Prof.Dr. A. Ahmed Finite element method 6

A coordinate system that is defined for a particular element and

A local coordinate system which permits identification /

Asst.Prof.Dr. A. Ahmed Finite element method 7

Asst.Prof.Dr. A. Ahmed Finite element method 8

It is possible to establish a mapping between Global coordinate

Asst.Prof.Dr. A. Ahmed Finite element method 9

Asst.Prof.Dr. A. Ahmed Finite element method 10

Asst.Prof.Dr. A. Ahmed Finite element method 11

Examples of element coordinates are triangular coordinates and

We will next describe the use of shape functions expressed in

Asst.Prof.Dr. A. Ahmed Finite element method 12

Elements for which geometry and displacements are given by

◮ This property may be generalized to arbitrary elements. Thus,

Elements which utilizes higher order shape function for

What is Iso-parametric element

Asst.Prof.Dr. A. Ahmed Finite element method 14

Asst.Prof.Dr. A. Ahmed Finite element method 15

Asst.Prof.Dr. A. Ahmed Finite element method 17

Jacobian matrix can be found by taking derivatives of

To compute the derivatives in global coordinates we invert the

Thus the computation of the derivatives with respect to x

Derivatives of shape functions are

Asst.Prof.Dr. A. Ahmed Finite element method 19

Integration and nodal forces

When internal forces are computed by integration over the

Asst.Prof.Dr. A. Ahmed Finite element method 20

Asst.Prof.Dr. A. Ahmed Finite element method 21

Isoparametric Q4 element - Geometry

Consider the quadrilateral element as shown

Isoparametric Q4 element - Shape functions

Asst.Prof.Dr. A. Ahmed Finite element method 24

a = M−1 û =⇒ u = pM−1 û =⇒ u = Nû

Asst.Prof.Dr. A. Ahmed Finite element method 25

Isoparametric Q4 element - Shape functions

Asst.Prof.Dr. A. Ahmed Finite element method 26

Asst.Prof.Dr. A. Ahmed Finite element method 27

Isoparametric Q4 element - Shape functions

Asst.Prof.Dr. A. Ahmed Finite element method 28