Can anyone explain the following code?

%% 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*-')
%----------------

[Walter Roberson]

The code is unusual in that it has comments describing the variables. It is not obvious what you were hoping to have explained. Are you hoping for explanation of MATLAB idioms such as assigning [] to something to delete it? Or are you hoping for explanation of the mathematics?

[Aparna Adhikari]

I would like to kow the mathematics behind it.It calulates first five natural frequency and plots mode shapes of a sandwich beam with honeyomb core.It uses first order shear deformation theory for that.Could not understand from line 24

[Walter Roberson]

Unfortunately I myself am not familiar with the mathematics; perhaps someone else will have some relevant experience.

Answer 1

[muthu g]

[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