THE GENERALIZED COORDINATE METHOD FOR DISCRETE SYSTEMS Revision F By Tom Irvine

THE GENERALIZED COORDINATE METHOD FOR DISCRETE SYSTEMS
Revision F
By Tom Irvine
Email: [email protected]
January 18, 2012
_______________________________________________________________________
Two-degree-of-freedom System
The method of generalized coordinates is demonstrated by an example. Consider the
system in Figure 1.
f2(t)
k3
x2
m2
k2
f1(t)
x1
m1
k1
Figure 1.
A free-body diagram of mass 1 is given in Figure 2. A free-body diagram of mass 2 is
given in Figure 3.
1
k2 (x2-x1)
f1(t)
x1
m1
k1x1
Figure 2.
Determine the equation of motion for mass 1.
 F  m1 x1
(1)
m1x1  f1 (t )  k 2 (x 2  x1 )  k1x1
(2)
m1x1  k 2 (x 2  x1 )  k1x1  f1 (t )
(3)
m1x1  k 2 (x 2  x1 )  k1x1  f1 (t )
(4)
m1x1  (k1  k 2 )x1  k 2 x 2  f1 (t )
(5)
2
f2(t)
k3 x2
x2
m2
k2 (x2-x1)
Figure 3.
Derive the equation of motion for mass 2.
 F  m 2 x 2
(6)
m 2 x 2  f 2 (t )  k 2 (x 2  x1 )  k 3x 2
(7)
m 2 x 2  k 2 (x 2  x1 )  k 3x 2  f 2 (t )
(8)
m 2 x 2  (k 2  k 3 )x 2  k 2 x1  f 2 (t )
(9)
Assemble the equations in matrix form.
m1 0   x1  k1  k 2
 0 m  x     k
2
2 2 

 k 2   x1   f1 ( t ) 

k 2  k 3  x 2  f 2 ( t )
(10)
Decoupling
Equation (10) is coupled via the stiffness matrix. An intermediate goal is to decouple the
equation.
3
Simplify,
M x  K x  F
(11)
where
0 
m
M 1

 0 m2 
(12)
k  k 2
K 1
  k2
(13)
 k2 
k 2  k 3 
x 
x   1
x 2 
(14)
 f (t) 
F 1 
f 2 ( t ) 
(15)
Consider the homogeneous form of equation (11).
M x  K x  0
(16)
Seek a solution of the form
x  q exp( jt )
(17)
The q vector is the generalized coordinate vector.
Note that
x  j q exp( jt )
(18)
x  2 q exp( jt )
(19)
4
Substitute equations (17) through (19) into equation (16).
 2 M q exp( jt )  K q exp( jt )  0
(20)
 2 M q  K q exp( jt)  0
(21)
 2 M q  K q  0
(22)
 2 M  Kq  0
K  2 Mq  0
(23)
(24)
Equation (24) is an example of a generalized eigenvalue problem. The eigenvalues can
be found by setting the determinant equal to zero.


det K  2 M  0
k  k 2
det  1
  k 2
(25)
 k2 
0 
m
 2  1

  0
k 2  k3 
 0 m 2 

k  k 2   2 m1
det  1
 k2




  0
2
k 2  k 3    m 2  
 k2
k1  k 2   2m1k 2  k3   2m2  k 22  0
(26)
(27)
(28)
5
k1  k 2 k 2  k 3   2 m1k 2  k 3   2 m 2 k1  k 2   4 m1m 2  k 2 2  0
(29)
m1m 2 4   m1 k 2  k 3   m 2 k1  k 2 2  k1k 3  k1  k 3 k 2  k 2 2  k 2 2  0
(30)
m1m 2 4   m1 k 2  k 3   m 2 k1  k 2 2  k1k 3  k1k 2  k 2 k 3  0
(31)
The eigenvalues are the roots of the polynomial.
1 
2
 b  b 2  4ac
2a
(32)
2
 b  b 2  4ac
2a
(33)
2 
where
a  m1m 2
(34)
b   m1 k 2  k 3   m 2 k1  k 2 
(35)
c  k1k 2  k1k 3  k 2 k 3
(36)
The eigenvectors are found via the following equations.
K  12 Mq1  0
(37)
K  22 Mq2  0
(38)
6
where
v 
q1   1 
v 2 
(39)
w 
q2   1 
w 2 
(40)
An eigenvector matrix Q can be formed. The eigenvectors are inserted in column
format.
Q  q1 q 2 
v
Q 1
v 2
w1 
w 2 
(41)
(42)
The eigenvectors represent orthogonal mode shapes.
Each eigenvector can be multiplied by an arbitrary scale factor. A mass-normalized
eigenvector matrix Q̂ can be obtained such that the following orthogonality relations
are obtained.
Q̂ T M Q̂  I
(43)
and
Q̂ T K Q̂  
(44)
where
I
is the identity matrix

is a diagonal matrix of eigenvalues
The superscript T represents transpose.
7
Note the mass-normalized forms
 v̂
Q̂   1
 v̂ 2
 v̂
Q̂ T   1
 ŵ1
ŵ1 
ŵ 2 
v̂ 2 
ŵ 2 
(45)
(46)
Rigorous proof of the orthogonality relationships is beyond the scope of this tutorial.
Further discussion is given in References 1 and 2.
Nevertheless, the orthogonality relationships are demonstrated by an example in this
tutorial.
Now define a generalize coordinate ( t ) such that
x  Q̂ 
(47)
Substitute equation (47) into the equation of motion, equation (11).
  K Q̂   F
M Q̂ 
(48)
Premultiply by the transpose of the normalized eigenvector matrix.
  Q̂ T K Q̂   Q̂ T F
Q̂ T M Q̂ 
(49)
The orthogonality relationships yield
     Q̂ T F
I
(50)
8
The equations of motion along with an added damping matrix become
 1  21 1
1 0  

0 1 

  2   0
 2
   1  1

2 2  2   2   0

0
0    1   v̂1


2     ŵ
2   2   1
v̂ 2   f1 ( t ) 
ŵ 2  f 2 ( t )
(51)
Note that the two equations are decoupled in terms of the generalized coordinate.
Equation (51) yields two equations
 1 21 1
 1  2 1 v̂1f1 (t )  v̂ 2 f 2 ( t )

1
(52)
 2  2 2 2 
 2   2  2  ŵ1f1 (t )  ŵ 2 f 2 (t )

2
(53)
The equations can be solved in terms of Laplace transforms, or some other differential
equation solution method.
Now consider the initial conditions. Recall
x  Q̂ 
(54)
Thus
x(0)  Q̂ (0)
(55)
Premultiply by Q̂ T M.
Q̂ T M x(0)  Q̂ T MQ̂ (0)
(56)
Recall
Q̂ T M Q̂  I
(57)
9
Q̂ T M x(0)  I (0)
(58)
Q̂ T M x(0)  (0)
(59)
Finally, the transformed initial displacement is
(0)  Q̂ T M x(0)
(60)
Similarly, the transformed initial velocity is
 (0)  Q̂ T M x (0)

(61)
A basis for a solution is thus derived.
Harmonic Force
Now consider the special case of harmonic forcing functions
f1 (t )  B1 sint 
(62)
f 2 (t )  B2 sint 
(63)
Thus,
 1 21 1
 1  2 1 v̂1B1 sint   v̂ 2 B2 sint 

1
(64)
 2  2 2 2 
 2   2  2  ŵ1B1 sin t   ŵ 2 B2 sin t 

2
(65)
Take the Laplace transform of equation (64).
 1 21 1
 1  2 1 v̂1B1 sint   v̂ 2 B2 sint 

1
(66)

(67)

 1 21 1
 1  1  Lv̂1B1 sin t   L{v̂ 2 B2 sin t }
L
1
2
10
s 2 ˆ 1(s)  s1 (0)   1(0)
 21 1sˆ 1(s)  21 11 (0)
  
  
2
 1 ˆ 1 (s)  v̂1B1 
 v̂ 2 B 2 


 s 2   2 
s 2   2 
(68)
s
2

2
 21 1s  1 ˆ 1 (s)  s  21 11(0)   1(0)
  
  
 v̂1B1 
 v̂ 2 B 2 


 s 2   2 
s 2   2 
(69)
s
2

2
 21 1s  1 ˆ 1 (s)

ˆ 1 (s) 


  
  

v̂1B1 
 v̂ 2 B 2 

  s  21 11 (0)   1 (0)
2
2
2
2


s


s   


(70)
  
v̂ 2 B 2



2
2
2
2
s 2  21 1s  1  s    s 2  21 1s  1
v̂1B1




  
 2

 s   2 
s  21 11 (0)   1 (0)
s
2
2
 21 1s  1

(71)
11
The solution is found via References 3 and 4. The inverse Laplace transform for the
first modal coordinate is
 1( t ) 

  2   cos t    2   2  sin t 


1 1
1 


 2

2 2
2
   1   21  1  



v̂1B1


 

v̂1B1
exp  11t 
 
 d,1 

 sin d,1t 

  2  cos  t   2   2  sin  t 


1 1
1 


v̂ 2B2
 2
2 2
2 




  21  1  

1





2


211 d,1  cos d,1t    2  1 1  212


 
 
 2
2 2
2 






2






1
1
1 






 

v̂ 2B2
exp  11t 
 
2


 d,1 
211 d,1  cos d,1t   2  1 1  212


 
 
 2
2 2
2 
   1   21  1  




 sin d,1t 



  (0)  11 1(0) 
 exp  11t 1(0) cos d,1t   1
sin

t

d,1 
d,1








(72)
12
Similarly,
 2( t ) 

 2
2 2
   2 




 

ŵ1B1
exp   22 t 



2

 d,2 
2 
 
  2

2 22 d,2  cos d,2 t     2 1  2 2  sin  d,2 t 

 
 

 2
 
2 2
2 
   2   2 2  2  






1 2


2
  2 22  cos t     2  sin t 




 2 2  2 2 

ŵ1B1


1 2
2
  2 22 cos t     2  sin t 




ŵ 2 B2
 2
2 2
2 
   2   2 2  2  








 

ŵ 2 B2
exp   22 t 



2

 d,2 
2
 
  2
2 22 d,2  cos d,2 t     2 1  2 2
2




 2


2
2 
   2   2 2  2  






  (0)   2 2  2(0) 
 exp   22 t  2(0) cos d,2 t   2
 sin d,2 t
d,2






 sin d,2t 




(73)
The physical displacements are found via
x  Q̂ 
(74)
13
References
1. Bathe, Finite Element Procedures in Engineering Analysis, Prentice-Hall, New
Jersey, 1982. Section 12.3.1.
2. Weaver and Johnston, Structural Dynamics by Finite Elements, Prentice-Hall,
New Jersey, 1987. Chapter 4.
3. T. Irvine, Table of Laplace Transforms, Vibrationdata, 2000.
4. T. Irvine, Partial Fraction Expansion, Rev F, Vibrationdata, 2010.
_______________________________________________________________________
APPENDIX A
Example
Consider the system in Figure 1 with the values in Table A-1.
Assume 5% damping for each mode. Assume zero initial conditions.
Table A-1. Parameters
Variable
Value
Unit
m1
3.0
lbf sec^2/in
m2
2.0
lbf sec^2/in
k1
400,000
lbf/in
k2
k3
300,000
lbf/in
100,000
lbf/in
B1
100
lbf
B2
200
lbf

55
Hz

100
Hz
The mass matrix is
0  3 0 
m
M 1


 0 m 2  0 2
(A-1)
14
The stiffness matrix is
k  k 2
K 1
  k2
 k 2   700,000  300,000

k 2  k 3   300,000 400,000 
(A-2)
The analysis is performed using a Matlab script.
>> twodof_sine_force
twodof_sine_force.m
ver 1.4 January 17, 2012
by Tom Irvine Email: [email protected]
This script calculates the response of a two-degree-of-freedom
system to sinusoidal force excitation.
Enter the units system
1=English 2=metric
1
Assume symmetric mass and stiffness matrices.
Select input mass unit
1=lbm 2=lbf sec^2/in
1
stiffness unit = lbf/in
Select file input method
1=file preloaded into Matlab
2=Excel file
1
Mass Matrix
Enter the matrix name:
Stiffness Matrix
Enter the matrix name:
Input Matrices
m_two
k_two
mass =
3
0
0
2
stiff =
700000
-300000
-300000
400000
Natural Frequencies
No.
f(Hz)
15
1.
2.
48.552
92.839
Modes Shapes (column format)
ModeShapes =
0.3797
0.5326
-0.4349
0.4651
Enter the damping ratio for mode 1 0.05
Enter the damping ratio for mode 2 0.05
Particpation Factors =
2.204
-0.3746
Enter the first force amplitude (lbf)
Enter the first force frequency (Hz)
Enter the second force amplitude (lbf)
Enter the second force frequency (Hz)
100
55
200
100
Enter the first initial displacement (in)
Enter the second initial displacement (in)
0
0
Enter the first initial velocity (in/sec)
Enter the second initial velocity (in/sec)
0
0
Enter the sample rate (samples/sec) 10000
Enter the duration (sec) 0.3
dof 1 displacement (in)
max= 0.001015
min= -0.001373
dof 2 displacement (in)
max= 0.002139
min= -0.00166
16
DISPLACEMENT
0.003
Mass 2
Mass 1
0.002
DISP (INCH)
0.001
0
-0.001
-0.002
-0.003
0
0.05
0.10
0.15
0.20
0.25
0.30
TIME (SEC)
Figure A-1.
17
RELATIVE DISPLACEMENT
(MASS 2 - MASS 1)
0.003
REL DISP (INCH)
0.002
0.001
0
-0.001
-0.002
-0.003
0
0.05
0.10
0.15
0.20
0.25
0.30
TIME (SEC)
Figure A-2.
18
VELOCITY
1.5
Mass 2
Mass 1
1.0
VEL (IN/SEC)
0.5
0
-0.5
-1.0
-1.5
0
0.05
0.10
0.15
0.20
0.25
0.30
TIME (SEC)
Figure A-3.
19
ACCELERATION
1.5
Mass 2
Mass 1
1.0
ACCEL (G)
0.5
0
-0.5
-1.0
-1.5
0
0.05
0.10
0.15
0.20
0.25
0.30
TIME (SEC)
Figure A-4.
20