%This script will attempt to solve the Quasi 1-d Euler equations % Cv=5; %Specific Heat in BTU/lb/deg R gamma=1.4; %Ratio of specific heats gamma_fac=(gamma-1)/(gamma+1); %Used coefficient Beta=gamma-1; %Used coefficient sigma=1; %Constant used in ROE finite_vol=1; %1=finite volume 2=finite difference flux_method=2; %1=Modified Staeger 2=Roe Method ghost_cell=2; %1=use ghost cell 2=no ghost cell flow_case=2; %1=No flow 2=pressure dropped at exit %The grid extends from -.5 to 1.5 and it's divided into 41 equally spaced points IL=41; IO=[1:IL]; IM=[2:IL-1]; %Surfaces %The grid encloses either volumes or represents the surfaces if finite_vol==1 IVL=IL-1; IVO=[1:IVL-1]; IVM=[2:IVL-2]; %IL-1 Volumes else IVL=IL; IVO=[1:IL]; IVM=[2:IL-1]; %IL Surfaces end xleft=-.5; xright=1.5; delx=(xright-xleft)/(IL-1); x=[xleft:delx:xright]; xvols=x(2:IL)-.5*delx; num_vols=IVL; %The leading edge of the chord is at 0 and the trailing edge is at 1 I_lead_edge=(IL-1)*.25+1; I_trail_edge=(IL-1)*.75+1; %Get indices of before, over, and after chord I_before_chord=[1:I_lead_edge]; I_over_chord=[I_lead_edge+1:I_trail_edge]; I_after_chord=[I_trail_edge+1:IL]; %Figure out the xvalues of the chord and the surface y value S(x) tau=.03; c=1; Surf=-.1*ones(1,IL); [xo,yo,R]=Circular_Chord_Quas_1d(tau,c); Surf(I_over_chord)=sqrt(R^2-(x(I_over_chord)-xo).^2)+yo-.1; %Now get the actual width of the channel and the corresponding volumes Surf=0-Surf; if ghost_cell==1 Ghost_cell=delx.*(Surf(2)+Surf(1))/2; Vol=[Ghost_cell,delx.*(Surf(2:IL)+Surf(1:IL-1))./2]; else Vol=delx.*(Surf(2:IL)+Surf(1:IL-1))./2; end %Set up the flux vectors. Initial conditions are determined using Newton's %method and isentropic relations %F is the flux vector---- F(1,:)=p F(2,:)=pu F(3,:)=e U=zeros(3,IVL); Vols=repmat(Vol,3,1); Surfs=repmat(Surf,3,1); Q=zeros(3,IVL); Q(1,IVL)=0; Q(3,IVL)=0; %Constants for flow cases if flow_case==1 %Flow_case=1-----NO FLOW po=2117; %Entrance Pressure in lb/ft^2 Mo=0; %Mach Number at Entrance=0 To=531.2; %Temperature at entrance p=po*ones(1,IVL); %Initial Pressure rho=p./((gamma-1)*Cv*To); %Initial density u=zeros(1,IVL); %Initial Velocity a=sqrt(2*gamma*gamma_fac*Cv*To); else po=2117; %Entrance Pressure in lb/ft^2 Mo=0; %Mach Number at Entrance=0 To=531.2; %Temperature at entrance p=po*ones(1,IVL); %Initial Pressure p(IVL)=1400; %Instantaneously drop pressure rho=p./((gamma-1)*Cv*To); %Initial density u=zeros(1,IVL); %Initial Velocity a=sqrt(2*gamma*gamma_fac*Cv*To); end %Set up conserved quantities U(1,:)=rho; U(2,:)=rho.*u; U(3,:)=rho.*(Cv.*To+.5.*(u.^2)); figure; set(gcf,'doublebuffer','on'); num_steps=500; for t=0:num_steps %Update Pressure first since it will be needed throughout p=(gamma-1)*U(1,:).*(U(3,:)./U(1,:)-.5*(U(2,:)./U(1,:)).^2); %pressure %-------------------------------------------------------------STEGER-WARMING METHOD if flux_method==1 %Exit Boundary Condition %State and Flux vectors and Q vector U_b=(U(:,1:IVL-1)+U(:,2:IVL))./2; %Averages F_msw=zeros(3,IVL-1); %Backward Flux Modified SW F_sw=zeros(3,IVL-1); %Forward Flux SW F=zeros(3,IVL-1); %Modified Flux %Jabcobians-----MSW A_b=zeros(3,3,IVL-1); A_b11=zeros(1,IVL-1); A_b12=ones(1,IVL-1); A_b13=zeros(1,IVL-1); A_b21=.5*(gamma-3).*U_b(2,:).^2./U_b(1,:).^2; A_b22=-(gamma-3).*U_b(2,:)./U_b(1,:); A_b23=Beta.*ones(1,IVL-1); A_b31=-gamma.*U_b(3,:).*U_b(2,:)./U_b(1,:).^2+Beta.*U_b(2,:).^3./U_b(1,:).^3; A_b32=gamma.*U_b(3,:)./U_b(1,:)-3*Beta.*U_b(2,:).^2./(2.*U_b(1,:).^2); A_b33=gamma.*U_b(2,:)./U_b(1,:); A_b(1,1,:)=A_b11; A_b(1,2,:)=A_b12; A_b(1,3,:)=A_b13; A_b(2,1,:)=A_b21; A_b(2,2,:)=A_b22; A_b(2,3,:)=A_b23; A_b(3,1,:)=A_b31; A_b(3,2,:)=A_b32; A_b(3,3,:)=A_b33; %Jacobians------SW A=zeros(3,3,IVL); A11=zeros(1,IVL); A12=ones(1,IVL); A13=zeros(1,IVL); A21=.5*(gamma-3).*U(2,:).^2./U(1,:).^2; A22=-(gamma-3).*U(2,:)./U(1,:); A23=Beta.*ones(1,IVL); A31=-gamma.*U(3,:).*U(2,:)./U(1,:).^2+Beta.*U(2,:).^3./U(1,:).^3; A32=gamma.*U(3,:)./U(1,:)-3*Beta.*U(2,:).^2./(2.*U(1,:).^2); A33=gamma.*U(2,:)./U(1,:); A(1,1,:)=A11; A(1,2,:)=A12; A(1,3,:)=A13; A(2,1,:)=A21; A(2,2,:)=A22; A(2,3,:)=A23; A(3,1,:)=A31; A(3,2,:)=A32; A(3,3,:)=A33; %Flux Split the Jacobians for vol_ind=2:num_vols %Get A+ A- from data from i and i-1 for MSW [Lam_v_b,Lam_e_b]=eig(A_b(:,:,vol_ind-1)); Lam_e_p_b=(Lam_e_b+abs(Lam_e_b))./2; Lam_e_m_b=(Lam_e_b-abs(Lam_e_b))./2; %Get A+ A- from data from i and i+1 for SW [Lam_v_i_1,Lam_e_i_1]=eig(A(:,:,vol_ind)); Lam_e_m=(Lam_e_i_1-abs(Lam_e_i_1))./2; [Lam_v_i,Lam_e_i]=eig(A(:,:,vol_ind-1)); Lam_e_p=(Lam_e_i+abs(Lam_e_i))./2; %Now get Generic Fluxes SW=Steger-warming and MSW=Modified SW F_sw(:,vol_ind-1)= Lam_v_i_1*Lam_e_m*Lam_v_i_1^-1*U(:,vol_ind)+... Lam_v_i*Lam_e_p*Lam_v_i^-1*U(:,vol_ind-1); F_msw(:,vol_ind-1)=Lam_v_b*Lam_e_p_b*Lam_v_b^-1*U(:,vol_ind-1)+... Lam_v_b*Lam_e_m_b*Lam_v_b^-1*U(:,vol_ind); %Use pressure gradient waiting to avoid going negative in pressure pg=(p(vol_ind)-p(vol_ind-1))/min(p(vol_ind),p(vol_ind-1)); wt=1/(1+pg^2); F(:,vol_ind-1)=wt*F_msw(:,vol_ind-1)+(1-wt)*F_sw(:,vol_ind-1); %Find the minimum delx/u+c (i.e. the greatest eigenvalue for the CFL) if vol_ind==2 lam_max=max(Lam_e_b(:)); uc_min=delx/lam_max; else if max(Lam_e_b(:)) > lam_max lam_max=max(Lam_e_b(:)); uc_min=delx/lam_max; end end end %----------------------------------ROE METHOD else %State and Flux vectors and Q vector U_b=(U(:,1:IVL-1)+U(:,2:IVL))./2; %Averages F_msw=zeros(3,IVL-1); %Backward Flux Modified SW F_sw=zeros(3,IVL-1); %Forward Flux SW F=zeros(3,IVL-1); %Modified Flux %Get the Roe Averaged Variables rho=U(1,:); u=U(2,:)./U(1,:); e=U(3,:); p=(gamma-1)*U(1,:).*(U(3,:)./U(1,:)-.5*(U(2,:)./U(1,:)).^2); %pressure rho_hat=sqrt(rho(1:IVL-1)).*sqrt(rho(2:IVL)); %RA density u_hat=(sqrt(rho(1:IVL-1)).*u(1:IVL-1)+sqrt(rho(2:IVL)).*u(2:IVL))./... %RA velocity (sqrt(rho(1:IVL-1))+sqrt(rho(2:IVL))); h_hat=(sqrt(rho(1:IVL-1)).*(e(1:IVL-1)+p(1:IVL-1))./rho(1:IVL-1)+... sqrt(rho(2:IVL)).*(e(2:IVL)+p(2:IVL))./rho(2:IVL))./... (sqrt(rho(1:IVL-1))+sqrt(rho(2:IVL))); c_hat=sqrt((gamma-1).*(h_hat-.5*u_hat.^2)); %Jabcobians-----MSW A_h=zeros(3,3,IVL-1); A_h11=zeros(1,IVL-1); A_h12=ones(1,IVL-1); A_h13=zeros(1,IVL-1); A_h21=.5*(gamma-3).*u_hat.^2; A_h22=-(gamma-3).*u_hat; A_h23=Beta.*ones(1,IVL-1); A_h31=(gamma/2-1).*u_hat.^3-c_hat.^2.*u_hat./(gamma-1); A_h32=(3/2-gamma).*u_hat.^2+c_hat.^2./(gamma-1); A_h33=gamma.*u_hat; A_h(1,1,:)=A_h11; A_h(1,2,:)=A_h12; A_h(1,3,:)=A_h13; A_h(2,1,:)=A_h21; A_h(2,2,:)=A_h22; A_h(2,3,:)=A_h23; A_h(3,1,:)=A_h31; A_h(3,2,:)=A_h32; A_h(3,3,:)=A_h33; %Jacobians------SW A=zeros(3,3,IVL); A11=zeros(1,IVL); A12=ones(1,IVL); A13=zeros(1,IVL); A21=.5*(gamma-3).*U(2,:).^2./U(1,:).^2; A22=-(gamma-3).*U(2,:)./U(1,:); A23=Beta.*ones(1,IVL); A31=-gamma.*U(3,:).*U(2,:)./U(1,:).^2+Beta.*U(2,:).^3./U(1,:).^3; A32=gamma.*U(3,:)./U(1,:)-3*Beta.*U(2,:).^2./(2.*U(1,:).^2); A33=gamma.*U(2,:)./U(1,:); A(1,1,:)=A11; A(1,2,:)=A12; A(1,3,:)=A13; A(2,1,:)=A21; A(2,2,:)=A22; A(2,3,:)=A23; A(3,1,:)=A31; A(3,2,:)=A32; A(3,3,:)=A33; %Split the Flux Difference vectors for vol_ind=2:num_vols %Get A_h+ A_h- from data from i and i-1 for ROE [Lam_v_h,Lam_e_h]=eig(A_h(:,:,vol_ind-1)); %Get eigenvalues of A to avoid shock-happiness Lam_e_i=eig(A(:,:,vol_ind)); Lam_e_i_1=eig(A(:,:,vol_ind-1)); %Check eigenvalues for e_ind=1:3 ei_half=sigma*max([0,... Lam_e_h(e_ind,e_ind)-Lam_e_i_1(e_ind),... Lam_e_i(e_ind)-Lam_e_h(e_ind,e_ind)]); if abs(Lam_e_h(e_ind,e_ind)) < ei_half Lam_e_h(e_ind,e_ind)=.5*((Lam_e_h(e_ind,e_ind))^2/ei_half+ei_half); end max_e_val=max(max([max(Lam_e_h)',Lam_e_i(:),Lam_e_i_1(:)])) end %Get the absolute value ROE averaged hat matrix Lam_e_p_h=(Lam_e_h+abs(Lam_e_h))./2; Lam_e_m_h=(Lam_e_h-abs(Lam_e_h))./2; abs_A_h=Lam_v_h*Lam_e_p_h*Lam_v_h^-1-... Lam_v_h*Lam_e_m_h*Lam_v_h^-1; %Now get the normal fluxes F_i=A(:,:,vol_ind)*U(:,vol_ind); F_i_1=A(:,:,vol_ind-1)*U(:,vol_ind-1); %And finally, get the ROE flux F(:,vol_ind-1)=(F_i+F_i_1)./2-.5*abs_A_h*(U(:,vol_ind)-U(:,vol_ind-1)); %Find the minimum delx/u+c (i.e. the greatest eigenvalue for the CFL) if vol_ind==2 lam_max=max_e_val; uc_min=delx/lam_max; else if max_e_val > lam_max lam_max=max_e_val; uc_min=delx/lam_max; end end end end %Calculate the time step using the delt=0.9*uc_min; Q(2,2:IVL-1)=(1./Vol(:,2:IVL-1)).*(Surf(:,3:IL-1)-Surf(:,2:IL-2)).*... p(2:IVL-1); U(:,2:IVL-1)=U(:,2:IVL-1)-... (delt./Vols(:,2:IVL-1)).*... (F(:,2:IVL-1).*Surfs(:,3:IL-1)-F(:,1:IVL-2).*Surfs(:,2:IL-2))+... delt.*Q(:,2:IVL-1); %Exit Boundary Conditions--------Works [U_IVL]=Quasi_1D_exit_SI_bc(delt,delx,Cv,a,To,po,gamma,... p(IVL),p(IVL-1),U(:,IVL),U(:,IVL-1),Surf(IVL),Surf(IVL-1)); U(:,IVL)=U_IVL; %Entrance Boundary Condition----- p(1) [delU]=Quasi_1D_entrance_SI_bc(delt,delx,Cv,a,To,po,gamma,... p(1),p(2),U(:,1),U(:,2),Surf(1),Surf(2)); U(:,1)=U(:,1)+delU; %Plot the conservative variables subplot(2,1,1); plot(xvols,U(2,:)./U(1,:)); %plot(xvols,U(2,:),xvols,U(1,:),xvols,U(3,:),'*') subplot(2,1,2); plot(xvols,(p/po)) set(gca,'ylim',[0 2]); pause(0.03) if t==100 pause; end end end