QUESTION 1
Consider three element model of fixed-free cantilever bar undergoing axial
vibration, find frequencies and mode shapes. 𝐿 = 1 m, 𝐴 = 30 × 10−6 m2, 𝐸 = 2 ×
1011 N/m2, 𝜌 = 7800 kg/m3.
FIGURE 1:
Solution
Characteristic stiffness of each
element
𝐴𝐸 3𝐴𝐸
𝑘=
= 𝐿
𝐿/3
So, element stiffness matrices are
[𝑘( 1)] = [𝑘( 2)] = [𝑘( 3)] 3𝐴𝐸
= ]
1 −1
[
𝐿 −1 1
Mass of each element is
𝜌𝐴𝐿
𝑚=
3
and the element consistent mass matrix is
[𝑚( 1)] = [𝑚( 2)] = [𝑚( 3)]𝜌𝐴𝐿
= ]
1 0
[
3 0 1
Using direct stiffness matrix
method, Global stiffness matrix,
1 −1 0 0
⎡ ⎤
3𝐴𝐸 ⎢−1 2 −1 0 ⎥
[𝑘] =
𝐿 ⎢ 0 −1 2 −1⎥⎥
⎢
⎣ 0 0 −1 1 ⎦
Global mass matrix,
1 0 0 0
⎡ ⎤
𝜌𝐴𝐿 ⎢0 2 0 0⎥
[𝑚] =
3 ⎢⎢0 0 2 0⎥⎥
⎣0 0 0 1⎦
Then we have characteristic equation,
[[𝑘] − 𝜔2[𝑚] ] {v̂} = {0}
1
Substituting and boundary condition, i.e. v̂1 = 0,
2 −1 0 2 0 0 ⎧ v̂ ⎫ ⎧0⎫
⎡ 3𝐴𝐸 ⎡ ⎤ 𝜔2𝜌𝐴𝐿 ⎡ ⎤⎤ { 2 } { }
⎢⎢ ⎢−1 2 −1 ⎥⎥ − ⎢⎢0 2 0⎥⎥⎥⎥ ⎨ v̂3 ⎬= ⎨ 0 ⎬
𝐿 ⎢ 3 { } { }
⎣ ⎣ 0 −1 1 ⎦ ⎣0 0 1⎦⎦ ⎩ v̂4 ⎭ ⎩ 0 ⎭
To obtain a solution to the set of homogeneous equation, we set the determinant of the
coefficient matrix to zero,
∣⎡ 2 −1 0 ⎡
2 0 0 ∣
⎤
∣⎢−1 2 −1 ⎤ − 𝛽 ⎢⎢0 2 0⎥⎥∣∣ = 0
∣⎢ ⎥⎥
∣ 0 −1 1 ∣
∣⎣ ⎦ ⎣0 0 1⎦∣
Where 𝛽 = (𝜔2𝜌𝐿2)/𝐸. On solving we get,
𝛽1 = 0.134,
𝛽2 = 1,
𝛽3 = 1.866.
Hence, the natural frequencies
will be
𝜔1 = 5560.850 rad/s,
𝜔2 = 15191.090 rad/
s, 𝜔3 = 20751.274
rad/s.
Mode shapes
Calculating the mode shapes based on the calculated natural
frequencies, TABLE 1: Mode shape for different modes.
Mode 1 Mode 2 Mode 3
𝜙2 0.500 -1.000 0.500
𝜙3 0.866 0 -0.866
𝜙4 1.000 1.000 1.000
2
QUESTION 2
Find natural frequencies and mode shapes for the truss shown in the Fig. 2.. 𝜌 =
2.6×10−4 lb − s2/in4, 𝐴 = 1.5 in2, 𝐸 = 10 × 106 psi.
FIGURE 2:
Solution
The element and global mass matrices for the bar element in two dimensions
are given by
FIGURE 3: Nodes, members, and displacement coordinates.
2 0 1 0
⎡ ⎤
𝜌𝐴𝐿 ⎢ 0 2 0 1⎥
[𝑚( 𝑒)] =
6 ⎢⎢1 0 2 0⎥⎥
⎣0 1 0 2⎦
3
As elements 1, 3, 4, 5, 7, and 8 have the same length, area, and density, we have
[𝑀 ( 1)] = [𝑀 ( 3)] = [𝑀 ( 4)]
= [𝑀 ( 5)] = [𝑀 ( 7)] = [𝑀 ( 8)]
2 0 1 0
⎡ ⎤
(2.6)(10)−4(1.5) ⎢0 2 0 1⎥
= ⎢1 0 2 0⎥
(40) 6 ⎢ ⎥
⎣0 1 0 2⎦
5.2 0 2.6 0
⎡ ⎤
⎢ 0 5.2 0 2.6⎥
=⎢ ⎥ (10)−3 lb ⋅ s2/in
⎢2.6 0 5.2 0 ⎥
⎣ 0 2.6 0 5.2⎦
while for elements 2 and 6
2 0 1 0
⎡ ⎤
2.6(10)−4(1.5) ⎢0 2 0 1⎥
[𝑀 ( 2)] = [𝑀 ( 6)] = ⎢1
(40√2) 6 ⎢ 0 2 0⎥⎥
⎣0 1 0 2⎦
7.36 0 3.68 0
⎡ ⎤
⎢ 0 7.36 0 3.68⎥
=⎢ ⎥ (10)−3 lb ⋅ s2/in
⎢3.68 0 7.36 0 ⎥
⎣ 0 3.68 0 7.36⎦
Using the direct assembly procedure, the global mass matrix is
12.56 0 0 0 2.6 0 3.68 0 0 0 0 0
⎡ ⎤
⎢ 0 12.56 0 0 2.6 0 3.68 0 0 0 0 0 ⎥
⎢ 0 0 5.2 0 0 0 2.6 0 0 0 0 0 ⎥⎥
⎢
⎢ 0 0 5.2 0 0 0 2.6 0 0 0 0 0 ⎥
⎢ 2.6 0 0 0 7.8 0 2.6 0 2.6 0 0 0 ⎥
⎢ ⎥
⎢ 0 2.6 0 0 7.8 0 2.6 0 2.6 0 0 0 ⎥
[𝑀] =⎢ (10)−3
⎢ 3.68 0 2.6 0 2.6 0 22.52 0 3.68 0 2.6 0 ⎥⎥
⎢ 0 3.68 0 2.6 0 2.6 22.52 0 3.68 0 2.6 0 ⎥
⎢ ⎥
⎢ 0 0 0 0 2.6 0 3.68 0 17.76 0 2.6 0 ⎥
⎢ 0 0 0 0 2.6 0 3.68 0 17.76 0 2.6 0 ⎥
⎢ ⎥
⎢ 0 0 0 0 2.6 0 2.6 0 10.4 0 10.4 0 ⎥
⎣ 0 0 0 0 0 0 0 0 2.6 0 2.6 10.4⎦
Applying the constraint conditions 𝑈1 = 𝑈2 = 𝑈3 = 𝑈4 = 0, the mass matrix for the active
degrees of freedom becomes
7.8 0 2.6 0 2.6 0 0 0
⎡ ⎤
⎢ 0 7.8 0 2.6 0 2.6 0 0 ⎥
⎢2.6 0 22.52 0 3.68 0 2.6 0 ⎥⎥
⎢
⎢ 0 2.6 0 22.52 0 3.68 0 2.6 ⎥
[𝑀𝑎] =⎢
2.6 0 3.68 0 17.76 0 2.6 0 ⎥ (10)−3 lb ⋅ s2/in
⎢ ⎥
⎢ 0 2.6 0 3.68 0 17.76 0 2.6 ⎥
⎢ ⎥
⎢0 0 2.6 0 2.6 0 10.4 0 ⎥
⎣0 0 0 2.6 0 2.6 0 10.4⎦
4
The stiffness matrix for the active degrees of freedom is
7.5 0 0 0 −3.75 0 0 0
⎡ ⎤
⎢ 0 3.75 0 −3.75 0 0 0 0 ⎥
⎢ 0 0 10.15 0 −1.325 1.325 −3.75 0 ⎥
⎢ ⎥
⎢ 0 −3.75 0 6.4 1.325 −1.325 0 0 ⎥
[𝐾𝑎] =⎢ ×105 lb/in.
⎢−3.7 0 −1.325 1.325 5.075 −1.325 0 0 ⎥⎥
⎢ 50 1.325 −1.325 −1.325 5.075 0 −3.75 ⎥
⎢ ⎥
⎢ 0 0 −3.75 0 0 0 0
3.75 0 ⎥
⎣ 0 0 0 0 0 −3.7 0 3.75⎦
5
The finite element model for the truss exhibits 8 degrees of freedom; hence, the
characteristic deter-minant
| − 𝜔2[𝑀] + [𝐾] | = 0
yields, theoretically, eight natural frequencies of oscillation and eight corresponding
mode shapes (modal amplitude vectors). The natural frequencies are shown in Table 2.
The corresponding modal amplitude vectors (normalized to the mass matrix as discussed
relative to orthogonality) are shown in Table 3.
TABLE 2: Natural frequencies for different modes.
Mode Frequency (rad/s)
1 767.1
2 2082.3
3 2958.7
4 4504.8
5 6790.9
6 7975.9
7 8664.5
8 8977.4
TABLE 3: Mode shapes for different modes.
Mode 1 2 3 4 5 6 7 8
𝜙5 -0.0385 -0.298 -0.1724 -0.4231 -0.193 -0.2092 0.5028 0.4733
𝜙6 -0.326 4 -0.4442 -0.252 0 -0.1085 0.1115 -0.078
𝜙7 6 0.4464
-0.097 -0.4105 0 0.7782
-0.0190 -0.444 -0.2290 8
𝜙8 0.1148
-0.3148 5 -0.2781 0.2789
-0.0509 -0.5709 6 -0.0360 0.0503
-0.0523
𝜙9 -0.0763 0.3653
-0.5243 -0.242 -0.463 -0.007 0.066 -0.2930 -0.428
𝜙10 -0.609 -0.3476 4 1 3 1 -0.349 8
𝜙11 5 -0.1321 0.2074
-0.594 0.133 0.121 0.230 7 0.477
𝜙12 0.1168
-0.6234 -0.3945 4 4 3 0
-0.3948 0.462 7
-0.5792
0.2683 0.578 0.098 0.189 7 0.130
7 4 1 0.508 6
0.331 0.080 0.705 7
6 6 6
5
QUESTION 3
For the quadrilateral element shown generate consistent mass matrix and calculate natural frequencies
and mode shapes. 𝜌 = 7.83 × 10−6 kg/mm3 thickness is 5 mm.
FIGURE 4:
We have shape functions,
1
𝑁1 (𝑟, 𝑠) = (1 − 𝑟)(1 − 𝑠)
4
1
𝑁2 (𝑟, 𝑠) = (1 + 𝑟)(1 − 𝑠)
4
1
𝑁3 (𝑟, 𝑠) = (1 + 𝑟)(1 + 𝑠)
4
1
𝑁1 (𝑟, 𝑠) = (1 − 𝑟)(1 + 𝑠)
4
Where,
𝑟 = (𝑥 − 25)/15 𝑑𝑥 = 15 𝑑𝑟,
𝑠 = (𝑦 − 20)/10 𝑑𝑦 = 10 𝑑𝑠.
1 1
𝑚11 = ∭ [𝑁]T [𝑁]𝜌 𝑑𝑉 (𝑒) = 𝜌𝑡 ∫ ∫ [𝑁]T [𝑁] 15𝑑𝑟 10𝑑𝑠
𝑉 (𝑒) −1 −1
150 × 5 1 1
= 𝜌 ∫ ∫ (1 − 𝑟)2 (1 − 𝑠)2 𝑑𝑟 𝑑𝑠
16 −1 −1
−3
= 2.6 × 10 kg
150 × 5 1 1
𝑚12 = 𝜌 ∫ ∫ 𝑁1 𝑁2 𝑑𝑟 𝑑𝑠
16 −1 −1
−3
= 1.3 × 10 kg
Similarly for other elements,
2.6 1.3 0.7 1.3 0 0 0 0
⎡ ⎤
⎢1.3 2.6 1.3 0.7 0 0 0 0⎥
⎢0.7 1.3 2.6 1.3 0 0 0 0 ⎥⎥
⎢
(𝑒) ⎢1.3 0.7 1.3 2.6 0 0 0 0⎥
[𝑚 ] = ⎢ × 10−3 kg
⎢0 0 0 0 2.6 1.3 0.7 1.3⎥⎥
⎢0 0 0 0 1.3 2.6 1.3 0.7⎥
⎢ ⎥
⎢0 0 0 0 0.7 1.3 2.6 1.3⎥
⎣0 0 0 0 1.3 0.7 1.3 2.6⎦
6
Calculated natural frequencies using Matlab,
𝜔1 = 0.0000 × 1000 + 1.7034 × 1004 𝑖 rad/s
𝜔1 = 0.0000 × 1000 + 9.0479 × 1005 𝑖 rad/s
𝜔1 = 2.6914 × 1004 rad/s
𝜔1 = 1.8043 × 1004 rad/s
𝜔1 = 2.1041 × 1004 rad/s
𝜔1 = 2.5880 × 1004 rad/s
𝜔1 = 2.9242 × 1004 rad/s
𝜔1 = 3.0880 × 1004 rad/s
Mode shapes,
TABLE 4: Mode shapes for different modes.
Mode 1 2 3 4 5 6 7 8
𝜙1 -0.193 0.5718 0.0339 0.4819 -0.416 -0.5 -0.1335 0
𝜙2 -0.092 -0.2018 -0.6288 -0.1335 -0.2774 0 -0.4819 0.5
𝜙3 -0.193 0.5718 0.0339 -0.4819 -0.416 0.5 0.1335 0
𝜙4 0.4151 0.2018 -0.1995 -0.1335 0.2774 0 -0.4819 -0.5
𝜙5 -0.531 0.3027 -0.2523 -0.4819 0.416 -0.5 0.1335 0
𝜙6 0.4151 0.2018 -0.1995 0.1335 0.2774 0 0.4819 0.5
𝜙7 -0.531 0.3027 -0.2523 0.4819 0.416 0.5 -0.1335 0
𝜙8 -0.092 -0.2018 -0.6288 0.1335 -0.2774 0 0.4819 -0.5
Source Code
%% Assignment 7, Q.3
% Given Data
h=5
rho=7.83e-6
E=2e5
NU=0.25
% Element Mass Matrix
disp('The Mass Matrix is:')
m=BilinearQuadElementMass(rho,h,10,10,40,10,40,30,10,30)
% Element Stifffness Matrix
disp('The Stifness Matrix is:')
k=BilinearQuadElementStiffness(E,NU,h,10,10,40,10,40,30,10,30,1)
% Calculation of Natural Frequencies and Mode Shapes
% Form the system matrix
A=m\k;
% Obtain eigenvalues and eigenvectors of A
[V,D]=eig(A);
% V and D above are matrices.
% V-matrix gives the eigenvectors and
% the diagonal of D-matrix gives the eigenvalues
% Sort eigen-values and eigen-vectors
[D_sorted, ind] = sort(diag(D),'ascend');
V_sorted = V(:,ind);
%Obtain natural frequencies and mode shapes
nat_freq_1 = sqrt(D_sorted(1))
7
nat_freq_2 = sqrt(D_sorted(2))
nat_freq_3 = sqrt(D_sorted(3))
nat_freq_4 = sqrt(D_sorted(4))
nat_freq_5 = sqrt(D_sorted(5))
nat_freq_6 = sqrt(D_sorted(6))
nat_freq_7 = sqrt(D_sorted(7))
nat_freq_8 = sqrt(D_sorted(8))
mode_shape_1 = V_sorted(:,1)
mode_shape_2 = V_sorted(:,2)
mode_shape_3 = V_sorted(:,3)
mode_shape_4 = V_sorted(:,4)
mode_shape_5 = V_sorted(:,5)
mode_shape_6 = V_sorted(:,6)
mode_shape_7 = V_sorted(:,7)
mode_shape_8 = V_sorted(:,8)
% Function for Defining Element Mass Matrix
function w = BilinearQuadElementMass(rho,h,x1,y1,x2,y2,x3,y3,x4,y4)
syms s t;
N1 = (1-s)*(1-t)/4;
N2 = (1+s)*(1-t)/4;
N3 = (1+s)*(1+t)/4;
N4 = (1-s)*(1+t)/4;
N = [N1 0 N2 0 N3 0 N4 0 ; 0 N1 0 N2 0 N3 0 N4];
Jfirst = [0 1-t t-s s-1 ; t-1 0 s+1 -s-t ;
s-t -s-1 0 t+1 ; 1-s s+t -t-1 0];
J = [x1 x2 x3 x4]*Jfirst*[y1 ; y2 ; y3 ; y4]/8;
NN = J*transpose(N)*N;
r = int(int(NN, t, -1, 1), s, -1, 1);
z = rho*h*r;
w = double(z);
end
% Function for Defining Element Stiffness Matrix
function w = BilinearQuadElementStiffness(E,NU,h,x1,y1,x2,y2,x3,y3,x4,y4,p)
syms s t;
a = (y1*(s-1)+y2*(-1-s)+y3*(1+s)+y4*(1-s))/4;
b = (y1*(t-1)+y2*(1-t)+y3*(1+t)+y4*(-1-t))/4;
c = (x1*(t-1)+x2*(1-t)+x3*(1+t)+x4*(-1-t))/4;
d = (x1*(s-1)+x2*(-1-s)+x3*(1+s)+x4*(1-s))/4;
B1 = [a*(t-1)/4-b*(s-1)/4 0; 0 c*(s-1)/4-d*(t-1)/4;
c*(s-1)/4-d*(t-1)/4 a*(t-1)/4-b*(s-1)/4];
B2 = [a*(1-t)/4-b*(-1-s)/4 0; 0 c*(-1-s)/4-d*(1-t)/4;
c*(-1-s)/4-d*(1-t)/4 a*(1-t)/4-b*(-1-s)/4];
B3 = [a*(t+1)/4-b*(s+1)/4 0; 0 c*(s+1)/4-d*(t+1)/4;
c*(s+1)/4-d*(t+1)/4 a*(t+1)/4-b*(s+1)/4];
B4 = [a*(-1-t)/4-b*(1-s)/4 0; 0 c*(1-s)/4-d*(-1-t)/4;
c*(1-s)/4-d*(-1-t)/4 a*(-1-t)/4-b*(1-s)/4];
Bfirst = [B1 B2 B3 B4];
Jfirst = [0 1-t t-s s-1 ; t-1 0 s+1 -s-t ; s-t -s-1 0 t+1 ; 1-s s+t -t-1 0];
J = [x1 x2 x3 x4]*Jfirst*[y1 ; y2 ; y3 ; y4]/8;
B = Bfirst/J;
if p == 1
D = (E/(1-NU*NU))*[1, NU, 0 ; NU, 1, 0 ; 0, 0, (1-NU)/2];
elseif p == 2
D = (E/(1+NU)/(1-2*NU))*[1-NU, NU, 0 ; NU, 1-NU, 0 ; 0, 0, (1-2*NU)/2];
end
BD = J*transpose(B)*D*B;
r = int(int(BD, t, -1, 1), s, -1, 1);
z = h*r;
8
w = double(z);
end