clear close all %************************************************************************** %****************** Início dos dados de entrada do usuário **************** %************************************************************************** % Mach desejado na saida do bocal (para resultado 3D) %M_e=input('Qual o Mach de saida desejado? '); answer=inputdlg({'Mach de exaustao desejado'}); M_e=str2num(answer{1}); % Número de características com origem na garganta n=100; % Ângulo inicial da expansão na garganta theta_min=0.375; %************************************************************************** %****************** Final dos dados de entrada do usuário ***************** %************************************************************************** % O Matlab trabalha com ângulos em radianos: theta_min=theta_min*pi/180; % Razão de calores específicos g=1.4; % Geração de tabela da função de Prandtl_Meyer e de Razão entre Áreas M_t=1:0.01:10; nu_t=sqrt((g+1)/(g-1))*atan(sqrt((g-1)/(g+1)*(M_t.^2-1)))-atan(sqrt(M_t.^2-1)); AR_t=sqrt(1./M_t.^2.*(2/(g+1)*(1+(g-1)/2*M_t.^2)).^((g+1)/(g-1))); % Cálculo do Mach na saida do bocal para a construção 2D AR_3D=interp1(M_t,AR_t,M_e,'spline'); AR_2D=sqrt(AR_3D); M_e=interp1(AR_t,M_t,AR_2D,'spline'); %M_e=2.4; % Cálculo do ângulo máximo do leque de expansão na garganta nu_g=interp1(M_t,nu_t,M_e,'spline'); theta_max=nu_g/2; %Cálculo do número de pontos count=n+1; j=0; k=0; while count>=2 for i=1:count j=j+1; end k=k+1; count=count-1; end %Pré-dimensionamento das variáveis (deixa o código mais rápido) num=zeros(j,1); theta=zeros(j,1); nu=zeros(j,1); Kn=zeros(j,1); Kp=zeros(j,1); M=zeros(j,1); mu=zeros(j,1); x=zeros(j,1); y=zeros(j,1); cn=zeros(j,1); cp=zeros(j,1); xp=zeros(k,1); yp=zeros(k,1); %************************************************************************** %M Mach local %theta ângulo local da direção do escoamento (em relação à horizontal) %nu ângulo de Prandtl Meyer local %Kp ângulo local da linha característica ao longo da qual theta+nu é constante %Kn ângulo local da linha característica ao longo da qual theta-nu é constante %mu ângulo local da linha de Mach %************************************************************************** %Ponto de início do perfil do bocal n_perfil=1; xp(n_perfil)=0; yp(n_perfil)=100; %Cálculo das linhas de Mach com origem na garganta for i=1:n num(i)=i; theta(i)=theta_min+(theta_max-theta_min)/(n-1)*(i-1); nu(i)=theta(i); Kn(i)=theta(i)+nu(i); Kp(i)=theta(i)-nu(i); M(i)=interp1(nu_t,M_t,nu(i),'spline'); mu(i)=asin(1/M(i)); cn(i)=tan(theta(i)-mu(i)); cp(i)=tan(theta(i)+mu(i)); if i==1 %Na linha de simetria x(i)=-100/cn(i); y(i)=0; else %Para os pontos na sequência até o limite do leque de expansão x(i)=(-x(i-1)*(cp(i-1)+cp(i))/2-(100-y(i-1)))/... (cn(i)-(cp(i-1)+cp(i))/2); y(i)=(x(i)-x(i-1))*(cp(i-1)+cp(i))/2+y(i-1); end end % Cálculo do primeiro ponto do perfil do bocal i=i+1; num(i)=i; theta(i)=theta(i-1); nu(i)=nu(i-1); Kn(i)=Kn(i-1); Kp(i)=Kp(i-1); M(i)=M(i-1); mu(i)=mu(i-1); x(i)=(100-y(i-1)+x(i-1)*cp(i-1))/(cp(i-1)-tan(theta(i))); y(i)=x(i)*tan(theta(i))+100; n_perfil=n_perfil+1; xp(n_perfil)=x(i); yp(n_perfil)=y(i); % Cálculo das demais linhas de Mach for j=1:n-1 for k=1:n-j i=i+1; num(i)=i; Kn(i)=Kn(i-(n-(j-1))); Kp(i)=-Kn(j+1); theta(i)=1/2*(Kn(i)+Kp(i)); nu(i)=1/2*(Kn(i)-Kp(i)); M(i)=interp1(nu_t,M_t,nu(i),'spline'); mu(i)=asin(1/M(i)); cn(i)=tan(theta(i)-mu(i)); cp(i)=tan(theta(i)+mu(i)); m=i-(n-(j-1)); if k==1 %Na linha de simetria x(i)=(x(m)*cn(m)-y(m))/cn(m); y(i)=0; else %Para os pontos na sequência até o limite do leque de expansão x(i)=(x(m)*(cn(m)+cn(i))/2-x(i-1)*(cp(i-1)+cp(i))/2-... (y(m)-y(i-1)))/((cn(m)+cn(i))/2-(cp(i-1)+cp(i))/2); y(i)=(x(i)-x(i-1))*(cp(i-1)+cp(i))/2+y(i-1); end end % Cálculo dos pontos do perfil do bocal i=i+1; m=i-(n-(j-1)); num(i)=i; Kn(i)=Kn(i-1); Kp(i)=Kp(i-1); theta(i)=theta(i-1); nu(i)=nu(i-1); M(i)=M(i-1); mu(i)=mu(i-1); x(i)=(y(m)-y(i-1)+x(i-1)*cp(i-1)-x(m)*tan(theta(m)))/... (cp(i-1)-tan(theta(m))); y(i)=(x(i)-x(m))*tan(theta(i))+y(m); n_perfil=n_perfil+1; xp(n_perfil)=x(i); yp(n_perfil)=y(i); end % dados=[num 180/pi*Kn 180/pi*Kp 180/pi*theta 180/pi*nu M 180/pi*mu]; % save('dados.txt','dados','-ascii') perfil=[xp yp zeros(size(xp,1),1)]; dlmwrite('perfil.txt',perfil,'precision',8) plot(xp,yp,'k') axis('equal') hold on plot(xp,-yp,'k') % figure % plot(x,y,'.') % axis('equal') formatSpec = 'Comprimento do bocal: %4.0f mm \n'; fprintf(formatSpec,round(max(xp),0)) formatSpec = 'Diâmetro de saída do bocal: %4.0f mm \n'; fprintf(formatSpec,2*round(max(yp),0)) AR_RES=pi*round(max(yp),0)^2/(pi*100^2); formatSpec = 'Relação entre áreas: %4.2f \n'; fprintf(formatSpec,AR_RES) M_RES=interp1(AR_t,M_t,AR_RES,'spline'); formatSpec = 'Mach resultante: %4.1f \n'; fprintf(formatSpec,M_RES)