%FINAL PROJECT #9--------------COMPUTATION OF LINEARIZED EULER EQUATIONS %Lance Simms 12/2005 % %u=velocity in x direction v=velocity in y direction %Flow is flowing in from left. Assumed to be uniform rho=1, u=M_inf, v=0 %Linearization about p=1, u=M_inf, v=0 % %Thin airfoil present having shape given by % ywall=Tx(1-x)/2 01 % %Euler Equations are written as % % DQ DQ DQ % -- + A -- + B -- = 0 % Dt Dx Dy % % where: | 0 1 0 | | rho | % A = | -M_inf^2 2 M_inf 0 |, Q = | rho u | % | 0 0 M_inf | | rho v | % % | 0 0 1 | % B = | 0 0 M_inf | % | 1 0 0 | % Inputs: % M_inf : Free Stream Mach number % nx : Number of points on airfoil (x=0,1) % cfl : dt = min(dx,dy)*cfl/(M_inf+1 + 1) % where Spectral radius of A is (M_inf+1) and B is 1 % nmax : Number of time steps to go % nout : Output frequency of plot of cp (0 for no plots) % meth : 1: Explicit Euler / 2: RK4 / 3: Implicit Euler % order : 0 : central 2nd / 1 : 1st order upwind / 2 : 2nd order upwind % ichar : 0 : simple bc / 1 : char bc / -1 simpler bc % % clear; %MACH NUMBER M_inf=.8; M_inf=input('Enter M_inf'); %Grid Generation xmin=-1.0; xmax=3.0; ymin=0.0; ymax=2.0; nx=10; nx=input('Enter number of points on body:'); dx=1.0/(nx-1); dy=dx; x=[xmin:dx:xmax]'; y=[ymin:dy:ymax]'; jmax=size(x,1); kmax=size(y,1); %Generate meshgrid having x values as rows and y values as columns [X,Y]=meshgrid(x,y); jle=nx; jte=2*nx-1; fprintf('Leading edge point= %g trailing edge point = %g \n', jle,jte); fprintf('X of Leading edge point =%g X of trailing edge point = %g \n', x(jle),x(jte)); %Surface Definition: tau is the thickness tau=0.1; %tau=input ('tau); s=[jle:jte]; ys=zeros(jmax,1); ys(s)=tau*0.5*x(s).*(1.0-x(s)); %Try to figure out indices of airfoil foil=0; for xind=1:length(x) newfoil=find(y2 residslope=resid(nstep)-resid(nstep-1); set(residslopetext,'string',sprintf('Slope=%2.8f',residslope)); end refresh; drawnow; %[flag]=figflag(gcf,0); pause(0.1); end end %%EULER EXPLICIT if (meth==1) for nstep=1:nmax_tot resid(nstep) = ... (norm(R(:,:,1)) + norm(R(:,:,2)) + norm(R(:,:,3)))/(3*jmax*kmax); resid(nstep) = resid(nstep)/resid_0; R=rhs_flux_split_so(Q,J,K,jmax,kmax,dx,dy,Ap,Am,Bp,Bm,epse); Q(JM,KM,:) = Q(JM,KM,:) + dt*R(JM,KM,:); %Boundary Conditions Q=bc8(Q,JM,KM,jmax,kmax,M_inf); %Position Airfoil in Plot PlotQ=Q(:,:,1)'; PlotQ(foil)=1; set(gcf,'CurrentObject',im); set(im,'cdata',flipud(PlotQ),'xdata',x,'ydata',y); set(maxp,'string',sprintf('p_{max}=%2.3f',max(max(Q(:,:,1))))); set(maxpu,'string',sprintf('pu_{max}=%2.3f',max(max(Q(:,:,2))))); set(maxpv,'string',sprintf('pv_{max}=%2.3f',max(max(Q(:,:,3))))); % set(gca,'title',text('string',sprintf('b=%2.2f',max(max(max(Q(:,:,1)))),... % 'color','r','fontsize',14))); % set(gcf,'CurrentObject',resid); set(residp,'XData',[1:nstep],'YData',resid); set(residtext,'string',sprintf('Residual=%2.6f',resid(nstep))); if nstep>2 residslope=resid(nstep)-resid(nstep-1); set(residslopetext,'string',sprintf('Slope=%2.8f',residslope)); end colorbar; refresh; drawnow; pause(0.1); end end %%EULER IMPLICIT if meth == 2 % Delta Form Implicit Approximate Factorization [a_x,b_x,c_x,d_x,e_x] = fill_x(JM,KM,A,dx,dt,jmax,kmax,epse); [a_y,b_y,c_y,d_y,e_y] = fill_y(JM,KM,B,dy,dt,jmax,kmax,epse); for nstep=1:nmax_tot resid(nstep) = ... (norm(R(:,:,1)) + norm(R(:,:,2)) + norm(R(:,:,3)))/(3*jmax*kmax); resid(nstep) = resid(nstep)/resid_0; % [I + dt A Del_cen_x ]^(-1) R % First Step of Implicit Factored Scheme : x block tridiagonal step R=rhs_flux_split_so(Q,J,K,jmax,kmax,dx,dy,Ap,Am,Bp,Bm,epse); R = dt*R; R = btrix_3x3(2,jmax-1,KM,b_x,c_x,d_x,R); % [I + dt B Del_cen_y ]^(-1) R % Second Step of Implicit Factored Scheme : y block tridiagonal step DelQ = btriy_3x3(2,kmax-1,JM,b_y,c_y,d_y,R); % Q^(n+1) = Q^n + Delta Q^n Q(JM,KM,:) = Q(JM,KM,:) + DelQ(JM,KM,:); %Boundary Conditions Q=bc8(Q,JM,KM,jmax,kmax,M_inf); %subplot(2,1,1); %quiver(X,Y,Q(:,:,2)',Q(:,:,3)'); %contourf(X,Y,Q(:,:,2)',50); %xlabel('rho'); PlotQ=Q(:,:,1)'; PlotQ(foil)=1; %subplot(2,1,2); %imagesc(x,y,flipud(PlotQ)); %imagesc(x,y(kmax:-1:1),Q(:,:,1)'); %set(gcf,'CurrentObject',im); %imagesc(x,y,flipud(PlotQ)); %imagesc(x,y(kmax:-1:1),Q(:,:,1)'); set(im,'cdata',flipud(PlotQ),'xdata',x,'ydata',y); set(maxp,'string',sprintf('p_{max}=%2.3f',max(max(Q(:,:,1))))); set(maxpu,'string',sprintf('pu_{max}=%2.3f',max(max(Q(:,:,2))))); set(maxpv,'string',sprintf('pv_{max}=%2.3f',max(max(Q(:,:,3))))); set(residtext,'string',sprintf('Residual=%2.9',resid(nstep))); %text(.7,1.1,sprintf('max_p=%2.2f',max(max(Q(:,:,1))))); streamslice(X,Y,Q(:,:,2)',Q(:,:,3)') % set(gcf,'CurrentObject',resid); set(residp,'XData',[1:nstep],'YData',resid); if nstep>2 residslope=resid(nstep)-resid(nstep-1); set(residslopetext,'string',sprintf('Slope=%2.9f',residslope)); end refresh; drawnow; pause(0.1); end end