Plane Stress Analysis by FEM

Forces on a differential element subjected to plane stress

Summation of Forces on Differential
Element

fx and fy are body forces in x and y directions, respectively.

PDE of Plane Stress
Constitutive Relation

Stress Vector

Strain Vector

Material Property Matrix

Plane strain
Kinematic Matrix
Fem Formulation
Applying the Galerkin’s method

Integrating by parts the first term

Writing in form of matrix

Fem Formulation
Writing the stress vector in terms of material property matrix and strain vector

Using kinematic relationships to express the strain vector the first term

Using FEM interpolation to express the displacement inside an element

Fem Formulation
Expressing the interpolation in form of a matrix

Expressing the kinematic matrix in terms of interpolation function

derivatives matrix and displacement vector
Fem Formulation
Hence the weak integral form can be expressed as

The element stiffness matrix is given by

The force vector due to the body force

For triangular element B is given as

Triangular Element
Linearly interpolating inside a triangular planer

Expressing in form of a matrix

Expressing nodal values as function of nodal coordinates

Triangular Element
Solving the coefficients of the equations

The area of the triang;eis given by

Triangular Element
Comparing the interpolation with the FEM interpolation the shape unction are
determined as follows
Triangular Element
The matrix [B] is detrmined by substituting for the derivatives of the shape
function
Fem Formulation
As [B] is constant, it can be taken out of the integral sign and the stiffness
matrix is obtained as follows
Axis Symmetric Problem
For an axis symmetric problem in cylindrical coordinate, the stress and
strain vector as given as

The material property matrix foraxis symmetric problems is given as

Axis Symmetric Problem
Kinematic Matrix for an axis symmetric problem

The strain in circumferential direction is given by

Axis Symmetric Problem

The strain vector is given by

Axis Symmetric Problem
Axis Symmetric Problem
Finite Element Method for
dynamic analysis
8.1 FEM Formulation

The Euler-Bernoulli equation for beam vibration is given by

Cont…
The integral form of weighted residual of the above differential equation
over a beam gives

L   2v  2   2v  
I   Wi   2  2  EI 2   q  dx  0 for i= 1 to n Eq.
 t x  x 
o

Where L is the length of the beam and W is the weighting function. The weak
formulation of the above equation is obtained from integration by parts twice for
the second term of the equation. In addition, discritization of the beam into a
number of finite elements gives .
L
L  2v  2 v  2Wi  Wi 
o Wi  t 2 dx  e EI x 2 x 2 dx   e qWi dx  VWi  M x  o  0
Cont ….

W here V   EI  3 v / x 3   
is the shear force, M  EI  2 v / x 2 is the bending
moment,  e is an element domain and n is the number of elements for the beam

v x   c o  c1 x  c 2 x 2  c3 x 3 Eq. 4

From the assumption for the Euler-Bernoulli beam, slope is computed from Eq.

  x   c1  2c 2 x  3c3 x 2 Eq. 5
Cont…
Evaluation of deflection and slope at both nodes yields

v 0   c 0  v1
  0   c1   1
Eq. 6
    c 0  c1  c 2 2  c3 3   2
    c1  2c 2  3c3 2   2

Solving for ci in terms of the nodal variables v i and  i from Eq. 6and substituting
the results into the Eq. 4. gives
v x   N 1  x  1  N 2  x   1  N 3  x  2  N 4  x   2 Eq. 7

W here

3x 2 2x3
N1  x   1  
l2 l3
2x 2 x3
N 2  x  x   2
l l
Eq. 8.
3x 2 2x 3
N 3  x  2  3
l l
x2 x3
N 4  x    2
l l
Cont…
The functions N i  x  are called Hermitian shape functions a which are C 1  type
which means both  and  / x must be continuous between two neighboring
elements..
Application of Hermitian shape functions and Galerkin’s method to Eq 3
L
L  2v  2v  2 Ni  N i 
 o
N i  2 dx   e EI 2
t  x x 2
dx   e
qN i dx 

 VN i  M 
x  o
 0 for

i= 1 to n Eq. 9

 2 v  2 N1 2 N2  2 N3 2N4
 v1  1  v2  2
x 2
x 2
x 2
x 2
x 2

 2v 2N T
{ } {d }  {B}{d } Eq. 10
x 2
x 2
Cont…
 v1 
 
 
{d }   1 
v 2 
 2 

 N1 
N 
 2
{N }   
N 3 
 N 4 

 2v  2 v1  21  2 v2  2 2
 N1  x  2  N 2  x  2  N 3  x  2  N 4  x  2
t 2
t t t t

v  {N }T {d} Eq. 11

Cont….
Substituting equation 10 and 11 in 9 gives

L
o
{N } {N }T dx {d}  
 e
{B} T EI {B} dx {d }    e {N } q{t}i dx Eq. 12

Which results in the stiffness matrix of the beam element

 K     B
e


0
T
EI  B  dx Eq. 13

Where

 B   N1'' N 2 '' N 3 '' N 4''  Eq. 14

The mass matrix is also

 M    {N} {N }
e
L

o
T
dx Eq. 15
Cont…
In Eq 14 double prime denotes the second derivative of the function and l is the
length of a beam element. Assuming the beam rigidity EI is Constant within the
elements, the element stiffness matrix is

12 6l  12 6l 
 
K 
e EI 6l
 3
4l 2
l   12  6l
 6l 2l 2 
12  6l 
Eq. 16.
 
6l 2l 2  6l 4l 2 

The mass matrix becomes

Eq 17
Cont….
The third term results in the element time dependent force vector. F or a generally
distributed pressure loading, we need to compute

 N1 
N 
F   
e 
0
 2
q  x, t  
N
dx Eq. 18.
 3

N 4 

If we have a uniform pressure load qo within the element, the elopement force
vector becomes

 N1  6l 
    2 
N 2  l
F 
e
 qo(t )   dx 
N3 
qo(t )


 Eq. 19.
o 
12 6l 

N 4 
  l 2 
 
The second order ordinary differential eqn obtained by discretizattion of dynamic
analysis of a beam is then

Eq. 20
Modal Analysis
For sinusoidal solution of displacement the above eqn becomes

Eq. 21

Eq. 22

For the above eqn to be satisfied the determinant of the matrix must be zero and
the square of the natural frequency is the eigen value