Finite element method for structural dynamic

and stability analyses

Module-1 Approximate methods and FEM

Lecture-3
Rayleigh Ritz method and method of weighted residuals

Prof C S Manohar
Department of Civil Engineering
IISc, Bangalore 560 012 India

1
 Recall

Exact solutions to PDE-s with inhomogeneous coefficients is in

general not possible even for linear systems.

Approximate solutions
 Natural frequencies and mode shapes
Forced response characteristics
A prelude to developing FEM

2
 Rayleigh’s quotient

Discrete MDOF system Axially vibrating bar

L
u t Ku
R u   t  AE  x    2
 x  dx
u Mu
R   x    0
 R  u  has units of (rad/s) 2 . L

 m  x   2
 x  dx
nt Kn
If u  n  R  u   t  n2 .
0

n M n Euler Bernoulli beam

  R  u   
2 2 L
1 N
 EI  x    2
 x  dx
u  n   y  R  u   n2  O   2  R   x    0
L

 m  x   2
 x  dx
Rayleigh’s principle 0

12  R   x   3
 Example x3
x1 x2
  
k 1.5k 3k k  105 N/m
2m m 1.5m m=1000 kg

For the system shown in the above figure the first eigen vector is obtained
as 1t  1.13 1.63 1.78 and the first natural frequency as 1  4.15 rad/s.
Due to a design modification the third mass is changed by 15%. Estimate the
change in the natural frequency.

Modified system
 2m 0 0   2.5k 1.5k 0 
M   0 m 0  ; K   1.5k
  4.5k 3k 
 0 0 1.5 1.15m   0 3k 3k 
4
Use  t  1.13 1.63 1.78  as the trial vector and estimate the
Rayleigh's quotient as
 t K
R    t  16.10 (rad/s) 2
 M
 MOD  4.0131
 Estimate of the change in natural frequency=
 4.15-4.01 / 4.15 *100  3.3%

Note: For the modified system the exact natural frequency

is 4.0121 rad/s
5
Example M
v  x, t 

L, EI  x  , m  x  K
a

L
1 1
V  t    EI v  x, t   dx  Kv 2  b, t 
2

20 2
L
1 1
T  t    mv  x, t  dx  Mv 2  a, t 
2

20 2
Normal mode oscillations: v  x, t     x  cos t   

6
 1 L 1  2
V  t     EI    x   dx  K  b   cos t   
2 2

2 0 2 

2 1
L
1  2
T  t      m  x  dx  M   a   sin t   
2 2

2 0 2 
1 L 1 
   EI    x   dx  K  b  
2 2
Vmax
2 0 2 

2 1
L
1 
Tmax     m  x  dx  M   a  
2 2

2 0 2 
L

     b 
2
EI x 
  x 
 dx  K  2

R   x     2  0
L

 m  x   2
 x  dx  M  2
a 7
0
Choice of trial function   x 
Mode shape of a simply supported beam with constant elastic and
x
mass properties:   x   sin
L
  x   x  L  x 
Deflected profile under an unit udl

Choice of   x  determines how good R   x   is as

an approximation to 12
12  R   x   8
 L

     b 
2
EI x 
  x 
 dx  K  2

R   x     2  0
L

    
   a
2 2
m x x dx M
0

x
Let EI  x   EI ; m  x   m;   x   sin
L
2
  L 2
2 b
EI  2   K sin
 
2  L  2 L
mL 2 a
 M sin
2 L
2
EI     1  EI
2
 K  0, M  0     2     9.87  2 
2
 OK 
mL  L  m 9
 L

      b 
2
EI x   x  dx  K  2

R   x     2  0
L

 m  x   2
 x  dx  M  2
a
0

Let EI  x   EI ; m  x   m;   x   x  L  x 
    x   2
4 EIL  Kb  L  b  2 2

  2

mL5
 Ma  L  a 
2 2

30
120 EI 1 EI
M  0, K  0    2  10.95 2
L m L m
10
Numerical illustration
L =3 m; E =210 GPa;  =7800 kg/m3 ; a  L / 3; b  3L / 4
Beam cross section (rectangular): width=0.1 m; depth=0.2 m
mL 2 EI
M ;K
2 L
Results on the estimate of the first natural frequency (rad / s)
x 
Beam parameters Exact   x   sin     x  x  L  x
 L 
M  0, K  0 328.51 328.51 364.63
mL 2 EI
M ;K  270.31 297.42
2 L

11
How to lower the value of R   x  ?  Rayleigh - Ritz Method
L

 EI  x    2
 x  dx
R   x    0
L

 m  x   2
 x  dx
0
N
  x    an n  x 
n 1

 n  x  , n  1, 2, , N : a set of known linearly independent functions

which satisfy all the boundary condtions.
an , n  1, 2, , N : a set of unknown constants which need
to be determined
Strategy: Select an , n  1, 2, , N such that R   x   is minimized.
12
 2
 
L N

0   
EI x an n    dx
  x
 A
R   x    R  a1 , a2 , , aN    
n 1
2
2

 
L N B
0 m  x  n 1
an n  x   dx

R
Condtion for optimality:  0 for i  1, 2, ,N
ai
1  A B   A B 
 2 B A  0B A 0
B  ai ai   ai ai 
N  N 
L L
B  2 EI  x    an n  x   i x  dx  A 2m  x    an n  x   i  x  dx  0
0  n 1  0  n 1 
   
L N L N
  EI  x    an n  x   i x  dx   2  m  x    an n  x   i  x  dx  0
0  n 1  0  n 1 
13
    
L N L N

0    n n   i  
          n n   i  x  dx  0

2
EI x a x x dx m x a x
n 1  0  n 1 
Denote
L L
K in   EI  x  n  x  i x  dx; M in   m  x  n  x  i  x  dx
0 0
N N
  K in an   2  M in an  0, i  1, 2, ,N
n 1 n 1

  K a   2  M a

 K  ,  M  are N  N matrices
a   a1 a2 aN  is a N 1 column
t

14
Remarks
 Ka   2 Ma represents an algebraic eigenvalue problem
The formulation leads to a N -dof model for the system
 K  Generalized stiffness matrix; M = Generalized mass matrix
K  K t ; M  M t
N
  x    an n  x   ak  0k  N  Increase in system stiffness
n 1

 Natural frequencies are overestimated.

L
If we select  n  x  such that  m  x  n  x  k  x  dx  nk  M  I
0

 Ka   2 a
The method also provides approximation to the k th eigenfunction as
N
k  x    akn n  x  15
n 1
Ka   2 Ma
1
  x    an n  x     1 1
n 1
2
  x    an n  x     2 1 , 2  2
n 1
3
  x    an n  x     3 1 , 3  2 , 3 3
n 1

N
  x    an n  x     N 1 , N  2 , N
N
n 1

Let n , n  1, 2, be the true (unknown) natural frequencies.

s
 r  r , r  1, 2, ,s
lim s  r  r2  Eigenvalues converge from above
s 
16
Example M
v  x, t 

L, EI  x  , m  x  K
a

b
Numerical illustration
L =3 m; E =210 GPa;  =7800 kg/m3 ; a  L / 3; b  3L / 4
Beam cross section (rectangular): width=0.1 m; depth=0.2 m
mL 2 EI
M ;K
2 L
N
 n x 
  x    an n  x ; n  x   sin  
n 1  L 
17
Ka   2 Ma
N
 n x 
  x    an n  x ; n  x   sin  
n 1  L 
 i   i b   j b 
4

K ij  EI    ij  K sin   sin  
 L  L   L 
mL  i a   j a 
M ij   ij  M sin   sin  
2  L   L 
Estimates of the natural frequencies in rad/s
N  1 270.309
N  2 266.024 1127.975
N  3 265.983 1127.917 2960.057
N  4 265.867 1124.862 2960.056 4624.599
N  5 265.846 1123.404 2960.0545 4549.745 7535.873
18
Approximate analysis of system response:
Method of Weighted Residues (MWR)

Field equation:  EI  x  v  x, t    m  x  v  x, t   c  x  v  x, t   f  x, t 

ICS: v  x, 0   v0  x  , v  x, 0   v0  x 
BCS: Appropriate geometric and natural BCS

We assume an approximate solution of the form

N
v  x, t    an  t  n  x 
n 1

an  t  , n  1, 2, , N : Unknown functions of time to be determined

(generalized coordinates)
n  x  , n  1, 2, , N : Known trial functions
19
 EI  x  v  x, t    m  x  v  x, t   c  x  v  x, t   f  x, t 
N
v  x, t    an  t  n  x 
n 1

 N
 N N
  EI  x   an  t  n  x    m  x   an  t  n  x   c  x   an  t  n  x 
 n 1  n 1 n 1

 f  x, t   e  x, t 
e  x, t   Error or the Residue
For the exact solution e  x, t   0t  0 & x   0, L  .
For an approximate solution this is not so.
Strategy:
Select the unknowns an  t  , n  1, 2, , N such that a specified measure
of total error is minimized in some sense. There is no unique sense of
"minimization". Depending upon the sense of minimization, we have
20
different methods.
Method of least squares
e  x, t  

 N
 N N

 EI  x   an  t  n  x    m  x   an  t  n  x   c  x   an  t  n  x   f  x, t 
 n 1  n 1 n 1
L
Total meansquare error   t    e 2  x, t  dx
0

  t 
Select an  t  , n  1, 2, , N such that  0 for n  1, 2, ,N
an
e  x, t 
L
 e  x, t  dx  0 for n  1, 2, ,N
0
an
L
e  x, t 
 w  x, t  e  x, t  dx  0 for n  1, 2,
0
n , N ; wn  x, t  
an
 Ma  Ca  Ka  P  t 
21
e  x, t  

 N
 N N

 EI  x   an  t  n  x    m  x   an  t  n  x   c  x   an  t  n  x   f  x, t 
 n 1  n 1 n 1

e  x, t    N
 

an

an 

EI  x  
n 1
an  t  n
  x   

  EI  x  n
  x  

L
   EI  x  n  x   e  x, t  dx  0 for n  1, 2, ,N
0

 Ma  Ca  Ka  P  t 
L
M ij    EI  x  i x   m  x  j  x  dx  M ji
0
L L
Cij    EI  x  i x   c  x  j  x  dx  C ji ; Pj  t     EI  x   j  x   f  x, t dx
0 0
L
K ij    EI  x  i x    EI  x   j  x  dx  K ji 22
0
Collocation Method
e  x, t  

 N
 N N

 EI  x   an  t  n  x    m  x   an  t  n  x   c  x   an  t  n  x   f  x, t 
 n 1  n 1 n 1

Select an  t  , n  1, 2, , N such that e  x, t   0 for x =xn ; n  1, 2, , N

L
    x  xn  e  x, t  dx  0 for n  1, 2, ,N
0
L

 w  x, t  e  x, t  dx  0 for n  1, 2,
0
n , N ; wn  x, t     x  xn 

 Ma  Ca  Ka  P  t 

     
x1 x2 xN
23
e  x, t  

 N
 N N

 EI  x   an  t  n  x    m  x   an  t  n  x   c  x   an  t  n  x   f  x, t 
 n 1  n 1 n 1
L

   x  x  e  x, t  dx  0 for n  1, 2,
0
n ,N

 Ma  Ca  Ka  P  t 
L
M ij     x  xi  m  x   j  x  dx m  xi   j  xi   M ji
0
L
Cij     x  xi  C  x   j  x  dx C  xi   j  xi   C ji
0
L
K ij     x  xi   EI  x   j  x   dx   EI  xi   j  xi    K ji
0

Pj  t   f  x j , t 
24
Subdomain Collocation Method
e  x, t  

 N
 N N

 EI  x   an  t  n  x    m  x   an  t  n  x   c  x   an  t  n  x   f  x, t 
 n 1  n 1 n 1

Select an  t  , n  1, 2, , N such that

xn

 e  x, t  dx  0 for n  1, 2,
xn1
, N ; x0  0; xN  L

 w  x, t  e  x, t  dx  0 for n  1, 2,
n , N ; wn  x, t   U  x  xn   U  x  xn 1 
0

 Ma  Ca  Ka  P  t 

x
       x
0 x1 x2 xN 1 N
25
Galerkin's method
e  x, t  

 N
 N N

 EI  x   an  t  n  x    m  x   an  t  n  x   c  x   an  t  n  x   f  x, t 
 n 1  n 1 n 1

Select an  t  , n  1, 2, , N such that

L
  n  x  e  x, t  dx  0 for n  1, 2, ,N
0

The error is orthogonal to the trial functions

L

 w  x, t  e  x, t  dx  0 for n  1, 2,
0
n , N ; wn  x, t   n  x 

 Ma  Ca  Ka  P  t 

26
e  x, t  

 N
 N N

 EI  x   an  t  n  x    m  x   an  t  n  x   c  x   an  t  n  x   f  x, t 
 n 1  n 1 n 1
L

 w  x, t  e  x, t  dx  0 for n  1, 2,
0
n , N ; wn  x, t   n  x 

 Ma  Ca  Ka  P  t 
L
M ij   m  x  i  x   j  x  dx  M ji
0
L
Cij   c  x  i  x   j  x  dx  C ji
0
L
K ij    EI  x  i x    j  x  dx
0
L
Pj  t    f  x, t   j  x  dx 27
0
 L
K ij    EI  x  i x    j  x  dx
0

 
L L
  EI  x  i x    j  x     EI  x  i x    j  x  dx
0 0

   
L L
  EI  x  i x       EI  x   x    x  dx
  EI  x  i x    j  x 
L
j x i j
0
0 0

If we take trial functions to satisfy all the boundary condtions, then

the first two terms in the above expression would be zero.
L
 K ij   EI  x  i x   j  x  dx  K ji
0

28
Initial condtions
N N
v  x, 0    an  0  n  x  & v  x, 0    an  0  n  x 
n 1 n 1

Collocation method
N N
v  xi , 0    an  0  n  xi  & v  xi , 0    an  0  n  xi  ; i  1, 2, ,N
n 1 n 1

 v  xi , 0   n  xi   an  0  & v  xi , 0   n  xi   an  0 

a  0     x  v  x , 0  & a  0     x  v  x , 0 

1 1
n n i i n n i i

Least squares and Galerkin

1
 L
 L 
an  0     k  x  n  x  dx   v  x, 0 k  x  dx 
0  0 
1
 L
 L 
an  0     k  x  n  x  dx   v  x, 0 k  x  dx 
0  0  29
 L
e  x, t 
Least squares  e  x, t  dx  0 for n  1, 2, ,N
0
an
L
Collocation    x  x  e  x, t  dx  0 for n  1, 2,
0
n ,N

L
Galerkin    x  e  x, t  dx  0 for n  1, 2,
0
n ,N

Subdomain
collocation
 U  x  x   U  x  x  e  x, t  dx  0
0
n 1 n

for n  1, 2, ,N
L
Petrov-Galerkin   x  e  x, t  dx  0 for n  1, 2,
0
n ,N
30
 L

 wn  x  e  x , t  dx  0 for n  1, 2, ,N
0 Weight Residue
Method of weighted resuidues

Choices to be made
A set of N trial functions 

A set of N weight functions 
Order to which the functions need to be differentiable
Issues 
Satisfaction of boundary condtions
Value of N
31
Field equation:  EI  x  v  x, t    m  x  v  x, t   c  x  v  x, t   f  x, t 
ICS: v  x, 0   v0  x  , v  x, 0   v0  x 
BCS: Appropriate geometric and natural BCS

N
v  x, t    an  t  n  x 
n 1


Ma  Ca  Ka  P  t  ; a  0 , a  0 

What MWR achieves?

A PDE governing the behavior of a continuous system
has been replaced by an equivalent set of ODE-s (IVP-s)
32
with a view to obtain an approximate solution.
  EI  x  v  x, t    m  x  v  x, t   0
Geometric (I)
Boundary condtions 
 Natural (II)
Normal mode oscillations: v  x, t     x  cos t   

  EI  x     x    m  x   2  x  (III)

Classification of trial functions

Geometric BCs (I) Natural BCs (II) Field Eqn. (III)
Admissible Satisfied Not satisfied Not satisfied
Comparison Satisfied Satisfied Not satisfied
Eigenfunction Satisfied Satisfied Satisfied

33
Rayleigh's quotient
L

    x  dx
  2
EI x
R   x    0
L

    x  dx
 2
m x
0

Trial functions must be at least admissble

and possess nonzero derivatives at least up to the 2nd order.

34
Rayleigh Ritz method
N   
L L N

0 EI  x   an n  x   i  x  dx    m  x    an n  x   i  x  dx  0

  2

n 1  0  n 1 
L L
K in   EI  x  n  x  i x  dx; M in   m  x  n  x  i  x  dx
0 0

 K a  2  M a
Trial functions must be at least admissible and
possess nonzero derivatives at least up to 2nd order.

Least squares, & Collocation

Trial functions need to be comparison functions
and possess nonzero derivatives up to 4 th order.

Galerkin
35
Trial functions need to be comparison functions.
 Assumed mode method and Lagrange's equation

Solving of  EI  x  y  m  x  y  0  appropriate BC-s and IC-s

 
t2 l
1
is equivalent to minimizing A    my  x, t   EI  y  x, t   dxdt
2 2

t1 0
2

Approximate solutions can be developed based on

  EI  x  y  m  x  y  0  appropriate BC-s and IC-s  MWR 

 
t2 l
1
 or by optimizing A    my  x, t   EI  y  x, t   dxdt
2 2

t1 0
2
 EI  x  , m  x  , L
N
Approximate solution: v  x, t    an  t n  x  ;0  x  l
n 1

an  t  ; n  1, 2, , N  unknown generalized coordinates

n  x  ; n  1, 2, , N  known trial functions; taken to satisfy BCS
n x
For example, n  x   sin
l
2
 
l N
1
T  t    m  x   an  t n  x   dx
20  n 1 
2
 
l N
1
V  t    EI  x   an  t n  x   dx

20  n 1 
 1 l
t2
 N

2
1
l
 N

2

A     m  x   an  t n  x   dx   EI  x   an  t n  x   dx dt
20
t1   n 1  20  n 1  

t2

  L  a1  t  , a2  t  , , aN  t  , a1  t  , a2  t  , , a N  t  dt
t1

d  L  L
Governing equations:    0; k  1, 2, ,N
dt  ak  ak
L
l N N
  m  x  an  t n  x  k  x  dx   an  t M nk ;
ak 0 n 1 n 1
l
M nk   m  x n  x  k  x  dx
0

d  L  N
     an  t M nk
dt  ak  n 1
L
l N N
   EI  x  an  t n  x  k  x  dx   an  t K nk ;
ak 0 n 1 n 1
l
K nk   EI  x n  x  k  x  dx
0
N N
  an  t M nk   an  t K nk  0; k  1, 2, ,N
n 1 n 1

Ma  Ka  0; a  0   a0 ; a  0   a0
Ma  Ka  0; a  0   a0 ; a  0   a0
Remarks
M and K are N  N matrices
a   a1  t  a2  t  aN  t   : generalized dof-s;N =dof.
t

M  M t ; K  K t
M is positive definite
 K is positive semidefinite
a   a1  t  a2  t  aN  t   does not have direct physical meaning
t

M and K depend upon choice of trial fucntions n  x  ; n  1, 2, ,N

Again
A PDE governing the behavior of a continuous system
has been replaced by an equivalent set of ODE-s (IVP-s)
with a view to obtain an approximate solution.
Strong (operational) form, Weighted residual form, and
weak (variational) form of governing equations
N
Objective : To find solutions in the form v  x    ann  x 
n 1

Strong (operational) form

Field equation:  EI  x  v  x    Pv  x   k  x  v  x   f  x 

BCS: Appropriate geometric and natural BCS
Weighted residual statement
L
   Pv x  k x v x  f x  dx  0
0 w  x  
  EI  x  v   x  
        

Select a set of n weight functions and obtain governing equations
for an , n  1, 2, , N . Trial functions n  x  n  1, 2, , N need to
be differentiable up to fourth order. No such restrictions on
41
weight functions.
Weak form: Carry out integration by parts in the WR form
L
 
  w  x    EI  x  v  x    Pv  x   k  x  v  x   f  x   dx  0
0  
L
   Pv  w x    EI x v x  w x
          
       
L
 EI x v  x  
  0 0

 w  x   EI  x  v  x   Pv  x  w  x   k  x  v  x  w  x   f  x  w  x  dx  0

0

 F2 w  L   F1w  0   M 2 w  L   M 1w  0  
L

 w  x   EI  x  v  x   Pv  x  w  x   k  x  v  x  w  x   f  x  w  x  dx  0

0

M 1  EI  0  v  0  ; M 2  EI  L  v  L 
   
F1    EI  x  v  x    Pv  @ x  0; F2    EI  x  v  x    Pv  @ x  L
   
42
Illustrations: Let the beam be simply supported at x  0 and x  L
and be subjected to end moments M 1 and M 2 .
F2 w  L   F1w  0   M 2 w  L   M 1w  0  

 w  x   EI  x  v  x   Pv  x  w  x   k  x  v  x  w  x   f  x  w  x  dx  0

L

M 1  EI  0  v  0  ; M 2  EI  L  v  L 
   
F1    EI  x  v  x    Pv  @ x  0; F2    EI  x  v  x    Pv  @ x  L
   

 w  x   EI  x  v  x   Pv  x  w  x   k  x  v  x  w  x   f  x  w  x  dx

L

 M 2 w  L   M 1w  0   0
Select   x  such that v  0   v  L   0 & w  L   w  0   0
43
Consider the term
L
       Pv   L  w  L  

  EI  x  v   x  
  Pv   w  x     EIv   
  0  
Primary
Secondary variable Variable=v  L 

 EIv   Pv  0 w 0
     
 
Primary
Secondary variable Variable=v  0 

Similarly

 EI  x  v  x  w  x    EI  L  v  L   w  L    EI  0  v  0   w  0 

L

0
Secondary variable Primary Secondary variable Primary
Variable=v 0  Variable=v ' 0 

The dependent variable expressed in the same form as the weight function
appearing in the boundary term= primary variables
Coefficients of weight functions: secondary variables
44
Remarks
In the weighted residual statement the trial functions n  x  , n  1, 2, ,N
must be differentiable up to the fourth order. w( x) can be any integrable
function. Equations for the undetermined coefficients can be obtained by
choosing a set of N weight functions.

In the weak form we distribute the requirements on differentiabilty evenly

between weight function and the trial function. Thus the continuity
requirements on the trial functions is "weakened" and hence the name
weak form. We have larger pool of functions to select from for the
trial functions.

45
Remarks (continued)

The trial functions need to satisfy only the geometric boundary condtions
as the natural boundary condtions are included in the weak form.

 Number of primary and secondary variables are the same. They appear
in pairs (e.g., translation and shear force; slope and bending moment).
Only one item in each pair can be specified at the boundary.

46
 Limitations
N
v  x, t    an  t  n  x 
n 1

n  x  ,n  1, 2, , N are global in nature, that is, they are

valid for all x   0, L .

Constructing these functions for simple geometries is

relatively easy. Not so when geometries become more complicated.

an  t  , n  1, 2, , N are generalized coordinates and they do not

have direct physical meaning.
47
48