#include #include #include #include #include #define NUM 1500 void files(void); void attach(void); void auto_spectrum(void); void viscosity(void); void speed(void); void frequency(void); void density(void); void franken(void); void mach(void); void choose_regime(void); void cc_plateau_tf(void); void cc_reattachment_tf(void); void cc_plateau_supersonic(void); void cc_separation_or_shockwave(void); void cc_auto(void); void exp_auto(void); void P2P1_ratio(void); void reynolds(void); void printdata(void); void exp_plateau(void); void exp_reattachment(void); void material(void); void maxcase(void); void maxpsdrms(void); void reverse_speed(void); void find_altitude(void); const double PI= atan2(0.,-1.); const double TPI=2.*atan2(0.,-1.); const double splref = (2.9e-09)*144.; const double ft_per_m = 3.28; const double ft_per_km = 3.28*1000.; const double ft_per_nmi = 6076.1; const double km_per_nmi = 1.852; const double conv = 0.00194; double al,alpha; double c,cexp,comp,constant,cmat,cref; double deltastar,df; double E; double F; double h; double k1,k2; double len; double M,M1,M2,mu; double omega; double poq,poqcomp,poqexp,prms,poqtbl,pshock,ptbl,P2P1; double q; double Rex,rho; double speed_of_sound=0.; double s2,scale,sfr,stationdiam,sum; double theta,thickness; double U; double va,vdb; double W; double x; double WT_ratio; double ar[30],fr[30],f[100],freq[100],G[100],Goct[100],psd[100],s[30],spl[NUM][100]; double altitude[NUM],macht[NUM],dp[NUM],velox[NUM]; double maxpsd[100],maxspl[100]; double maxsplcase; double maxq; double maxMach; double maxaltitude; long maxikk; double maxoaspl; int icc,iex,iflow,ifm,ikk,ilast,imat,j,jk,last; int i=0; int imax=0; int iunit,i_vel_unit; int i_first; int i_second; char filename[8][FILENAME_MAX]; FILE *pFile[8]; char string[30]; void main() { strcpy(string,"Flowtable, ver 2.5, August 27, 2007"); printf("\n\n\n\n %s",string); printf("\n\n by Tom Irvine \n\n"); files(); fprintf(pFile[3],"\n\n\n\n %s",string); fprintf(pFile[3],"\n\n by Tom Irvine \n"); maxsplcase =0.; maxq =0.; maxMach =0.; maxaltitude =0.; maxikk =0; maxoaspl =0.; for(i=0; i<100; i++) { maxpsd[i]=0.; maxspl[i]=0.; } frequency(); choose_regime(); printf("\n Enter the station x (inch) relative to front end.\n"); scanf("%lf",&x); fprintf(pFile[3],"\n Station = %8.3lf in\n",x); printf("\n\n"); x/=12.; for(ikk=0; ikk<=last; ikk++) { imax=0; h=altitude[ikk]; M=macht[ikk]; if(M<1.0) { printf("\n mach error \n"); printf("\n ikk=%d h=%12.4g mach=%12.4g \n",ikk,h,M); exit(1); } mach(); speed(); density(); if(i_first !=3 && i_second !=3) { q = 0.5*rho*pow(U,2); } else { q=dp[ikk]; } reynolds(); if(iflow==1) // Attached Flow (Type I) { if(ikk==0) { printf("\n Attached Flow \n via the Laganelli and Wolfe method, AIAA paper 89-1064 \n\n"); fprintf(pFile[3],"\n Attached Flow (Type I) \n via the Laganelli and Wolfe method, AIAA paper 89-1064 \n\n"); } attach(); } if(iflow==2) // Separated Flow & Shockwaves: Compression Corner (IV or V) { if(icc==1) { if(ikk==0) { fprintf(pFile[3],"\n Plateau Region, Transonic Flow (Type IV) "); } cc_plateau_tf(); } if(icc==2) { if(ikk==0) { fprintf(pFile[3],"\n Reattachment Region, Transonic Flow (Type V)"); } cc_reattachment_tf(); } if(icc==3) { if(ikk==0) { fprintf(pFile[3],"\n Plateau Region, Supersonic Flow (Type IV)"); } cc_plateau_supersonic(); } if(icc==4) { if(ikk==0) { fprintf(pFile[3],"\n Separation or Reattachment Shockwave (Type V)\n"); } cc_separation_or_shockwave(); } cc_auto(); } if( iflow==3) // Separated Flow & Shockwaves: Expansion Corner (II or III) { if(iex==1) { if(ikk==0) { fprintf(pFile[3],"\n Plateau Region, Transonic and Supersonic (Type II)"); } exp_plateau(); } if(iex==2) { if(ikk==0) { fprintf(pFile[3],"\n Reattachment Region, Transonic (Type III)"); } exp_reattachment(); } exp_auto(); } auto_spectrum(); if(sum > maxsplcase ) { maxsplcase=sum; maxq = q; maxMach=M; maxaltitude=altitude[ikk]; maxikk = ikk; maxoaspl=sum; } if(ikk==0) { fprintf(pFile[3],"\n\n n alt(nmi) Mach q(psf) OASPL(dB) \n"); } fprintf(pFile[3]," %ld %8.4g %8.4g %8.4g %8.4g \n",ikk,altitude[ikk]/ft_per_nmi,M,q,sum); printf(" %ld alt=%8.4g nmi Mach=%8.4g q=%8.4g psf OASPL=%8.4g dB\n",ikk,altitude[ikk]/ft_per_nmi,M,q,sum); } printf("\n worst case OASPL\n\n %ld alt=%8.4g nmi maxMach=%8.4g q=%8.4g psf OASPL=%8.4g dB\n",maxikk,maxaltitude/ft_per_nmi,maxMach,maxq,maxoaspl); fprintf(pFile[3],"\n worst case OASPL\n\n %ld alt=%8.4g nmi maxMach=%8.4g q=%8.4g psf OASPL=%8.4g dB\n",maxikk,maxaltitude/ft_per_nmi,maxMach,maxq,maxoaspl); printf("\n\n Press any key to continue.\n"); getch(); q=maxq; M=maxMach; h=maxaltitude; // printf("\n Reference 1: U=%8.4g M=%8.4g maxMach=%8.4g \n",U,M,maxMach); maxcase(); // calls mach, speed printf("\n\n Calculate resulting vibration via the Franken method? "); printf("\n 1=yes 2=no \n"); scanf("%d",&ifm); if(ifm==1) { printf("\n Franken upper limit curve is used per NASA CR-1302 \n"); material(); for(ikk=0; ikk<=last; ikk++) { franken(); } fprintf(pFile[3],"\n Maximum envelope PSD of all flow cases \n"); fprintf(pFile[3],"\n\n f(Hz) \t PSD(G^2/Hz) \n"); for(i=0; i<=ilast; i++) { fprintf(pFile[3]," %8.3lf \t %9.3e\n",f[i],maxpsd[i]); fprintf(pFile[4]," %12.4e \t %12.4e\n",f[i],maxpsd[i]); fprintf(pFile[6]," %12.4e \t %12.4e\n",f[i],2.*maxpsd[i]); } maxpsdrms(); } printf("\n close files \n"); for(i=0;i<=6;i++) { fclose(pFile[i]); } printf("\n\n ************ Output files ************** "); printf("\n\n \t aeroflow.txt - text file "); printf("\n\n Worst case individual SPL \n"); printf("\n \t wc_spl.dat - freq(Hz) & SPL(dB)"); printf("\n\n Maximum SPL envelope of all flow cases\n"); printf("\n \t maxe_spl.dat - freq(Hz) & SPL(dB)"); if(ifm==1) { printf("\n\n Maximum PSD envelope of all flow cases: \n"); printf("\n \t maxe_psd.dat - freq(Hz) & PSD(G^2/Hz) \n"); printf("\n\n Maximum PSD + 3 dB: \n"); printf("\n \t maxe_psd3dB.dat - freq(Hz) & PSD(G^2/Hz) \n"); } printf("\n\n Press any key to exit.\n"); getch(); exit(1); } void frequency(void) { freq[0]=20.; freq[1]=25.; freq[2]=31.5; freq[3]=40.; freq[4]=50.; freq[5]=63.; freq[6]=80.; freq[7]=100.; freq[8]=125.; freq[9]=160.; freq[10]=200.; freq[11]=250.; freq[12]=315.; freq[13]=400.; freq[14]=500.; freq[15]=630.; freq[16]=800.; freq[17]=1000.; freq[18]=1250.; freq[19]=1600.; freq[20]=2000.; freq[21]=2500.; freq[22]=3150.; freq[23]=4000.; freq[24]=5000.; freq[25]=6300.; freq[26]=8000.; freq[27]=10000.; ilast = 27; for(i=0; i<=27; i++) { f[i]=freq[i]; } } void viscosity(void) { // kinematic viscosity ft^2/sec double ht=h/1000.; if(h<360000.) { mu = 1.560401 + 0.0395762*ht + 3.009485e-04*pow(ht,2.)+ 1.675145e-05*pow(ht,3.); } else { mu = -4.48529 + 0.4546508*ht - 1.090837e-02*pow(ht,2.)+ 1.373416e-04*pow(ht,3.); } mu*=0.0001; if(imax==1) { printf("\n Kinematic viscosity = %12.4g ft^2/sec\n",mu); fprintf(pFile[3],"\n Kinematic viscosity = %12.4g ft^2/sec\n",mu); } } void attach(void) { if(imax==1) { // printf("\n deltastar = %8.3lf ft \n",deltastar); printf("\n F = %8.3lf \n",F); } // G is the pressure autospectrum prms = (0.01/F)*q; poq = prms/q; if(imax==1) { // printf("\n factor = %8.3lf \n",factor); printf("\n prms=%12.4g psf q=%12.4g psf poq=%12.4g\n",prms,q,poq); } constant = 4.*pow(poq,2.)*pow(F,1.433); constant*= pow(q,2.)*(deltastar/U); constant/= pow(splref,2.); for(i=0; i<=ilast; i++) { omega= TPI*freq[i]*deltastar/U; G[i] = constant/(1. + pow(F,2.867)*pow(omega,2.)); // printf(" freq=%12.4g G=%12.4g U=%12.4g\n",fff,G[i],U); } } void auto_spectrum(void) { if(imax==1) { fprintf(pFile[3],"\n\n Worst Case SPL \n"); fprintf(pFile[3],"\n freq (Hz) \t SPL(dB) \n"); } sum=0.; for(jk=0; jk<=27; jk++) { Goct[jk]= G[jk]*(freq[jk]*0.2301) ; // printf("%12.4g %12.4g \n",freq[jk],Goct[jk]); spl[ikk][jk]=10.*log10(Goct[jk]); if( spl[ikk][jk] > maxspl[jk] ){ maxspl[jk] = spl[ikk][jk] ;} if(imax==1) { fprintf(pFile[1],"%12.4e \t %12.4e \n",freq[jk],Goct[jk]); fprintf(pFile[2],"%12.4e \t %12.4e \n",freq[jk],spl[ikk][jk]); fprintf(pFile[3],"%8.1lf \t %8.1lf \n",freq[jk],spl[ikk][jk]); fprintf(pFile[5],"%8.1lf \t %8.1lf \n",freq[jk],maxspl[jk]); } sum+=Goct[jk]; } // sum=sqrt(sum); sum=10.*log10(sum); if(imax==1) { printf("\n Overall SPL = %8.2lf dB \n",sum); fprintf(pFile[3],"\n\n Overall SPL = %8.2lf dB \n",sum); } } void speed(void) { double LH1=11. *ft_per_km; double LH2=20. *ft_per_km; if( h < LH1) // troposphere { double Tz=288.; double L=6.5; double gamma=1.402; double ROM = (8314.3/28.97); if( h < 0.) { printf("\n Error \n"); printf("\n\n Press any key to exit.\n"); getch(); exit(1); } c=sqrt(fabs(gamma*ROM*(Tz-L*(h/(ft_per_km)))))*3.28; // printf("\n gamma=%12.4e ROM=%12.4e Tz=%12.4e L=%12.4e h=%12.4e\n",gamma,ROM,Tz,L,h); } if( h >= LH1 && h < LH2 ) // lower stratosphere { c = 295.*3.28; } if( h >= LH2 ) { printf("\n Warning: altitude is too high for speed of sound calculation. \n\n"); } if(speed_of_sound > 0. ) { c=speed_of_sound; } U= M*c; } void density(void) { double alt[30],dens[30]; for(i=0; i<=20; i++) { alt[i]=double(i)*ft_per_km; } dens[0]=1.226; dens[1]=1.112; dens[2]=1.007; dens[3]=0.9096; dens[4]=0.8195; dens[5]=0.7365; dens[6]=0.6600; dens[7]=0.5898; dens[8]=0.5254; dens[9]=0.4666; dens[10]=0.4129; dens[11]=0.3641; dens[12]=0.3104; dens[13]=0.2652; dens[14]=0.2266; dens[15]=0.1936; dens[16]=0.1654; dens[17]=0.1413; dens[18]=0.1207; dens[19]=0.1032; dens[20]=0.0881; rho=0.; for(i=0; i<=20; i++) { dens[i]*= conv; //slugs/ft^3 } double len,c1,c2; for(i=0; i< 20; i++) { if( h >= alt[i] && h < alt[i+1] ) { len = alt[i+1]-alt[i]; c2 = (h-alt[i])/len; c1 = 1. -c2; rho = c1*dens[i] + c2*dens[i+1]; break; } } } void files(void) { printf("\n\n The input file must have two of the following three parameters:\n"); printf("\n altitude, velocity, dynamic pressure. \n"); while(1) { printf("\n\n Select first column metric: "); printf("\n 1=altitude 2=velocity 3=dynamic pressure \n"); scanf("%d",&i_first); printf("\n Select second column metric: "); printf("\n 1=altitude 2=velocity 3=dynamic pressure \n"); scanf("%d",&i_second); if(i_first==1 || i_first==2 || i_first==3) { if(i_second==1 || i_second==2 || i_second==3) { if(i_first != i_second) { break; } } } } if(i_first==2 || i_second==2) { printf("\n\n Select input file velocity unit "); printf("\n 1=Mach 2=ft/sec \n"); scanf("%d",&i_vel_unit); if(i_vel_unit !=1 && i_vel_unit !=2 ) { printf("\n input error \n"); exit(1); } } if(i_first==1 || i_second==1) { printf("\n Select input file altitude unit "); printf("\n 1=nmi 2=feet 3=km \n"); scanf("%d",&iunit); if(iunit !=1 && iunit !=2 && iunit !=3 ) { printf("\n input error \n"); exit(1); } } if(i_first==1 && i_second==2) { if(iunit==1 && i_vel_unit ==1) { printf("\n\n The input file must have two columns: altitude(nmi) and Mach \n"); } if(iunit==2 && i_vel_unit ==1) { printf("\n\n The input file must have two columns: altitude(feet) and Mach \n"); } if(iunit==3 && i_vel_unit ==1) { printf("\n\n The input file must have two columns: altitude(km) and Mach \n"); } if(iunit==1 && i_vel_unit ==2) { printf("\n\n The input file must have two columns: altitude(nmi) and vel(ft/sec) \n"); } if(iunit==2 && i_vel_unit ==2) { printf("\n\n The input file must have two columns: altitude(feet) and vel(ft/sec) \n"); } if(iunit==3 && i_vel_unit ==2) { printf("\n\n The input file must have two columns: altitude(km) and vel(ft/sec) \n"); } } if(i_first==2 && i_second==1) { if(iunit==1 && i_vel_unit ==1) { printf("\n\n The input file must have two columns: Mach and altitude(nmi) \n"); } if(iunit==2 && i_vel_unit ==1) { printf("\n\n The input file must have two columns: Mach and altitude(feet) \n"); } if(iunit==3 && i_vel_unit ==1) { printf("\n\n The input file must have two columns: Mach and altitude(km) \n"); } if(iunit==1 && i_vel_unit ==2) { printf("\n\n The input file must have two columns: vel(ft/sec) and altitude(nmi) \n"); } if(iunit==2 && i_vel_unit ==2) { printf("\n\n The input file must have two columns: vel(ft/sec) and altitude(feet) \n"); } if(iunit==3 && i_vel_unit ==2) { printf("\n\n The input file must have two columns: vel(ft/sec) and altitude(km) \n"); } } if(i_first==1 && i_second==3) { if(iunit==1) { printf("\n\n The input file must have two columns: altitude(nmi) and q(psf) \n"); } if(iunit==2) { printf("\n\n The input file must have two columns: altitude(feet) and q(psf) \n"); } if(iunit==3) { printf("\n\n The input file must have two columns: altitude(km) and q(psf) \n"); } } if(i_first==3 && i_second==1) { if(iunit==1) { printf("\n\n The input file must have two columns: q(psf) and altitude(nmi) \n"); } if(iunit==2) { printf("\n\n The input file must have two columns: q(psf) and altitude(feet)\n"); } if(iunit==3) { printf("\n\n The input file must have two columns: q(psf) and altitude(km) \n"); } } if(i_first==2 && i_second==3) { if(i_vel_unit ==1) { printf("\n\n The input file must have two columns: Mach and q(psf) \n"); } if(i_vel_unit ==2) { printf("\n\n The input file must have two columns: vel(ft/sec) and q(psf) \n"); } } if(i_first==3 && i_second==3) { if(i_vel_unit ==1) { printf("\n\n The input file must have two columns: q(psf) and Mach \n"); } if(i_vel_unit ==2) { printf("\n\n The input file must have two columns: q(psf) and vel(ft/sec) \n"); } } printf("\n\n Enter input filename:\n"); scanf("%s",&filename[0]); pFile[0]=fopen(filename[0], "rb"); if(pFile[0] == NULL ) { printf(" Failed to open file: %s \n", filename[0]); } // printf("\n ref 1 \n"); i=0; // char string[50]; int numBytes=50; // while( fgets(string,numBytes,pFile[0] ) >= 0 ) double aaa,bbb; // printf("\n i_first=%d i_second=%d \n",i_first,i_second); while( fscanf(pFile[0],"%lf %lf", &aaa,&bbb ) >=0 ) { // printf("aaa=%lf bbb=%lf \n",aaa,bbb); if(i_first==1 && i_second==2) { altitude[i]=aaa; macht[i]=bbb; } if(i_first==1 && i_second==3) // altitude & q { altitude[i]=aaa; dp[i]=bbb; } if(i_first==2 && i_second==1) { altitude[i]=bbb; macht[i]=aaa; } if(i_first==2 && i_second==3) // velocity & q { macht[i]=aaa; dp[i]=bbb; } if(i_first==3 && i_second==1) { altitude[i]=bbb; dp[i]=aaa; } if(i_first==3 && i_second==2) { macht[i]=bbb; dp[i]=aaa; } // printf("%ld %s \n", i,string); // sscanf(string,"%lf %lf \n",&altitude[i],&macht[i] ); // printf(" %lf \n",alt[i]); if(i_first==1 || i_second==1 ) { if(iunit==1){ altitude[i]*= ft_per_nmi;} if(iunit==3){ altitude[i]*= ft_per_km;} } if(altitude[i] > 65000.) { printf("\n Warning: altitude is too high for speed of sound calculation. \n"); printf("\n The upper limit is: \n\t 20 km \n\t 10.8 nmi \n\t 65000 ft \n"); printf("\n Calculation will continue using data below the altitude limit. \n"); break; } // if( i > 2 && altitude[i]==altitude[i-1] && macht[i]==macht[i-1] ){break;} if(i>NUM){break;} i++; } last=i-1; printf("\n lines read = %ld \n",last+1); if(i_first==3 || i_second==3) { printf("\n\n Enter the speed of sound(ft/sec) \n"); scanf("%lf",&speed_of_sound); } if(i_vel_unit==1) // Mach { if(speed_of_sound<=0.) { printf("\n\n Enter the speed of sound(ft/sec) \n"); scanf("%lf",&speed_of_sound); } for(i=0;i<=last;i++) { velox[i]=macht[i]*speed_of_sound; if(velox[i] < 1.0e-20) { printf("\n R1: Error: velox[%ld]=%8.3g macht[%ld]=%8.3g speed_of_sound=%8.3g\n",i,velox[i],i,macht[i],speed_of_sound); exit(1); } } } else // feet/sec { for(i=0;i<=last;i++) { velox[i]=macht[i]; if(velox[i] < 1.0e-20) { printf("\n R2: Error: velox[%ld]= %g macht[%ld]= %g\n",i,velox[i],i,macht[i]); exit(1); } } reverse_speed(); } if(i_first==3 || i_second==3) { printf("\n find altitude \n"); find_altitude(); } strcpy(filename[1],"press1.ps"); pFile[1]=fopen(filename[1], "w"); strcpy(filename[2],"wc_spl.dat"); pFile[2]=fopen(filename[2], "w"); strcpy(filename[3],"aeroflow.txt"); pFile[3]=fopen(filename[3], "w"); strcpy(filename[4],"maxe_psd.dat"); pFile[4]=fopen(filename[4], "w"); strcpy(filename[5],"maxe_spl.dat"); pFile[5]=fopen(filename[5], "w"); strcpy(filename[6],"maxe_psd3db.dat"); pFile[6]=fopen(filename[6], "w"); } void cc_plateau_tf(void) { comp = 3.; poqcomp = 0.025/(1.+M2); prms = poqcomp*q; } void cc_reattachment_tf(void) { // // Note: reattachment is four times greater than plateau. // Wilby 131 document has decimal error in equation (34) on page 93. // // comp = 9.; poqcomp = 0.10/(1.+M2); prms = poqcomp*q; } void cc_plateau_supersonic(void) { if(ikk==0){ P2P1_ratio(); } poqtbl = (0.006/F); // turbulent boundary layer poq = poqtbl*P2P1; // plateau prms = poq*q; comp = 10.; poqcomp = poqtbl*P2P1; } void cc_separation_or_shockwave(void) { if(ikk==0){ P2P1_ratio(); } poqtbl = (0.006/F); // turbulent boundary layer // ptbl = poqtbl*q; ptbl = (0.006/F)*q; if( P2P1 < 0.593 ) { printf("\n (p2/p1) must be > 0.593 \n"); printf("\n\n Press any key to exit.\n"); getch(); exit(1); } // pshock= -1.181 + 1.713*P2P1 + 0.468*pow( P2P1, 2.); // pshock *=ptbl; pshock= ( -1.181 + 1.713*P2P1 + 0.468*pow( P2P1, 2.) )*(0.006/F)*q; if( pshock < .0 ) { printf("\n pshock must be >= 0. \n"); printf("\n\n Press any key to exit.\n"); getch(); exit(1); } // poq = poqtbl*(pshock/ptbl); // poq = poqtbl * ( -1.181 + 1.713*P2P1 + 0.468*pow( P2P1, 2.); poq = ( -1.181 + 1.713*P2P1 + 0.468*pow( P2P1, 2.) )*(0.006/F); comp = 30.; // poqcomp = poqtbl*(pshock/ptbl); poqcomp = ( -1.181 + 1.713*P2P1 + 0.468*pow( P2P1, 2.) )*(0.006/F); // printf("\n P2P1=%12.4g \n",P2P1); printf("\n poqcomp=%12.4g ",poqcomp); printf("\n poqtbl =%12.4g \n",poqtbl); printf("\n pshock =%12.4g ",pshock); printf("\n ptbl =%12.4g \n",ptbl); } void cc_auto(void) { if(imax==1) { printf("\n U=%12.4g F=%12.4g deltastar=%12.4g \n",U,F,deltastar); } constant = 4.*pow(poqcomp,2.)*comp*pow(F,1.433); constant*= pow(q,2.)*(deltastar/U); constant/= pow(splref,2.); for(i=0; i<=ilast; i++) { omega= TPI*freq[i]*deltastar/U; G[i] = constant/(1. + pow(F,2.867)*pow((comp*omega),2.)); } } void exp_plateau(void) { poqexp = 0.040/(1. + M2); prms=poqexp*q; cexp = 3.; } void exp_reattachment(void) { poqexp = 0.16/(1. + M2); prms=poqexp*q; cexp = 9.; } void exp_auto(void) { if(ikk==0) { printf("\n poqexp=%12.4g \n",poqexp);} if(imax==1) { printf("\n U=%12.4g F=%12.4g deltastar=%12.4g \n",U,F,deltastar); } constant = 4.*pow(poqexp,2.)*cexp*pow(F,1.433); constant*= pow(q,2.)*(deltastar/U); constant/= pow(splref,2.); for(i=0; i<=ilast; i++) { omega= TPI*freq[i]*deltastar/U; G[i] = constant/(1. + pow(F,2.867)*pow((cexp*omega),2.)); } } void P2P1_ratio(void) { /* P1 = static pressure upstream of shockwave at the separation point P2 = static pressure downstream of shockwave at the separation point */ printf("\n Enter frustrum angle (deg)\n"); scanf("%lf",&alpha); fprintf(pFile[3],"\n Frustrum Angle = %12.4g deg \n",alpha); alpha*=PI/180.; printf("\n\n Enter upstream Mach number (must be >= 1.0)\n"); scanf("%lf",&M1); fprintf(pFile[3],"\n Upstream Mach number = %12.4g \n",M1); if( M1 < 1.0 ) { printf("\n Upstream Mach number must be >= 1.0 \n"); printf("\n\n Press any key to exit.\n"); getch(); exit(1); } theta = alpha + asin(1./M1); s2 = pow( (sin(theta)),2.); P2P1 = ( (2.8*pow(M1,2)*s2) - 0.4 )/2.4; /* printf("\n Upstream Mach number = %12.4g \n",M1); printf("\n asin(1./M1) = %12.4g \n",asin(1./M1)); printf("\n alpha = %12.4g rad ",alpha); printf("\n theta = %12.4g rad ",theta); printf("\n sin(theta) = %12.4g ",sin(theta)); printf("\n sin2(theta) = %12.4g \n",s2); */ printf("\n P2/P1 = %12.4g \n",P2P1); if(P2P1 < 0 ) { printf("\n P2/P1 Error \n"); printf("\n\n Press any key to exit.\n"); getch(); exit(1); } } void choose_regime(void) { printf("\n Choose flow regime \n"); printf("\n 1=Attached Flow (Type I) "); printf("\n 2=Separated Flow & Shockwaves: Compression Corner (IV or V) "); printf("\n 3=Separated Flow & Shockwaves: Expansion Corner (II or III) \n"); scanf("%d",&iflow); if(iflow != 1 && iflow !=2 && iflow !=3 ) { printf("\n Input Error \n"); printf("\n\n Press any key to exit.\n"); getch(); exit(1); } if(iflow==2) // Separated Flow & Shockwaves: Compression Corner (IV or V) { printf("\n Choose Separated Flow, Compression Corner option \n"); printf("\n 1=Plateau Region, Transonic Flow (IV) "); printf("\n 2=Reattachment Region, Transonic Flow (V) "); printf("\n 3=Plateau Region, Supersonic Flow (IV) "); printf("\n 4=Separation or Reattachment Shockwave (Local Supersonic Flow)(V) \n"); scanf("%d",&icc); fprintf(pFile[3],"\n Compression Corner "); if(icc != 1 && icc !=2 && icc !=3 && icc !=4) { printf("\n Input Error \n"); printf("\n\n Press any key to exit.\n"); getch(); exit(1); } } if(iflow==3) // Separated Flow & Shockwaves: Expansion Corner (II or III) { printf("\n Choose Separated Flow, Expansion Corner option \n"); printf("\n 1=Plateau Region, Transonic and Supersonic (II)"); printf("\n 2=Reattachment Region, Transonic (III)\n"); scanf("%d",&iex); fprintf(pFile[3],"\n Expansion Corner "); if(iex != 1 && iex !=2 ) { printf("\n Input Error \n"); printf("\n\n Press any key to exit.\n"); getch(); exit(1); } } } void reynolds(void) { viscosity(); Rex = U*x/mu; // Equation (6.23) in heat transfer text if(imax==1) { printf("\n Reynolds Number = %12.4g \n",Rex); fprintf(pFile[3],"\n Reynolds Number = %12.4g \n",Rex); } deltastar = x*0.0371*pow(Rex,-0.2)*((9./7.)+0.475*M2)/pow((1.+0.13*M2),0.64); if(imax==1) { printf("\n deltastar = %12.4g ft U=%12.4g ft/sec\n",deltastar,U); fprintf(pFile[3],"\n deltastar = %12.4g ft \n",deltastar); } } void printdata(void) { printf("\n Speed of Sound = %12.4g ft/sec ",c); printf("\n Freestream velocity = %12.4g ft/sec \n",U); printf("\n Air Density = %12.3g slugs/ft^3 ",rho); printf("\n = %12.3g kg/m^3 \n",rho/conv ); printf("\n Dynamic Pressure = %12.4g psf\n",q); fprintf(pFile[3],"\n Altitude = %8.3lf ft ",h); fprintf(pFile[3],"\n = %8.3lf nmi ",h/ft_per_nmi); fprintf(pFile[3],"\n = %8.3lf km \n",h/ft_per_km); fprintf(pFile[3],"\n Speed of Sound = %12.4g ft/sec ",c); fprintf(pFile[3],"\n Freestream velocity = %12.4g ft/sec \n",U); fprintf(pFile[3],"\n Air Density = %12.3g slugs/ft^3 ",rho); fprintf(pFile[3],"\n = %12.3g kg/m^3 \n",rho/conv ); fprintf(pFile[3],"\n Dynamic Pressure = %12.4g psf\n",q); } void mach(void) { // printf("\n Enter the Mach number \n"); // scanf("%lf",&M); // fprintf(pFile[3],"\n Mach number = %8.3lf \n",M); WT_ratio=1.; M2 = pow(M,2.); F= 0.5 + WT_ratio*(0.5 + 0.09*M2) + 0.04*M2; } void material(void) { cref=200000./12.; printf("\n Enter material: \n 1=aluminum 2=steel 3=titanium 4=magnesium 5=other \n"); scanf("%d",&imat); if(imat !=1 && imat !=2 && imat != 3 && imat != 4 && imat != 5) { printf("\n Input Error \n"); printf("\n\n Press any key to exit.\n"); getch(); } if(imat==1) // Aluminum { fprintf(pFile[3],"\n Aluminum \n"); W=0.1; cmat=cref; } if(imat==2) // Steel { fprintf(pFile[3],"\n Steel \n"); W=0.29; cmat=cref; } if(imat==3) // Titanium { fprintf(pFile[3],"\n Titanium \n"); W=0.16; cmat=cref; } if(imat==4) // Magnesium { fprintf(pFile[3],"\n Magnesium \n"); W=0.063; cmat=cref; } if(imat==5) // Other { fprintf(pFile[3],"\n Other Material \n"); printf("\n Enter cylindrical skin weight density (lbm/in^3) \n"); char sW[16]; scanf("%s",&sW); sscanf(sW,"%lf",&W); printf("\n Enter elastic modulus (lbf/in^2) \n"); char sE[16]; scanf("%s",&sE); sscanf(sE,"%lf",&E); cmat = sqrt( (E/W)*(32.2/12.) ); fprintf(pFile[1],"\n\n Elastic Modulus = %10.3g lbf/in^2 \n",E); } fprintf(pFile[3],"\n\n Speed of Sound in Material = %10.3g ft/sec ",cmat); printf("\n\n Speed of Sound in Material = %10.3g ft/sec \n\n",cmat); scale = cmat/cref; } void franken() { if(ikk==0) { printf("\n Enter diameter (in)\n"); char sdia[12]; scanf("%s",&sdia); sscanf(sdia,"%lf",&stationdiam); fprintf(pFile[3],"\n Diameter = %7.3g in", stationdiam); stationdiam/=12.; printf("\n Enter skin thickness (in) \n"); char sthick[16]; scanf("%s",&sthick); sscanf(sthick,"%lf",&thickness); fprintf(pFile[3],"\n Skin Thickness = %12.3f inches \n",thickness); thickness/=12.; fprintf(pFile[3],"\n\n Weight Density = %10.3g lbm/in^3 ",W); W*=pow(12.,3.); W*=thickness; fprintf(pFile[3],"\n\n Weight Density (per area) = %12.3f lbm/ft^2 \n",W); } //************************************************************** if(ikk==0) { fr[0]=1.0; fr[1]=2.0; fr[2]=2.7; fr[3]=3.0; fr[4]=3.3; fr[5]=3.6; fr[6]=3.7; fr[7]=3.78; fr[8]=3.9; fr[9]=4.0; fr[10]=4.1; fr[11]=4.3; fr[12]=4.7; ar[0]=-158; ar[1]=-146; ar[2]=-140.5; ar[3]=-136.; ar[4]=-130.; ar[5]=-121.; ar[6]=-119.; ar[7]=-118.8; ar[8]=-119.; ar[9]=-120.; ar[10]=-121.7; ar[11]=-123.; ar[12]=-124.5; for(i=0; i<=12; i++) { fr[i]= log10( (pow(10.,fr[i])*scale) ); } } //************************************************************** int jlast = 12; for(i=0; i<= ilast; i++) { sfr = log10( f[i]*stationdiam ); // printf("\n %12.4e %12.4e\n",sfr,f[i]); if( sfr <= fr[0]) { va= ar[0]; } if(sfr >= fr[jlast] ) { va= ar[jlast]; } if( sfr > fr[0] && sfr < fr[jlast] ) { for(j=0; j<=11; j++) { if(sfr >= fr[j] && sfr <= fr[j+1] ) { len = fr[j+1] -fr[j]; k2 = (sfr -fr[j])/len; k1=1.-k2; va = k1*ar[j] + k2*ar[j+1]; break; } } } vdb = va - 20.*log10(W) + spl[ikk][i]; al = 1.*pow(10.,(vdb/20.)); df = (pow(2.,(1./6.))-1./pow(2.,(1./6.)))*f[i]; // psd[i] = 0.5*pow(al,2.)/df; // assume vdb is in terms of RMS psd[i] = pow(al,2.)/df; /* if(f[i] > 98. && f[i] < 102.) { printf(" f=%8.3lf psd=%9.3e df=%9.3lf va=%9.3lf 20logW=%8.4lf spl=%8.3lf\n",f[i],psd[i],df,va,20.*log10(W),spl[i]); } */ if(maxpsd[i] -0.9999 ) { rb+= ( maxpsd[i+1] * f[i+1]- maxpsd[i]*f[i])/( s[i]+1.); } else { rb+= maxpsd[i]*f[i]*log( f[i+1]/f[i]); } if(f[i+1] < 2050.){ra=rb;} } grms=sqrt(ra); grmsb=sqrt(rb); // printf("\n cmat=%18.4g PI=%8.4g stationdiam=%8.4g \n",cmat,PI,stationdiam); printf("\n Ring Frequency = %12.4g Hz \n",cmat/(PI*stationdiam)); fprintf(pFile[3],"\n Ring Frequency = %12.4g Hz \n",cmat/(PI*stationdiam)); fprintf(pFile[3],"\n\n Max Envelope PSD Overall Skin Vibration Level (up to 2000 Hz) = %12.4g GRMS\n",grms ); printf("\n Max Envelope PSD Overall Skin Vibe Level (up to 2000 Hz) = %12.4g GRMS\n",grms ); fprintf(pFile[3],"\n Max Envelope PSD Overall Skin Vibration Level (up to 10,000 Hz) = %12.4g GRMS\n",grmsb ); printf("\n Max Envelope PSD Overall Skin Vibe Level (up to 10,000 Hz)= %12.4g GRMS\n",grmsb ); } void reverse_speed(void) { for(ikk=0; ikk<=last; ikk++) { h=altitude[ikk]; double LH1=11. *ft_per_km; double LH2=20. *ft_per_km; if( h < LH1) // troposphere { double Tz=288.; double L=6.5; double gamma=1.402; double ROM = (8314.3/28.97); if( h < 0.) { printf("\n Error \n"); printf("\n\n Press any key to exit.\n"); getch(); exit(1); } c=sqrt(fabs(gamma*ROM*(Tz-L*(h/(ft_per_km)))))*3.28; // printf("\n gamma=%12.4e ROM=%12.4e Tz=%12.4e L=%12.4e h=%12.4e\n",gamma,ROM,Tz,L,h); } if( h >= LH1 && h < LH2 ) // lower stratosphere { c = 295.*3.28; } if( h >= LH2 ) { printf("\n Error: altitude is too high for speed of sound calculation. \n\n"); printf(" Maximum allowable altitude = %12.6g ft \n",LH2); printf(" = %12.5g nmi \n",LH2/ft_per_nmi); printf(" = %12.5g km \n\n",LH2/ft_per_km); printf("\n\n Press any key to exit.\n"); getch(); exit(1); } if(speed_of_sound > 0.) { c=speed_of_sound; } macht[ikk]/=c; if(macht[ikk]<1.0) { printf("\n mach error \n"); printf("\n ikk=%d h=%12.4g mach=%12.4g c=%12.4g \n",ikk,h,macht[ikk],c); exit(1); } } } void find_altitude(void) { // have velox ft/sec want altitude[i] ft... double alt[30],dens[30]; double c1,c2,x,len; for(i=0; i<=20; i++) { alt[i]=double(i)*ft_per_km; } dens[0]=1.226; dens[1]=1.112; dens[2]=1.007; dens[3]=0.9096; dens[4]=0.8195; dens[5]=0.7365; dens[6]=0.6600; dens[7]=0.5898; dens[8]=0.5254; dens[9]=0.4666; dens[10]=0.4129; dens[11]=0.3641; dens[12]=0.3104; dens[13]=0.2652; dens[14]=0.2266; dens[15]=0.1936; dens[16]=0.1654; dens[17]=0.1413; dens[18]=0.1207; dens[19]=0.1032; dens[20]=0.0881; for(i=0;i<=last;i++) { if(velox[i] < 1.0e-20) { printf("\n Error: velox[%ld]= %g \n",i,velox[i]); exit(1); } rho = 2.*dp[i]/pow(velox[i],2.); // printf("\n rho=%8.3g dp=%8.3g velox=%8.3g \n",rho,dp[i],velox[i]); rho/=conv; for(j=0;j<20;j++) { if(rho <= dens[j] && rho >= dens[j+1] ) { x=rho-dens[j]; len=dens[j+1]-dens[j]; c2=x/len; c1=1.-c2; altitude[i]=c1*alt[j]+c2*alt[j+1]; // printf("\n rho=%8.3g alt=%8.3g \n",rho,altitude[i]); break; } } } }