disp('  ');
disp(' circular_elastic.m   ver 1.0   May 23, 2005  ');
disp('  ');
disp(' by Tom Irvine');
disp(' Email:  tomirvine@aol.com ');
disp('  ');
disp(' This program calculates the fundamental frequency of a ');
disp(' circular 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 beta;
clear D;
clear t;
%
disp(' ');
disp(' Enter diameter (inch)')
diam=input(' ');
a=diam/2.;
%
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*(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);
%
num=((D*pi^3)/(a^2))*(0.2978+0.25*mu);
den=(rho*thickness*a^2)*((pi/4)-(1/pi));
%
fssp=sqrt(num/den)/(2*pi);
%
out1= sprintf('\n plate simply-supported fn = %9.4g Hz\n',fssp);
disp(out1);
%
minimum=1.0e+90;
min_beta=0.;
%
for( i=0:10000)
    beta=i/500;
    c1=pi/2;
    c2=4*beta*(1-(2/pi));
    c3=beta^2*((pi/4)-(1/pi));
    den=c1+c2+c3;
    den=den*rho*thickness*(a^2);
    num=(beta^2*D*pi^3/(a^2))*(0.2978 +0.25*mu) + 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(' Eigenvector = [ 1 %7.3g] \n',(1+min_beta));
disp(out1);