Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
This problem is solved for two forms of governing equations that are Conservative and Non conservative form for two types of grids that is 31 & 61:
Theory:
Maccormack Method is used to solve the Non-Conservative and Conservative form of governing equations which is as shown below:
1.Non-Conservative:
Code for Non-Conservative form is given as:
%1D Super-Sonic Nozzle Flow in Non-Conservative Form
close all
clc
%inputs
L=input(\'Enter the length of the domain: \');
n=input(\'Enter the number of grid points: \');
gamma=input(\'Enter the value of gamma: \');
%Specifying grid
x=linspace(0,L,n);
dx=x(2)-x(1);
%initial Conditions
rho=1-(0.3146*x);
T=1-(0.2314*x);
v=(0.1+(1.09*x)).*sqrt(T);
a=1+2.2*((x-1.5).^2);
n_throat= (n+1)/2;
c=0.5;
%time step
nt=1400;
dt=min(c.*(dx./((T.^0.5)+v)));
tic
%initial Condition
for k=1:nt
iter=0;
rho_old=rho;
T_old=T;
v_old=v;
massflow_old=rho_old.*v_old.*a;
%predictor step
for i=2:n-1
drho_dx=(rho(i+1)-rho(i))/dx;
dv_dx=(v(i+1)-v(i))/dx;
dT_dx=(T(i+1)-T(i))/dx;
dlog_a_dx=(log(a(i+1))-log(a(i)))/dx;
%continuity equation
drho_dt_p(i)= -rho(i)*dv_dx-rho(i)*v(i)*dlog_a_dx-v(i)*drho_dx;
%energy equation
dT_dt_p(i)=-v(i)*dT_dx-(gamma-1)*T(i)*(dv_dx+v(i)*dlog_a_dx);
%momentum equation
dV_dt_p(i)=-v(i)*dv_dx-(dT_dx+(T(i)/rho(i))*drho_dx)/gamma;
%solution update
rho(i)=rho_old(i)+ (drho_dt_p(i)*dt);
T(i)=T_old(i)+ dT_dt_p(i)*dt;
v(i)=v_old(i)+dV_dt_p(i)*dt;
end
%Corrector step
for i=2:n-1
drho_dx=(rho(i)-rho(i-1))/dx;
dv_dx=(v(i)-v(i-1))/dx;
dT_dx=(T(i)-T(i-1))/dx;
dlog_a_dx=(log(a(i))-log(a(i-1)))/dx;
%continuity equation
drho_dt_c(i)= -rho(i)*dv_dx-rho(i)*v(i)*dlog_a_dx-v(i)*drho_dx;
%energy equation
dT_dt_c(i)=-v(i)*dT_dx-(gamma-1)*T(i)*(dv_dx+v(i)*dlog_a_dx);
%momentum equation
dV_dt_c(i)=-v(i)*dv_dx-(dT_dx+(T(i)/rho(i))*drho_dx)/gamma;
end
drho_dt_avg=0.5*(drho_dt_c+ drho_dt_p);
dT_dt_avg=0.5*(dT_dt_c+dT_dt_p);
dV_dt_avg=0.5*(dV_dt_c+dV_dt_p);
%final Solution updates
for i= 2:n-1
rho(i)=rho_old(i)+drho_dt_avg(i)*dt;
T(i)=T_old(i)+dT_dt_avg(i)*dt;
v(i)=v_old(i)+dV_dt_avg(i)*dt;
end
%applying Boundary conditions;
%inlet Condition
v(1)=2*v(2)-v(3);
%outlet Condition
rho(n)=2*rho(n-1)-rho(n-2);
T(n)=2*T(n-1)-T(n-2);
v(n)=2*v(n-1)-v(n-2);
%mach number
M=v./sqrt(T);
%mass flow rate
massflow=rho.*v.*a;
p=rho.*T;
%throat values
rho_throat=rho(n_throat);
T_throat=T(n_throat);
p_throat=p(n_throat);
M_throat=M(n_throat);
v_throat=v(n_throat);
error=max(abs(massflow-massflow_old));
figure(1);
if k==1
plot(x,massflow,\'--k\')
hold on
elseif k==50
plot(x,massflow,\'c\')
hold on
elseif k==100
plot(x,massflow,\'g\')
hold on
elseif k==200
plot(x,massflow,\'b\')
hold on
elseif k==800
plot(x,massflow,\'r\')
hold on
elseif k==1400
plot(x,massflow,\'y\')
hold on
title(\'Mass flow Rates(Non Conservative)\')
xlabel(\'X/L\')
ylabel(\'\\rhoVA/\\rho_{o}v{o}a{o}\')
legend(\'50\\Deltat\',\'100\\Deltat\',\'150\\Deltat\',\'200\\Deltat\',\'800\\Deltat\',\'1400\\Deltat\')
saveas(gcf,\'massflow rate non conservative\',\'jpeg\')
iter=iter+1;
end
end
figure(2)
plot(x,rho)
title(\'Plot of Density Ratio vs Distance\')
xlabel(\'X/L\')
ylabel(\'\\rho/\\rho_{o}\')
saveas(gcf,\'Plot of Density vs Distance non conservative\',\'jpeg\')
hold on
figure(3)
plot(x,v,\'g\')
title(\'Plot of Velocity Ratio vs Distance\')
xlabel(\'X/L\')
ylabel(\'v/v{o}\')
saveas(gcf,\'Plot of Velocity vs Distance non conservative\',\'jpeg\')
hold on
figure(4)
plot(x,T,\'r\')
title(\'Plot of Temperature Ratio vs Distance\')
xlabel(\'X/L\')
ylabel(\'T/T{o}\')
saveas(gcf,\'Plot of Temperature vs Distance non conservative\',\'jpeg\')
hold on
figure(5)
plot(x,real(M),\'y\')
title(\'Plot of Mach number vs Distance\')
xlabel(\'X/L\')
ylabel(\'M\')
saveas(gcf,\'Plot of Mach number vs Distance non conservative\',\'jpeg\')
hold on
figure(6)
plot(x,p,\'b\')
title(\'Plot of Pressure Ratio vs Distance\')
xlabel(\'X/L\')
ylabel(\'P\')
saveas(gcf,\'Plot of Pressure Ratio vs Distance non conservative\',\'jpeg\')
simulation_time=toc;
sprintf(\'Computational Results from Non-Conservative Method are: \\n Computation Time is = %d s\\n Density Ratio From Results = %d \\n Temperature Ratio= %d \\n Velocity Ratio=%d \\n Mach No=%d \\n Pressure Ratio=%d\',simulation_time,rho_throat,T_throat,v_throat,M_throat,p_throat)
For Grid Point:31
Inputs:
Plots
A.Plot of Mass Flow Rate vs Distance of the Nozzle:
B.Plot of Density Ratio vs Distance of the Nozzle:
C.Plot of Velocity Ratio vs Distance of the Nozzle:
D.Plot of Temperature Ratio vs Distance of the Nozzle:
E.Plot of Pressure Ratio vs Distance of the Nozzle:
F.Plot of MACH NO vs Distance of the Nozzle:
Outputs:
For Grid = 61
2.Conservative:
Code for Conservative form is given as:
For Grid Point:31
%1D Super-Sonic Nozzle Flow in Conservative Form
close all
clc
%inputs
L=input(\'Enter the length of the domain: \');
n=input(\'Enter the number of grid points: \');
gamma=input(\'Enter the value of gamma: \');
%Specifying grid
x=linspace(0,L,n);
dx=x(2)-x(1);
nt=1400;
%initial Conditions
a=1+2.2*((x-1.5).^2);
for i=1:n
if (x(i)>=0 && x(i)<0.5)
rho_i(i)=1;
t_i(i)=1;
else
if (x(i)>=0.5 && x(i)<=1.5)
rho_i(i)=1-(0.366*(x(i)-0.5));
t_i(i)=1-(0.167*(x(i)-0.5));
else
if (x(i)>=1.5 && x(i)<=3.5)
rho_i(i)=0.634-(0.3879*(x(i)-1.5));
t_i(i)=0.833-(0.3507*(x(i)-1.5));
end
end
end
end
%velocity
v=(0.59)./((rho_i).*(a));
throat= (n+1)/2;
c=0.5;
%time step
for i=1:n
del_t(i)=min(c*(dx/((t_i(i)^0.5)+v(i))));
end
dt=min(del_t(i));
%variables
p=rho_i.*t_i;
U1=rho_i.*a;
U2=rho_i.*a.*v;
U3=rho_i.*((t_i./(gamma-1))+((gamma/2).*(v.^2))).*a;
tic
for k=1:nt
U1_old=U1;
U2_old=U2;
U3_old=U3;
%flux Variables
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));
%predictor step
for j=2:n-1
dF1dx=(F1(j+1)-F1(j))/dx;
dF2dx=(F2(j+1)-F2(j))/dx;
dF3dx=(F3(j+1)-F3(j))/dx;
dlogadx=(log(a(j+1))-log(a(j)))/dx;
dadx=(a(j+1)-a(j))/dx;
J2=(1/gamma).*rho_i(j).*t_i(j).*dadx;
%continuity equation
dU1dt_p(j)=-dF1dx;
%momentum equation
dU2dt_p(j)=-dF2dx+J2;
%energy equation
dU3dt_p(j)=-dF3dx;
%solution update
U1(j)= U1(j)+dU1dt_p(j)*dt;
U2(j)= U2(j)+dU2dt_p(j)*dt;
U3(j)= U3(j)+dU3dt_p(j)*dt;
end
%corrector step
rho_i=U1./a;
v=U2./U1;
t_i=(gamma-1)*((U3./U1)-((gamma/2)*((U2./U1).^2)));
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));
for j=2:n-1
dF1dx=(F1(j)-F1(j-1))/dx;
dF2dx=(F2(j)-F2(j-1))/dx;
dF3dx=(F3(j)-F3(j-1))/dx;
dlogadx=(log(a(j))-log(a(j-1)))/dx;
dadx=(a(j)-a(j-1))/dx;
J2=(1/gamma).*rho_i(j).*t_i(j).*dadx;
%continuity equation
dU1dt_c(j)=-dF1dx;
%momentum equation
dU2dt_c(j)=-dF2dx+J2;
%energy equation
dU3dt_c(j)=-dF3dx;
end
%avg derivative
dU1dt=0.5*(dU1dt_c+dU1dt_p);
dU2dt=0.5*(dU2dt_c+dU2dt_p);
dU3dt=0.5*(dU3dt_c+dU3dt_p);
%final solution update
for i=2:n-1
U1(i)=U1_old(i)+dU1dt(i)*dt;
U2(i)=U2_old(i)+dU2dt(i)*dt;
U3(i)=U3_old(i)+dU3dt(i)*dt;
end
%applying boundary conditions
%inlet
U1(1)=rho_i(1)*a(1);
U2(2)=rho_i(2)*a(2)*v(2);
U3(1)=U1(1)*((t_i(1)/(gamma-1))+((gamma/2)*(v(1)^2)));
%outlet
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);
rho_i=U1./a;
v=U2./U1;
t_i=(gamma-1)*((U3./U1)-((gamma/2)*((U2./U1).^2)));
%calculating solutions
%mass flow rate
m=rho_i.*a.*v;
%pressure
p=rho_i.*t_i;
%Mach no
mach=v./sqrt(t_i);
%throat values
rho_t=rho_i(throat);
t_t=t_i(throat);
p_t=p(throat);
mach_t=mach(throat);
m_t=m(throat);
v_t=v(throat);
figure(1);
if k==1
plot(x,m,\'--k\')
hold on
elseif k==50
plot(x,m,\'c\')
hold on
elseif k==100
plot(x,m,\'g\')
hold on
elseif k==200
plot(x,m,\'b\')
hold on
elseif k==800
plot(x,m,\'r\')
hold on
elseif k==1400
plot(x,m,\'y\')
hold on
title(\'Mass flow Rates(Conservative)\')
xlabel(\'X/L\')
ylabel(\'\\rhoVA/\\rho_{o}v{o}a{o}\')
legend(\'50\\Deltat\',\'100\\Deltat\',\'150\\Deltat\',\'200\\Deltat\',\'800\\Deltat\',\'1400\\Deltat\')
saveas(gcf,\'massflow rate non conservative\',\'jpeg\')
end
end
figure(2)
plot(x,rho_i)
title(\'Plot of Density vs Distance\')
xlabel(\'X/L\')
ylabel(\'\\rho/\\rho_{o}\')
saveas(gcf,\'Plot of Density vs Distance conservative\',\'jpeg\')
hold on
figure(3)
plot(x,v,\'g\')
title(\'Plot of Velocity vs Distance\')
xlabel(\'X/L\')
ylabel(\'v/v{o}\')
saveas(gcf,\'Plot of Velocity vs Distance conservative\',\'jpeg\')
hold on
figure(4)
plot(x,t_i,\'r\')
title(\'Plot of Temperature vs Distance\')
xlabel(\'X/L\')
ylabel(\'T/T{o}\')
saveas(gcf,\'Plot of Temperature vs Distance conservative\',\'jpeg\')
hold on
figure(5)
plot(x,mach,\'y\')
title(\'Plot of Mach number vs Distance\')
xlabel(\'X/L\')
ylabel(\'M\')
saveas(gcf,\'Plot of Mach number vs Distance conservative\',\'jpeg\')
hold on
figure(6)
plot(x,p,\'b\')
title(\'Plot of Pressure ratio vs Distance\')
xlabel(\'X/L\')
ylabel(\'P/P{o}\')
saveas(gcf,\'Plot of Pressure vs Distance conservative\',\'jpeg\')
simulation_time=toc;
sprintf(\'Computational Results from Conservative Method are: \\n Computation Time is = %d s\\n Density Ratio From Results = %d \\n Temperature Ratio= %d \\n Velocity Ratio=%d \\n Mach No=%d \\n Pressure Ratio=%d\',simulation_time,rho_t,t_t,v_t,mach_t,p_t)
Inputs:
Plots
A.Plot of Mass Flow Rate vs Distance of the Nozzle:
OUTPUT:
FOR GRID POINT:61
CONCLUSION FROM THE RESULTS:
-
- From the Results it can be seen that the values coming from the computational results are matching with the exact results in the Anderson Book.
- Grid Independence check is performed in the Both Conservative and Non Conservative part of the domain and results are not varying.
- Moreover the fluctuations in mass flow rate plot for conservative part is lesser than the non conservative part. Hence Conservative Part gives much satisfying results comparing to the non conservative part.
- The accuracy of the solution increases as we are increasing Grid Points.
- The rate of decay of residuals is much faster in non Conservative form.
- For same number of grid points Non Conservative solution yeilds better results compare to conservative results.