Adding more rows to the existing matrix(more element)Think of how to get the connectivity matrix in a general case to the existing matlab code.clc;clear all;close all;E=2e11; %(N/m^2) Young's...

2 answer below »
Adding more rows to the existing matrix(more element)
Think of how to get the connectivity matrix in a general case to the existing matlab code.




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=[...






3454 922*b -922*a 1226 398*b 548*a




922*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




1226 398*b -548*a 3454 922*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=[...






3454 -922*b 922*a 1226 -398*b -548*a




-922*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




1226 -398*b 548*a 3454 -922*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];




m21_p=[...






394 232*b -232*a 1226 548*b 398*a




-232*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




1226 548*b -398*a 394 232*b 232*a




-548*b -240*b*b 168*a*b -232*b -120*b*b -112*a*b




-398*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^2)+((1+4*v)/5));


k11_31=2*a*((-2*(beta^2))-((1+4*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^2)-((1+4*v)/5));


k21_13=2*a*((2*(beta^2))+((1-v)/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=[...






-1 0 0




0 1 0




0 0 1];




I2=[...




1 0 0




0 -1 0




0 0 1];




I3=[...




1 0 0




0 1 0




0 0 -1];




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.');









py : D Adding were wows £0 his maka wo move, eles YY) . ) A made an @ Fest Cao. tly) gl [NEw - 3 (x ce] Pp poe | \ 2+ \3 In % How tr webs - dlc ee ax nodes ——7N las ) OA hort chee * ley cf We Conn echo msdn vw will he Lv a Nec+1+41 pa] il > lise c { loot - i Se Loops 45 (ook re whl poof Jer hae - Giese a Conmechrasty sabi foc Joist ent Foe of ons —=nd Mebe & p J Loep 0 Sf “te ret becawsy atl sto | VORA SAM lap UA Cae A Coedr a Loop = a just cece The usy ——— Ceealp doe Ee Covad Ne sn odbet $ 3 Gers X Loot Ctl, Mk x $ Copp af Mak wet pec $2 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=[... 3454 922*b -922*a 1226 398*b 548*a 922*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 1226 398*b -548*a 3454 922*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=[... 3454 -922*b 922*a 1226 -398*b -548*a -922*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 1226 -398*b 548*a 3454 -922*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]; m21_p=[... 394 232*b -232*a 1226 548*b 398*a -232*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 1226 548*b -398*a 394 232*b 232*a -548*b -240*b*b 168*a*b -232*b -120*b*b -112*a*b -398*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^2)+((1+4*v)/5)); k11_31=2*a*((-2*(beta^2))-((1+4*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^2)-((1+4*v)/5)); k21_13=2*a*((2*(beta^2))+((1-v)/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=[... -1 0 0 0 1 0 0 0 1]; I2=[... 1 0 0 0 -1 0 0 0 1]; I3=[... 1 0 0 0 1 0 0 0 -1]; 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.');
Answered 12 days AfterJun 06, 2023

Answer To: Adding more rows to the existing matrix(more element)Think of how to get the connectivity matrix in...

Kshitij answered on Jun 13 2023
32 Votes
SOLUTION.PDF

Answer To This Question Is Available To Download

Related Questions & Answers

More Questions »

Submit New Assignment

Copy and Paste Your Assignment Here