%%Lance Simms 1/30/06 %%This script will attempt to solve the Transonic Disturbance equations for %%a Circular Chord Thin Airfoil. plot_density=0; M_inf=0; M_inf=input('Enter M_inf (Mach_numer) '); V_inf=M_inf; num_steps=200 num_steps=input('Enter Number of Steps '); %Ratio of Specific Heats and quantities gamma=1.4; rho_inf=1; a_inf=V_inf/M_inf; p_inf=rho_inf*a_inf^2/gamma; gamma_exp=(gamma/(gamma-1)); %%GRID INDEX DEFINITIONS I=2; J=2; %Beginning Indices IL=51; JL=51; %Maximum Indices %Lead edge occurs at 3/10 of x and trailing edge at 7/10+1 I_lead_edge=(IL-1)*.3+1; I_trail_edge=(IL-1)*.7+2; %SEPERATE X INDICES FOR 3 REGIONS (SUBSONIC--SONIC---SUBSONIC) I_before_foil=[2:I_lead_edge]; I_over_foil=[I_lead_edge:I_trail_edge-1]; I_after_foil=[I_trail_edge:IL]; IO=[1:IL]; JO=[1:JL]; %Full Range Indices IM=[2:IL-1]; JM=[2:JL-1]; %Interior Indices %DISTANCES c=1; %Chord length Dy=50*c; Dx=50*c; %Max DIistance del_y_min=.06*c/10; del_x_min=c/20; %The minimum separation in y near the foil %The Grid was solved with the exponential mesh stretching formula discussed %in the text on page 11 of Chapter 5 %The value of k was determined [kx,ky]=Stretching_Function(I,J,I_lead_edge,JL,Dx,Dy,del_y_min,del_x_min); J_fac=(JO-1)/(JL-1); I_fac=(I_before_foil-1)/(I_lead_edge-1); %GENERATE THE Y SPACING-----Spaced with the stretchin function y=Dy*(exp(ky*J_fac)-1)/(exp(ky)-1); %GENERATE THE X SPACING-----Spaced with the stretchin function I_before_foil=[2:I_lead_edge]; I_over_foil=[I_lead_edge:I_trail_edge-1]; I_after_foil=[I_trail_edge:IL]; X_trail=Dx*(exp(kx*I_fac)-1)/(exp(kx)-1)+1; X_over=(0:c/20:1); X_lead=-fliplr(X_trail)+1; x=[X_lead,X_over,X_trail]; %GENERATE THE MESH [X,Y]=meshgrid(x,y); %DRAW IT IF YOU"D LIKE %figure; set(gca,'XLIM',[min(x),max(x)]); %set(gca,'XTICK',x); set(gca,'YTICK',y); grid; %figure; %PUT THE AIRFOIL AT Y=0 tau=.03; c=1; [xo,yo,R]=Circular_Chord(tau,c); ys=zeros(1,IL); ys(I_over_foil)=sqrt(R^2-(x(I_over_foil)-xo).^2)+yo; dys=zeros(1,IL); dys(I_over_foil)=-(x(I_over_foil)-xo)./... sqrt(R^2-(x(I_over_foil)-xo).^2); %SET UP THE INITIAL VALUES OF THE POTENTIAL PHI=ones(JL,IL); %%%SET UP DENOMINATOR TERMS FOR Y----ONLY NEED ONE COLUMN HERE %%INDIVIDUAL yj1=y(3:JL)'; yjm1=y(1:JL-2)'; yj0=y(2:JL-1)'; Yj1=Y(3:JL,3:IL-1); Yjm1=Y(1:JL-2,3:IL-1); Yj0=Y(2:JL-1,3:IL-1); %%DENOMINATORS y_j1_j0=(yj1-yj0).*(yj1-yjm1); y_j0_jm1=(yj0-yjm1).*(yj1-yjm1); Y_j1_j0=(Yj1-Yj0).*(Yj1-Yjm1); Y_j0_jm1=(Yj0-Yjm1).*(Yj1-Yjm1); %%SET UP DENOMINATOR TERMS FOR X-----ONLY NEED ONE ROW HERE %%INDIVIDUAL xi1=X(2:JL-1,4:IL); xim2=X(2:JL-1,1:IL-3); xi0=X(2:JL-1,3:IL-1); xim1=X(2:JL-1,2:IL-2); %%DENOMINATORS x_i1_i0_c=(xi1-xi0).*(xi1-xim1); x_im1_im2_b=(xim1-xim2).*(xi0-xim2); x_i0_im1_c=(xi0-xim1).*(xi1-xim1); x_i0_im1_b= (xi0-xim1).*(xi0-xim2); %%Y difference prefactors---These do not change-----b*phi_i_j-1--c*phi_i_j+1 b(1:JL-2)=(1./y_j0_jm1); b(JL-1)=0; c(2:JL-1)=(1./y_j1_j0); c(1)=-1; %%Plotting set (gcf,'doublebuffer','on'); if plot_density==1 subplot(3,1,3); im=imagesc(x,y,PHI); end cp_ind=[min(I_over_foil)-3:max(I_over_foil)+3]; subplot(3,1,1); cpplot=plot(x(cp_ind),ones(size(cp_ind))) subplot(3,1,2); Residplot=semilogy([1],[1]) R_max_old=1; %%Begin the loop for n=1: num_steps for i=2: IL-1 %%Case where i=2, we cannot use backward differencing. No A_im1_j %%term if i==2 Aij=1-M_inf^2-(gamma+1)*M_inf^2*(PHI(JM,3)-PHI(JM,1))./... (X(JM,3)-X(JM,1))./V_inf; mu=(Aij<0); x_i1_i0_c_first=(X(JM,3)-X(JM,2)).*(X(JM,3)-X(JM,1)); x_i0_im1_c_first=(X(JM,2)-X(JM,1)).*(X(JM,3)-X(JM,1)); %%LHS OF EQN---TERMS WITH PHI_i--bijPhi_i_j-1+aijPhi_i_j+cijPhi_i_j+1 a(2:JL-1)=(1-mu).*Aij.*(-1./x_i1_i0_c_first-1./x_i0_im1_c_first)... -1./y_j1_j0-1./y_j0_jm1; a(1)=1; a(JL)=1; %%RHS OF EQUATION----TERMS INVOLVING PHI_i+1 and PHI_i-1 f(2:JL-1)=-(1-mu).*Aij.*PHI(JM,3)./x_i1_i0_c_first... -(1-mu).*Aij.*PHI(JM,1)./x_i0_im1_c_first; f(1)=-del_y_min*V_inf*dys(i); f(JL)=1; %%Do the solver [L,U,alph,gamm]=LU_decomposition(a,b,c,JL); PHI(:,i)=LU_solve(alph,b,gamm,f,JL); if plot_density==1 set(im,'cdata',flipud(PHI),'xdata',x,'ydata',y); drawnow; end else %%Now we can use backward differencing----Get on it Aij=1-M_inf^2-(gamma+1)*M_inf^2*(PHI(JM,i+1)-PHI(JM,i-1))./... (X(JM,i+1)-X(JM,i-1))./V_inf; Aim1j=1-M_inf^2-(gamma+1)*M_inf^2*(PHI(JM,i)-PHI(JM,i-2))./... (X(JM,i)-X(JM,i-2))./V_inf; muij=(Aij<0); muim1j=(Aim1j<0); %%LHS %%LHS OF EQN---TERMS WITH PHI_i--bijPhi_i_j-1+aijPhi_i_j+cijPhi_i_j+1 a(2:JL-1)= (1-muij).*Aij.*(-1./x_i1_i0_c(:,i-2)-1./x_i0_im1_c(:,i-2))... +(muim1j).*Aim1j.*(1./x_i0_im1_b(:,i-2))+... (-1./y_j1_j0-1./y_j0_jm1); a(1)=1; a(JL)=1; %%RHS OF EQUATION----TERMS INVOLVING PHI_i+1 and PHI_i-1 f(2:JL-1)=-(1-muij).*Aij.*PHI(JM,i+1)./x_i1_i0_c(:,i-2)... -(1-muij).*Aij.*PHI(JM,i-1)./x_i0_im1_c(:,i-2)... +muim1j.*Aim1j.*(1./x_i0_im1_b(:,i-2)+1./x_im1_im2_b(:,i-2)).*PHI(JM,i-1)... -muim1j.*Aim1j.*(1./x_im1_im2_b(:,i-2)).*PHI(JM,i-2); f(1)=-del_y_min*V_inf*dys(i); f(JL)=1; %%Do the solver [L,U,alph,gamm]=LU_decomposition(a,b,c,JL); PHI(:,i)=LU_solve(alph,b,gamm,f,JL); %%PLOT if plot_density==1 set(im,'cdata',flipud(PHI),'xdata',x,'ydata',y); drawnow; end end end %The span in x has been covered. Plot cp vs x u=V_inf+(PHI(1,cp_ind+1)-PHI(1,cp_ind-1))./... (x(cp_ind+1)-x(cp_ind-1)); v=(PHI(2,cp_ind)-PHI(1,cp_ind))./... (del_y_min); vel=(u.^2+v.^2)./V_inf^2; p=p_inf*(1-(gamma-1)*(M_inf^2.*(vel-1))/2).^gamma_exp; neg_cp=-(p-p_inf)/(.5*rho_inf*V_inf^2); set(cpplot,'ydata',neg_cp); drawnow; %Get the maximum Residual Rij Aij=1-M_inf^2-(gamma+1)*(M_inf^2).*(PHI(JM,4:IL)-PHI(JM,2:IL-2))./... (X(JM,4:IL)-X(JM,2:JL-2))./V_inf; Aim1j=1-M_inf^2-(gamma+1)*(M_inf^2).*(PHI(JM,3:IL-1)-PHI(JM,1:IL-3))./... (X(JM,3:IL-1)-X(JM,1:IL-3))./V_inf; muij=(Aij<0); muim1j=(Aim1j<0); R=(1-muij).*Aij.*((PHI(JM,4:IL)-PHI(JM,3:IL-1))./(x_i1_i0_c)-... (PHI(JM,3:IL-1)-PHI(JM,2:IL-2))./(x_i0_im1_c))... +muim1j.*Aim1j.*((PHI(JM,3:IL-1)-PHI(JM,2:IL-2))./(x_i0_im1_b)-... (PHI(JM,2:IL-2)-PHI(JM,1:IL-3))./(x_im1_im2_b))... +((PHI(3:JL,3:IL-1)-PHI(2:JL-1,3:IL-1))./(Y_j1_j0)-... (PHI(2:JL-1,3:IL-1)-PHI(1:JL-2,3:IL-1))./(Y_j0_jm1)); R_max=max(max(abs(R))) if n==1 R_max_arr=[R_max]; n_arr=[n] else R_max_arr=[R_max_arr,R_max]; n_arr=[n_arr,n]; end set(Residplot,'xdata',n_arr,'ydata',R_max_arr); end %%Generate Final Figures to turn in %%CP PLOT figure; subplot(2,1,1); cpplot=plot(x(cp_ind),neg_cp) xlabel('x/c'); ylabel('-Cp'); title(sprintf('Cp vs x for flow with M_{inf}=%2.3f',M_inf)); %%Residual PLOT subplot(2,1,2); semilogy(n_arr,R_max_arr); xlabel('Step Number'); ylabel('Residual'); Title(sprintf('Maximum residual for M_{inf}=%2.3f',M_inf));