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

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 ${\displaystyle kg/m^{3}}$. Contruct the diagonal mass matrix and find the eigenpairs. With this, design a matlab and calfem program to determine the deformation shape.

Matlab Solution

```%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

CALFEM Solution

```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
```