Finite Element Method

MEE710

Presented by
Kallol Khan

Mechanical Engineering Department

National Institute of Technology Durgapur
1
 Approximation Methods for solving Differential
Equations, weak form of differential equation
 One-dimensional FE formulation
 FE formulation of truss and frames
 Two dimensional FE formulation, Plane stress/ plane
strain problem, Axisymmetric problem.

 FE formulation for bending of beam 2

Text Books:
1. Text book of Finite Element Analysis by P. Sesu (PHI)
2. Introduction to Finite Elements in Engineering by T. R. Chandrupatla, A. D.
Belegundu ( Prentice- Hall)
3. An Introduction to the Finite Element Method by J. N. Reddy (Tata McGraw Hill)

Reference Books:
1. Finite Element Procedures by K. J. Bathe ( Prentice Hall)
2. Finite Element analysis Theory and Programming by C. S. Krishnamoorthy (Tata
McGraw Hill)
3. Concepts and applications of finite element analysis by R. D. Cook, D. S. Malkus
etc. (Wiley)

3
https://www.youtube.com/watch?v=KR74TQesUoQ
https://www.youtube.com/watch?v=LCTp7H6Tb8w
https://www.youtube.com/watch?v=MldJ6WHCsvQ
:
:
:

Dr. R. Krishnakumar,Department of Mechanical Engineering,IIT Madras

4
 Mathematical Modeling
SYSTEM ( A set of ODE/PDE)

d 2 g
2
 sin   0
dt L
 g 
 (t )   0 cos  t 
 L  Exact Solution Approximate Solution

d 
at t  0,   0 and   0
dt

Analytical Solution Numerical Solution

(closed form) Finite Element Method(FEM)
Finite Difference Method(FDM)
Finite Volume Method(FVM)
Boundary Element Method(BEM)
Differential Quadrature Method(DQM)
5
 Approximate solution of differential equation
( Weighted residual methods)
1. Bubnov Galerkin Method (Galerkin Method)
2. Petrov Galerkin Method
3. Method of Least square
4. Collocation Method
5. Subdomain Collocation Method
The governing equation of a bar subjected to uniformly distributed load

qxN/m

L
Δx
dF 1  d 2F 
F+ΔF  F  x   2   x   ......
2
F dx 2!  dx 
Δx
F x 0
dF
F F Δx  qx Δx  0
dx
dF
  qx  0..................( 1) 6
dx
 dF
 q x  0............(1)
dx
du σ
Now, F  σ xx A and ε xx  and xx  E
dx ε xx
du
F  EA ..........(2)
dx
From Equations (1)and (2)
d  du 
EA  q x  0.............(3)
dx  dx 

Boundary conditions :
du du
u  0 at x  0; and EA 0  0 at x  L
dx dx

 du 
 where EA is force, and at the free 
dx
 
 end no force is acting in the bar 
7
 The governing equation of a beam subjected to uniformly distributed load
z
z
qz

h y
x
A B
L

Cross-sectional area of beam: A b

Young’s modulus of beam material: E
Second moment of area (about neutral axis): I
qz M B 0
 Fz  0 Δx
 M   M  ΔM   v x   qz Δx  0
2
 v   v  Δv   q z Δx  0 M+ΔM ΔM qz Δx
M v  v ......(2)
v+Δv Δx 2
v
  q z .........(1) For Δx  0
x
A B v dv ΔM dM
 and 
Δx x dx Δx dx 8
 dv dM
  q z ........(3) and   v........(4)
dx dx

d 2 w( x) d  d 2 w( x) 
If the transverse deflection of the beam is w( x), then M  EI 2
, EI 2 
v
dx dx  dx 
d2  d 2 w( x) 
EI 2 
 qz  0
dx 2  dx 
d 4 w( x)
EI 4
 qz  0
x
dx
w(x) Fixed Support: dw(x)
Deflection : w(x)  0 , and Slope : 0
dx
x Free edge:
w(x) d3w d3w d2w d2w
z Shear Force : EI 3  0  3  0 , Moment : EI 2  0  2  0
dx dx dx dx
x z Roller Support/Hinge Support:
d2w d2w
x Deflection : w( x)  0 , Moment : EI 2  0  2  0
dx dx
9
 1.Galerkin’s Method
Solution of bar problem
qxN/m
d  du  d 2u
EA  qx  0  EA 2  q0  0
x 
dx  
dx  dx
L
du du
Boundary Conditions: 1: at x  0, u  0 2: at x  L, EA  0  dx  0
dx

Assume u  c 0  c1 x  c 2 x 2  c 3 x 3  ..................
c 0 , c1 , c 2 , c 3 .......... ..... are unknown constants
Let u  u  c 0  c1 x  c 2 x 2  c3 x 3
Applying boundary condition 1: c 0  0
Applying boundary condition 2:
du
 c1  2 c 2 x  3 c 3 x 2  0  c1  2 c 2 L  3 c3 L  0  c1   2 c 2 L  3 c3 L
2 2

dx at x  L

 
u   2c 2 L  3 c3 L2 x  c 2 x 2  c 3 x 3 = (x2-2Lx) c2 + (x3 -3L2x) c3
         
1st approximation 2nd approximation
function function
10
 u  (x2-2Lx) c2 + (x3 -3L2x) c3
d 2u
  2c 2  6 x c 3 
dx 2
d 2u
EA 2  q x  0
dx
d 2u
EA 2  q x  R d R d : Residual
dx
x 2  4 x  4  0 The roots are 2
(1.9) 2  4(1.9)  4  R d
(1.99) 2  4(1.99)  4  R d
Rd  EA  2 c 2  6 x c3   q x

For Galerkin’s method

L

R
0
d Wi dx  0 Wi are called weight function

For Galerkin’s method the weight functions are same as approximation functions
u (x2-2Lx) c + (x3 -3L2x) c3
   2      11
1st weight function  W1 2nd weight function  W2
 L

R
0
d Wi dx  0

Rd  EA  2 c 2  6 x c3   q x
W1 = (x2-2Lx) W2= (x3 -3L2x)
L

  EA  2 c 2  6 x c3   qx 
0
(x2-2Lx) dx = 0………….(1)
L

  EA 2c 2  6 x c3   q x   x 3  3 L2 x  dx  0...........(2)
0

Assignment 1A: Calculate the values of c2 and c3

Assignment 1B: Calculate the values of the constants for 5 terms solution
u  u  c 0  c1 x  c 2 x 2  c3 x 3  c 4 x 4

12
 2.Petrov Galerkin’s Method
u  (x2-2Lx) c2 + (x3 -3L2x) c3
d 2u
  2c 2  6 x c 3 
dx 2
Rd  EA  2 c 2  6 x c3   q x
For Petrov Galerkin’s method the weight functions are not same as approximation functions
L

R
0
d Wi dx  0

 EA 2c 2  6 x c3  qx  (1)

0 
dx = 0…….(1)
W1

L
Rd

  EA 2c  6 xc q  x  dx  0...........(2)

0
2 3 x
W2
Rd

Assignment 1C: Calculate the values of c2 and c3

Assignment 1D: Calculate the values of the constants for 5 terms solution
u  u  c 0  c1 x  c 2 x 2  c3 x 3  c 4 x 4 13
 3.Method of Least square
In least square method the parameters (in this case c 2 and c3) are determined by
minimizing the integral of the square of the residual.
L

R
2
d dx
0

  2 
L

ci  0
 R d dx   0

Rd  EA  2 c 2  6 x c3   q x
  
L
   
c 2  0
2
 EA 2c 2  6 xc 3  q x dx   0..........(1)

u  (x2-2Lx) c2+ (x3 -3L2x) c3
  
L
   
c 3  0
2
 EA 2c 2  6 xc 3  q x dx   0.......... ...( 2)

L
 Rd
 Rd
0
ci
dx  0

Assignment 1E: Calculate the values of c2 and c3

Assignment 1F: Calculate the values of the constants for 5 terms solution 14
 4. Collocation Method
.
In this method the residual (Rd) is made zero at some selected points called
collocation point
L 2L
Let for the present case the collocation points are and
3 3
Rd  EA  2 c 2  6 x c3   q x
L  L
R d  0 at x   EA 2a2  6a3   q x  0..........(1)
3  3
u  (x2-2Lx) c2+ (x3 -3L2x) c3
2L  2L 
R d  0 at x   EA 2a2  6a3   q x  0.........( 2)
3  3 
The selection of collocation points is crucial in obtaining a well conditioned system of
equations and a convergent solution. The collocation points should be located as evenly as
possible to avoid ill-conditioning of the resulting equations.

Assignment 1G: Calculate the values of c2 and c3

Assignment 1H: Calculate the values of the constants for 5 terms solution
Assignment 1I: What do you mean by ill conditioned equations explain with example15
 5. Subdomain Collocation Method
Here the domain is first subdivided into n subdomains and the integral of the
residual over each subdomain is then required to be zero.
 L L 
If the subdomains are  0 to  and  to L 
 2 2 
Rd  EA  2 c 2  6 x c3   q x

L
2

  EA 2c 2  6 xc 3   q x  dx..........(1)
0

L
u  (x2-2Lx) c2+ (x3 -3L2x) c3
  EA 2c
L
2  6 xc3   q x  dx..........(2)
2

Assignment 1J: Calculate the values of c2 and c3

Assignment 1K: Calculate the values of the constants for 5 terms solution

16
 Example: 1
Consider a 1mm diameter, 50mm long aluminium pin-fin as shown in Fig. used to
enhance heat transfer from a surface wall maintained at 300 0C. The governing
differential equation and the boundary conditions are given by:
h Tα=300C
d 2T p h
dia 1mm k 2   T  T 
wall x dx A
3000C Boundary conditions :
L=50mm tip is insulated dT
1 : T(0)  Tw  3000 C 2: 0
dx x L
k = coefficient of thermal conductivity = 200W/m/ C 0

p = perimeter = π d = 0.001 π m
A = cross-sectional area = π d2/4=(π/4)(0.001)(0.001) m2
h = convection heat transfer coefficient = 20W/m 2 0C
Tw = wall temperature ; T∞= ambient temperature = 300C

d 2T
2
 400 T  30
dx

17
Tˆ ( x )
d 2T
 400 T  30  0
dx 2
Assume a trial solution
T(x) = c0 +c1 x + c2x2

BC1 : at x  0 , T  300  c 0  300

dT
BC 2 : 0  c1   2c 2 L
dx x L


T  300  c 2 x 2  2 Lx 
weight function W1  x 2  2 L x  
Residual = Rd = 2c2-400[270+c2(x2-2Lx)]

Apply Galerkin’s method we can write:

 x  x   
L 0.05
2
 2 L x Rd dx   2
 2 L x  2c2  400 270  c2 ( x 2  2 Lx)  dx  0
0 0

solving for c2 we get c2= 38751.43

T= 300 + 38751.43(x2-2Lx 18
 T= 300 + 38751.43(x2-2Lx)

The approximate solution obtained just now is quadratic. We can improve our
approximation by taking higher degree terms in the polynomial trial function.
Let,
T= co + c1x + c2 x2+c3x3

from the boundary conditions we get c0=300, c1= -2c2 L-3c3L2

T=300 + c2(x2-2Lx) + c3(x3-3L2x)

The weight functions are W1= (x2-2Lx) and W2 = (x3-3L2x)

Using Galerkin method we can write

 x  2 L x  Rd dx  0 and  x  3L2 x  Rd dx  0
L L
2 3
0 0

T=300 + 48860.4 (x2-2Lx) -109229 (x3-3L2x)

19
 Trial function Solution
T = c0 +c1 x + c2x2 T= 300 + 38751.43(x2-2Lx)

T= co + c1x + c2 x2+c3x3 T=300 + 48860.4 (x2-2Lx) -109229 (x3-3L2x)

T= co + c1x + c2 x2+c3 x3+c4x4 T= 300 + 53702.8 (x2-2Lx) - 2555668 (x3-3L2x) + 1.32x106(x4-4L3x)

 T T  hp
Exact Solution T  T   w   cosh  m ( L  x) where m 2 
 cosh (m L)  kA
Axial Quadratic Cubic Quartic Exact
Location solution solution Solution Solution
Tα=300C (mm) (0C) (0C) (0C) (0C)

0 300.00 300.00 300.00 300.00

wall x
3000C
10 265.12 264.11 264.00 264.02
L=50mm 20 238.00 237.33 237.39 237.43
40 207.00 208.51 208.43 208.49
50 203.12 205.16 204.91 204.97
20
 Weak form of weighted residual statement

d  du 
We are considering the differential equation EA  qx  0
dx  dx 
L

The weighted residual statement is:  W R d dx  0

0
L
d  du  
0  dx  dx   dx  0
W  EA   q x
where u is the approximate solution and u  u
    
Rd
L L
d  du  
0  dx  dx   dx  0 W qx dx  0
W  EA

dV d
d
 W V   W dV  V dW W   W V   V dW
dx dx dx dx dx dx

Integrating over the interval (0, L)

L
L
dV
L
d dW
L
dW
  V
L
0 dx
W dx  0 dx W V dx  0 dx dx  W V
V 0
dx
dx
0

L L
dV dW
0 dx 0 dx dx
L
W dx  W V 0
 V

21
 L L
dV dW
0 dx 0 dx dx
L
W dx  W V 0
 V

L L
d  du  
0 W  dx  EA dx   dx  0 W qx dx  0
V
L L L
du dW du
 W EA  EA dx   W qx dx  0
dx 0 0
dx dx 0

Boundary terms

Primary variables: Displacements u (x-direction), v (y-direction), w(z-direction)

 du   d3w 
Secondary variables: Axial Force EA  in bar, Shear Force EI 3  in beam,
 dx   dx 
 d w2

Moment  EI dx 2  in beam.
 

Essential/Geometric/Dirichlet: Specification of primary variables at boundary

Natural/Forced/Neumann/dynamic: Specification of secondary variables at boundary

22
 Example 1:
Derive the weak form of the differential equation of a beam under uniformly
distributed load.
z qz
Cross-sectional area of beam: A
Young’s modulus of beam material: E
x
Second moment of area (about neutral axis): I
L
d2  d2w 
The governing differential equation of beam is: 2  EI 2   q z  0
dx  dx 
Weighted residual statement is:
L
 d2  d2w   L
d2  d2w 
L

0 W  dx 2  EI dx 2   qz  dx  0  0 W dx 2  EI dx 2  dx  0 Wqz dx  0

 

d  d  d 2 w 
L L L L
dV dW
0 dx 0 dx dx
L
0 W dx  dx  EI dx 2  dx  0 Wqz dx  0 W dx  W V 0
 V
 

VL
d  d w dW  d  d2w  
2L L
W  EI 2      EI 2   dx -  W q z dx  0
dx  dx  0 0 dx  dx  dx   0
W V
L L
d  d2w 
L L
dW d 2 w d2W d2w
W  EI 2   EI 2   EI 2 dx   Wq z dx  0
dx  dx  0 dx dx 0 0 dx 2 dx 0
23
 L L
d  d2w 
L L
dW d 2 w d2W d2w
W  EI 2   EI 2   EI 2 dx   Wq z dx  0
dx  dx  0 dx dx 0 0 dx 2 dx 0

1 2 3 4
In general , a 2mth order differential equation require m integration by parts to transfer m
derivatives from w to W and therefore there will be m boundary terms involving m primary
variables and m secondary variables.

The weight functions (W1, W2 etc.) can be considered as variation of primary variable (or
virtual displacement, ie. W = δw) so if v is specified at some point, in that point the
variation of w (δw) that is W = 0.

=0

=0
L
d  d2w  d  d2w  d  d2w 
1: W  EI 2    w  EI 2 
dx  dx  0 dx  dx  x  L
 w  EI 2 
dx  dx  x  0
z
qz
Shear force Shear force
L
dW d 2 w dw d 2 w dw d 2 w
x
2: dx
EI 2
dx 0

dx
EI 2
dx x L

dx
EI 2
dx x 0

=0
L
=0
 dw  d2w  dw  d2w
SFD    EI 2    EI 2
X=0 X=L  dx  dx x L  dx  dx x 0
Bending Bending
moment
BMD moment
24
 Approximate solution using weak form

z qz
Cross-sectional area of beam: A
x Young’s modulus of beam material: E
Second moment of area (about neutral axis): I
L
d2  d2w 
Governing differential equation of beam is:  EI   qz  0
dx 2  dx 2 
Boundary conditions are: w  0   0.....(1) w  L   0......(3)

dw dw
 0.....(2)  0......(4)
dx x 0 dx x  L

In weak form the trial function should satisfy only the essential boundary conditions
Assume the approximate solution as : w  a0  a1 x  a2 x 2  a3 x 3  a4 x 4
dw
 a1  2 a2 x  3 a3 x 2  4 a4 x 4
dx
BC1 : a0  0 BC3 : a2 L2  a3 L3  a4 L4  0 a3  2 a4 L
BC 2 : a1  0 BC4 : 2 a2 L  3 a3 L  4 a4 L  0
2 3 a2  a4 L2

Approximation function becomes: w  x   a4 L x  2Lx  x

2 2 3 4
  25
 Approximation function becomes: w  x   a4 L x2 2
 2Lx 3  x 4 
The weak form of the governing equation is: weight function W
L L
d  d2w 
L L
dW d 2 w d2W d2w
W  EI 2   EI 2   EI 2 dx   Wq z dx  0
dx  dx  0 dx dx 0 0 dx 2 dx 0
            
Boundary term Boundary term

L L L L
d  d2w 
L L
2
 dw  d w d2W d2w d2W d2w
w  EI 2     EI 2   2 EI 2 dx   Wq z dx  0  2
EI 2 dx   Wq z dx  0
dx  dx   dx  dx dx dx 0
dx dx 0
       0        0 0 0

Boundary term Boundary term

d2W

W  L x  2L x  x
2 2 3 4
 
dW
dx
 2 L x  6 L x  4 x  2  2 L  12 L x  12 x
2 2 3
dx
2 2

Substituting in weak form we get:

L L

 2 L  12 L x  12 x  EIa4 (2 L  12 L x  12 x ) dx   q z  L2 x 2  2 L x 3  x 4  dx  0
2 2 2 2

0 0

qz
a4 
24EI

w x 
qz
24EI

L2 x 2  2Lx 3  x 4 
26
 Assignment 1L: Derive the week form of the differential equation of the
following beam.

Approximate solution is w x   a0  a1 x  a2 x  a3 x  .......

z 2 3
qz

x Send your assignment to this mail address

A, E, I [email protected]
L

Assignment 1K :Derive the weak form and find the solution. Calculate the deflection at the
free end and middle of the beam.
E=210GPa, A=15cmx15cm
10N
z 20N/m

x 15Nm

L=1m

27