Simulating 1D Quasi Supersonic Nozzle Flow using Macormack Method
OBJECTIVE:
To obtain the numerical solution of a 1D Quasi Supersonic Nozzle Flow using Macormack Method, compare the flow parametres in conservative and non-conservative forms of equation and perform grid dependence test.
INTRODUCTION:
Macormack Method
Macormach method is an explicit finite difference technique which provides second order accuracy in both space and time.It follows a similar time marching solution like the Lax-wendroff method. It includes two steps i.e the predictor step and the corrector step. In the predictor step forward differencing of the space terms are done while in the corrector step backward differencing of the space terms are done.
We assume the flow field variables are known at a time say \'t\' and then we proceed to evaulate the flow field variables for same grid at time \'t+dt\'. After done with both steps we compute the average of the predictor and the corrector step and finally update the solutions for next time loop.
For example a representative mean value of density``((delrho)/(delt))_(avg` between times t and t+dt for same grid point (i,j) would be :-
``((delrho)/(delt))_(avg)=1/2[((delrho_(1))/(delt))^t+((delrho_(2))/(delt))^(t+dt)]`
Flow Problem Statement
Consider a steady isentrophic flow in a convergent-divergent nozzle. let the flow to the inlet of the nozzle come from a reservoir where the stagnation pressure and temperature are `p_o` and `T_o`. The flow expands isentrophically to supersonic speeds at the nozzle exit. The flow is locally subsonic in the convergent section, sonic in the throat and supersonic in the divergent section. A mach number of 1 at the throat indicates that the local velocity of this location is equal to velocity of sound. We will now compute the flow field variables for both conservative and non conservative forms for this profile.
GOVERNING EQUATIONS (CONSERVATIVE AND NON-CONSERVATIVE FORMS) :
The Continuity, Momentum and Energy equations forms the basis of all fluid problems in real life. The conservative and non-conservative forms of governing equations are derived using the integral forms of continuity, momentum and energy equations. The conservative and non-conservative forms of equation are mathematically equivalent i.e can be converted from one form to another. They only differ in the way they are numerically discreatized. Other than that analytic results of both forms would be same.
For the derivations/equations for conservative and non-conservative forms one can refer to the book -Computational Fluid Dynamics by JHON D ANDERSON2 JR.
SCRIPTS
Non-conservative Solver
function [rho_con,v_con,t_con,m_con,mach_con,p_con,v_t_con,rho_t_con,t_t_con,m_t_con,mach_t_con,p_t_con]=conservative_form(n,L,c,nt,throat)
%variable to be input
x=linspace(0,L,n);
dx=x(2)-x(1);
%initial conditions
rho_con=1-0.3146*x;
t_con=1-0.2314*x;
v_con=(0.1+1.09*x).*t_con.^0.5;
%Area
a =1+2.2*(x-1.5).^2;
gamma=1.4;
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-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)
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);
p_con=rho_con.*t_con;
%sol vectors defination
u1 = rho_con.*a;
u2 = rho_con.*a.*v_con;
term1= t_con./(gamma-1);
term2 = (0.5*gamma*(v_con.^2));
u3 = u1.*(term1 + term2);
%%
for k = 1:nt
u1_old = u1;
u2_old = u2;
u3_old = u3;
%flux terms
f1 = u2;
term3 = ((gamma-1)/gamma)*(u3-(gamma*u2.^2)./(2*u1));
f2 = (u2.^2./u1)+ term3;
f3 =((gamma*u2.*u3)./u1)-(gamma*(gamma-1)/2)*(u2.^3./u1.^2);
a_con = sqrt(t_con);
time_step= (c*dx)./(a_con+v_con);
dt =min(time_step);
for i=2:n-1
j2(i)=(1/gamma).*rho_con(i).*t_con(i).*((a(i+1)-a(i))/dx);
df1_dx(i)= (f1(i+1)-f1(i))/dx;
df2_dx(i)= (f2(i+1)-f2(i))/dx;
df3_dx(i) =(f3(i+1)-f3(i))/dx;
%continuity equation
du1dt_p(i) = -df1_dx(i);
%momentum
du2dt_p(i) = -df2_dx(i)+j2(i);
% energy equation
du3dt_p(i) = -df3_dx(i);
u1(i)= u1(i)+du1dt_p(i)*dt;
u2(i)=u2(i)+du2dt_p(i)*dt;
u3(i)=u3(i)+du3dt_p(i)*dt;
end
rho_con = u1./a;
v_con = u2./u1;
t_con = (gamma-1)*((u3./u1)-((gamma*v_con.^2)/2));
p_con = rho_con.*t_con;
term3 = ((gamma-1)/gamma)*(u3-(gamma*u2.^2)./(2*u1));
f1 = u2;
f2 = ((u2.^2)./u1)+ term3;
f3 =(gamma*u2.*u3./u1)-(gamma*(gamma-1)/2)*(u2.^3./u1.^2);
% corrector step
for i = 2:n-1
j2(i)=(1/gamma)*rho_con(i).*t_con(i)*((a(i)-a(i-1))/dx);
% df1_dx(i)= (f1(i)-f1(i-1))/dx;
%df2_dx(i)= (f2(i)-f2(i-1))/dx;
%df3_dx(i) =(f3(i)-f3(i-1))/dx;
%continuity equation
du1dt_c(i) = -(f1(i)-f1(i-1))/dx;
%momentum
du2dt_c(i) = -(f2(i)-f2(i-1))/dx + j2(i);
% energy equation
du3dt_c(i) = -(f3(i)-f3(i-1))/dx;
end
%average derivative
du1dt = 0.5*(du1dt_p+du1dt_c);
du2dt = 0.5*(du2dt_p+du2dt_c);
du3dt = 0.5*(du3dt_p+du3dt_c);
%final sol update
for j=2:n-1
u1(j)=u1_old(j)+du1dt(j)*dt;
u2(j)=u2_old(j)+du2dt(j)*dt;
u3(j)=u3_old(j)+du3dt(j)*dt;
end
% apply BC's
%inlet
u1(1)=rho_con(1)*a(1);
u2(1)= 2*u2(2)-u2(3);
v_con(1)=u2(1)/u1(1);
t_con(1)=1;
%update u3
u3(1) = u1(1)*((t_con(1)/(gamma-1))+((gamma/2)*v_con(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-1);
rho_con = u1./a;
v_con=u2./u1;
t_con = (gamma-1)*((u3./u1)-((gamma*v_con.^2)/2));
%calculate mass flow rates and other variables
m_con=rho_con.*v_con.*a;
mach_con = v_con./((t_con).^0.5);
p_con=rho_con.*t_con;
%throat values input
rho_t_con(k)=rho_con(throat);
t_t_con(k)=t_con(throat);
p_t_con(k)= p_con(throat);
mach_t_con(k)= mach_con(throat);
v_t_con(k)=v_con(throat);
m_t_con(k)=m_con(throat);
figure(2);
xlabel('Non dimensional distance in nozzle (x/L)','linewidth',1.5);
ylabel(' \rhova/\rho_{o}v_{o}a_{o}');
title('Non dimensional Mass flow rate vs non dimensional distance (Conservative) ');
legend('50 time steps','100 time steps','150 time steps','200 time steps','500 time steps','700 time steps','1400 time steps') ;
if k==50
plot(x,m_con,'r','linewidth',1.5);
hold on
elseif k==100
plot(x,m_con,'b','linewidth',1.5);
hold on
elseif k==150
plot(x,m_con,'g','linewidth',1.5);
hold on
elseif k==200
plot(x,m_con,'c','linewidth',1.5);
hold on
elseif k==500
plot(x,m_con,'k','linewidth',1.5);
hold on
elseif k==700
plot(x,m_con,'y','linewidth',1.5);
hold on
elseif k==1400
plot(x,m_con,'m','linewidth',1.5);
hold on
end
end
end
Conservative Solver
function [rho_con,v_con,t_con,m_con,mach_con,p_con,v_t_con,rho_t_con,t_t_con,m_t_con,mach_t_con,p_t_con]=conservative_form(n,L,c,nt,throat)
%variable to be input
x=linspace(0,L,n);
dx=x(2)-x(1);
%initial conditions
rho_con=1-0.3146*x;
t_con=1-0.2314*x;
v_con=(0.1+1.09*x).*t_con.^0.5;
%Area
a =1+2.2*(x-1.5).^2;
gamma=1.4;
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-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)
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);
p_con=rho_con.*t_con;
%sol vectors defination
u1 = rho_con.*a;
u2 = rho_con.*a.*v_con;
term1= t_con./(gamma-1);
term2 = (0.5*gamma*(v_con.^2));
u3 = u1.*(term1 + term2);
%%
for k = 1:nt
u1_old = u1;
u2_old = u2;
u3_old = u3;
%flux terms
f1 = u2;
term3 = ((gamma-1)/gamma)*(u3-(gamma*u2.^2)./(2*u1));
f2 = (u2.^2./u1)+ term3;
f3 =((gamma*u2.*u3)./u1)-(gamma*(gamma-1)/2)*(u2.^3./u1.^2);
a_con = sqrt(t_con);
time_step= (c*dx)./(a_con+v_con);
dt =min(time_step);
for i=2:n-1
j2(i)=(1/gamma).*rho_con(i).*t_con(i).*((a(i+1)-a(i))/dx);
df1_dx(i)= (f1(i+1)-f1(i))/dx;
df2_dx(i)= (f2(i+1)-f2(i))/dx;
df3_dx(i) =(f3(i+1)-f3(i))/dx;
%continuity equation
du1dt_p(i) = -df1_dx(i);
%momentum
du2dt_p(i) = -df2_dx(i)+j2(i);
% energy equation
du3dt_p(i) = -df3_dx(i);
u1(i)= u1(i)+du1dt_p(i)*dt;
u2(i)=u2(i)+du2dt_p(i)*dt;
u3(i)=u3(i)+du3dt_p(i)*dt;
end
rho_con = u1./a;
v_con = u2./u1;
t_con = (gamma-1)*((u3./u1)-((gamma*v_con.^2)/2));
p_con = rho_con.*t_con;
term3 = ((gamma-1)/gamma)*(u3-(gamma*u2.^2)./(2*u1));
f1 = u2;
f2 = ((u2.^2)./u1)+ term3;
f3 =(gamma*u2.*u3./u1)-(gamma*(gamma-1)/2)*(u2.^3./u1.^2);
% corrector step
for i = 2:n-1
j2(i)=(1/gamma)*rho_con(i).*t_con(i)*((a(i)-a(i-1))/dx);
% df1_dx(i)= (f1(i)-f1(i-1))/dx;
%df2_dx(i)= (f2(i)-f2(i-1))/dx;
%df3_dx(i) =(f3(i)-f3(i-1))/dx;
%continuity equation
du1dt_c(i) = -(f1(i)-f1(i-1))/dx;
%momentum
du2dt_c(i) = -(f2(i)-f2(i-1))/dx + j2(i);
% energy equation
du3dt_c(i) = -(f3(i)-f3(i-1))/dx;
end
%average derivative
du1dt = 0.5*(du1dt_p+du1dt_c);
du2dt = 0.5*(du2dt_p+du2dt_c);
du3dt = 0.5*(du3dt_p+du3dt_c);
%final sol update
for j=2:n-1
u1(j)=u1_old(j)+du1dt(j)*dt;
u2(j)=u2_old(j)+du2dt(j)*dt;
u3(j)=u3_old(j)+du3dt(j)*dt;
end
% apply BC's
%inlet
u1(1)=rho_con(1)*a(1);
u2(1)= 2*u2(2)-u2(3);
v_con(1)=u2(1)/u1(1);
t_con(1)=1;
%update u3
u3(1) = u1(1)*((t_con(1)/(gamma-1))+((gamma/2)*v_con(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-1);
rho_con = u1./a;
v_con=u2./u1;
t_con = (gamma-1)*((u3./u1)-((gamma*v_con.^2)/2));
%calculate mass flow rates and other variables
m_con=rho_con.*v_con.*a;
mach_con = v_con./((t_con).^0.5);
p_con=rho_con.*t_con;
%throat values input
rho_t_con(k)=rho_con(throat);
t_t_con(k)=t_con(throat);
p_t_con(k)= p_con(throat);
mach_t_con(k)= mach_con(throat);
v_t_con(k)=v_con(throat);
m_t_con(k)=m_con(throat);
figure(2);
xlabel('Non dimensional distance in nozzle (x/L)','linewidth',1.5);
ylabel(' \rhova/\rho_{o}v_{o}a_{o}');
title('Non dimensional Mass flow rate vs non dimensional distance (Conservative) ');
legend('50 time steps','100 time steps','150 time steps','200 time steps','500 time steps','700 time steps','1400 time steps') ;
if k==50
plot(x,m_con,'r','linewidth',1.5);
hold on
elseif k==100
plot(x,m_con,'b','linewidth',1.5);
hold on
elseif k==150
plot(x,m_con,'g','linewidth',1.5);
hold on
elseif k==200
plot(x,m_con,'c','linewidth',1.5);
hold on
elseif k==500
plot(x,m_con,'k','linewidth',1.5);
hold on
elseif k==700
plot(x,m_con,'y','linewidth',1.5);
hold on
elseif k==1400
plot(x,m_con,'m','linewidth',1.5);
hold on
end
end
end
Main Script
clear all
close all
clc
n=61;
L=3;
x=linspace(0,L,n);
a =1+2.2*(x-1.5).^2; %area
c=0.5; %courants number
throat=((n-1)/2)+1;
nt=1450;% graph has been plotted for nt =1400
%% calculating time required to solve non conservative equation
tic
[rho,v,t,m,mach,p,v_t,rho_t,t_t,m_t,mach_t,p_t]=nonconservative_eq(n,L,c,nt,throat);
Time_required_nonconservative = toc;
%% variation of transient properties at throat non conservative
figure(3)
subplot(4,1,1);
plot(rho_t,'r','linewidth',1.5);
title('Variation of properties at nozzle with number of time step (Non-Conservative Form )');
legend('Density ratio');
ylabel('\rho/\rho_{o}');
axis([0 nt 0.5 1]);
grid minor;
subplot(4,1,2);
plot(t_t,'b');
legend('Temperature');
ylabel('T/T_{o}');
axis([0 nt 0.6 1.2]);
grid minor;
subplot(4,1,3);
plot(p_t,'g');
legend('Pressure');
ylabel('P/P_{o}');
axis([0 nt 0.4 1]);
grid minor;
subplot(4,1,4);
plot(mach_t,'m');
legend('Mach Number');
ylabel('M');
axis([0 nt 0.8 1.3]);
grid minor;
%% variation of steady properties at throat Non-conservative
figure(4)
subplot(4,1,1);
plot(x,rho,'r','linewidth',1.5);
title('Variation of profiles in 1D Supersonic Flow (Non-Conservative Form )');
legend('Density ratio');
ylabel('\rho/\rho_{o}');
grid minor;
subplot(4,1,2);
plot(x,t,'b');
legend('Temperature ratio');
ylabel('T/T_{o}');
grid minor;
subplot(4,1,3);
plot(x,p,'g');
legend('Pressure');
ylabel('P/P_{o}');
grid minor;
subplot(4,1,4);
plot(x,mach,'m');
legend('Mach Number');
ylabel('M');
grid minor;
%% calculating time required to solve conservative equation
tic
[rho_con,v_con,t_con,m_con,mach_con,p_con,v_t_con,rho_t_con,t_t_con,m_t_con,mach_t_con,p_t_con]=conservative_form(n,L,c,nt,throat);
Time_required_conservative=toc;
%% variation of transient properties at throat %conservative
figure(5)
subplot(4,1,1);
plot(rho_t_con,'r','linewidth', 1.5);
title('Variation of properties at nozzle with number of time step (Conservative Form )');
legend('Density ratio');
ylabel('\rho/\rho_{o}');
axis([0 nt 0.5 1]);
grid minor;
subplot(4,1,2);
plot(t_t_con,'b');
legend('Temperature');
ylabel('T/T_{o}');
axis([0 nt 0.6 1.2]);
grid minor;
subplot(4,1,3);
plot(p_t_con,'g');
legend('Pressure');
ylabel('P/P_{o}');
axis([0 nt 0.4 1]);
grid minor;
subplot(4,1,4);
plot(mach_t_con,'m');
legend('Mach Number');
ylabel('M');
axis([0 nt 0.8 1.3]);
grid minor;
%% variation of steady properties at throat conservative
figure(6)
subplot(4,1,1);
plot(x,rho_con,'r','linewidth',1.5);
title('Variation of profiles in 1D Supersonic Flow (Conservative Form )');
legend('Density ratio');
ylabel('\rho/\rho_{o}');
grid minor;
subplot(4,1,2);
plot(x,t_con,'b');
legend('Temperature ratio');
ylabel('T/T_{o}');
grid minor;
subplot(4,1,3);
plot(x,p_con,'g');
legend('Pressure');
ylabel('P/P_{o}');
grid minor;
subplot(4,1,4);
plot(x,mach_con,'m');
legend('Mach Number');
ylabel('M');
grid minor;
%% steady mass flow rate profile
figure(7)
plot(x,m,'r','linewidth',1.5);
hold on
plot(x,m_con,'b','linewidth',1.5);
hold on
plot(x,0.579,'*','linewidth',3.5);
title('Steady state mass flow rate comparison')
legend('non-conservative','conservative','Exact');
xlabel('Non-dimensional Mass flow rate');
ylabel('Non-dimensional distance across nozzle');
PLOTTING GRAPHS
1. Non-Conservative Form
a.) Variation of flow field variables at throat of nozzles with number of time teps
In the first few iterations there were slight disturbances causing a instable solution but after a time step of 200 the solution converges and becomes stable. Notice that the mach number reaches the value of 1 indicating that the flow is sonic in the throat. This result was expected.
b.) Steady State Profile Variables
There is a drop of density ratio, temperature ratio and pressure ratio at the exit of nozzle. Mach number value increases indicating the value of the fluid velocity is now more than velocity of sound i.e Super sonic flow condition.
c.) Non-Dimensional Transient Mass Flow Rate
Here we can notice that after about 500 iterations there is almost a constant mass flow rate 0.579 which is approximately equal to the analytical mass-flow rate . We will compare the time required and number of time steps to converge for conservative and non-conservative forms later.
2. Conservative Form
a.) Variation of flow field variables at throat of nozzles with number of time teps
Compared to the non conservative here the fluid flow variabes reach a stable state sooner. Hence convergence is faster and lesser variations in values are seen.
b.) Steady State Profile Variables
Here it is clearly understood that the conservative and non-conservative steady state profiles are alike.
c.) Non-Dimensional Transient Mass Flow Rate
Here we get to see that compared to non-conservative form, variation of mass flow rates in conservative form is lesser and more closer to the exact analytical solution 0.579.
STEADY STATE MASS FLOW RATE
From the above graph we can clearly conclude that the conservative form gives us a closer insight to the actual analytical mass flow rate value. On the other hand there are some un desirable oscillations in the non-conservative form in the inlet of nozzle. Hence conservative form is more accurate in this case.
TIME REQUIRED FOR CONVERGENCE
From the above picture we see that computational time for non-conservative solver is a bit less compared to conservative solver, but could a faster computational time compete over accuracy and stability? Accuracy and stability should be considered before computational time
GRID DEPENDENCY TEST
let us further dig into the concept of stability and accuracy by performing a grid dependency test by increasing the number of grid points from 31 to 61 and compute the same. We will try to figure out if the results obtained previously are dependent or independent of the grid size.
Density ratio
Exact Analytical 0.634
Non-Conservative (n=31) 0.6425
Non-Conservative (n=61) 0.6315
Conservative (n=31) 0.648
Conservative (n=61) 0.6412
There is no such significant change in the primitive variables hence independent of grid size. However the non-conservative form is seen to give better accuracy in terms of calculating primitive variables.