%This script will attempt to solve the Quasi 1-d Euler equations % Cv=.17; %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 finite_vol=1; %1=finite volume 2=finite difference flux_method=1; %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 press=po*ones(1,IVL); %Initial Pressure rho=press./((gamma-1)*Cv*To); %Initial density u=zeros(1,IVL); %Initial Velocity a=sqrt(2*gamma*(gamma-1)/(gamma-1)*Cv*To); else po=2117; %Entrance Pressure in lb/ft^2 Mo=0; %Mach Number at Entrance=0 To=531.2; %Temperature at entrance press=po*ones(1,IVL); %Initial Pressure press(IVL)=1400; %Instantaneously drop pressure rho=press./((gamma-1)*Cv*To); %Initial density u=zeros(1,IVL); %Initial Velocity a=sqrt(2*gamma*(gamma-1)/(gamma-1)*Cv*To); end %Set up conserved quantities U(1,:)=rho; U(2,:)=rho.*u; U(3,:)=rho.*(Cv.*To+.5.*(u.^2)); plot(press) pause; num_steps=300; for t=0:num_steps %Boundary Conditions 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_p=zeros(3,IVL-2); %Forward Flux F_m=zeros(3,IVL-2); %Backward Flux %Jabcobiabn 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; end %Flux Split the Jacobians for vol_ind=2:num_vols-1 %Get A+ A- from data from i and i-1 [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 [Lam_v_f,Lam_e_f]=eig(A_b(:,:,vol_ind)); Lam_e_p_f=(Lam_e_f+abs(Lam_e_f))./2; Lam_e_m_f=(Lam_e_f-abs(Lam_e_f))./2; %Now get Generic Fluxes F_p=F_i+1/2 F_m=F_i-1/2 F_p(:,vol_ind-1)=Lam_v_f*Lam_e_p_f*Lam_v_f^-1*U(:,vol_ind)+... Lam_v_f*Lam_e_m_f*Lam_v_f^-1*U(:,vol_ind+1); F_m(:,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); %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_f(:)) > lam_max lam_max=max(Lam_e_f(:)); uc_min=delx/lam_max; end end end %Calculate the time step using the delt=0.9*uc_min; p=(gamma-1)*U(1,:).*(U(3,:)./U(1,:)-.5*U(2,:).^2); %pressure 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_p.*Surfs(:,3:IL-1)-F_m.*Surfs(:,2:IL-2))+... delt.*Q(:,2:IVL-1); p(IVL) pause; %Exit Boundary Conditions--------Works [U_IVL]=Quasi_1D_exit_SI_bc(delt,delx,Cv,a,To,po,gamma,... U(2,IVL),U(2,IVL-1),U(1,IVL),U(1,IVL-1),Surf(IVL),Surf(IVL-1)); U(:,IVL)=U_IVL; %Entrance Boundary Condition----- [delU]=Quasi_1D_entrance_SI_bc(delt,delx,Cv,a,To,po,gamma,... U(2,1),U(2,2),U(1,1),U(1,2),Surf(1),Surf(2)); U(:,1)=U(:,1)+delU; %Plot the conservative variables %plot(xvols,U(2,:),xvols,U(1,:),xvols,U(3,:),'*') plot(p) end