1D Super-sonic nozzle flow simulation using Macormark method technique
Mach number:
It is defined as the ratio of the speed of the body to speed of the soud
If Mach number<1,it is sub sonic flow
Mach number=1 ,it is sonic flow
Mach number>1, it is supersonic flow
The main program of the 1D Supersonic Nozzle flow byusing Macormack method as follows
The program for the conservative method
clear all
close all
clc
n=31;
x=linspace(0,3,n);
dx=x(2)-x(1);
gamma=1.4;
% area
a=1+(2.2*(x-1.5).^2);
%calculate the initial profiles
for i=1:n
if (x(i)>=0 && x(i)<=0.5)
rho(i)=1;
t(i)=1;
elseif (x(i)>0.5 && x(i)<1.5)
rho(i)=1-0.366.*(x(i)-0.5);
t(i)=1.0-0.167*(x(i)-0.5);
elseif (x(i)>=1.5 && x(i)<=3.5)
rho(i)=0.634-0.3879*(x(i)-1.5)
t(i)=0.833-0.3507*(x(i)-1.5)
end
end
v=0.59./(rho.*a);
%time step
nt=1400;
dt=0.0201;
% the solution vector
u1=rho.*a;
u2=rho.*a.*v;
u3= rho.*a.*((t/(gamma-1))+((gamma/2).*v.^2));
%outer loop
for k=1:nt
u1_old=u1;
u2_old=u2;
u3_old=u3;
% predicator method
for j=1:n-1
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));
df1dx(j)=(f1(j+1)-f1(j))/dx;
df2dx(j)=(f2(j+1)-f2(j))/dx;
df3dx(j)=(f3(j+1)-f3(j))/dx;
dlogadx(j)=((a(j+1)-a(j)))/dx;
J2=(1/gamma)*rho(j).*t(j).*dlogadx(j);
% continuity equation
du1dt_p(j)=-df1dx(j);
%Momentum equation
du2dt_p(j)=-df2dx(j)+J2;
% energy equation
du3dt_p(j)=-df3dx(j);
% 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
rho=u1./a;
v=u2./u1;
t=(gamma-1).*((u3/u1)-(0.5.*gamma*v.^2));
% corrector method
for j=2:n-1
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));
df1dx(j)=(f1(j)-f1(j-1))/dx;
df2dx(j)=(f2(j)-f2(j-1))/dx;
df3dx(j)=(f3(j)-f3(j-1))/dx;
dlogadx(j)=((a(j)-a(j-1)))/dx;
J2=(1/gamma)*rho(j).*t(j).*dlogadx(j);
% continuity equation
du1dt_c(j)=-df1dx(j);
%Momentum equation
du2dt_c(j)=-df2dx(j) + J2;
% energy equation
du3dt_c(j)=-df3dx(j);
end
% compute the average period
du1dt=0.5*(du1dt_p+du1dt_c);
du2dt=0.5*(du2dt_p+du2dt_c);
du3dt=0.5*(du3dt_p+du3dt_c);
% 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
% apply boundary conditions
% inlet boundary
u2(1)=2*u2(2)-u2(3);
% outlet boundary
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 the flow parameters
rho=u1./a;
v=u2./u1;
t=(gamma-1).*((u3/u1)-(0.5.*gamma*v.^2));
% calculation of the variables
m=v./sqrt(t); % Mach number
massflow=rho.*v.*a; % mass flow
p=rho.*t; % pressure
% throat parameters
rho_throat(k)=rho(16); % density at throat
p_throat(k)=p(16); % pressure at throat
t_throat(k)=t(16); %temperature at throat
m_throat(k)=m(16); % mach number at throat
% Mass flow rate profile for Consevation
figure(1)
hold on
if k==100
plot(x,massflow,\'color\',\'r\')
hold on
elseif k==150
plot(x,massflow,\'color\',\'y\')
hold on
elseif k==200
plot(x,massflow,\'color\',\'g\')
hold on
elseif k==700
plot(x,massflow,\'color\',\'k\')
xlabel(\'N0N Dimensional Mass flow\')
ylabel(\'NON Dimensional flow Through Nozzle\')
title(\'Mass Dimensional flow Transient state\')
legend(\'100Deltat\',\'150Deltat\',\'200Deltat\',\'700Deltat\')
end
%axis([0 3 0.5 0.8])
% Flow parameters
% Qualitative parameters of 1D Quasi Nozzle
figure(2)
subplot(4,1,1)
plot(rho_throat)
title(\'Time wise variations of Flowparameters at the Nozzle throat\')
ylabel(\'rho/rho_{o}\')
xlabel(\'Number of time steps\')
legend(\'Density Ratio\')
axis([0 k 0.5 1])
grid on;
subplot(4,1,2)
plot(p_throat)
ylabel(\'p/p_{o}\')
xlabel(\'Number of time steps\')
legend(\'Pressure Ratio\')
axis([0 k 0 1]);
grid on
subplot(4,1,3)
plot(t_throat)
ylabel(\'t/t_{o]\')
xlabel(\'Number of time steps\')
legend(\'Temperature Ratio\')
axis([0 k 0.5 1]);
grid on
subplot(4,1,4)
plot(m_throat)
ylabel(\'Mach number\')
xlabel(\'Number of time steps\')
legend(\'Mach Number\')
axis([0 k 0.5 4]);
figure(3)
subplot(4,1,1)
plot(x,m)
ylabel(\'Mach number\')
legend(\'Mach number\')
title(\'Qualitative aspects of 1D Quasi Nozzle\')
grid on
subplot(4,1,2)
plot(x,t)
ylabel(\'t/t_{o}\')
legend(\'Temperature Ratio\')
grid on
subplot(4,1,3)
plot(x,p)
ylabel(\'p/p_{o}\')
legend(\'Pressure Ratio\')
grid on
subplot(4,1,4)
plot(x,rho)
ylabel(\'rho/rho_{o}\')
legend(\'Density ratio\')
end The graphs for the conservative equations are obserevd below:
The Non Dimensional mass flow for the conservative is plotted with the different time values along the Domain of the nozzle.It looks like stable at 700 time steps.
The time wise variations of flowparameters at the nozzle throat are observed that unstable at the initial iterations and it becomes stable almost for maximum iterations.
The staedy state profile for the Mach number,temperature ratio,pressure ratio,Density ratio are plotted along the length of the nozzle.FRom the figure we can conclude that the Density,pressure and temperature ratio are dropping as the fluid accelerating along the domain.
Initially the Mach number is low because it is sub sonic where as in middle it is sonic at the throat and at the end the mach number is gradually increase because it is supersonic.
Conservative N=61
The program is same as above.just change the value of n=61
The following graphs are observed
The variations of the properties at the throat with the change of iterations are shown in the plot.It can be observed that in between 500 and 600 the flow simulation has achieved stability.whereas,in initial it obseved that a peak changes occur .
The staedy state profile for the Mach number,temperature ratio,pressure ratio,Density ratio are plotted along the length of the nozzle.From the figure we can conclude that the Density,pressure and temperature ratio are dropping as the fluid accelerating along the domain.
Initially the Mach number is low because it is sub sonic where as in middle it is sonic at the throat and at the end the mach number is gradually increase because it is supersonic.
2)Non conservative
% non conservative equation of 1d supersonic Nozzle
clear all
close all
clc
n=31;
x=linspace(0,3,n);
dx=x(2)-x(1);
gamma=1.4;
%calculate the initial profiles
rho=1-(0.3146*x); % density
t=1-(0.2314*x); %temperature
v=(0.1+(1.09*x)).*t.^0.5; % velocity
% area
a=1+(2.2*(x-1.5).^2);
%Time step
nt=1000;
dt=0.01;
% outer time loop
for k=1:nt
rho_old=rho;
v_old=v;
t_old=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;
%continuity equation
drhodt_p(j)=-(rho(j)*dvdx)-rho(j)*v(j)*dlogadx-v(j)*drhodx;
% Momentum equation
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
rho(j)=rho_old(j)+drhodt_p(j)*dt;
v(j)=v_old(j)+dvdt_p(j)*dt;
t(j)=t_old(j)+dtdt_p(j)*dt;
end
% corrector method
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;
%continuity equation
drhodt_c(j)=-(rho(j)*dvdx)-rho(j)*v(j)*dlogadx-v(j)*drhodx;
% Momentum equation
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);
% solution update
rho(j)=rho_old(j)+drhodt_c(j)*dt;
v(j)=v_old(j)+dvdt_c(j)*dt;
t(j)=t_old(j)+dtdt_c(j)*dt;
end
%compute the average time method
drhodt=0.5*(drhodt_p+drhodt_c);
dvdt=0.5*(dvdt_p+dvdt_c);
dtdt=0.5*(dtdt_p+dtdt_c);
%final solution update
for i=2:n-1
rho(i)=rho_old(i)+drhodt(i)*dt;
v(i)=v_old(i)+dvdt(i)*dt;
t(i)=t_old(i)+dtdt(i)*dt;
end
% apply boundary conditions
% inlet
rho(1)=1;
v(1)=2*v(2)-v(3);
t(1)=1;
% outlet
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);
% calculating the solutions
% Mass flow rate
massflow=rho.*a.*v;
% pressure
p=rho.*t;
%mach number
m=v./((t).^0.5);
% throat parameters
rho_throat(k)=rho(16); % density at throat
p_throat(k)=p(16); % pressure at throat
t_throat(k)=t(16); %temperature at throat
m_throat(k)=m(16); % mach number at throat
% to check the mass flow rate
figure(1)
hold on
if k==1
plot(x,massflow,\'color\',\'y\')
hold on
elseif k==50
plot(x,massflow,\'color\',\'g\')
hold on
elseif k==100
plot(x,massflow,\'color\',\'k\')
hold on
elseif k==150
plot(x,massflow,\'color\',\'b\')
hold on
elseif k==200
plot(x,massflow,\'color\',\'r\')
elseif k==700
plot(x,massflow,\'color\',\'m\')
plot(x(end),massflow(end),\'color\',\'c\')
end
xlabel(\'N0N Dimensional Mass flow\')
ylabel(\'NON Dimensional flow Through Nozzle\')
title(\'NON Dimensional flow Transient state\')
legend(\'0Deltat\',\'50Deltat\',\'100Deltat\',\'150Deltat\',\'200Deltat\',\'700Deltat\')
axis([0 3 0 2]);
% Flow parameters
figure(2)
subplot(4,1,1)
plot(rho_throat)
title(\'Time wise variations of Flowparameters at the Nozzle throat\')
ylabel(\'rho/rho_{o}\')
xlabel(\'Number of time steps\')
legend(\'Density Ratio\')
axis([0 k 0.5 1])
subplot(4,1,2)
plot(p_throat)
ylabel(\'p/p_{o}\')
xlabel(\'Number of time steps\')
legend(\'Pressure Ratio\')
axis([0 k 0.5 1]);
subplot(4,1,3)
plot(t_throat)
ylabel(\'t/t_{o]\')
xlabel(\'Number of time steps\')
legend(\'Temperature Ratio\')
axis([0 k 0.5 1]);
subplot(4,1,4)
plot(m_throat)
ylabel(\'Mach number\')
xlabel(\'Number of time steps\')
legend(\'Mach Number\')
axis([0 k 0.9 2]);
end
% Qualitative parameters of 1D Quasi Nozzle
figure(3)
subplot(4,1,1)
plot(x,m)
ylabel(\'Mach number\')
legend(\'Mach number\')
title(\'Qualitative aspects of 1D Quasi Nozzle\')
grid on
subplot(4,1,2)
plot(x,t)
ylabel(\'t/t_{o}\')
legend(\'Temperature Ratio\')
grid on
subplot(4,1,3)
plot(x,p)
ylabel(\'p/p_{o}\')
legend(\'Pressure Ratio\')
grid on
subplot(4,1,4)
plot(x,rho)
ylabel(\'rho/rho_{o}\')
legend(\'Density ratio\')The Non Dimensional mass flow for the Non-conservative is plotted with the different time values along the Domain of the nozzle.It looks like stable at 700 time steps.
Comparing the both Conservation form of the equations and the Non conservativon form of the equations we can observe that the stability occurs after the time step 700.
The variations of the properties at the throat with the change of iterations are shown in the plot.It can be observed that in between 500 and 600 the flow simulation has achieved stability.whereas,in initial it obseved that a peak changes occur .
The staedy state profile for the Mach number,temperature ratio,pressure ratio,Density ratio are plotted along the length of the nozzle.From the figure we can conclude that the Density,pressure and temperature ratio are dropping as the fluid accelerating along the domain.
Initially the Mach number is low because it is sub sonic where as in middle it is sonic at the throat and at the end the mach number is gradually increase because it is supersonic.
NON-Conservative N=61:
For the non conservative n=61 the program is same as above.just change the value of n=61
The following graphs are observed
The Non Dimensional mass flow for the Non-conservative is plotted with the different time values along the Domain of the nozzle.It looks like stable at 700 time steps.
The variations of the properties at the throat with the change of iterations are shown in the plot.It can be observed that at 500 time step the flow simulation has achieved stability.whereas,in initial it obseved that a peak changes occur .
The staedy state profile for the Mach number,temperature ratio,pressure ratio,Density ratio are plotted along the length of the nozzle.From the figure we can conclude that the Density,pressure and temperature ratio are dropping as the fluid accelerating along the domain.
Initially the Mach number is low because it is sub sonic where as in middle it is sonic at the throat and at the end the mach number is gradually increase because it is supersonic.
3)The main program snippet for the Mass flow of both Conservative and Non conservative by using the functions
Main program:
clear all
close all
clc
n=31;
x=linspace(0,3,n)
dx=x(2)-x(1);
gamma=1.4;
a=1+(2.2*(x-1.5).^2);
% calling the function of conservative and Non conservative
tic
[rho,v,t,p,m,massflow]=Non_conservative_massflow(n,x,dx)
time=toc;
fprintf(\'The time elapsed for the Non conservative form is %f\',time);
%calling the function of conservative
tic
[rho_con,v_con,t_con,p_con,m_con,massflow_con]=conservative_massflow(n,x,dx)
time_elapsed=toc;
fprintf(\'The time elapsed for the conservative form is %f\',time_elapsed);
figure(1)
axis([0 3 0.58 0.6])
plot(x,massflow,\'color\',\'r\')
hold on
plot(x,massflow_con,\'color\',\'b\')
xlabel(\'x/L\')
ylabel(\'Mass dimensional flow\')
title(\'comparison of conservative and Non conservative Mass flows\')
legend(\'conservative\',\'Non conservative\')The conservative function:
function[rho_con,v_con,t_con,p_con,m_con,massflow_con]=conservative_massflow(n,x,dx)
% area
gamma=1.4;
a=1+(2.2*(x-1.5).^2);
%calculate the initial profiles
for i=1:n
if (x(i)>=0 && x(i)<=0.5)
rho_con(i)=1;
t_con(i)=1;
elseif (x(i)>0.5 && x(i)<1.5)
rho_con(i)=1-0.366.*(x(i)-0.5);
t_con(i)=1.0-0.167*(x(i)-0.5);
elseif (x(i)>=1.5 && x(i)<=3.5)
rho_con(i)=0.634-0.3879*(x(i)-1.5)
t_con(i)=0.833-0.3507*(x(i)-1.5)
end
end
v_con=0.59./(rho_con.*a);
%time step
nt=1400;
dt=0.0201;
% the solution vector
u1=rho_con.*a;
u2=rho_con.*a.*v_con;
u3= rho_con.*a.*((t_con/(gamma-1))+((gamma/2).*v_con.^2));
%outer loop
for k=1:nt
u1_old=u1;
u2_old=u2;
u3_old=u3;
% predicator method
for j=1:n-1
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));
df1dx(j)=(f1(j+1)-f1(j))/dx;
df2dx(j)=(f2(j+1)-f2(j))/dx;
df3dx(j)=(f3(j+1)-f3(j))/dx;
dlogadx(j)=((a(j+1)-a(j)))/dx;
J2=(1/gamma)*rho_con(j).*t_con(j).*dlogadx(j);
% continuity equation
du1dt_p(j)=-df1dx(j);
%Momentum equation
du2dt_p(j)=-df2dx(j)+J2;
% energy equation
du3dt_p(j)=-df3dx(j);
% 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
rho_con=u1./a;
v_con=u2./u1;
t_con=(gamma-1).*((u3/u1)-(0.5.*gamma*v_con.^2));
% corrector method
for j=2:n-1
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));
df1dx(j)=(f1(j)-f1(j-1))/dx;
df2dx(j)=(f2(j)-f2(j-1))/dx;
df3dx(j)=(f3(j)-f3(j-1))/dx;
dlogadx(j)=((a(j)-a(j-1)))/dx;
J2=(1/gamma)*rho_con(j).*t_con(j).*dlogadx(j);
% continuity equation
du1dt_c(j)=-df1dx(j);
%Momentum equation
du2dt_c(j)=-df2dx(j) + J2;
% energy equation
du3dt_c(j)=-df3dx(j);
end
% compute the average period
du1dt=0.5*(du1dt_p+du1dt_c);
du2dt=0.5*(du2dt_p+du2dt_c);
du3dt=0.5*(du3dt_p+du3dt_c);
% 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
% apply boundary conditions
% inlet boundary
u2(1)=2*u2(2)-u2(3);
% outlet boundary
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 the flow parameters
rho_con=u1./a;
v_con=u2./u1;
t_con=(gamma-1).*((u3/u1)-(0.5.*gamma*v_con.^2));
% calculation of the variables
m_con=v_con./sqrt(t_con); % Mach number
massflow_con=rho_con.*v_con.*a; % mass flow
p_con=rho_con.*t_con; % pressure
end
end
Non Conservative function:
function[rho,v,t,p,m,massflow]=Non_conservative_massflow(n,x,dx,nt)
gamma=1.4;
%calculate the initial profiles
rho=1-(0.3146*x); % density
t=1-(0.2314*x); %temperature
v=(0.1+(1.09*x)).*t.^0.5; % velocity
% area
a=1+(2.2*(x-1.5).^2);
%Time step
nt=1000;
dt=0.001;
% outer time loop
for k=1:nt
rho_old=rho;
v_old=v;
t_old=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;
%continuity equation
drhodt_p(j)=-(rho(j)*dvdx)-rho(j)*v(j)*dlogadx-v(j)*drhodx;
% Momentum equation
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
rho(j)=rho_old(j)+drhodt_p(j)*dt;
v(j)=v_old(j)+dvdt_p(j)*dt;
t(j)=t_old(j)+dtdt_p(j)*dt;
end
% corrector method
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;
%continuity equation
drhodt_c(j)=-(rho(j)*dvdx)-rho(j)*v(j)*dlogadx-v(j)*drhodx;
% Momentum equation
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);
% solution update
rho(j)=rho_old(j)+drhodt_c(j)*dt;
v(j)=v_old(j)+dvdt_c(j)*dt;
t(j)=t_old(j)+dtdt_c(j)*dt;
end
%compute the average time method
drhodt=0.5*(drhodt_p+drhodt_c);
dvdt=0.5*(dvdt_p+dvdt_c);
dtdt=0.5*(dtdt_p+dtdt_c);
%final solution update
for i=2:n-1
rho(i)=rho_old(i)+drhodt(i)*dt;
v(i)=v_old(i)+dvdt(i)*dt;
t(i)=t_old(i)+dtdt(i)*dt;
end
% apply boundary conditions
% inlet
rho(1)=1;
v(1)=2*v(2)-v(3);
t(1)=1;
% outlet
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);
% calculating the solutions
% Mass flow rate
massflow=rho.*a.*v;
% pressure
p=rho.*t;
%mach number
m=v./((t).^0.5);
end
end
Time dependence:
1)N=31
The time elapsed for the Non-conservative form is 0.105953
The time elapsed for the conservative form is 1.073888
2)N=61
The time elapsed for the Non-conservative form is 0.147297
The time elapsed for the conservative form is 3.348838
Hence from above results we can conclude that the Non_conservative form of equations execution time is less compared to the conservative form of equations.