User:Egm4507.s13.team4.colocar/FEA.s13.R4.4.cc

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



For animation using smartgifcreator, for j=1:28 clear plot ed = extract(Edof,X(:,j)); P = eldisp2(Ex,Ey,ed,plotpar,scale); W(j) = getframe; drawnow pause end