Can anyone explain the following code?
5 views (last 30 days)
Show older comments
%% Vibration of nomex honeycomb sandwich structure
%=============================================================================
clear all % To clear all variables
close all % To close all windows
clc; % To clear the screan
%=============================================================================
%% Pre-Processor
%=============================================================================
% Geometrical Properties
L=0.25; % Length of the beam in meters
B=0.03; % Breadth of the beam in meters
t=0.25e-3; % Thickness of each lamina in meters
tc=2e-3; % Thickness of nomex honeycomb
%==============================================================================
% Material Properties (All units are in SI)
% =======================
% Face Sheet properties
%========================
E1=70e9; % Youngs modulus in direction 1
E2=6e9; % Youngs modulus in direction 2
G12=3.2e9; % shear modulus in the plane 12
v12=0.22; % poisson ratio in the plane 12
rho=1627;
theta=[0 90 90 0]*(pi/180); % Ply angle of face sheet
N=size(theta,2);
d=tc+N*t;
e=d/tc;
% ===================================
% Nomex honeycomb material properties
%===================================
Ecm=0.9e9;
Gc=0.32e9; % shear modulus in the plane 12
rhocm=724;
wt=0.063;
wl=3.2;
% ===================================
% Nomex honeycomb properties
%====================================
Gc13=(5/(3*sqrt(3)))*(wt/wl)*Gc;
Ec1=(4/sqrt(3))*(wt/wl)*Ecm;
rhoc=(2/sqrt(3))*(wt/wl)*rhocm;
%================================
%ABD Matrix
%--------------------------------
% Compliance matrix
%--------------------------------
s = [1/E1 -v12/E1 0 ;
-v12/E1 1/E2 0 ;
0 0 1/G12;];
%-----------------------------------
% ABD and I matrix
%-----------------------------------
Q1=s\eye(3);
A1=0;
B1=0;
D1=0;
I0=0;
I1=0;
I2=0;
Zk =t*(0-N/2);
for i1=1:N
teta=theta(i1);
T = [ cos(teta)^2 sin(teta)^2 -2*sin(teta)*cos(teta) ;
sin(teta)^2 cos(teta)^2 2*sin(teta)*cos(teta) ;
cos(teta)*sin(teta) -cos(teta)*sin(teta) cos(teta)^2-sin(teta)^2;];
Q = T*Q1*transpose(T) ;
Zk1 = t*(i1-N/2);
A1 = A1+Q(1,1)*(Zk1-Zk);
B1 = B1+(1/2)*Q(1,1)*(Zk1^2-Zk^2);
D1 = D1+(1/3)*Q(1,1)*(Zk1^3-Zk^3);
I0 = I0 + (rho)*(Zk1-Zk);
I1 = I1 + (1/2)*(rho)*(Zk1^2-Zk^2);
I2 = I2 + (1/3)*(rho)*(Zk1^3-Zk^3);
Zk = Zk1;
end
[ABD]=[A1 0 B1;
0 A1 B1;
B1 B1 2*D1;];
[I] =[I0 0 I1;
0 I0 I1;
I1 I1 2*I2;];
I0c = rhoc*tc;
%==============================================================================
% Meshing
%==============================================================================
Nel=100; % No. of Elements
Le=L/Nel; % Length of the element
nodes=Nel+1; % Total no. of nodes
dof=4; % degrees of freedom per node
sdof=nodes*dof; % System degrees of freedom
nnodes=1:nodes; % Numbering of nodes
xcoord=Le*(nnodes-1);
% ===========================================================================
% boundary conditions
% ===========================================================================
bc=[1 2 3 4]; % Cantilever beam
% bc=[1 2 3 4 4*nnodes(end)-3 4*nnodes(end)-2 4*nnodes(end)-1 4*nnodes(end)]; % Fixed beam
% bc=[1 2 3 4*nnodes(end)-3 4*nnodes(end)-2 4*nnodes(end)-1]; % Simply supported beam
%=============================================================================
%% Processor
%=============================================================================
% initilization of Global Stiffness and Mass matrix
K=zeros(sdof); % Global stiffness Matrix;
M=zeros(sdof); % Global Mass Matrix;
%=================================================================
% Gauss points for numerical Integration
zeta=[-1/sqrt(3) 1/sqrt(3)];
gw=[1 1];
[N1dbar, N1dbeam]=shapeFn1db(zeta,Le); % Shape function values at gauss point
[dN1dbar, dN1dbeam]=dshapeFn1db(zeta,Le); % differentiation of shape function values at gauss point
[ddN1dbeam]=ddshapeFn1db(zeta,Le); % double differentiation of shape function values at gauss point
%=================================================================
Kel=zeros(2*dof); % Element stiffness matrix
Mel=zeros(2*dof); % Element mass matrix
for i2=1:Nel
conn=[nnodes(i2) nnodes(i2+1)]; %connectivity Matrix
index=[4*(conn(1)-1)+1 4*(conn(1)-1)+2 4*(conn(1)-1)+3 4*(conn(1)-1)+4 4*(conn(2)-1)+1 4*(conn(2)-1)+2 4*(conn(2)-1)+3 4*(conn(2)-1)+4];
xgp(i2,:)=xcoord(conn)*N1dbar;
for i3=1:size(zeta,2)
J=dN1dbar(:,i3).'*xcoord(conn).';
dN1dbar=(1/J)*dN1dbar;
dN1dbeam=(1/J)*dN1dbeam;
ddN1dbeam=(1/J^2)*ddN1dbeam;
Bb1=[dN1dbar(1,i3) 0 0 0 dN1dbar(2,i3) 0 0 0;
0 dN1dbar(1,i3) 0 0 0 dN1dbar(2,i3) 0 0;
0 0 ddN1dbeam(1,i3) ddN1dbeam(2,i3) 0 0 ddN1dbeam(3,i3) ddN1dbeam(4,i3);];
S1=[N1dbar(1,i3) 0 0 0 N1dbar(2,i3) 0 0 0;
0 N1dbar(1,i3) 0 0 0 N1dbar(2,i3) 0 0;
0 0 dN1dbeam(1,i3) dN1dbeam(2,i3) 0 0 dN1dbeam(3,i3) dN1dbeam(4,i3);];
S2=[0 0 N1dbeam(1,i3) N1dbeam(2,i3) 0 0 N1dbeam(3,i3) N1dbeam(4,i3)];
S3=e*[N1dbar(1,i3)/d -N1dbar(1,i3)/d dN1dbeam(1,i3) dN1dbeam(2,i3) N1dbar(2,i3)/d -N1dbar(2,i3)/d dN1dbeam(3,i3) dN1dbeam(4,i3);];
B2=e*[N1dbar(1,i3)/d -N1dbar(1,i3)/d dN1dbeam(1,i3) dN1dbeam(2,i3) N1dbar(2,i3)/d -N1dbar(2,i3)/d dN1dbeam(3,i3) dN1dbeam(4,i3);];
B3=(1/2)*[dN1dbar(1,i3) dN1dbar(1,i3) 0 0 ddN1dbeam(2,i3) dN1dbar(2,i3) 0 0;];
% Kel=Kel+B*(Bb1.'*ABD*Bb1+B2.'*Gc13*tc*B2+B3.'*Ec1*tc*B3)*J*gw(i3);
Kel=Kel+B*(Bb1.'*ABD*Bb1+B2.'*Gc13*tc*B2)*J*gw(i3);
Mel=Mel+B*(S1.'*I*S1+S2.'*(2*I0+I0c)*S2+S3.'*(rhoc*tc^3/3)*S3)*J*gw(i3);
end
%------------------------------------------------------------------
%Assemble stiffness and Mass vector
%------------------------------------------------------------------
K(index,index)=K(index,index)+Kel;
M(index,index)=M(index,index)+Mel;
end
%% apply boundary conditions
K(bc,:)=[];
K(:,bc)=[];
M(bc,:)=[];
M(:,bc)=[];
%% Natural frequiencies
[eigvec,eigval]=eig(K,M);
nfrq=sort(sqrt(diag(eigval))); % natural frequencies in rad/s
natfre=(nfrq)/(2*pi); % natural frequencies in Hz/s
[lamda1,kn1]=sort(diag(eigval)); % to sort eigvector
eigvector=eigvec(:,kn1); % mode shpe
disp('First Five Natural Frequencies')
disp(natfre(1:5))
% mode shapes of cantelever beam
%--------------------------------------------------------------------------
K1=4;
ms=eigvector(3:4:end,K1);
ms=[0; ms;];
ms=(1/(max(abs(ms))))*ms;
x1=linspace(0,1,Nel+1);
plot(x1,ms,'r*-')
%--------------------------------------------------------------------------
% mode shapes of Fixed beam
%--------------------------------------------------------------------------
% K1=1;
% ms=eigvector(3:4:end,K1);
% ms=[0; ms; 0;];
% ms=(1/(max(abs(ms))))*ms;
% x1=linspace(0,1,Nel+1);
% plot(x1,ms,'r*-')
%--------------------------------------------------------------------------
% mode shapes of cantelever beam
%--------------------------------------------------------------------------
% K1=4;
% ms=eigvector(4:4:end,K1);
% ms=[0; ms; 0;];
% ms=(1/(max(abs(ms))))*ms;
% x1=linspace(0,1,Nel+1);
% plot(x1,ms,'r*-')
%----------------
3 Comments
Walter Roberson
on 8 May 2021
Unfortunately I myself am not familiar with the mathematics; perhaps someone else will have some relevant experience.
Answers (1)
muthu g
on 2 Dec 2021
[N1dbar, N1dbeam]=shapeFn1db(zeta,Le); % Shape function values at gauss point
[dN1dbar, dN1dbeam]=dshapeFn1db(zeta,Le); % differentiation of shape function values at gauss point
[ddN1dbeam]=ddshapeFn1db(zeta,Le); % double differentiation of shape function values at gauss point
function file is missing
0 Comments
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!