disp('  ');
disp(' annular_elastic.m   ver 1.3   May 16, 2007  ');
disp('  ');
disp(' by Tom Irvine');
disp(' revised by Matthew Fifield 5/16/07');
disp(' Email:  tomirvine@aol.com ');
disp('  ');
disp(' This program calculates the fundamental frequency of a ');
disp(' annular honeycomb bulkhead supported by isolators');
disp(' mounted along the outside rim.');
disp('  ');
disp(' Assume that the face sheets have the same thickness.');
disp('  ');
%
clear a;
clear b;
clear beta;
clear D;
clear t;
clear rho;
clear mass;
clear volume;
%
iflag=0;
while(iflag==0)
    disp(' ');
    disp(' Enter outside diameter (inch)')
    diam=input(' ');
    b=diam/2.;
%
    disp(' ');
    disp(' Enter inside diameter (inch)')
    diam=input(' ');
    a=diam/2.;
%
    if(a>=b)
        disp(' ')
        disp(' Input data error.  O.D. must be > I.D. ');
    else
        iflag=1;
    end
%
   if(a<0)
       disp(' Input data error. ');    
       iflag=0;
   end
end
%
disp(' ');
disp(' Enter core thickness (inch)')
c=input(' ');
%
disp(' ');
disp(' Enter individual face sheet thickness (inch)')
t=input(' ');
%
disp(' ');
disp(' Enter total mass (lbm)');
mass=input(' ');
mass=mass/386.;
%
disp(' ');
disp(' Select face sheet material ');
disp(' 1=aluminum  2=other ');
imat = input(' ');
%
if(imat==1)
    E=1.0e+07;
else
    disp(' ');
    disp(' Enter elastic modulus of face sheets (lbf/in^2)');
    E=input(' ');
end
%
disp(' ');
disp(' Enter the total isolator stiffness (lbf/in) ');
k=input(' ');
%
mu=0.3;
%
H=c+t;
D=E/(1-mu^2);
D=(D*t*H^2)/2;
%
thickness = c+2*t;
volume = pi*( b^2 - a^2 )*thickness;
rho = mass/volume;
%
out1= sprintf('\n Bending Stiffness/Elastic Modulus = %9.4g in^3',D/E);
disp(out1);
%
out1= sprintf('\n Bending Stiffness = %9.4g lbf in^2 ',D);
disp(out1);
%
out1= sprintf('\n Total Thickness = %9.4g in ',thickness);
disp(out1);
%
out1= sprintf('\n Volume = %9.4g in^3 \n',volume);
disp(out1);
%
fsm=sqrt(k/(rho*volume))/(2*pi);
%
out1= sprintf('\n spring mass fn = %9.4g Hz ',fsm);
disp(out1);
%
adb=a/b;
padb=pi*adb;
%
scale= 0.00194*exp(9.79*adb) +4.82 -0.0479*adb;
num=D;
den=(rho*thickness);
fssp=scale*(1/b^2)*sqrt(num/den)/(2*pi);
%
out1= sprintf('\n annular plate simply-supported outside, free inside \n\n fn = %9.4g Hz\n',fssp);
disp(out1);
%
fdu=(1/fssp^2)+(1/fsm^2);
fdu=1./fdu;
fdu=sqrt(fdu);
out1= sprintf(' Dunkerley fn = %9.4g Hz\n',fdu);
disp(out1);
%
minimum=1.0e+90;
min_beta=0.;
%
for( i=0:10000)
    beta=i/500;   
%   
    d1=(pi/4)*(b^2-a^2)*(2+beta^2);
%    
    d2=b-a*sin((padb)/2.)-((2*b)/pi)*cos((padb)/2.);
    d2=d2*4*b*beta;
%    
    d3=-b-a*pi*sin(padb)-b*cos((padb));
    d3=d3*beta^2*b/(2*pi);
%    
    den=d1+d2+d3;
    den=den*rho*thickness;
% 
    n11=b^2*((1/2)-(1/pi^2));
    n12=(((a^2)/2)+(a*b/pi)*sin(padb)+(b^2/pi^2)*cos(padb));
    n1=n11+n12;
    n1=n1*(0.5*(pi/(2*b))^2);
    n2=0.688*adb^3 -1.57*adb^2 +0.0552*adb +0.823;
    n3=mu*(-1-cos(padb));
%   
    num=n1+n2+n3;
    num=num*beta^2*D*pi^3/(4*b^2);
    num=num+0.5*k;
%    
    om2=num/den;
%
    if(om2 < minimum )
        minimum=om2;
        min_beta=beta;
    end
%
end
%
fn=sqrt(minimum)/(2*pi);
%
out1= sprintf('\n Rayleigh fn = %9.4g Hz\n',fn);
disp(out1);
%
disp('');
%
out1= sprintf(' Rayleigh Eigenvector = [ 1 %7.3g] \n',(1+min_beta));
disp(out1);
%
if(fdu<fn)
    out1= sprintf(' Final Frequency Estimate: %7.4g Hz < fn < %7.4g Hz \n',fdu,fn);
else
    out1= sprintf(' Final Frequency Estimate: %7.4g Hz < fn < %7.4g Hz \n',fn,fdu);  
end
disp(out1);