A Three day workshop on Finite

Element Applications in Engineering

Analysis
(Application to 2D problems &
Higher order Elements)

Dr. K. Gopikrishna
Department of Civil Engineering
National Institute of Technology, Warangal
kgopi@nitw.ac.in 1
Dr. K. Gopikrishna, kgopi@nitw.ac.in 2
 Introduction
Two-dimensional (planar) elements are thin-plate elements such that two
coordinates define a position on the element surface.The elements are
connected at common nodes and/or along common edges to form
continuous structures

Nodal compatibility is then enforced during the formulation of the nodal

equilibrium equations for two-dimensional elements. If proper
displacement functions are chosen, compatibility along common edges is
also obtained.

The two-dimensional element is extremely important problems such

as plates with holes, fillets, or other changes in geometry that are loaded in their
plane resulting in local stress concentrations.

Dr. K. Gopikrishna, kgopi@nitw.ac.in 3

Dr. K. Gopikrishna, kgopi@nitw.ac.in 4
Dr. K. Gopikrishna, kgopi@nitw.ac.in 5
Dr. K. Gopikrishna, kgopi@nitw.ac.in 6
Dr. K. Gopikrishna, kgopi@nitw.ac.in 7
Dr. K. Gopikrishna, kgopi@nitw.ac.in 8
Dr. K. Gopikrishna, kgopi@nitw.ac.in 9
Finite Elements

Dr. K. Gopikrishna, kgopi@nitw.ac.in 10

Dr. K. Gopikrishna, kgopi@nitw.ac.in 11
Element Formulation

The Constant Strain Triangle element:

• 2D element used in plane

stress and plane strain
problems
• Nominal thickness = h
(small; can be variable)
• Three corner nodes with
coordinates (xi, yi)
• Two degrees of freedom
per node: (ui, vi)
Element Formulation (cont.)

Two approaches for generating shape

functions:
 “Interpolation approach”:
◦ Matrix-based method
◦ Works best for small numbers of d.o.f.
 “Direct approach”:
◦ More geometric method
◦ Works best for higher-order elements
 Element Formulation (cont.)
Some basic ideas (for both approaches):
 6 d.o.f. total  6 shape functions.
 Fundamental unknowns are horizontal displacement u(x,y) and vertical
displacement v(x,y).
  Each displacement expected to use 3 shape functions.

 Rule of thumb: simple shape functions = better shape functions.

(easier to integrate, more widely applicable, …)
 For 2D elements, polynomials in x and y are most common choice for
shape functions.
 If you have derivatives of order n in your variational principle, it is best
to choose your shape functions so that they can form a complete
polynomial of order n.
(Gives control over errors, faster convergence, …)
 Element Formulation (cont.)

Pascal’s triangle

(Row n+1 based upon expansion of (x+y)n. )

Element Formulation (cont.)
“Interpolation approach”:
◦ Approximate u(x,y) and v(x,y) by complete 1st order polynomials:
u  x, y   a1  a 2 x  a 3 y; v  x, y   a 4  a 5 x  a 6 y.

◦ At each node, require u(xi,yi) = ui and v(xi,yi) = vi :

u1  a1  a 2 x1  a 3 y1 ; v1  a 4  a 5 x1  a 6 y1.
u2  a1  a 2 x2  a 3 y2 ; v2  a 4  a 5 x2  a 6 y2 .
u3  a1  a 2 x3  a 3 y3 ; v3  a 4  a 5 x3  a 6 y3 .

6 equations for the 6 unknowns!

Element Formulation (cont.)
“Interpolation approach”:
◦ Write this in matrix form:

 u1  1 x1 y1 0 0 0   a1 
 v  0 0 0 1 x1 y1   a 2 
 1   
 u2  1 x2 y2 0 0 0   a3 
       d   C  a  .
 v2  0 0 0 1 x2 y2   a 4 
 u3  1 x3 y3 0 0 0   a5 
    
 v3   0 0 0 1 x3 y3   a 6 

  d .
◦ Solution (in symbolic form) is a   C
1
Element Formulation (cont.)
“Interpolation approach”:
◦ Now, rewrite interpolation functions in matrix/vector form:

 a1 
a 
 2
 u  x, y    1 x y 0 0 0   a 3 
       u  x     P  x    a  .
 v  x, y    0 0 0 1 x y   a 4 
 a5 
 
 a6 
 
◦ Substitute previous result: u  x    P  x   C
    d .
1

 N x  

Matrix of shape functions!

Element Formulation (cont.)
“Interpolation approach”:
◦ For CST, can show that
 u1 
v 
 1
 u  x, y    N1  x, y  0 N 3  x, y  0 N 5  x, y  0   u2 
 
    ,
 v  x , y    0 N 2  x , y  0 N 4  x , y  0 N 6  x , y    v2 
 u3 
 
 v3 
N1  x, y   N 2  x, y    x2 y3  x3 y2   x  y2  y3   y  x3  x2   2 A;
N 3  x, y   N 4  x, y    x3 y1  x1 y3   x  y3  y1   y  x1  x3   2 A;
N 5  x, y   N 6  x, y    x1 y2  x2 y1   x  y1  y2   y  x2  x1   2 A;
1 x1 y1 
A  area of triangle = 12 det 1 x2 y2  .
 
1 x3 y3 
Element Formulation (cont.)
Notes on “Interpolation approach”:
◦ This approach generalizes to different shapes, different node
locations, and different numbers of d.o.f
◦ However, the matrix [C] is not always invertible for general
choices of nodal locations.
◦ As number of d.o.f. increases, matrix inversion becomes
more difficult, and thus exact functions become harder to
determine.
Element Formulation (cont.)
 Visually, this looks like:
Element Formulation (cont.)
Check the required properties:
◦ Well-posedness: 6 d.o.f and 6 shape functions.
◦ Differentiability: Will show that only 1st derivatives show up
in variational principle, so need continuous shape functions.
◦ Completeness –
 Rigid body motion: Suppose u1  u2  u3  uˆ; v1  v2  v3  0.
 uˆ 
0
 
 u  x, y    N1  x, y  0 N 3  x, y  0 N 5  x, y  0   uˆ   uˆ * ( N1  N 3  N 5 ) 
 
    
 v  x , y    0 N 2  x, y  0 N 4  x, y  0 N 6  x, y    0   0 .

 uˆ 
 
0
N1  x, y   N 3  x, y   N 5  x, y    x2 y3  x3 y2    x3 y1  x1 y3    x1 y2  x2 y1   2 A  1. 
Element Formulation (cont.)
Check the required properties:
◦ Completeness –
 Constant strain: Can u
x   constant?
x
u N N N
x   u1 1  u2 2  u3 3  u1  y2  y3   u2  y3  y1   u3  y1  y2  2 A  constant. 
x x x x
In fact, strain must be a constant !
◦ Continuity:
Neither N3(x,y) nor N4(x,y) can
influence what happens on the line
between pt. 1 and pt. 2.
Only d.o.f. at pt. 1 and pt. 2
matter! 
Implementation of FEA – the CST

Overview: What is Finite Element Analysis?

 Assume a given problem is to be solved.
FEA proceeds as follows:
1. Discretize the problem domain into a finite number of
“local” approximate problems on regions called
elements.
2. Set up each of the local approximate problems.
3. Assemble the local problems into a global problem.
4. Solve the global problem and check convergence.
(If not converged, repeat steps #2 and #3.)
5. Once converged, evaluate the solution.
 Concept of Assembly
Concept of assembly:
◦ Each degree of freedom di in a given
element corresponds to a unique
degree of freedom Dk in the overall
object.
◦  Each local stiffness matrix [k]e
contributes to part of the global
stiffness matrix of the object, [K].
◦ Also, each local force vector (f)e
contributes to part of the global
force vector of the object, (F).
ne

K   D    F   "  "k e  d    f e .
e1
Global Local
 Assembly (cont.)
Concept of Assembly:
 Assembly is not straight addition –
 Assembly (cont.)

Pseudocode for Assembly:

For e = 1, numel ← sum over all elements
For i = 1, numdof(e) ← sum over all local d.o.f.
For j = i, numdof(e) ← local d.o.f.
ii = map(e,i); jj = map(e,j); ← get global d.o.f.
K(ii,jj) = K(ii,jj) + k(e,i,j); ← assemble [K]
Continue
F(ii) = F(ii) + f(e,i); ← assemble (F)
Continue
Continue
 Assembly (cont.)
Assembly by hand:
 Assembly (cont.)
 Create a set of “new” local coordinates that are parallel –
 d1  1   d1 
d   1 d 
 2   2 
 d3   cos  sin    d3   d3   cos  sin    d3 
  
 d     sin  cos    d   d      ,
 4   4   4    sin  cos    d4 
 d5   1   d5 
    
 6 
d 1   d6 
     d. or  d  T  d  .

T
d T

 Will also need to transform element force vector:

 f     T  f  .
 Assembly (cont.)
 Stiffness matrix then transforms –
k   d    f   Tk   d   T  f   f  .
But  d    T  d    Tk  T  d    f   .
T T

 k  These can now be assembled!

 Can show that same approach works for other types of

“transformations” (e.g., renumbering d.o.f., linking d.o.f, …)
Boundary Conditions (cont.)
◦ Three general techniques for enforcing
displacement boundary conditions:
 Condensation
 Penalty Method
 Lagrange Multipliers
◦ In all methods, must discretize the boundary
conditions to apply only at the nodes.
Boundary Conditions (cont.)
Condensation –
◦ Idea: formally remove constrained d.o.f. from the
calculation, but keep their effects on other d.o.f.

D1  0;
D2  0;
D6  0.
3.4: Boundary Conditions (cont.)
Condensation –
◦ For each constrained d.o.f., remove (“condense
out”) corresponding row and column from global
equation.
 K11 K12 K13 K14 K15 K16 K17 K18 K19 K110 K111 K112   0   F1 
K    
 21 K 22 K 23 K 24 K 25 K 26 K 27 K 28 K 29 K 210 K 211 K 212   0   F2 
 K 31 K 32 K 33 K 34 K 35 K 36 K 37 K 38 K 39 K 310 K 311 K 312   D3   F3 
    
    
 K 51 K 52 K 53 K 54 K 55 K 56 K 57 K 58 K 59 K 510 K 511 K 512   D5   F5 
    
New (F) vector
 K 61 K 62 K 63 K 64 K 65 K 66 K 67 K 68 K 69 K 610 K 611 K 612   0   F6 

K K 72 K 73 K 74 K 75 K 76 K 77 K 78 K 79 K 710 K 711 K 712   D7   F7 
 71    
    
    
    
    
    
    
 K121 K122 K123 K124 K125 K126 K127 K128 K129 K1210 K1211 K1212   D12   F12 
 

New [K] matrix

What if Di ≠ 0?
Boundary Conditions (cont.)
Notes on condensation –
◦ Very powerful method; can handle a variety of displacement
boundary conditions exactly.
◦ Reduces bandwidth by eliminating d.o.f.
◦ If all boundary conditions are single-point, this method is equivalen
to simply “condensing out” the constrained d.o.f.
and adding in “extra” forces due to the imposed displacements.
◦ If there are many multi-point constraints, the process of
re-numbering, transforming, and condensing can be very
time-consuming.
◦ Condensation can sometimes be done at element level.
Example

kgopi@nitw.ac.in 35
Example Problem

◦ Given: Cantilevered beam with dimensions shown;

rigidly fixed at x = 0; applied traction
t  2.4 j ksi
at x = 18. E = 27,000 ksi;  = 0.25 .
◦ Required: Using CST plane stress elements, find
the approximate deflection of the free end at y =
0.
Example (cont.)
◦ Mesh the beam as shown:

6 elements;
16 total d.o.f.;
4 constraints.

◦ Exact solution is known:

 2.4 ksi   0.5 in  18 in 

2 3
 AL
3
v  x  18, y  0      0.5184 in.
3EI 3  27000 ksi   0.1667 in  4
 Formulate the elements –
◦ Shape functions:

Element #1: Element #2:

N1  x, y   1  16 x  12  y  1 . N1  x, y   1  16  x  6   12  y  1 .
N 2  x, y   16 x. N 2  x, y    16  x  6  .
N 3  x, y   12  y  1 . N 3  x, y    12  y  1 .
Example (cont.)
Formulate the elements –
◦ Shape functions for other elements:

◦ Element #3 is simply Element #1 shifted from x

= 0 to x = 6; thus, can “shift” each shape function:
N1  x, y   1  16  x  6   12  y  1 .
N 2  x, y   16  x  6  .
N 3  x, y   1
2  y  1 .
◦ Element #4 = Element #2 shifted from x = 6 to x
= 12:N1  x, y   1  16  x  12   12  y  1 ,
N 2  x, y    16  x  12  , N 3  x, y    12  y  1 .

◦ Element #5 = Element #1 shifted from x = 0 to x

= 12: N1  x, y   1  16  x  12   12  y  1 ,
N 2  x, y   16  x  12  , N 3  x, y   1
2  y  1 .
◦ Element #6 = Element #2 shifted from x = 6 to x
= 18:
N1  x, y   1  6  x  18   2  y  1 ,
1 1

N 2  x, y    16  x  18  , N 3  x, y    12  y  1 .
 Example (cont.)
Formulate the elements –
◦ [B] matrix for Element #1:
 
 0 
 x 
 1
   x y 0
1 1
x 0
1
y1 1
0 
B  x       N  x   0
2 6 2 6 2 2

#1 #1
y  0 1
 x y 0
1 1
x 0
1 1
y  
1

  2 6 2 6 2 2

 
 y x 
 
1
6  0
1
6 0 0 0
 0  0 0 0  .
1 1
2
 2

1
2
1
6    0
1
6
1
2 0 
 Example (cont.)
Formulate the elements –
◦ [B] matrix for Element #2:
 
 0 
 x 
    
1 1
x y 0 1
1 x 0
1 1
 y 0
1

B  x    0
2 6 2 6 2 2

#2
y  0   x y 0
1 1 1
1 x 0
1 1
 y 
1

  2 6 2 6 2 2

 
 y x 
 
1
6  0  0 1
6 0 0 
 01
0 0 0     B  x   .
1

 2 2
 #1
1
2 1
6 0  1
6  1
2 0 
 Example
Other [B] matrices:
◦ Since Elements #3 and #5 are both simply shifts of
Element #1, coefficients of x and y do not change.
 1
6 0 1
6 0 0 0
 B  x    B  x    0  1
0 0 0 . 1
#3 #5  2
 2

  1
2  1
6 0 0 
1
6
1
2

◦ Likewise, Elements #4 and #6 are shifts of Element

#2, so have same [B] matrices.
 0 
1
6
1
6 0 0 0 
 B  x    B  x    0 0 1
0 0  .
1
#4 #6  2 2

1
2 0 1
6  1
6  1
2 0 
 Example..
Element stiffness matrices –
◦ Elastic matrix is:
1  0   28800 7200 0 
E   
C    1 0   7200 28800 0  ksi.
1 2  
 0 0 1
2 1    0 0 10800 

◦ Since [B] and [C] are constant matrices, the

volume integral for [k] reduces to
       
     CB * hAtriangle .
T
 
T
k B C B * hdA k B
triangle
◦ For Element #1:
  16 0  12 
 0  1  1
 2 6
 28800 7200 0    16 0 16 0 0 0 
 6 0 0 
1
 0  1 0 0 0 1 
 #1  #1   #1 triangle 
  
T
k B C B * hA 0.25*6 1   7200 28800 0  2 2
 0 0  0
6
 0 10800    12  16 0 16 12 0 
0 0 2  1
 
 0 12 0 
 5250 2250 1200 1350 4050 900 
 2250 11250 900 450 1350 10800 
 
 1200 900 1200 0 0 900 
  kips/in.
 1350  450 0 450 1350 0 
 4050 1350 0 1350 4050 0 
 
  900 10800 900 0 0 10800 
 Example…
◦ For the other elements, we notice the following:
B #n #1 #n 
   B   k     B 
T

#1 C B  * hA
#1 triangle
 k  .
#1

All elements have the same stiffness matrix!

 ◦ Only element force vector for Element #6:
  16 x  12 y  52 0     12 y  12 0    0 
 5   1   
6 x  2 y  2 2 y  2  
1 1 1
 0   0   1.2 y 1.2 
  3  16 x 0   0     0 0   0    0 
 N  x   #6  t    
T
     
  .
 0 3  16 x   2.4    0 0   2.4    0 
  12  12 y 0     12  12 y 0    0 
       
  0 1
2  1
2 y   x18   0
1
2  2 
1
y   1.2 y  1.2 

 0   0 
 1.2 y  1.2   0.6 
   
y 1
 0   0 
 f #6    N  x    t  dA  h *  
T
dy    kips
Ae , y 1 
0   0 
 0   0 
   
 1.2 y  1.2   0.6 
 Assembling K
matrix
◦ Element #1:
1 2 3 4 5 6
1  5250 2250 1200 1350 4050 900  3

2  2250 11250 900 450 1350 10800  4
3  1200 900 1200 0 0 900  7
 
4  1350 450 0 450 1350 0 8
5  4050 1350 0 1350 4050 0 1
 
6  900 10800 900 0 0 10800  2
3 4 7 8 1 2
◦ Element #2:
1 2 3 4 5 6
1  5250 2250 1200 1350 4050 900  5
2  2250 11250 900 450 1350 10800  6
3  1200 900 1200 0 0 900  1
 
4  1350 450 0 450 1350 0 2
5  4050 1350 0 1350 4050 0 7
 
6  900 10800 900 0 0 10800  8
5 6 1 2 7 8
 Assembly…

◦ Stiffness matrix for #1 + #2 only:

1 2 3 4 5 6 7 8
1 5250 0 4050 1350 1200 900 0 
2250
2 0 11250 900 10800 1350 450 2250 0 

3  4050 900 5250 2250 0 0 1200 1350 
 
4  1350 10800 2250 11250 0 0 900 450 
kips/in.
5  1200 1350 0 0 5250 2250 4050 900 
 
6  900 450 0 0 2250 11250 1350 10800 
7 0 2250 1200 900 4050 1350 5250 0 
 
8  2250 0 1350 450 900 10800 0 11250 

◦ What about stiffness matrix for Element #3 + #4?

 Only real difference between #3+#4 and #1+#2 is the
numbering of the d.o.f. – all shifted by 4
 Have same matrix for different d.o.f.
◦ Stiffness matrix for #3 + #4 only:
5 6 7 8 9 10 11 12
5 5250 0 4050 1350 1200 900  0 2250
6 0 11250 900 10800 1350 450 2250 0 

7  4050 900 5250 2250 0 0 1200 1350 
 
8  1350 10800 2250 11250 0 0 900 450 
kips/in.
9  1200 1350 0 0 5250 2250 4050 900 
 
10  900 450 0 0 2250 11250 1350 10800 
11 0 2250 1200 900 4050 1350 5250 0 
 
12  2250 0 1350 450 900 10800 0 11250 

◦ Stiffness matrix for #5 + #6 only:

9 10 11 12 13 14 15 16
9 5250 0 4050 1350 1200 900  0 2250
10 0 11250 900 10800 1350 450 2250 0 

11  4050 900 5250 2250 0 0 1200 1350 
 
12  1350 10800 2250 11250 0 0 900 450 
kips/in.
13  1200 1350 0 0 5250 2250 4050 900 
 
14  900 450 0 0 2250 11250 1350 10800 
15 0 2250 1200 900 4050 1350 5250 0 
 
16  2250 0 1350 450 900 10800 0 11250 
 Assembled stiffness matrix
◦ Stiffness matrix for entire structure:
5250 0 4050 1350 1200 900 0 2250 0 0 0 0 0 0 0 0 
0 11250 900 10800 1350 450 2250 0 0 0 0 0 0 0 0 0 
 
 4050 900 5250 2250 0 0 1200 1350 0 0 0 0 0 0 0 0 
 
 1350 10800 2250 11250 0 0 900 450 0 0 0 0 0 0 0 0 
 1200 1350 0 0 10500 2250 8100 2250 1200 900 0 2250 0 0 0 0 
 
 900 450 0 0 2250 22500 2250 21600 1350 450 2250 0 0 0 0 0 
0 2250 1200 900 8100 2250 10500 2250 0 0 1200 1350 0 0 0 0 
 
 2250 0 1350 450 2250 21600 2250 22500 0 0 900 450 0 0 0 0 
K     kips/in
0 0 0 0 1200 1350 0 0 10500 2250 8100 2250 1200 900 0 2250 
0 0 0 0 900 450 0 0 2250 22500 2250 21600 1350 450 2250 0 
 
0 0 0 0 0 2250 1200 900 8100 2250 10500 2250 0 0 1200 1350 
0 0 0 0 2250 0 1350 450 2250 21600 2250 22500 0 0 900 450 
 
0 0 0 0 0 0 0 0 1200 1350 0 0 5250 2250 4050 900 
0 0 0 0 0 0 0 0 900 450 0 0 2250 11250 1350 10800 
 
0 0 0 0 0 0 0 0 0 2250 1200 900 4050 1350 5250 0 
 
0 0 0 0 0 0 0 0 2250 0 1350 450 900 10800 0 11250 
 Force vector & BCS
◦ Force vector for entire structure:
0 
 
0 
0 
 
0 
0 
 
0 
0 
 
0 
 F     kips.
0 
0 
 
0 
0 
 
0 
 -0.6 
 
0 
 
 -0.6 
 Boundary conditions
Enforce constraints –
◦ Using condensation, you must simply eliminate the
first 4 rows and columns of [K] and first 4 rows of
(F):
10500 2250 8100 2250 1200 900 0 2250 0 0 0  0 0 
 2250   
 22500 2250 21600 1350 450 2250 0 0 0 0 0  0 
 8100 2250 10500 2250 0 0 1200 1350 0 0 0 0  0 
   
 2250 21600 2250 22500 0 0 900 450 0 0 0 0  0 
 1200 1350 0 0 10500 2250 8100 2250 1200 900 0 2250  0 
   
900 450 22500 2250 21600 1350 450 0 
K    
0 0 2250 2250 0
kips/in;  F   =  kips.
0 2250 1200 900 8100 2250 10500 2250 0 0 1200 1350  0 
   
 2250 0 1350 450 2250 21600 2250 22500 0 0 900 450  0 
   
0 0 0 0 1200 1350 0 0 5250 2250 4050 900  0 
0 0 0 0 900 450 0 0 2250 11250 1350 10800   -0.6 
   
0 0 0 0 0 2250 1200 900 4050 1350 5250 0  0 
0 1350 450 900 10800 0 11250   
 0 0 0 2250 0  -0.6 
 ◦ This can now be solved for the remaining d.o.f:
1.759 
 
 -6.563 
 -1.702 
 
 -6.487 
 2.845 
 
 -21.37 
 D  =   10-3in.
-2.674 
 
 -21.29 
 
 3.245 
 -40.28 
 
 -2.959 
 -40.21 
 

◦ Use this to interpolate the requested displacement:

v  x  18, y  0   D14 N1(6)  x  18, y  0   D10 N 2(6)  x  18, y  0   D16 N 3(6)  x  18, y  0 
=  -40.28  10-3in   0.5    -21.37  10-3in   0    -40.21  10-3in   0.5   40.24  10-3in.
Very poor approximation!
• Properties of shape functions
• Higher order elements in 1D
• Higher order triangular elements (using
area coordinates)
• Higher order rectangular elements
• Lagrange family
• Serendipity family

Dr. K. Gopikrishna, kgopi@nitw.ac.in 55

 Recall that the finite element shape functions need to satisfy
the following properties

1. Kronecker delta property

 1 at node i
Ni  
0 at all other nodes

Inside an element
u  N1u1  N 2 u 2  ...
At node 1, N1=1, N2=N3=…=0, hence

u node 1  u1

Facilitates the imposition of boundary conditions

 2. Polynomial completeness

If u  1   2 x   3 y
Then
N
i
i 1

N x
i
i i x

N y
i
i i y
 Higher order elements in 1D

2-noded (linear) element:

x  x2
N1 
x1  x 2
x1 x2 x  x1
x N2 
x 2  x1
1 2

In “local” coordinate system (shifted to center of element)

ax
N1 
x 2a
ax
2 N2 
1 2a
a a
 3-noded (quadratic) element:
x  x 2 x  x3 
N1 
x1 x3 x2 x1  x 2 x1  x3 
x  x1 x  x3 
x N2 
1 3 2 x 2  x1 x 2  x3 
x  x1 x  x 2 
N3 
x3  x1 x3  x 2 
In “local” coordinate system (shifted to center of element)
x1  a ; x2  0 ; x3  a
x
xa  x 
N1  
1 3 2 2a 2
a a xa  x 
N2 
2a 2
a2  x2
N3 
a2
 x  x 2 x  x3 x  x 4 
N1 
4-noded (cubic) element: x1  x 2 x1  x3 x1  x4 
x  x1 x  x3 x  x4 
N2 
x1 x3 x4 x2 x 2  x1 x 2  x3 x 2  x 4 
x x  x1 x  x2 x  x 4 
N3 
1 3 4 2 x3  x1 x3  x 2 x3  x4 
x  x1 x  x 2 x  x3 
N4 
x 4  x1 x 4  x 2 x 4  x3 
In “local” coordinate system (shifted to center of element)
9
N1   ( x  a )( x  a / 3)( x  a / 3)
2a/3 2a/3 2a/3 x 16a 3

9
N2  ( x  a )( x  a / 3)( x  a / 3)
1 3 4 2 16a 3
a a
27
N3  3
( a  x)(a / 3  x)(a  x)
16a
27
N4   3
( x  a )( x  a / 3)( x  a / 3)
16a
 Polynomial completeness

1 Convergence
rate (displacement)
x 2 node; k=1; p=2
x2
3 node; k=2; p=3
x3
4 node; k=3; p=4
x4

Recall that the convergence in displacements
u  uh 0  Ch p ; p  k  1
k=order of complete polynomial
Triangular elements

Area coordinates (L1, L2, L3)

1
Total area of the triangle A=A1+A2+A3
At any point P(x,y) inside
y P A2 the triangle, we define
A3 A1
A1 L1 
3 A
A2
L2 
2 A
x A3
L3 
A
Note: Only 2 of the three area coordinates are independent, since
L1+L2+L3=1
 ai  bi x  ci y
Li 
2A

1 x 1 y1 
1 
A  area of triangle  det 1 x 2 y2 
2
1 x 3 y 3 
a1  x 2 y 3  x3 y 2 b1  y 2  y 3 c1  x3  x 2
a 2  x3 y1  x1 y 3 b2  y 3  y1 c 2  x1  x3
a3  x1 y 2  x 2 y1 b3  y1  y 2 c3  x 2  x1
Check that

L1  L2  L3  1
L1 x1  L2 x2  L3 x3  x
L1 y1  L2 y 2  L3 y3  y
 1
L1= constant
P’
y P
A1
3
2
x

Lines parallel to the base of the triangle are lines of constant ‘L’
 L1= 1
1

y P L1= 0
A1
3
2 L2= 0
L3= 1
x L =1
2
L3= 0

We will develop the shape functions of triangular elements in

terms of the area coordinates
 For a 3-noded triangle

N1  L1
N 2  L2
N 3  L3
 Higher order elements are useful if the boundary of the geometry is curve
in nature.
 For curved case, higher order triangular element can be suited effectively
while generating the finite element mesh.
 Moreover, in case of flexural action in the member, higher order elements
can produce more accurate results compare to those using linear element
 For a 6-noded triangle

L2= 0 L1= 1
1
L2= 1/2 L1= 1/2
6
4
y L1= 0
L2= 1
3
5
2
x
L3= 0 L3= 1/2 L3= 1
How to write down the expression for N1?

Realize the N1 must be zero along edge 2-3 (i.e., L1=0) and at
nodes 4&6 (which lie on L1=1/2)
N1  cL1  0L1  1 / 2
Determine the constant ‘c’ from the condition that N1=1 at
node 1 (i.e., L1=1)

N1 ( at L1  1)  c11  1 / 2   1
c2
 N1  2 L1 L1  1 / 2 
N1  2 L1 ( L1  1 / 2)
N 2  2 L2 ( L2  1 / 2)
N 3  2 L3 ( L3  1 / 2)
N 4  4 L1 L2
N 5  4 L3 L2
N 6  4 L3 L1

Such six node triangular element is commonly known as linear strain

triangle (LST) as its strain is assumed to vary linearly inside the element
The main advantage of this element is that it can capture the variation of strains
and therefore stresses of the element
 For a 10-noded triangle

L2= 0 L1= 1
L2= 1/3 1 L1= 2/3
9
L2= 2/3 4 L1= 1/3
8
y 10 L1= 0
L2= 1 5
3
7
2 6
x
L3= 0 L3= 1/3 L3= 2/3 L3= 1
 9
N1  L1 ( L1  1 / 3)( L1  2 / 3)
2
9
N2  L2 ( L2  1 / 3)( L2  2 / 3)
2
9
N3  L3 ( L3  1 / 3)( L3  2 / 3)
2
27
N4  L1 L2 ( L1  1 / 3)
2
27
N5  L1 L2 ( L2  1 / 3)
2
27
N6  L2 L3 ( L2  1 / 3)
2
27
N7  L2 L3 ( L3  1 / 3)
2
:
N 10  27 L1 L2 L3
 NOTES:
1. Polynomial completeness
1 Convergence
rate (displacement)
x y
3 node; k=1; p=2
2 2
x xy y
3 2 2 3
6 node; k=2; p=3
x x y xy y
10 node; k=3; p=4

 Shape Function by Degrading
Technique
 Sometimes, the geometry of the structure or its loading and boundary
conditions are such that the stresses developed in few locations are
quite high.
 On the other hand, variations of stresses are less in some areas and as
a result, refinement of finite element mesh is not necessary.
 It would be economical in terms of computation, if higher order elements
are chosen where stress concentration is high and lower order elements
at area away from the critical area

Type 1 has only three nodes, type 2

element has five nodes, type 3 has four
nodes and type 4 has six nodes

The shape functions of 6-node element can

suitably be degraded to derive shape functions
of other two types of triangular elements

kgopi@nitw.ac.in 74
 5 noded triangular Element
 u = N u1 +N2u2 +N3u3 +N4u4 +N5u5 + N6u6
If there is no node between 1 and 3, the
displacement along line 1-3is considered to
vary linearly
Thus the displacement at an assumed node6´
may be written as

75
Thus, for a five node triangular element, the above shape
function can be used for finite element analysis.

kgopi@nitw.ac.in 76
 2. Integration on triangular
domain
k! m! n!
 L1 L2 L3 dA  2 A
k m n
1.
1 A (2  k  m  n)!
k! m!
2.  L1 L2 dS  l12
k m
l1- 1 2 edge (1  k  m)!
y 2

3
2
x
 3. Computation of derivatives of shape functions: use chain rule
e.g.,
N i N i L1 N i L2 N i L3
  
x L1 x L2 x L3 x
L1 b1 L2 b2 L3 b3
But  ;  ; 
x 2 A x 2 A x 2 A

e.g., for the 6-noded triangle

N 4  4 L1 L2
N 4 b b
  4 L1 2  4 L2 1
x 2A 2A
The derivation of interpolation function in terms of Cartesian coordinate
system is algebraically complex as seen from earlier section. However, the
complexity can be reduced by the use of natural coordinate system, where the
natural coordinates will vary from -1 to +1 in place of –a to +a or –b to +b.
Rectangular elements
Lagrange family
Serendipity family

Lagrange family

4-noded rectangle In local coordinate system

y
2 1 ( a  x )(b  y)
a a N1 
b 4ab
( a  x )(b  y)
x N2 
b 4ab
( a  x )(b  y)
3 4 N3 
4ab
( a  x )(b  y)
N4 
4ab
9-noded
Corner
quadratic
nodes  x(a  x)   y(b  y)   x(a  x)   y (b  y ) 
y 5 N1     2b 2  N 2   2a 2   2b 2 
     
2
2 a
2 a 1
a  x(a  x)   y (b  y )   x(a  x)   y (b  y ) 
N 3      N 4   2a 2   2b 2 
b      
2 2
2 a 2b
9 8
6 x Midside
b 7 nodes
 a 2  x 2   y (b  y )   x(a  x)   b  y 
2 2
3 4 N5     N 6    
 a   2b 
2 2
 2a 2   b 2 
 a 2  x 2   y (b  y )   x(a  x)   b  y 
2 2
N7      N 8   2a 2   b 2 
   
2 2
 a 2b 
Center
node
 a2  x2  b2  y 2 
N9   2   2 
 a  b 
 NOTES:
1. Polynomial completeness
Convergence
1 rate (displacement)
x y 4 node; p=2
x2 xy y2 9 node; p=3
x3 x2 y xy 2 y 3
x4 x3 y x 2 y 2 xy 3 y 4
x 5 x 4 y x 3 y 2 x 2 y 3 xy 4 y 5

Lagrange shape functions contain higher order terms but miss

out lower order terms
Serendipity
family
4-noded same as
Lagrange
8-noded rectangle: how to generate the shape
functions? First generate the shape functions
y 5 of the midside nodes as appropriate
2 a a 1
products of 1D shape functions,
b e.g., 2 2
8
x N   a  x   (b  y)  ; N   (a  x)   b  y 
2 2
6
5  2   8  2a   b 2 
b 7  a   2 b    
3 4 Then go to the corner nodes. At each
corner node, first assume a bilinear
shape function as in a 4-noded element
and then modify: ˆ (a  x)(b  y )
“bilinear” shape fn at N1 
4ab
node 1:
ˆ N5 N8
actual shape fn at N1  N1  
2 2
8-noded
rectangle Midside
y 5 nodes
2 a a 1 a  x
2 2
  (b  y )   (a  x)   b  y 
2 2
N5   2   N6    2 
b  a   2b   2a   b 
8
6 x  a 2  x 2   (b  y )   (a  x)   b  y 
2 2
N7   2   N8    2 
b 7  a   2b   2a   b 
3 4

Corner
nodes
(a  x)(b  y ) N 5 N 8 (a  x)(b  y ) N 5 N 6
N1    N2   
4ab 2 2 4ab 2 2
(a  x)(b  y ) N6 N 7 (a  x)(b  y ) N8 N 7
N3    N4   
4ab 2 2 4ab 2 2
 NOTES:
1. Polynomial completeness
Convergence
1 rate (displacement)
x y 4 node; p=2
x2 xy y2 8 node; p=3
x3 x2 y xy 2 y 3 12 node; p=4
x4 x3 y x 2 y 2 xy 3 y 4
5 4 3 2 2 3 4 5
16 node; p=4
x x y x y x y xy y

More even distribution of polynomial terms than Lagrange

shape functions but ‘p’ cannot exceed 4!
 Using natural coordinates
It may be noted that the cubic terms
are omitted and geometric invariance
is ensured

The nodal field variables can be obtained from the above expression
after putting the coordinates at nodes

kgopi@nitw.ac.in 85
kgopi@nitw.ac.in 86
The shape functions of rectangular elements with higher nodes can be
derived in similar manner using appropriate polynomial satisfying all
necessary criteria. However, difficulty arises due to the inversion of large size
of the matrix because of higher degree of polynomial chosen.
kgopi@nitw.ac.in 87
Lagrange Interpolation Function
 An alternate and simpler way to derive shape functions
is to use Lagrange interpolation polynomials. This
method is suitable to derive shape function for elements
having higher order of nodes

kgopi@nitw.ac.in 88
3 noded element

kgopi@nitw.ac.in 89
 2D elements
Those elements whose shape functions are derived from the products
of one dimensional Lagrange interpolation functions are called
Lagrange elements. The Lagrange interpolation function for a
rectangular element can be obtained from the product of appropriate
interpolation functions in the ξ direction [fi(ξ )] and η direction [fi(η)].

kgopi@nitw.ac.in 90
9 noded Element

kgopi@nitw.ac.in 91
Thus, it is observed that the two dimensional Lagrange element contains
internal nodes which are not connected to other nodes, do not contribute
to the inter element connectivity
However, these can be eliminated by condensation procedure which needs
extra computation. The elimination of these internal nodes results in reduction
in size of the element matrices. Alternatively, one can develop shape functions
of two dimensional elements which contain nodes only on the boundaries.
These elements are called serendipity elements
kgopi@nitw.ac.in 92
Serendipity Elements
 interpolation functions can be derived by inspection

kgopi@nitw.ac.in 93
kgopi@nitw.ac.in 94
It may be observed that the Lagrange elements have a better degree of
completeness in polynomial function compare to serendipity elements.
Therefore, Lagrange elements produce comparatively faster and better
accuracy

kgopi@nitw.ac.in 95
Concept of Jacobian Matrix
 As the both element geometry and variation of the shape functions
are represented in terms of the natural coordinates of the parent
element, some additional mathematical obstacle arises.
 The relationship between two coordinate systems may be computed
by using the chain rule of partial differentiation as

One of the indices

for distortion and
important parameter
for quality of shape
of an element kgopi@nitw.ac.in 96
 Resources
 The material in this academic
presentation has been compiled from
following sources.
 http://www.colorado.edu/engineering/CA
S/courses.d/IFEM.d
 http://nptel.ac.in/courses/105105041/
 P.Seshu, “Textbook of Finite Element
Analysis”,PHI Learning Private limited,
seventh Printing,2009.
Thank you for your attention

Dr. K. Gopikrishna, Design of Compression

members, kgopi@nitw.ac.in 98