View source for University of Florida/Eml4507/s13.team4ever.Wulterkens.R5.4

{{Big3|Problem 5.4 Eigenvalues and Eigenvectors for Mode Shape}}

On my honor, I have neither given nor recieved unauthorized aid in doing this assignment.
==Problem Statement==
Assume the mass density to be 5,000 <math>kg/m^{3}</math>. Contruct the diagonal mass matrix and find the eigenpairs. With this, design a matlab and calfem program to determine the deformation shape. 

==Matlab Solution==
<syntaxhighlight lang="matlab">

%Constants
p = 5000;
E = 100000000000;
A = 0.0001;
L = 0.3;
%Degree of Freedom
dof = zeros(25,5);
for i = 1:6
    j = 2*i-1;
    dof(i,:) = [i,j,j+1,j+2,j+3];
end
for i =7:12
    j = i*2+1;
    dof(i,:) = [i,j,j+1,j+2,j+3];
end
for i = 13:19
    j = (i-12)*2-1;
    dof(i,:) = [i,j,j+1,j+14,j+15];
end
for i = 20:25
    j = (i-19)*2-1;
    dof(i,:) = [i,j,j+1,j+16,j+17];
end
%coordinates
pos = zeros(2,14);
for i=1:7
    pos(:,i) = [(i-1)*L;0];
    pos(:,i+7) = [(i-1)*L;L];
end
%Connections
conn = zeros(2,25);
for i = 1:6
    conn(1,i)= i;
    conn(2,i)= i+1;
end
for i = 7:12
    conn(1,i) = i+1;
    conn(2,i) = i+2;
end
for i = 13:19
    conn(1,i) = i-12;
    conn(2,i) = i-5;
end
for i = 20:25;
    conn(1,i) = i-19;
    conn(2,i) = i-11;
end
%seperates into x and y coordinates
x = zeros(14,2);
y = zeros(14,2);
for i = 1:25
    x(i,:) = [pos(1,conn(1,i)),pos(1,conn(2,i))];
    y(i,:) = [pos(2,conn(1,i)),pos(2,conn(2,i))];
end
%Set up K and M matrix
K = zeros(28);
M = zeros(28);
for i = 1:25
    xm = x(i,2)-x(i,1);
    ym = y(i,2)-y(i,1);
    L = sqrt(xm^2+ym^2);
    l = xm/L;
    m = ym/L;
    k = E*A/L*[l^2 l*m -l^2 -l*m;l*m m^2 -l*m -m^2;-l^2 -l*m l^2 l*m;-l*m -m^2 l*m m^2;];
    m = L*A*p;
    m = [m/2 0 0 0;0 m/2 0 0;0 0 m/2 0;0 0 0 m/2];
    Dof = dof(i,2:5);
    K(Dof,Dof) = K(Dof,Dof)+k;
    M(Dof,Dof) = M(Dof,Dof)+m;
end

%Makes an alternate K and M matrix to calncel rows and columns for boundary
%conditions.
Kr = K;
Mr = M;
Kr([1 15 16],:) = [];
Kr(:,[1 15 16]) = [];
Mr([1 15 16],:) = [];
Mr(:,[1 15 16]) = [];
%Find Eigenvalue and Eigenvector
[Vec Val] = eig(Kr,Mr);
 
Vec1 = zeros(25,1);
Vec2 = zeros(25,1);
Vec3 = zeros(25,1);
for i = 1:25
    Vec1(i) = Vec(i,1);
    Vec2(i) = Vec(i,2);
    Vec3(i) = Vec(i,3);
end
 
%Insert back in the values for initial conditions.
Vec1 = [0;Vec1(1:13);0;0;Vec1(14:25)];
Vec2 = [0;Vec2(1:13);0;0;Vec2(14:25)];
Vec3 = [0;Vec3(1:13);0;0;Vec3(14:25)];

%Add the Eigenvectors to the nodes

X1 = zeros(14,1);
Y1 = zeros(14,1);
X2 = zeros(14,1);
X3 = zeros(14,1);
Y2 = zeros(14,1);
Y3 = zeros(14,1);
X = zeros(14,1);
Y = zeros(14,1);

for i = 1:14
X1(i) = Vec1((i*2)-1,1);
Y1(i) = Vec1((i*2),1);
X(i) = pos(1,i);
Y(i) = pos(2,i);
X2(i) = Vec2((i*2)-1,1);
Y2(i) = Vec2((i*2),1);
X3(i) = Vec3((i*2)-1,1);
Y3(i) = Vec3((i*2),1);
end

%plot origional shape
for i = 1:25
xc = [X(conn(1,i)),X(conn(2,i))];
yc = [Y(conn(1,i)),Y(conn(2,i))];
axis([0 2 -1 1])
plot(xc,yc)
hold on
end

X1 = X+-0.1*X1;
Y1 = Y+-0.1*Y1;

%plot mode 1
for i = 1:25
xc = [X1(conn(1,i)),X1(conn(2,i))];
yc = [Y1(conn(1,i)),Y1(conn(2,i))];
axis([0 2 -1 1])
plot(xc,yc,'r')
hold on
end


X2 = X+X2;
Y2 = Y+Y2;
X3 = X-0.9*X3;
Y3 = Y-0.9*Y3;

%plot mode 2
for i = 1:25
xc = [X2(conn(1,i)),X2(conn(2,i))];
yc = [Y2(conn(1,i)),Y2(conn(2,i))];
axis([0 2 -1 1])
plot(xc,yc,'r')
hold on
end

%plot mode 3
for i = 1:25
xc = [X3(conn(1,i)),X3(conn(2,i))];
yc = [Y3(conn(1,i)),Y3(conn(2,i))];
axis([0 2 -1 1])
plot(xc,yc,'r')
hold on
end
</syntaxhighlight>
References for Matlab code:
[http://en.wikiversity.org/wiki/Team_Negative_Damping_(3):_Report_4 Team 3 Report 4]
[http://en.wikiversity.org/wiki/Report_4_Team_7 Team 7 Report 4]

===Mode Shapes===
[[file:EML4507.s13.team4ever.Wulterkens.R5.4.1.gif|600px|Mode 1]]
[[file:EML4507.s13.team4ever.Wulterkens.R5.4.2.gif|600px|Mode 1]]
[[file:EML4507.s13.team4ever.Wulterkens.R5.4.3.gif|600px|Mode 1]]

== CALFEM Solution ==

 function R4p4
 <syntaxhighlight lang="matlab">
 p = 5000;  %kg/m3
 A = 0.01;  %m2
 L = 0.3;   %m
 
 h = p*A*L;
 v = h;
 d = (sqrt((L^2)+(L^2)))*A*p;
 
 m = [(h/2)+(d/2)+(v/2);(h/2)+(d/2)+(v/2);
     (h/2)+(v/2);(h/2)+(v/2);
     (h/2)+(h/2)+(d/2)+(v/2);(h/2)+(h/2)+(d/2)+(v/2);
     (h/2)+(h/2)+(d/2)+(v/2);(h/2)+(h/2)+(d/2)+(v/2);
     (h/2)+(h/2)+(d/2)+(v/2);(h/2)+(h/2)+(d/2)+(v/2);
     (h/2)+(h/2)+(d/2)+(v/2);(h/2)+(h/2)+(d/2)+(v/2);
     (h/2)+(h/2)+(d/2)+(v/2);(h/2)+(h/2)+(d/2)+(v/2);
     (h/2)+(h/2)+(d/2)+(v/2);(h/2)+(h/2)+(d/2)+(v/2);
     (h/2)+(h/2)+(d/2)+(v/2);(h/2)+(h/2)+(d/2)+(v/2);
     (h/2)+(h/2)+(d/2)+(v/2);(h/2)+(h/2)+(d/2)+(v/2);
     (h/2)+(h/2)+(d/2)+(v/2);(h/2)+(h/2)+(d/2)+(v/2);
     (h/2)+(h/2)+(d/2)+(v/2);(h/2)+(h/2)+(d/2)+(v/2);
     (h/2)+(v/2);(h/2)+(v/2);
     (h/2)+(d/2)+(v/2);(h/2)+(d/2)+(v/2)];
 
 M = diag(m);
 
 Coord = [0 0;0 0.3;0.3 0;0.3 0.3;0.6 0;
     0.6 0.3;0.9 0;0.9 0.3;1.2 0;1.2 0.3;
     1.5 0;1.5 0.3;1.8 0;1.8 0.3];
 Dof = [1 2;3 4;5 6;7 8;9 10;1 12;13 14;15 16;
     17 18;19 20;21 22;23 24;25 26;27 28];
 Edof = [1 1 2 5 6;2 5 6 9 10;3 9 10 13 14;4 13 14 17 18;5 17 18 21 22;
     6 21 22 25 26;7 3 4 7 8;8 7 8 11 12;9 11 12 15 16;10 15 16 19 20;
     11 19 20 23 24;12 23 24 27 28;13 1 2 3 4;14 5 6 7 8;15 9 10 11 12;
     16 13 14 15 16;17 17 18 19 20;18 21 22 23 24;19 25 26 27 28;20 1 2 7 8;
     21 5 6 11 12;22 9 10 15 16;23 13 14 19 20;24 17 18 23 24;25 21 22 27 28];
 
 [Ex,Ey] = coordxtr(Edof,Coord,Dof,2);
 K = zeros(28);
 ep = [100000000000 0.01];
 
 for i=1:25
     Ke = bar2e(Ex(i,:),Ey(i,:),ep);
     K = assem(Edof(i,:),K,Ke);
 end
 
 [L,X] = eigen(K,M);
 
 eigval = L;
 eigvect = X;

</syntaxhighlight>

Three lowest eigenvalues and eigenvectors


 j1eig = eigval(1)
 j1eigvx = eigvect(:,1)
 j1eigvy = eigvect(:,2)
 
 j2eig = eigval(2)
 j2eigvx = eigvect(:,3)
 j2eigvy = eigvect(:,4)
 
 j3eig = eigval(3)
 j3eigvx = eigvect(:,5)
 j3eigvy = eigvect(:,6)
 
 
 j1eig =
  -4.1672e-08
 
 j1eigvx =
   -0.0018
    0.0853
    0.0136
    0.0853
   -0.0018
    0.0699
    0.0136
    0.0699
   -0.0018
    0.0545
    0.0136
    0.0545
   -0.0018
    0.0391
    0.0136
    0.0391
   -0.0018
    0.0238
    0.0136
    0.0238
   -0.0018
    0.0084
    0.0136
    0.0084
   -0.0018
   -0.0070
    0.0136
   -0.0070
 
 j1eigvy =
   -0.0485
    0.0068
   -0.0488
    0.0068
   -0.0485
    0.0071
   -0.0488
    0.0071
   -0.0485
    0.0074
   -0.0488
    0.0074
   -0.0485
    0.0077
   -0.0488
    0.0077
   -0.0485
    0.0080
   -0.0488
    0.0080
   -0.0485
    0.0083
   -0.0488
    0.0083
   -0.0485
    0.0086
   -0.0488
    0.0086
 
 j2eig =
  -1.0740e-08
 
 j2eigvx =
    0.0155
   -0.0334
   -0.0053
   -0.0334
    0.0155
   -0.0127
   -0.0053
   -0.0127
    0.0155
    0.0081
   -0.0053
    0.0081
    0.0155
    0.0289
   -0.0053
    0.0289
    0.0155
    0.0496
   -0.0053
    0.0496
    0.0155
    0.0704
   -0.0053
    0.0704
    0.0155
    0.0912
   -0.0053
    0.0912
 
 j2eigvy =
   -0.0242
    0.0747
    0.0200
    0.0761
   -0.0197
    0.0156
    0.0197
    0.0202
   -0.0092
   -0.0374
    0.0148
   -0.0335
    0.0040
   -0.0572
    0.0040
   -0.0572
    0.0148
   -0.0335
   -0.0092
   -0.0374
    0.0197
    0.0202
   -0.0197
    0.0156
    0.0200
    0.0761
   -0.0242
    0.0747
 
 j3eig =
   1.2133e-07
 
 j3eigvx =
    0.0206
   -0.0518
   -0.0370
   -0.0552
    0.0097
    0.0333
   -0.0347
    0.0284
   -0.0031
    0.0552
   -0.0190
    0.0633
    0.0002
   -0.0083
   -0.0002
    0.0083
    0.0190
   -0.0633
    0.0031
   -0.0552
    0.0347
   -0.0284
   -0.0097
   -0.0333
    0.0370
    0.0552
   -0.0206
    0.0518
 
 j3eigvy =
    0.0616
   -0.0146
    0.0727
   -0.0164
    0.0453
    0.0010
    0.0647
   -0.0047
    0.0125
    0.0180
    0.0453
    0.0165
   -0.0204
   -0.0017
    0.0204
    0.0017
   -0.0453
   -0.0165
   -0.0125
   -0.0180
   -0.0647
    0.0047
   -0.0453
   -0.0010
   -0.0727
    0.0164
   -0.0616
    0.0146


Plotting the three lowest mode shapes

<syntaxhighlight lang="matlab">
 for j=1:3
     figure
     plot(eigvect(:,j));
     title(['Mode Shape ',num2str(j)])
 end

</syntaxhighlight>

{{CourseCat}}