University of Florida/Eml4507/s13.team4ever.Wulterkens.R5.4

From Wikiversity
Jump to navigation Jump to search

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[edit]

Assume the mass density to be 5,000 . Contruct the diagonal mass matrix and find the eigenpairs. With this, design a matlab and calfem program to determine the deformation shape.

Matlab Solution[edit]

%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

References for Matlab code: Team 3 Report 4 Team 7 Report 4

Mode Shapes[edit]

Mode 1 Mode 1 Mode 1

CALFEM Solution[edit]

function R4p4
 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;

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

 for j=1:3
     figure
     plot(eigvect(:,j));
     title(['Mode Shape ',num2str(j)])
 end