1 D Supersonic Nozzle Flow Simulation - Macormack Method
1D SUPERSONIC NOZZLE FLOW SIMULATION - MACORMACK METHOD
Objective
- To solve and compare the numerical solution of conservation and non conservation forms of governing equation.
- Creating user defined functions for calculating the flow parameters in both the forms.
- To demonstrate time wise variation of density , temperature , velocity and mach number at the Nozzle Throat.
- To obtain the minimum number of cycle for which the simulation should be run in order for convergence.
- To compare the normalized mass flow rate between the conservation forms and non conservation forms.
- To determine which form is the fastest.
- Perform Grid dependency test for various flow quantities
Introduction
The Conservation and Non Conservation equations are identically identical . Both are derived by applying the integral conservation method to control volume represented by the nozzle and therefore they are identical .
When the Control Volume is fixed in space with the fluid is moving then the equations governing the process shall be of the Conservative Form. And when the control volume is moving with the fluid in a sense that fluid is within the control volume then the Governing Equation shall be of Non Conservative Form.
Assumptions
Steady, Isentropic flow through the Convergent Divergent Nozzle .
At Inlet the flow comes from a reservoir where Pressure and Temperature are at Stagnation.
Cross Sectional Area of the reservoir is very large therefore the velocity is almost “ `V~=0` ”.
The flow expands adiabatically to Super Sonic Speeds at the Nozzle exit .
The Flow is Sub-Sonic in the Convergent section ie (M < 1) , Sonic at the Throat Section (M = 1) , Super - Sonic at the divergent (M > 1)
We define the necessary Non Dimensional Variables :
Temperature :`T' = T/T_0`
Density : `rho' = rho/rho_0`
Non Dimensional Distance :`x'=x/L`
Non - Dimensional Velocity :`V' = V/a_0` where ‘ `a_0`’ is the local sonic velocity(`a_0 = sqrt(gammaRT_0)`)
Non - Dimensional Time :`t' = t/(L/a_0)` where `L/a_0` is characterstic time scale
Non - Dimensional Area :`A'=A/(A^star)`
Following are the Governing Equations :
CONSERVATIVE FORM :
- a) Continuity Equation :`(del(rho'A'))/(delt') + (del(rho'A'V'))/(delx') = 0`
- b) Momentum Equation : `(del(rho'A'V'))/(delt') + (del(rho'A'V'^2+1/gammap'A'))/(delx')=1/gammap'(delA')/(delx')`
- c) Energy Equation :`(del[rho'((T')/(gamma-1)+gamma/2V'^2)A'])/(delt') + (del[rho'((T')/(gamma-1) +gamma/2V'^2)V'A' + p'A'V'])/(delx') = 0`
NON CONSERVATIVE FORM:
- a) Continuity Equation : `(delrho')/(delt') = rho'(delV')/(delx') - rho'V'(del(lnA))/(delx') - V'(delp')/(delx')`
- b) Momentum Equation :`(delV')/(delt') = V'(delV')/(delx') - 1/gamma((delT')/(delx') + (T')/(rho')(delrho')/(delx'))`
- c) Energy Equation :`(delT')/(delt') = V'(delT')/(delx') - (gamma-1)T'[(delV')/(delx') + V'(del(lnA'))/(delx')]`
The following equations are to be solved in the Generic form which is depicted in the Mat-lab Coding.
Main Script
clear all
close all
clc
%Inputs
n=31;
x=linspace(0,3,n);
dx=x(2)-x(1);
gamma=1.4;
%CFL number for the problem
C=0.5;
% Convergence control
error1=1e9;
error2=1e9;
tol=1e-5;
%Area
a=1+2.2*(x-1.5).^2;
n_throat=(n-1)/2;
tic
%% Solution By Non - Conservative form:
[rho , v , t , p , Mach ,m ,rho_throat , t_throat,p_throat,Mach_throat,k]=NonConservative_Form(n,x,dx,gamma,a,n_throat,C,error1,tol);
time_nc = toc;
fprintf(' Minimum Iteration (NonConservative Form): %f\n',k);
fprintf('Computation time (NonConservative Form): %f s\n',time_nc);
%Transient variation of properties at the Throat
Transient_Throat=figure(2);
subplot(4,1,1);
plot(rho_throat,'m');
title('Variations of the properties at the nozzle throat (Non-Conservative Form)');
legend('Density ratio');
ylabel('\rho/\rho_{o}');
axis([0 k 0.5 1]);
grid minor;
subplot(4,1,2);
plot(t_throat,'r');
legend('Temperature');
ylabel('T/T_{o}');
axis([0 k 0.6 1]);
grid minor;
subplot(4,1,3);
plot(p_throat,'g');
legend('Pressure');
ylabel('p/p_{o}');
axis([0 k 0.4 1]);
grid minor;
subplot(4,1,4);
plot(Mach_throat,'b');
legend('Mach Number');
ylabel('M');
axis([0 k 0.8 1]);
grid minor;
xlabel('Number of iterations');
%Steady profile
figure(3)
subplot(4,1,1);
plot(x,rho,'m');
title('Profiles of the properties at the Nozzle (Non-Conservative Form)');
ylabel('\rho/\rho_{o}');
legend('Density ratio','location','southwest');
grid minor;
subplot(4,1,2);
plot(x,t,'r');
ylabel('T/T_{o}');
legend('Temperature ratio','location','southwest');
grid minor;
subplot(4,1,3);
plot(x,p,'g');
ylabel('p/p_{o}');
legend('Pressure ratio','location','southwest');
grid minor;
subplot(4,1,4);
plot(x,Mach,'b');
ylabel('M');
legend('Mach Number','location','northwest');
grid minor;
%% Solution By Conservative form:
tic
[rho_c , v_c , t_c , p_c , Mach_c , m_c , rho_throat_c , t_throat_c , p_throat_c , Mach_throat_c ,k2]=Conservative_Form_(n,x,dx,gamma,a,n_throat,C,error1,tol);
time_c = toc;
fprintf('Minimum Iterations for Convergence (Conservative Form): %f\n',k2);
fprintf('Computation time (Conservative Form): %f s\n',time_c);
%Transient variation of properties at the throat
Transient_Throat_con=figure(5);
subplot(4,1,1);
plot(rho_throat_c,'m');
title('Timewise variations of the properties at the nozzle throat (Conservative Form)');
legend('Density ratio');
ylabel('\rho/\rho_{o}');
axis([0 k2 0.5 1]);
grid minor;
subplot(4,1,2);
plot(t_throat_c,'r');
legend('Temperature');
ylabel('T/T_{o}');
axis([0 k2 0.6 1]);
grid minor;
subplot(4,1,3);
plot(p_throat_c,'g');
legend('Pressure');
ylabel('p/p_{o}');
axis([0 k2 0.4 1]);
grid minor;
subplot(4,1,4);
plot(Mach_throat_c,'b');
legend('Mach Number');
ylabel('M');
axis([0 k2 0.8 1]);
grid minor;
xlabel('Number of iterations');
%Steady profile of the properties at the nozzle
figure(6)
subplot(4,1,1);
plot(x,rho_c,'m');
title('Profiles of the properties at the Nozzle (Conservative Form)');
ylabel('\rho/\rho_{o}');
legend('Density ratio','location','southwest');
grid minor;
subplot(4,1,2);
plot(x,t_c,'r');
ylabel('T/T_{o}');
legend('Temperature ratio','location','southwest');
grid minor;
subplot(4,1,3);
plot(x,p_c,'g');
ylabel('p/p_{o}');
legend('Pressure ratio','location','southwest');
grid minor;
subplot(4,1,4);
plot(x,Mach_c,'b');
ylabel('M');
legend('Mach Number','location','northwest');
grid minor;
%Steady massflow rate profile (Comparisson Non-Conservative vs Conservative)
figure(7)
plot(x,m)
hold on
plot(x,m_c)
title('Steady Mass-Flow Rate Profile')
legend('Non-Conservative','Conservative')
grid minorUser Defined Non Conservative Form
function [rho , v , t , p , Mach , m , rho_throat , t_throat , p_throat , Mach_throat ,k]=NonConservative_Form(n,x,dx,gamma,a,n_throat,C,error,tol)
%Initial the Variable
rho=1-0.3146*x;
t=1-0.2314*x;
v=(0.1+1.09*x).*t.^0.5;
%Initialization for the time loop
k=0;
%Outer time loop
while error>tol
k=k+1;
%saving initial data for flow comparisson
rho_old=rho;
v_old=v;
t_old=t;
a_old=sqrt(t_old);
m_old=rho_old.*v_old.*a;
%CFL Time Loop Control
tt=C*dx./(a_old+v_old);
dt=min(tt);
%Predictor Step
for j=2:n-1
%drho_dt_p =
%-rho(j)*(v(j+1)-v(j))/dx-rho(j)*v(j)*(log(a(j+1))-log(a(j)))/dx-v(j)*(rho(j+1)-rho(j))/dx;
dvdx=(v(j+1)-v(j))/dx;
dlogadx=(log(a(j+1))-log(a(j)))/dx;
drhodx=(rho(j+1)-rho(j))/dx;
dtdx=(t(j+1)-t(j))/dx;
%Continuity Equation
drho_dt_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(j)+drho_dt_p(j)*dt;
v(j)=v(j)+dvdt_p(j)*dt;
t(j)=t(j)+dtdt_p(j)*dt;
end
%Corrector Step
for j=2:n-1
%drho_dt_p= -rho(j)*(v(j+1)-v(j))/dx -rho(j)*v(j)*(log(A(j+1))-log(A(j)))/dx - v(j)*(rho(j+1)-rho(j))/dx
dvdx=(v(j)-v(j-1))/dx;
dlogadx=(log(a(j))-log(a(j-1)))/dx;
drhodx=(rho(j)-rho(j-1))/dx;
dtdx=(t(j)-t(j-1))/dx;
%Continuity Equation
drho_dt_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);
end
%Compute the Average Time derivative
drhodt_avg=0.5*(drho_dt_p+drho_dt_c);
dvdt_avg=0.5*(dvdt_p+dvdt_c);
dTdt_avg=0.5*(dtdt_p+dtdt_c);
for i=2:n-1
%Final Solution Update
rho(i)=rho_old(i)+drhodt_avg(i)*dt;
v(i)=v_old(i)+dvdt_avg(i)*dt;
t(i)=t_old(i)+dTdt_avg(i)*dt;
end
%Apply the boundary conditions
%Inlet
rho(1)=1;
v(1)=2*v(2)-v(3);
t(1)=1;
%Outlet
rho(n)=2*rho(n-1)-rho(n-2);
v(n)=2*v(n-1)-v(n-2);
t(n)=2*t(n-1)-t(n-2);
% Defining mass flow rate , pressure and mach number
Mach=v./sqrt(t); %Mach Number
m=rho.*v.*a; %Non-dimensional Mass Flow Rate
p=rho.*t; %Non-dimensional pressure
%Transient Throat Values
rho_throat(k)=rho(n_throat);
t_throat(k)=t(n_throat);
p_throat(k)=p(n_throat);
Mach_throat(k)=Mach(n_throat);
%Check time convergence
error=max(abs(m-m_old));
% massflow rate variation before achieving steady state
TransientMassFlow=figure(1);
hold on;
if k== 1
plot(x,m,'--k');
else if k == 50
plot(x,m,'r');
else if k == 100
plot(x,m,'g');
else if k == 150
plot(x,m,'c');
else if k == 200
plot(x,m,'b');
else if k == 700
plot(x,m,'-.b','linewidth',1.25);
plot(x(end),m(end),'k*');
end
end
end
end
end
end
end
title('Non-dimensional Mass Flow Rates(Non-Conservative Form)');
xlabel('x/L');
ylabel('\rhoVA/\rho_{o}V_{o}A_{o}');
legend('0\Deltat','50\Deltat','100\Deltat','150\Deltat','200\Deltat','700\Deltat');
axis([0 3 0.3 2]);
grid minor;
end User Defined Conservative Form
function [rho_c , v_c , t_c , p_c , Mach_c , m_c , rho_throat_c , t_throat_c , p_throat_c , Mach_throat_c , k2]=Conservative_Form_(n,x,dx,gamma,a,n_throat,C,error2,tol)
%Calculate Initial Profiles
for i = 1:length(x)
if x(i)>=0 && x(i)<0.5
rho_c(i)=1;
t_c(i)=1;
else if x(i)>=0.5 && x(i)< 1.5
rho_c(i)=1.0-0.366*(x(i)-0.5);
t_c(i)=1.0-0.167*(x(i)-0.5);
else if x(i)>=1.5 && x(i)<=3.0
rho_c(i)=0.634-0.3879*(x(i)-1.5);
t_c(i)=0.833-0.3507*(x(i)-1.5);
end
end
end
end
v_c =0.59./(rho_c.*a);
p_c=rho_c.*t_c;
%Calculate the initial solution vectors
%Solution Vector = U
%Flux Vector = F
%Source Term = J
% U1=ρ′A′
% U2=ρ′A′V′
% U3=ρ′(T′/γ-1+(γ/x)V′^2)A′
U1=rho_c.*a;
U2=rho_c.*v_c.*a;
U3=rho_c.*a.*((t_c./(gamma - 1)) + (gamma/2)*v_c.^2);
%Initialization for the time loop
k2=0;
%Outer time loop with a convergence criteria
while error2>tol
%Convergence variables
k2=k2+1;
massflow_old_con=rho_c.*v_c.*a;
%We store the Old Solution Vectors
U1_old=U1;
U2_old=U2;
U3_old=U3;
%Defining the flux vectors:
%F1=ρ′A′V′
% F1= U2;
F1= U2;
%F2=ρ′A′V′^2+1\γ*(p′A′)
%F2=(U2.^2./U1)+((gamma-1)/gamma)*(U3-(gamma*U2.^2)./(2*U1));
F2=(U2.^2./U1)+((gamma-1)/gamma)*(U3-(gamma*U2.^2)./(2*U1));
%F3=ρ′(T′γ-1+γ2V′2)V′A′+p′A′V′
% F3=(gamma*U2.*U3./U1)-(gamma*(gamma-1)/2)*(U2.^3./U1.^2);
F3=(gamma*U2.*U3./U1)-(gamma*(gamma-1)/2)*(U2.^3./U1.^2);
%CFL Adaptative Time Loop Control
a_con=sqrt(t_c);
TimeStep=C*dx./(a_con+v_c);
dt=min(TimeStep);
%Predictor Step
for j=2:n-1
% Defining the source term:
%J2=(1/γ)p(′∂A′/∂x′)
J2(j)=(1/gamma).*rho_c(j).*t_c(j).*((a(j+1)-a(j))/dx);
%Continuity Equation
dU1dt_p(j)=-(F1(j+1)-F1(j))/dx;
%Momentum Equation
dU2dt_p(j)=-(F2(j+1)-F2(j))/dx+J2(j);
%Energy Equation
dU3dt_p(j)=-(F3(j+1)-F3(j))/dx;
%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
% Predicted values of the primitive variables and flux vectors
rho_c=U1./a;
v_c=U2./U1;
t_c=(gamma-1)*(U3./U1-(gamma/2)*v_c.^2);
p_c=rho_c.*t_c;
%F1=U2
F1=U2;
%F2=(U2)^2U1+γ-1γ(U3-γ\2(U2)^2U1)
F2=(U2.^2./U1)+((gamma-1)/gamma)*(U3-(gamma*U2.^2)./(2*U1));
%F3=γU2U3U1-γγ-12(U2)3(U1)2
F3=(gamma*U2.*U3./U1)-(gamma*(gamma-1)/2)*(U2.^3./U1.^2);
%Corrector Step
for j=2:n-1
%Updating the source term:
%J2=γ-1γ(U3-γ2(U2)2U1)∂(lnA′)∂x′
J2(j)=(1/gamma)*rho_c(j).*t_c(j)*((a(j)-a(j-1))/dx);
%Continuity Equation
dU1dt_c(j)=-(F1(j)-F1(j-1))/dx;
%Momentum Equation
dU2dt_c(j)=-(F2(j)-F2(j-1))/dx+J2(j);
%Energy Equation
dU3dt_c(j)=-(F3(j)-F3(j-1))/dx;
end
%Compute the Average Time derivative
dU1dt_avg=0.5*(dU1dt_p+dU1dt_c);
dU2dt_avg=0.5*(dU2dt_p+dU2dt_c);
dU3dt_avg=0.5*(dU3dt_p+dU3dt_c);
%Final Solution Update
for i=2:n-1
U1(i)=U1_old(i)+dU1dt_avg(i)*dt;
U2(i)=U2_old(i)+dU2dt_avg(i)*dt;
U3(i)=U3_old(i)+dU3dt_avg(i)*dt;
end
%Apply the boundary conditions
%Inlet
rho_c(1)=1;
t_c(1)=1;
U1(1)=rho_c(1)*a(1);
U2(1)=2*U2(2)-U2(3);
U3(1)=U1(1)*(t_c(1)/(gamma-1)+(gamma/2)*v_c(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);
%Calculation of the new primitive variables
% Density:ρ=U1/A′
rho_c=U1./a;
% Velocity : V′=U2U1
v_c=U2./U1;
% Temperature:(T′=(γ-1)(U3U1-γ2V′2)
t_c=(gamma -1)*(U3./U1-(gamma/2)*v_c.^2);
% Pressure: p′=ρ′T′
p_c=rho_c.*t_c;
% Defining mass flow rate , pressure and mach number
Mach_c=v_c./sqrt(t_c); %Mach Number
m_c=rho_c.*v_c.*a; %Non-dimensional Mass Flow Rate
p_c=rho_c.*t_c; %Non-dimensional pressure
%Transient Throat Values
rho_throat_c(k2)=rho_c(n_throat);
t_throat_c(k2)=t_c(n_throat);
p_throat_c(k2)=p_c(n_throat);
Mach_throat_c(k2)=Mach_c(n_throat);
%Check time convergence
error2=max(abs(m_c-massflow_old_con));
% massflow rate variation before achieving steady state
TransientMassFlow_Con=figure(4);
hold on;
if k2 == 100
plot(x,m_c,'g');
else if k2 == 200
plot(x,m_c,'r');
else if k2 == 700
plot(x,m_c,'-.b','linewidth',1.25);
end
end
end
end
title('Non-dimensional Mass Flow Rates(Non-Conservative Form)');
xlabel('x/L');
ylabel('\rhoVA/\rho_{o}V_{o}A_{o}');
legend('100\Deltat','200\Deltat','700\Deltat');
axis([0 3 0.5 0.8]);
grid minor;
end UNDERSTANDING GRAPHS :
SOLVED FOR GRID POINTS : 31
Time Wise Variation of Pressure Temperature Velocity and Mach Number at Throat
The graph depicts variation of the main properties at the THROAT during Iterations.
Non – Conservative Form
We can observe the stability of the scheme and notice high/low results in the values. There is slight disturbances during the initial iterations and then the solution becomes stable. Mach number M = 1 at Throat.
Conservative Form
The solution becomes stable sooner to the non conservative form . The number of Iterations is more as we have stated the convergence condition in the mass flow rate
Mass Flow Profile at Different TimeStep
Mass flow rate is also used to check the solution convergence.
Non Conservative Form
At 0 timestep the curve obtained is variation of profile at initial condition , the curve is the product of initial density , velocity and nozzle area ratio. At 700 timestep we achieve the constant mass flow rate profile.
Conservative Form
The solution converges smoothly than the Non Conservative Form .
Steady Solution
The final steady profiles :
Non Conservative Form
Conservative Form
The density , temperature and pressure are dropping through the exit of the nozzle with the increase in Mach number. The result is identical to the analytical solution in the reference book.
Time and Iteration to Converge
The computational time and number of Iteration for each solver is mentioned.
Non Conservative Form
Minimum Iteration (Non Conservative Form): 767.000000
Computation time (Non Conservative Form): 1.155475 s
Conservative Form
Minimum Iterations for Convergence (Conservative Form): 951.000000
Computation time (Conservative Form): 1.464537 s
Since the number of iterations is based on mass flow rate criterion instead of the primitive variable therefore the number of iterations are more in the case of Conservative form.
Steady Mass Flow Comparison (Conservative & Non – Conservative)
Both the schemes are stable But the conservative form seems to be better than non conservative form as in non conservative form one can notice some form of oscillations .
Grid Dependency Test
Grid dependency is an important factor to determine if the problem can be solved faster .
Grid points are changed from 31 to 61 and result are displayed at index 15.
The results:
Conclusion
More accurate results were obtained in the Non Conservative solver.