Week 7 - Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
Objectives of the project –
- Create MATLAB script to solve the 1D Quasi supersonic nozzle flow governing equations to obtain primitive parameters like pressure, velocity, temperature, Mach number and mass flowrate using the Macormack Method
- Derive the differential equations in both conservative and non-conservative approach and compare the accuracy of the results of both approaches.
- Obtain the steady state flow condition and observe and report the variation of primitive variables and residuals as the solution converges.
- Perform a grid dependence test to observe and report the effect of the grid on the accuracy of the result as well as the simulation run time.
Theory
The macormack method is a popular method to solve the governing equations of 1D supersonic nozzle flow. It is an explicit time marching finite difference method which employs two steps to calculate the space derivative terms – predictor and corrector. Both the space and time derivatives are first order approximated.
As per the method the value of a variable 'v' at the next time step at a particular non-boundry grid point is given as-
And the average term is calculated as..
Predictor step-
The space terms are discretized and forward differenced, and the value of time derivative is calculated. The time derivative is then stored as the predictor value.
The value in the next time step is calculated and updated as shown. These values will be passed on the next step the corrector step.
Corrector step-
In this step the variable ‘v’ is updated from the previous step and the time derivative is once again calculated , the space terms are backward differenced .
After the two steps the final time derivative is taken as the average of the two-time derivatives from both the steps.
The value of the variable at the grid point is then updated with the average value .
Governing equations of quasi 1D supersonic nozzle flow
The governing equations have been derived through the integral form of the continuity(mass conservation) , momentum conservation and energy conservation.
The derivation has two different approaches non-conservative (control volume moving) , conservative (control volume not moving).
Non-conservative :
These quantities are all in non-dimension form. This makes the comparisons and graph plotting easier because all these values vary between 0 and 1.
The three values that will be solver for through macromack method would be density (rho) , temperature(T) , velocity(u).
The other values are obtained through the relations.
Boundary Conditions-
The values of the primitive variables in the inlet and outlet are not used in the macromack procedure. Their values are updated after the bulk interior has reached a convergent solution.
Outlet conditions
Inlet conditions
Initial conditions
The starting value of the solution at t=0, are given according to a nozzle distance-based distribution.
The primitive variable will change their values according to the time step , however the area will remain constant .
Calculation of time step
Courant flow number is used for calculating a time step that allows the solution to be stable, for a linear hyperbolic equation , CFL <1 is the ideal stable solution
Conservative form –
The simplified governing equations for the conservative forms do not have single variables inside the time derivative , therefore a solution vector U , flux vector F ,and source term J are used to simplify the equations that can be fed into the macromack solver .
Before each time step the values of U1,U2,U3 are calculated from the primitive variables , and after the macromack solver , the updated primitive values are obtained from the solution vectors .
The relation of the flux vectors are shown below –
Boundary conditions-
Inlet-
U2 , U3 change with time however U1 remains constant .
Outlet-
Initial conditions-
, initial velocity
All these are at t=0.
Results for Non-Conservative .
Primitive variables of 1D supersonic nozzle flow
These plots were obtained, they show the variation of qualitative aspect of nozzle flow along the nozzle length. These graphs agree with the Anderson Textbook.
The pressure decreases along the nozzle and the density along with the temperature also decrease, the Mach number increases along the nozzle. This is because it is a convergent -divergent nozzle.
The following graphs depict the variation of the qualitative aspect’s pf the flow at the nozzle throat as the time simulation progresses. A comparison with the exact solution is presented.
As it is clear from the graphs, the solution approaches steady state over about 700 iterations
The above plot depicts the mass flow rate distribution in the nozzle, at different time intervals , as it can be seen by the 700th iteration , the mass flow rate converges and becomes constant .
Results of conservative:
Primitive varibles for the 1D supersonic flow
The results are similar to the non-conservative . Mach number increases from inlet to outlet since the nozzle is convergent - divergent.
The following four plots depict the values at the throat as the simulation progresses , exact solution added in for comparison.
These graphs show that unlike the non-conservative , the conservative solution approaches steady state faster.
The following graphs depicts the residuals at different time intervals and shows how they converge . The residuals in this case is the average time derivative of each of the three solution vectors.
The following graph depicts the mass flow rate at different time steps as the solution progresses. By the 700th iteration , the mass flow rate approaches a steady value. Mass conservation is more accurately solved in the conservative form.
Grid dependence test
Simulations were run on the non-conservative and conservative forms for the primitive variables for two different number of grids (n)= 31 and 61.
The values were compared to the exact value.
| Density | Temperature | Pressure | Mach number | Mass flow rate | |
| Non-con , 31 | 0.6392 | 0.8367 | 0.5289 | 0.9991 | 0.5842 |
| Con ,31 | 0.6569 | 0.8406 | 0.5522 | 1.1274 | 0.59 |
| Non-con,61 | 0.6376 | 0.8354 | 0.5284 | 0.9997 | 0.5826 |
| Con,61 | 0.6388 | 0.8348 | 0.5333 | 1.1454 | 0.5841 |
| Exact | 0.634 | 0.833 | 0.528 | 1 | 0.579 |
Conclusion
1.The increase in the number of grid points gives a more accurate result in both forms. However, the run time also increases, therefore a careful compromise needs to be arrived at.
2.Non-conservative has a more accurate result for a given number of grid points.
3.Non conservative is easier to code and also runtime of the simulation is lesser.
- The conservative form approaches steady state faster (prefer above graphs). Lesser number of iterations are required.
MATLAB Script-
Code description-
'Project_2_Timemarching'- Main program for fixed time step Non conservative.
'Project_2_non_cons_steady_state'-Main program for steady state Non conservative.
'Project_2_Conser_timemarching'-Main program for fixed time step Conservative.
'Project_2_conser_Steady_state'-Main program for steady state Conservative.
'plot_mass_flow_dtvalues'- Function for plotting variation of mass flow rates for both forms
'plot_throat_time_conser'-Function for plotting primitive variables in the throat for conservative.
'plot_throat_time_non_conser' -Function for plotting primitive variables in the throat for non-conservative .
'plot_residuals'- Function for plotting the three residuals.
'plot_primitive_non_conser' - Function used to plot the primitive variables along the nozzle for both forms , initially only for non conservative , but used for conservative as well.
1.Project_2_Timemarching
clear all
clc
close all
%inputs
n=31;
x = linspace(0,3,n);
gamma=1.4;
dx=x(2)-x(1);
%initial conditions
rho = 1-0.3146*x;
k=1;
t = 1- 0.2314*x;
v=(0.1 + 1.09*x).*t.^2;
a= 1+2.2*(x-1.5).^2;
nt=1400;
th=(n+1)/2;
for i =1:n
intidt(i)= 0.5*dx/((t(i)^0.5)+v(i));
end
actualdt=min(intidt);
dt=actualdt;
index_for_massflow=1;
for k=1:nt
rhoold=rho
vold=v;
told=t;
%predictor method
for j=2:n-1
dvdx = (v(j+1)-v(j))/dx;
drhodx=(rho(j+1)-rho(j))/dx;
dtdx=(t(j+1)-t(j))/dx;
dlogadx=(log(a(j+1))-log(a(j)))/dx;
%continiuity
drhodt_p(j) = -rho(j)*dvdx - rho(j)*v(j)*dlogadx - v(j)*drhodx;
%momentum
dvdt_p(j) = -v(j)*dvdx - (1/gamma)*(dtdx + (t(j)/rho(j))*drhodx);
%energy equation
dtdt_p(j) = -v(j)*dtdx -(gamma-1)*t(j)*(dvdx + v(j)*dlogadx);
%solution update
v(j) =v(j) + dvdt_p(j)*dt;
rho(j) = rho(j) + drhodt_p(j)*dt;
t(j)= t(j) + dtdt_p(j)*dt;
end
for j=2:n-1
dvdx = (v(j)-v(j-1))/dx;
drhodx=(rho(j)-rho(j-1))/dx;
dtdx=(t(j)-t(j-1))/dx;
dlogadx=(log(a(j))-log(a(j-1)))/dx;
%continiuity
drhodt_c(j) = -rho(j)*dvdx - rho(j)*v(j)*dlogadx - v(j)*drhodx;
%momentum
dvdt_c(j) = -v(j)*dvdx - (1/gamma)*(dtdx + (t(j)/rho(j))*drhodx);
%energy equation
dtdt_c(j) = -v(j)*dtdx -(gamma-1)*t(j)*(dvdx + v(j)*dlogadx);
end
dvdt=0.5*(dvdt_p + dvdt_c);
drhodt=0.5*(drhodt_p+drhodt_c);
dtdt=0.5*(dtdt_c + dtdt_p);
for i=2:n-1
v(i)=vold(i) + dvdt(i)*dt;
rho(i) = rhoold(i) + drhodt(i)*dt;
t(i)= told(i) + dtdt(i)*dt;
M(i) = v(i)/(t(i)^0.5);
p(i) = ( 1 + ((gamma-1)/2*M(i)*M(i)))^(-gamma/(gamma-1));
m(i)=rho(i)*a(i)*v(i);
end
v(1)=2*v(2)-v(3);
v(n)=2*v(n-1)-v(n-2);
rho(n)=2*rho(n-1)-rho(n-2);
t(n)=2*t(n-1)-t(n-2);
m(n)=rho(n)*a(n)*v(n);
M(n) = v(n)/(t(n)^0.5);
p(n) = ( 1 + ((gamma-1)/2*M(n)*M(n)))^(-gamma/(gamma-1));
rho_throat(k)=rho(16);
time_iter(k)=k;
end
2.Project_2_non_cons_steady_state
clear all
clc
close all
%inputs
n=31;
x = linspace(0,3,n);
gamma=1.4;
dx=x(2)-x(1);
%initial conditions
rho = 1-0.3146*x;
k=1;
t = 1- 0.2314*x;
v=(0.1 + 1.09*x).*t.^2;
a= 1+2.2*(x-1.5).^2;
nt=1400;
th=(n+1)/2;
for i =1:n
intidt(i)= 0.5*dx/((t(i)^0.5)+v(i));
end
actualdt=min(intidt);
dt=actualdt;
index_for_massflow=1;
for k=1:nt
rhoold=rho
vold=v;
told=t;
%predictor method
for j=2:n-1
dvdx = (v(j+1)-v(j))/dx;
drhodx=(rho(j+1)-rho(j))/dx;
dtdx=(t(j+1)-t(j))/dx;
dlogadx=(log(a(j+1))-log(a(j)))/dx;
%continiuity
drhodt_p(j) = -rho(j)*dvdx - rho(j)*v(j)*dlogadx - v(j)*drhodx;
%momentum
dvdt_p(j) = -v(j)*dvdx - (1/gamma)*(dtdx + (t(j)/rho(j))*drhodx);
%energy equation
dtdt_p(j) = -v(j)*dtdx -(gamma-1)*t(j)*(dvdx + v(j)*dlogadx);
%solution update
v(j) =v(j) + dvdt_p(j)*dt;
rho(j) = rho(j) + drhodt_p(j)*dt;
t(j)= t(j) + dtdt_p(j)*dt;
end
for j=2:n-1
dvdx = (v(j)-v(j-1))/dx;
drhodx=(rho(j)-rho(j-1))/dx;
dtdx=(t(j)-t(j-1))/dx;
dlogadx=(log(a(j))-log(a(j-1)))/dx;
%continiuity
drhodt_c(j) = -rho(j)*dvdx - rho(j)*v(j)*dlogadx - v(j)*drhodx;
%momentum
dvdt_c(j) = -v(j)*dvdx - (1/gamma)*(dtdx + (t(j)/rho(j))*drhodx);
%energy equation
dtdt_c(j) = -v(j)*dtdx -(gamma-1)*t(j)*(dvdx + v(j)*dlogadx);
end
dvdt=0.5*(dvdt_p + dvdt_c);
drhodt=0.5*(drhodt_p+drhodt_c);
dtdt=0.5*(dtdt_c + dtdt_p);
for i=2:n-1
v(i)=vold(i) + dvdt(i)*dt;
rho(i) = rhoold(i) + drhodt(i)*dt;
t(i)= told(i) + dtdt(i)*dt;
M(i) = v(i)/(t(i)^0.5);
p(i) = ( 1 + ((gamma-1)/2*M(i)*M(i)))^(-gamma/(gamma-1));
m(i)=rho(i)*a(i)*v(i);
end
v(1)=2*v(2)-v(3);
v(n)=2*v(n-1)-v(n-2);
rho(n)=2*rho(n-1)-rho(n-2);
t(n)=2*t(n-1)-t(n-2);
m(n)=rho(n)*a(n)*v(n);
M(n) = v(n)/(t(n)^0.5);
p(n) = ( 1 + ((gamma-1)/2*M(n)*M(n)))^(-gamma/(gamma-1));
rho_throat(k)=rho(16);
time_iter(k)=k;
if(k==1)
mass_flow_throat_dt(index_for_massflow)=m(
end
plot_primitive_non_conser(rho,t,Mach,v,x,pressure);
plot_throat_time_non_conser(rho_throat,temp_throat,pressure_throat,Mach_throat,time_iter);
3.Project_2_Conser_timemarching
clear all
close all
clc
L=3;
n=31;
gamma=1.4;
nt=1400;
dt=0.0206;
x=linspace(0,L,n);
dx=x(2)-x(1);
for i=1:n
a(i)=1 + 2.2*(x(i)-1.5)^(2);
end
for i=1:length(x)
if((x(i)>=0)&&(x(i)<=0.5))
rho(i)=1.0;
t(i)=1;
else if((x(i)>=0.5)&&(x(i)<=1.5))
rho(i)= (1 - 0.366*(x(i)-0.5));
t(i) =(1 - 0.167*(x(i)-0.5));
else
rho(i)=(0.634 - 0.3879*(x(i)-1.5));
t(i)=(0.833 - 0.3507*(x(i)-1.5));
end
end
end
v= 0.59./(rho.*a);
u1=rho.*a;
num=1;
u2=rho.*a.*v;
u3=rho.*a.*((t/(gamma-1)) + (0.5*gamma.*v.*v));
rhoi=rho
ai=a;
ti=t;
vi=v;
th=(n+1)/2;
for j=1:nt
u1old=u1;
u2old=u2;
u3old=u3;
rhoold=rho;
told=t;
vold=v;
%flux vairables
f1=u2;
f2=((u2.^2)./u1)+((gamma-1)/gamma).*(u3 - ((gamma/2)*(u2.^2./u1)));
f3=((gamma.*u2.*u3)./u1) - ((gamma*(gamma-1))/2)*(u2.^3)./(u1.^2);
%predictormethod
for k=2:n-1
dadax=(a(k+1)-a(k))/dx;
dlogadx=(log(a(k+1))-log(a(k)))/dx;
j2=dlogadx*((gamma-1)/gamma)*(u3(k)-(0.5*gamma*u2(k)^2)/(u1(k)));
%predictor step
du1dt_p(k) = -((f1(k+1)-f1(k))/dx);
du2dt_p(k)=-((f2(k+1)-f2(k))/dx) + j2;
du3dt_p(k)=-((f3(k+1)-f3(k))/dx);
%solution update
u1(k) = (du1dt_p(k)*dt)+u1old(k);
u2(k) = (du2dt_p(k)*dt)+u2old(k);
u3(k) = (du3dt_p(k)*dt)+u3old(k);
end
%boundry conditions
u1(1)=a(1);
u2(1)=2*u2(2)-u2(3);
u3(1)=u1(1)*((t(1)/(gamma-1)+((gamma/2)*v(1)^2)));
u1(n)=(2*u1(n-1))-u1(n-2);
u2(n)=(2*u2(n-1))-u2(n-2);
u3(n)=(2*u3(n-1))-u3(n-2);
%calculating primitive
rho = u1./a;
t=(gamma-1).*((u3./u1)-((gamma/2).*(u2./u1).^2));
v=u2./u1;
%updating flux variable
f1=u2;
f2=((u2.^2)./u1)+((gamma-1)/gamma).*(u3 - ((gamma/2)*(u2.^2./u1)));
f3=((gamma.*u2.*u3)./u1) - ((gamma*(gamma-1))/2)*(u2.^3)./(u1.^2);
%corrector step
for k=2:n-1
dadx=(a(k)-a(k-1))/dx;
dlogadx=(log(a(k))-log(a(k-1)))/dx;
j2=dlogadx*((gamma-1)/gamma)*(u3(k)-(0.5*gamma*u2(k)^2)/(u1(k)));
%corrector
du1dt_c(k)=-(f1(k)-f1(k-1))/dx;
du2dt_c(k)=-(f2(k)-f2(k-1))/dx + j2;
du3dt_c(k)=-(f3(k)-f3(k-1))/dx;
end
for p=2:n-1
du1dt_avg(p)=0.5*(du1dt_c(p) + du1dt_p(p));
du2dt_avg(p)=0.5*(du2dt_c(p) + du2dt_p(p));
du3dt_avg(p)=0.5*(du3dt_c(p) + du3dt_p(p));
%final for u1 u2 u3
u1(p)=(du1dt_avg(p)*dt)+u1old(p);
u2(p)=(du2dt_avg(p)*dt)+u2old(p);
u3(p)=(du3dt_avg(p)*dt)+u3old(p);
end
u1(1)=a(1);
u2(1)=2*u2(2)-u2(3);
u3(1)=u1(1)*((t(1)/(gamma-1)+((gamma/2)*v(1)^2)));
u1(n)=(2*u1(n-1))-u1(n-2);
u2(n)=(2*u2(n-1))-u2(n-2);
u3(n)=(2*u3(n-1))-u3(n-2);
%recalculating
rho=u1./a;
v=u2./u1;
t=(gamma-1).*((u3./u1)-((gamma/2).*(u2./u1).^2));
for i =1:n
mach(i)=(v(i)^0.5)/t(i);
p(i)=rho(i)*t(i);
massflow(i)=rho(i)*v(i)*a(i)
end
temp_throat(j)=t(th);
mach_throat(j)=mach(th);
rho_throat(j)=rho(th);
pressure_throat(j)=p(th);
u1_res(j)=du1dt_avg(th);
u2_res(j)=du2dt_avg(th);
u3_res(j)=du3dt_avg(th);
if(j==10)
massflow_10dt=massflow;
end
if(j==50)
massflow_50dt=massflow;
end
if(j==100)
massflow_100dt=massflow;
end
if(j==150)
massflow_150dt=massflow;
end
if(j==200)
massflow_200dt=massflow;
end
if(j==700)
massflow_700dt=massflow;
end
if(j==1000)
massflow_1000dt=massflow;
end
end
time_array=(1:nt);
plot_throat_time_conser(rho_throat,temp_throat,pressure_throat,mach_throat,time_array);
plot_residuals(u1_res,u2_res,u3_res,time_array)
plot_mass_flow_dtvalues(massflow_10dt,massflow_50dt,massflow_100dt,massflow_150dt,massflow_200dt,massflow_700dt,massflow_1000dt,x)
rho_31=rho(th);
t_31=t(th);
p_31=p(th);
mach_31=mach(th);
m_31=massflow(th);
4.Project_2_conser_steady_state
clear all
close all
clc
L=3;
n=31;
gamma=1.4;
nt=1400;
dt=0.0206;
x=linspace(0,L,n);
dx=x(2)-x(1);
th=(n+1)/2;
for i=1:n
a(i)=1 + 2.2*(x(i)-1.5)^(2);
end
for i=1:length(x)
if((x(i)>=0)&&(x(i)<=0.5))
rho(i)=1.0;
t(i)=1;
else if((x(i)>=0.5)&&(x(i)<=1.5))
rho(i)= (1 - 0.366*(x(i)-0.5));
t(i) =(1 - 0.167*(x(i)-0.5));
else
rho(i)=(0.634 - 0.3879*(x(i)-1.5));
t(i)=(0.833 - 0.3507*(x(i)-1.5));
end
end
end
v= 0.59./(rho.*a);
u1=rho.*a;
num=1;
u2=rho.*a.*v;
u3=rho.*a.*((t/(gamma-1)) + (0.5*gamma.*v.*v));
rhoi=rho
ai=a;
ti=t;
vi=v;
errort=9e9;
errorv=9e9;
errorrho=9e9;
error=9e9;
temp=9e9;
num=1;
tol=1e-5;
while(error>tol)
num=num+1;
u1old=u1;
u2old=u2;
u3old=u3;
rhoold=rho;
told=t;
vold=v;
%flux vairables
f1=u2;
f2=((u2.^2)./u1)+((gamma-1)/gamma).*(u3 - ((gamma/2)*(u2.^2./u1)));
f3=((gamma.*u2.*u3)./u1) - ((gamma*(gamma-1))/2)*(u2.^3)./(u1.^2);
%predictormethod
for k=2:n-1
dadax=(a(k+1)-a(k))/dx;
dlogadx=(log(a(k+1))-log(a(k)))/dx;
j2=dlogadx*((gamma-1)/gamma)*(u3(k)-(0.5*gamma*u2(k)^2)/(u1(k)));
%predictor step
du1dt_p(k) = -((f1(k+1)-f1(k))/dx);
du2dt_p(k)=-((f2(k+1)-f2(k))/dx) + j2;
du3dt_p(k)=-((f3(k+1)-f3(k))/dx);
%solution update
u1(k) = (du1dt_p(k)*dt)+u1old(k);
u2(k) = (du2dt_p(k)*dt)+u2old(k);
u3(k) = (du3dt_p(k)*dt)+u3old(k);
end
%boundry conditions
u1(1)=a(1);
u2(1)=2*u2(2)-u2(3);
u3(1)=u1(1)*((t(1)/(gamma-1)+((gamma/2)*v(1)^2)));
u1(n)=(2*u1(n-1))-u1(n-2);
u2(n)=(2*u2(n-1))-u2(n-2);
u3(n)=(2*u3(n-1))-u3(n-2);
%calculating primitive
rho = u1./a;
t=(gamma-1).*((u3./u1)-((gamma/2).*(u2./u1).^2));
v=u2./u1;
%updating flux variable
f1=u2;
f2=((u2.^2)./u1)+((gamma-1)/gamma).*(u3 - ((gamma/2)*(u2.^2./u1)));
f3=((gamma.*u2.*u3)./u1) - ((gamma*(gamma-1))/2)*(u2.^3)./(u1.^2);
%corrector step
for k=2:n-1
dadx=(a(k)-a(k-1))/dx;
dlogadx=(log(a(k))-log(a(k-1)))/dx;
j2=dlogadx*((gamma-1)/gamma)*(u3(k)-(0.5*gamma*u2(k)^2)/(u1(k)));
%corrector
du1dt_c(k)=-(f1(k)-f1(k-1))/dx;
du2dt_c(k)=-(f2(k)-f2(k-1))/dx + j2;
du3dt_c(k)=-(f3(k)-f3(k-1))/dx;
end
for p=2:n-1
du1dt_avg(p)=0.5*(du1dt_c(p) + du1dt_p(p));
du2dt_avg(p)=0.5*(du2dt_c(p) + du2dt_p(p));
du3dt_avg(p)=0.5*(du3dt_c(p) + du3dt_p(p));
%final for u1 u2 u3
u1(p)=(du1dt_avg(p)*dt)+u1old(p);
u2(p)=(du2dt_avg(p)*dt)+u2old(p);
u3(p)=(du3dt_avg(p)*dt)+u3old(p);
end
u1(1)=a(1);
u2(1)=2*u2(2)-u2(3);
u3(1)=u1(1)*((t(1)/(gamma-1)+((gamma/2)*v(1)^2)));
u1(n)=(2*u1(n-1))-u1(n-2);
u2(n)=(2*u2(n-1))-u2(n-2);
u3(n)=(2*u3(n-1))-u3(n-2);
%recalculating
rho=u1./a;
v=u2./u1;
t=(gamma-1).*((u3./u1)-((gamma/2).*(u2./u1).^2));
errort=max(abs(told-t));
errorv=max(abs(vold-v));
errorrho=max(abs(rhoold-rho));
temp=max(errort,errorv);
error=max(temp,errorrho);
end
5.plot_mass_flow_dtvalues
function out= plot_mass_flow_dtvalues(m10,m50,m100,m150,m200,m700,m1000,x)
figure(1)
plot(x,m10)
hold on
plot(x,m50)
hold on
plot(x,m100)
hold on
plot(x,m150)
hold on
plot(x,m200)
hold on
plot(x,m700)
hold on
plot(x,m1000)
hold on
title('Mass flow rates at different time intervals');
xlabel('Nozzle length');
ylabel('Mass flow');
set(gcf,'position',[10,10,2000,400])
legend('10dt','50dt','100dt','150dt','200dt','700dt','1000dt');
end
6.plot_throat_time_conser
function out=plot_throat_time_conser(rho,temp,pressure,mach,time)
rho_ex=0.634*ones(1,length(time));
temp_ex=0.8338*ones(1,length(time));
pressure_ex=0.528*ones(1,length(time));
mach_ex=ones(1,length(time));
figure(1)
plot(time,rho,'r');
hold on
plot(time,rho_ex,'b');
title('Density at the throat ');
axis([0 1400 0 1]);
xlabel('Time steps');
ylabel('Density Ratio');
set(gcf,'position',[10,10,2000,400])
legend('Numerical solution','Analytical silution');
figure(2)
plot(time,temp,'r');
hold on
plot(time,temp_ex,'b');
axis([0 1400 0.7 1]);
set(gcf,'position',[10,10,2000,400])
title('Temperature at the throat');
xlabel('Time steps');
ylabel('Temperature ratio');
legend('Numerical solution','Analytical solution');
figure(3)
plot(time,pressure,'r');
hold on
plot(time,pressure_ex,'b');
axis([0 1400 0.2 0.7]);
title('Pressure ratio at the throat');
set(gcf,'position',[10,10,2000,400])
xlabel('Time steps');
ylabel('Pressure ratio');
legend('Numerical solution','Analytical solution');
figure(4)
plot(time,mach,'r');
hold on
plot(time,mach_ex,'b');
axis([0 1400 0 2]);
title('Mach number at the throat');
set(gcf,'position',[10,10,2000,400]);
xlabel('Time steps');
ylabel('Mach number');
legend('Numerical solution','Analytical solution');
end
7.plot_throat_time_non_conser
function out = plot_throat_time_non_conser(rho,temp,pressure,mach,time)
rho_ex=0.634*ones(1,length(time));
temp_ex=0.8338*ones(1,length(time));
pressure_ex=0.528*ones(1,length(time));
mach_ex=ones(1,length(time));
figure(1)
xlabel('Time steps');
ylabel('Density ratio');
plot(time,rho,'r')
hold on
plot(time,rho_ex,'b')
grid on
set(gcf,'position',[10,10,2000,400])
title('Variation of Density ratio at throat');
legend('Numerical Solution','Exact solution');
figure(2)
xlabel('Time steps');
ylabel('Temperature ratio');
plot(time,temp,'r')
hold on
plot(time,temp_ex,'b')
grid on
set(gcf,'position',[10,10,2000,400])
title('Variation of Temperature ratio at throat');
legend('Numerical Solution','Exact solution')
figure(3)
xlabel('Time steps');
ylabel('Pressure ratio');
plot(time,pressure,'r')
hold on
plot(time,pressure_ex,'b')
grid on
set(gcf,'position',[10,10,2000,400])
title('Variation of Pressure ratio at throat');
legend('Numerical Solution','Exact solution')
figure(4)
xlabel('Time steps');
ylabel('Mach number');
plot(time,mach,'r')
hold on
plot(time,mach_ex,'b')
grid on
set(gcf,'position',[10,10,2000,400])
title('Variation of Mach number at throat');
legend('Numerical Solution','Exact solution')
end
8.plot_residuals
function out = plot_residuals(u1,u2,u3,time)
figure(1)
plot(time,u1,'r');
hold on
plot(time,u2,'b');
hold on
plot(time,u3,'g');
title('Residuals at the throat section');
xlabel('Number of iterations');
ylabel('Residual');
set(gcf,'position',[10,10,2000,400])
legend('U1-residual','U2-residual','U3-residual');
end
9.plot_primitive_non_conser
function out=plot_primitive_non_conser(rho,temp,mach,velocity,x,p)
figure(1)
plot(x,rho)
grid on
set(gcf,'position',[10,10,2000,400])
title('Density Ratio over the grid')
xlabel('Nozzle length');
ylabel('Density ratio');
figure(2)
plot(x,temp)
grid on
set(gcf,'position',[10,10,2000,400])
title('Temperature Ratio over the grid')
xlabel('Nozzle length');
ylabel('Temperature ratio');
figure(3)
plot(x,mach)
grid on
set(gcf,'position',[10,10,2000,400])
title('Mach number over the grid')
xlabel('Nozzle length');
ylabel('Mach number');
figure(4)
plot(x,p)
grid on
set(gcf,'position',[10,10,2000,400])
title('Pressure Ratio over the grid')
xlabel('Nozzle length');
ylabel('Pressure ratio');
end