clc;
clear all;
close all;
E=2e11; %(N/m^2) Young's Modulus
v=1/6; %Poisson Ratio
rho=7850;%(kg/m^2) Density
a=1; %(m)Length of one side
b=2; %(m)Length of the other side
h=0.02;%(m)
Nnodes=9; %number of nodes
Connectivity= 1;
alpha=a/b;
beta=b/a;
%coordinates=[...
% ];
%Coefficients of the Mass Matrix
m11_p=[...
XXXXXXXXXX*b -922*a XXXXXXXXXX*b 548*a
922*b 320*b*b -252*a*b 398*b 160*b*b 168*a*b
XXXXXXXXXX*a -252*a*b 320*a*a -548*a -168*a*b -240*a*a
XXXXXXXXXX*b -548*a XXXXXXXXXX*b 922*a
398*b 160*b*b -168*a*b 922*b 320*b*b 252*a*b
548*a 168*a*b -240*a*a 922*a 252*a*b 320*a*a];
m22_p=[...
XXXXXXXXXX*b 922*a XXXXXXXXXX*b -548*a
XXXXXXXXXX*b 320*b*b -252*a*b -398*b 160*b*b 168*a*b
922*a -252*a*b 320*a*a 548*a -168*a*b -240*a*a
XXXXXXXXXX*b 548*a XXXXXXXXXX*b -922*a
XXXXXXXXXX*b 160*b*b -168*a*b -922*b 320*b*b 252*a*b
XXXXXXXXXX*a 168*a*b -240*a*a -922*a 252*a*b 320*a*a];
m21_p=[...
XXXXXXXXXX*b -232*a XXXXXXXXXX*b 398*a
XXXXXXXXXX*b -120*b*b 112*a*b -548*b -240*b*b -168*a*b
232*a 112*a*b -120*a*a 398*a 168*a*b 160*a*a
XXXXXXXXXX*b -398*a XXXXXXXXXX*b 232*a
XXXXXXXXXX*b -240*b*b 168*a*b -232*b -120*b*b -112*a*b
XXXXXXXXXX*a -168*a*b 160*a*a -232*b -112*a*b -120*a*a];
% me=(rho*h*a*b/6300)*[...
% m11 m21'
% m21 m22];
m11=m11_p(1:3,1:3);
m21=m11_p(4:6,1:3);
m22=m11_p(4:6,4:6);
m12=m21';
m31=m21_p(1:3,1:3);
m41=m21_p(4:6,1:3);
m13=m31';
m14=m41';
m32=m21_p(1:3,4:6);
m42=m21_p(4:6,4:6);
m23=m32';
m24=m42';
m33=m22_p(1:3,1:3);
m43=m22_p(4:6,1:3);
m34=m43';
m44=m22_p(4:6,4:6);
%Elemental Mass Matrix;
me=[...
m11 m12 m13 m14
m21 m22 m23 m24
m31 m32 m33 m34
m41 m42 m43 m44];
% T3=[...
% lx mx nx
% ly my ny
% lz mz nz];
%
% T=diag(T3);
%
% M=T'*me*T;
% Coefficients of the stiffness matrix
k11_11=4*(beta^2+alpha^2)+2*(7-2*v)/5;
k11_21=2*b*(2*(alpha^ XXXXXXXXXX*v)/5));
k11_31=2*a*((-2*(beta^ XXXXXXXXXX*v)/5));
k11_22=4*(b^2)*((4*(alpha^2)/3)+(4*(1-v)/15));
k11_32=-4*v*a*b;
k11_33=4*(a^2)*((4*(beta^2)/3)+(4*(1-v)/15));
k11=[...
k11_11 k11_21 k11_31;
k11_21 k11_22 k11_32;
k11_31 k11_32 k11_33];
k21_11=-2*(2*beta^2-alpha^2)+2*(7-2*v)/5;
k21_12=2*b*((alpha^ XXXXXXXXXX*v)/5));
k21_13=2*a*((2*(beta^ XXXXXXXXXXv)/5));
k21_21=k21_12;
k21_31=-k21_13;
k21_22=4*(b^2)*((2*(alpha^2)/3)-(4*(1-v)/15));
k21_33=4*(a^2)*((2*(beta^2)/3)-(1*(1-v)/15));
k21=[...
k21_11 k21_12 k21_13
k21_21 k21_22 0
k21_31 0 k21_33];
k31_11=-2*(beta^2+alpha^2)-2*(7-2*v)/5;
k31_12=2*b*(-(alpha^2)+((1-v)/5));
k31_13=2*a*((beta^2)-(1-v)/5);
k31_21=2*b*((alpha^2)-((1-v)/5));
k31_22=4*(b^2)*(((alpha^2)/3)+(1*(1-v)/15));
k31_31=2*a*(-(beta^2)+(1-v)/5);
k31_33=4*(a^2)*(((beta^2)/3)+(1*(1-v)/15));
k31=[...
k31_11 k31_12 k31_13
k31_21 k31_22 0
k31_31 0 k31_33];
k41_11=2*((beta^2)-2*(alpha^2))-2*(7-2*v)/5;
k41_12=2*b*(-2*(alpha^2)-(1-v)/5);
k41_13=2*a*(-(beta^2)+(1+4*v)/5);
k41_21=2*b*(2*(alpha^2)+(1-v)/5);
k41_22=4*(b^2)*((2*(alpha^2)/3)-(1*(1-v)/15));
k41_31=2*a*(-(beta^2)+(1+4*v)/5);
k41_33=4*(a^2)*((2*(beta^2)/3)-(4*(1-v)/15));
k41=[...
k41_11 k41_12 k41_13
k41_21 k41_22 0
k41_31 0 k41_33];
I1=[...
XXXXXXXXXX
XXXXXXXXXX
XXXXXXXXXX];
I2=[...
XXXXXXXXXX
XXXXXXXXXX
XXXXXXXXXX];
I3=[...
XXXXXXXXXX
XXXXXXXXXX
XXXXXXXXXX];
k22=I3'*k11*I3;
k32=I3'*k41*I3;
k42=I3'*k31*I3;
k33=I1'*k11*I1;
k43=I1'*k21*I1;
k44=I2'*k11*I2;
% exp_k=1;
kn=E*(h^3)/(48*(1-(v^2))*a*b);
k=kn*[...
k11 k21' k31' k41'
k21 k22 k32' k42'
k31 k32 k33 k43'
k41 k42 k43 k44];
% eqn1=k-(omega^2)*me;
% solve(det(k-(omega^2).*me==0),omega)
% solu=solve(eqn1=0,omega);
z=zeros(3,3);
% Node coordinates and connectivity matrix
% GNodeCoords = [0 0; Lx 0; Lx Ly; 0 Ly];
B = [1 2 5 4;2 3 6 5;6 9 8 5;8 7 4 5];
%K_global=assemble(K_local,B)
%K_global is global stiffness matrix
%k is local stiffness matrices of all elements
% B = connectivity matrix
% when B=ElConnMatrix 1
[NE,j] = size(B); %number of elements and number of nodes per element
Nnodes = max(max(B));
K_global = zeros(3*Nnodes,3*Nnodes);
M_global = zeros(3*Nnodes,3*Nnodes);
for i=1:NE
aux=[3*B(i,:)-2;3*B(i,:)-1;3*B(i,:)];
aux=aux(:);
K_global(aux,aux)=K_global(aux,aux)+k;
M_global(aux,aux)=M_global(aux,aux)+me;
end
%Global Stiffness Matrix
K=kn*[...
k11 k21' z k41' k31' z z z z
k21 k22+k11 k21' k42' k32'+k41' k31' z z z
z k21 k22 z k42' k32' z z z
k41 k42 z k44+k33 k43+k43' z k32 k31 z
k31 k32+k41 k42 k43'+k43 k33+3*k44 k43+k41 k42 k43+k41 k42
z k31 k32 z k43'+k41' k33+k11 z k31' k21'
z z z k32' k42' z k22 k21 z
z z z k31' k43'+k41' k31 k21' k33+k11 k32
z z z z k42' k21 z k32' k22];
%Global Mass Matrix
% M=(rho*h*a*b/6300)*[...
M=[...
m11 m21' z m41' m31' z z z z
m21 m22+m11 m21' m42' m32'+m41' m31' z z z
z m21 m22 z m42' m32' z z z
m41 m42 z m44+m33 m43+m43' z m32 m31 z
m31 m32+m41 m42 m43'+m43 m33+3*m44 m43+m41 m42 m43+m41 m42
z m31 m32 z m43'+m41' m33+m11 z m31' m21'
z z z m32' m42' z m22 m21 z
z z z m31' m43'+m41' m31 m21' m33+m11 m32
z z z z m42' m21 z m32' m22];
[vecfreq,freq]=eig(K,M);
% freq = diag(freq) ;
% freq=sqrt(freq); % UNITS :rad per sec
% freqHz = abs(freq/(2*pi)) ; % UNITS : Hertz
% freqs=sort(freqHz);
[freq_sorted, ind] = sort(diag(freq),'ascend');
V_sorted = vecfreq(:,ind);
req_freq=sqrt(freq_sorted(1:10)); %Top 10 natural frequencies
freqHz = abs(req_freq/(2*pi)); % UNITS : Hertz
for jj=1:10
V = zeros(size(V_sorted,2),size(V_sorted,2)) ;
x_a=linspace(0,a*2,27);
y_b=linspace(0,b*2,27);
[s,t] = meshgrid(x_a,y_b);
for ii = 1:size(V_sorted,2)
V(:,ii) = V_sorted(:,jj);
% figure(ii)
% plot(V(:,ii));
end
% figure(jj)
% surf (s,t,abs(V));
end
% N1 = 1/4*(1-t).*(1-s);
% subplot(2,2,1);
% surf(s,t,abs(vecfreq));
% xlabel('a');
% ylabel('b');
%
% N2 = @(s,t)1/4*(1-t).*(1+s);
% subplot(2,2,2);
% surf(s,t,N2(s,t));
% xlabel('a');
% ylabel('b');
%
% N3 = @(s,t)1/4*(1+t).*(1+s);
% subplot(2,2,3);
% surf(s,t,N3(s,t));
% xlabel('a');
% ylabel('b');
%
% N4 = @(s,t)1/4*(1+t).*(1-s);
% subplot(2,2,4);
% surf(s,t,N4(s,t));
% xlabel('a');
% ylabel('b');
SymmetryM=(M-M.');
SymmetryK=(K-K.');
SymmetryMe=(me-me.');