Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
Objective:
To write a code to solve the 1D supersonic nozzle flow equations using the Macormack Method
- Compare the results obtained using conservative and non-conservative forms of the governing equations
- Perform grid dependence test
Problem Setup:
Consider a steady , isentropic flow through a convergent -divergent nozzle as shown in the Fig.1.
The flow at the inlet to the nozzle comes from a reservoir where the pressure and temperature is constant.The cross-sectional area of the reservoir is so large that the velocity is almost Zero (`V~=0`).The flow expands isentropically to supersonic speeds at the nozzle exit .The flow is locally subsonic in the convergent section of the nozzle (`M<1`),sonic at the throat section (`M=1`) and supersonic at the divergent section (`M>1`).
Though in reality the flow field varies in the two-dimensional (xy plane), an assumption is made that the flow properties varies only along x axis ; the flow is defined as quasi-one-dimensional flow.
Fig.1
Governing Equations
A.Conservative Form
Non-Dimensional Form:
1.Continuity Equation:
`(delta(rho'.A'))/(deltat')+(delta(rho'.A'.v'))/(deltax')=0`
where:
`rho'=rho/rho_0` , `A'=A/(A_ast)` , `t'=t/((L/a_0))`
`v'=v/a_0` , `x'=x/L`
`rho_0rarr`Density of the reservoir `Lrarr`Length of the nozzle
`a_0rarr`Speed of sound in reservoir `A_astrarr`Sonic throat area
`rhorarr`Density of the fluid `vrarr`Velocity of the fluid
`Ararr`Area of the Nozzle
2.Momentum Equation:
`(delta(rho'.A'.v'))/(deltat')+(delta(rho'.A'.v' ^2+(1/gamma)rho'.A'))/(deltax')=(1/gamma).p'.(deltaA')/(deltax')`
where:
`p'=p/p_0`
`p_0rarr`Pressure of the reservoir `prarr`Pressure of the fluid
`gammararr`Specific heat ratio
3.Energy Equation:
`delta((rho'.((e')/(gamma-1)+gamma/2.v'^2).A'))/(deltat')+(delta(rho'.((e')/(gamma-1)+gamma/2.v'^2).v'.A'+p'.A'.v'))/(deltax')=0`
where:
`e'=e/e_0`
`e_0rarr`Energy of the reservoir
`e=C_V.T` `Trarr`Temperature
`e'=T'` `C_vrarr`Specific heat at constant volume
B.Non-Conservative Form
Non-Dimensional Form:
1.Continuity Equation:
`(deltarho')/(deltat')=-rho'.(deltav')/(deltax')-rho'.v'(delta(lnA'))/(deltax')-v'.(deltarho')/(deltax')`
2.Momentum Equation:
`(deltav')/(deltat')=-v'.(deltav')/(deltax')-(1/gamma)[(deltaT')/(deltax')+((T')/(rho')).(deltap')/(deltax')]`
3.Energy Equation:
`(deltaT')/(deltat')=-v'.(deltaT')/(deltax')-(gamma-1).T'[(deltav')/(deltax')+v'.(delta(ln(A')))/(deltax')]`
Macormack Method is used to solve the above partial differential equations .The derivations of the above equations and the background required to understand this project is explicated in the Chapter 7 of the book: John D. Anderson Jr. Computational Fluid Dynamics: The Basics with Applications.
Results:
1. Qualitative aspects of Quasi 1D nozzle flow at steady state:
Conservative form vs. Non-Conservative form(Grid-size n = 31)
It is clear that both the conservative and the non-conservative form return identical profiles at the steady state and the profiles matches exactly with the analytical values.
2. Non-Dimensional mass flow rate distribution at various time-steps:
a. Conservative Form: (Grid-size n = 31)
b. Non-Conservative Form: (Grid-size n = 31)
Although the non-conservation form and the conservation form have different initial conditions , it can be observed that variations in the mass flow rates are far less in the conservative form and closer to the exact analytical value of 0.579 .
3.Density, Temperature,Pressure,Mach number variation at the nozzle throat(Grid-size n = 31)
From the above graph , we can clearly observe that the conservative form converges quicker (approximately 400 iterations) and has lot less variations initially compared to non -conservative form .
Therefore, the minimum number of cycles for which the simulation should be run in order for convergence is around 400 for conservative form and around 1200 for non-conservative form.
After convergence , it can be observed that non-conservative form returns more accurate results compared to conservative form for the given grid size of n=31.
4. Comparison of the non-dimensional mass flow rates at steady state:
As evident from the plot the conservative form gives a mass flow rate variation that is close to a constant value and more accurate results at the steady state. The non-conservative form tends to show some oscillations at the inflow and outflow boundaries and is less satisfactory compared to the conservation form.
5.Grid dependency
For n=61
Non-dimensional mass flow rate comparision plot at steady state
Raise in the number of grid points to 61 has definitely increased the accuracy of the results. But the change is not so significant compared to results obtained from 31 grind points .Also at 61 grid points , both conservative and non-conservative results were very close to analytical solution.
6.Simulation time
For n=31
Execution time for Non-Conservative form is 2.46 Sec.
Execution time for Conservative form is 1.57 Sec. For n=61
Execution time for Non-Conservative form is 3.02 Sec.
Execution time for Conservative form is 1.87 Sec. Although convervative form requires more effort to code compared to non-conservative form. Once developed ,computational time for conservative form is much lesser than the non-conservative form.
Matlab Program:
Main Program
%********************* Nozzle Main program *******************************
clc
clear all
close all
%inputs
n=31; %number of nodes
x=linspace(0,3,n);
dx=x(2)-x(1);
gamma=1.4;
throat= ((n-1)/2)+1;% Throat point
C = 0.5; % Courant Number (specified for the problem)
% Non -Dimensional Area
A=1+(2.2*(x-1.5).^2);
%Time step
nt=1400;
dt=0.01;
tic
[rho,T,v,P,M,m,rho_throat,T_throat,v_throat,P_throat,M_throat,m_throat]= non_conservative_form(n,x,dx,gamma,throat,nt,dt,A);
Time_non_conservative=toc;
fprintf('Execution time for Non-Conservative form is %0.3g Sec. n' ,Time_non_conservative)
tic
[rho_c,T_c,v_c,P_c,M_c,m_c,rho_throat_c,T_throat_c,v_throat_c,P_throat_c,M_throat_c,m_throat_c]= conservative_form(n,x,dx,gamma,throat,nt,dt,A);
Time_conservative=toc;
fprintf('Execution time for Conservative form is %0.3g Sec. n ',Time_conservative)
figure(3)
subplot(2,3,1)
plot(x,rho,'r',x,rho_c,'b*')
xlabel('Length of nozzle(m)','FontWeight', 'bold')
ylabel('Non-Dimensional density','FontWeight', 'bold')
title('Density Profile','FontWeight', 'bold')
grid minor
legend('Non-Conservative','Conservative')
subplot(2,3,2)
plot(x,T,'r',x,T_c,'b*')
xlabel('Length of nozzle(m)','FontWeight', 'bold')
ylabel('Non-Dimensional Temperature','FontWeight', 'bold')
title('Temperature Profile','FontWeight', 'bold')
grid minor
legend('Non-Conservative','Conservative')
subplot(2,3,3)
plot(x,v,'r',x,v_c,'b*')
xlabel('Length of nozzle(m)','FontWeight', 'bold')
ylabel('Non-Dimensional Velocity','FontWeight', 'bold')
title('Velocity Profile','FontWeight', 'bold')
grid minor
legend('Non-Conservative','Conservative')
subplot(2,3,4)
plot(x,P,'r',x,P_c,'b*')
xlabel('Length of nozzle(m)','FontWeight', 'bold')
ylabel('Non-Dimensional Pressure','FontWeight', 'bold')
title('Pressure Profile','FontWeight', 'bold')
grid minor
legend('Non-Conservative','Conservative')
subplot(2,3,5)
plot(x,M,'r',x,M_c,'b*')
xlabel('Length of nozzle(m)','FontWeight', 'bold')
ylabel('Mach number','FontWeight', 'bold')
title('Mach number Profile','FontWeight', 'bold')
grid minor
legend('Non-Conservative','Conservative')
subplot(2,3,6)
plot(x,m,'r',x,m_c,'b*')
xlabel('Length of nozzle(m)','FontWeight', 'bold')
ylabel('Non-Dimensional Mass Flow Rate','FontWeight', 'bold')
title('Mass Flow Rate Profile','FontWeight', 'bold')
grid minor
legend('Non-Conservative','Conservative')
figure(4);
subplot(2,3,1)
plot(1:nt,rho_throat,'r',1:nt,rho_throat_c,'b')
hold on
line([0 nt],[0.634 0.634], 'color', 'm')
xlabel('Time(s)','FontWeight', 'bold')
ylabel('Non-Dimensional Density','FontWeight', 'bold')
title('Variation of Density at the Nozzle Throat','FontWeight', 'bold')
grid minor
legend('Non-Conservative','Conservative','Exact Solution')
subplot(2,3,2)
plot(1:nt,T_throat,'r',1:nt,T_throat_c,'b')
hold on
line([0 nt],[0.833 0.833], 'color', 'm')
xlabel('Time(s)','FontWeight', 'bold')
ylabel('Non-Dimensional Temperature','FontWeight', 'bold')
title('Variation of Temperature at the Nozzle Throat','FontWeight', 'bold')
grid minor
legend('Non-Conservative','Conservative','Exact Solution')
subplot(2,3,3)
plot(1:nt,v_throat,'r',1:nt,v_throat_c,'b')
xlabel('Time(s)','FontWeight', 'bold')
ylabel('Non-Dimensional Velocity','FontWeight', 'bold')
title('Variation of Velocity at the Nozzle Throat','FontWeight', 'bold')
grid minor
legend('Non-Conservative','Conservative')
subplot(2,3,4)
plot(1:nt,P_throat,'r',1:nt,P_throat_c,'b')
hold on
line([0 nt],[0.528 0.528], 'color', 'm')
xlabel('Time(s)','FontWeight', 'bold')
ylabel('Non-Dimensional Pressure','FontWeight', 'bold')
title('Variation of Pressure at the Nozzle Throat','FontWeight', 'bold')
grid minor
legend('Non-Conservative','Conservative','Exact Solution')
subplot(2,3,5)
plot(1:nt,M_throat,'r',1:nt,M_throat_c,'b')
hold on
line([0 nt],[1 1], 'color', 'm')
xlabel('Time(s)','FontWeight', 'bold')
ylabel(' Mach number','FontWeight', 'bold')
title('Variation of Mach number at the Nozzle Throat','FontWeight', 'bold')
grid minor
legend('Non-Conservative','Conservative','Exact Solution')
subplot(2,3,6)
plot(1:nt,m_throat,'r',1:nt,m_throat_c,'b')
hold on
line([0 nt],[0.579 0.579], 'color', 'm')
xlabel('Time(s)','FontWeight', 'bold')
ylabel('Non-Dimensional Mass Flow Rate','FontWeight', 'bold')
title('Variation of Mass Flow Rate at the Nozzle Throat','FontWeight', 'bold')
grid minor
legend('Non-Conservative','Conservative','Exact Solution')
Functions:
1.Conservatie_form
function [rho,T,v,P,M,m,rho_throat,T_throat,v_throat,P_throat,M_throat,m_throat]= conservative_form(n,x,dx,gamma,throat,nt,C,A)
%Calculating initial profiles
for i=1:n
if x(i)>=0 && x(i)<0.5
rho(i)=1;
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 if 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
end
end
% Defining non-dimensional velocity and non-dimensional pressure
v=0.59./(rho.*A);
P=rho.*T;
% Defining solution vectors U1, U2 and U3
U1=rho.*A;
U2=rho.*A.*v;
U3=rho.*A.*((T/(gamma-1))+(0.5*gamma.*v.*v));
%Outer time loop
for k=1:nt
% Copies of solution vectors
U1_ini=U1;
U2_ini=U2;
U3_ini=U3;
% Adaptable time step
dt=min(0.5.*dx./(sqrt(T) + v));
% Defining Flux vectors
F1=U2;
F2=((U2.^2)./(U1))+ ((gamma-1)/gamma).*(U3-(0.5*gamma.*U2.^2)./U1);
F3=((gamma.*U2.*U3)./U1)-(0.5*gamma*(gamma-1).*U2.^3)./(U1.^2);
% Predictor step
for i=2:n-1
dF1_dx=(F1(i+1)-F1(i))/dx;
dF2_dx=(F2(i+1)-F2(i))/dx;
dF3_dx=(F3(i+1)-F3(i))/dx;
% Continuity equation
dU1_dt_p(i)=-dF1_dx;
% Momentum equation
J2(i) = (1/gamma).*rho(i).*T(i).*((A(i+1)-A(i))/dx);
dU2_dt_p(i)=-dF2_dx+J2(i);
% Energy equation
dU3_dt_p(i)=-dF3_dx;
%solution update
U1(i)=U1(i)+dU1_dt_p(i)*dt;
U2(i)=U2(i)+dU2_dt_p(i)*dt;
U3(i)=U3(i)+dU3_dt_p(i)*dt;
end
% Updation of primitive variables and flux vectors
rho = U1./A;
v = U2./U1;
T = (gamma-1).*((U3./U1)-((gamma*v.^2)/2));
P=rho.*T;
F1=U2;
F2=((U2.^2)./(U1))+ ((gamma-1)/gamma).*(U3-(0.5*gamma.*U2.^2)./U1);
F3=((gamma.*U2.*U3)./U1)-(0.5.*gamma.*(gamma-1).*U2.^3)./(U1.^2);
% Corrector step
for j=2:n-1
dF1_dx=(F1(j)-F1(j-1))/dx;
dF2_dx=(F2(j)-F2(j-1))/dx;
dF3_dx=(F3(j)-F3(j-1))/dx;
% Continuity equation
dU1_dt_c(j)=-dF1_dx;
% Momentum equation
J2(j) = (1/gamma).*rho(j).*T(j).*((A(j)-A(j-1))/dx);
dU2_dt_c(j)=-dF2_dx+J2(j);
% Energy equation
dU3_dt_c(j)=-dF3_dx;
end
% Compute the average time derivative
dU1_dt_avg=0.5.*(dU1_dt_p+dU1_dt_c);
dU2_dt_avg=0.5.*(dU2_dt_p+dU2_dt_c);
dU3_dt_avg=0.5.*(dU3_dt_p+dU3_dt_c);
% Final solution update
for l=2:n-1
U1(l)=U1_ini(l)+dU1_dt_avg(l)*dt;
U2(l)=U2_ini(l)+dU2_dt_avg(l)*dt;
U3(l)=U3_ini(l)+dU3_dt_avg(l)*dt;
end
% Apply boundary conditions
% Inlet
U2(1) = 2*U2(2) - U2(3);
% 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);
% Updation of flow field variables
rho = U1./A;
v = U2./U1;
T = (gamma-1).*((U3./U1)-((gamma*v.^2)/2));
P=rho.*T;
% Calculating Mach number(M) and mass flow rate(m)
M= v./((T).^0.5);
m = rho.*v.*A;
% Capturing the values at throat for each iterations
rho_throat(k)=rho(throat);
v_throat(k)=v(throat);
T_throat(k)=T(throat);
P_throat(k)=P(throat);
M_throat(k)=M(throat);
m_throat(k)=m(throat);
% Ploting mass flow rate at different times
b=figure(2);
hold on;
if k ==50
plot(x,m,'b','linewidth',2);
else if k == 100
plot(x,m,'g','linewidth',2);
else if k == 150
plot(x,m,'r','linewidth',2);
else if k == 200
plot(x,m,'m','linewidth',2);
else if k == 700
plot(x,m,'-.','linewidth',2);
legend('50Deltat','100Deltat','150Deltat','200Deltat','700Deltat');
title('Variation of Mass flow rate at different time step(Conservative form)','FontWeight', 'bold', 'FontSize',12);
xlabel('Non-Dimensional Length of Nozzle (x/L)','FontWeight', 'bold', 'FontSize', 11 );
ylabel('Non -Dimensional Mass flow rate', 'FontWeight', 'bold', 'FontSize', 11)
grid minor
saveas(b, 'con_mass.png')
end
end
end
end
end
end
2. Non_conservative_form
function [rho,T,v,P,M,m,rho_throat,T_throat,v_throat,P_throat,M_throat,m_throat]= non_conservative_form(n,x,dx,gamma,throat,nt,C,A)
%Calculating initial profiles
rho=1-(0.3146*x);
T=1-(0.2314*x);
v=(0.1+(1.09*x)).*T.^(0.5);
%Outer time loop
for k=1:nt
rho_ini=rho;
T_ini=T;
v_ini=v;
% Time-step:
dt = 0.01;
% Predictor step
for i=2:n-1
dv_dx=(v(i+1)-v(i))/dx;
drho_dx=(rho(i+1)-rho(i))/dx;
dlogA_dx=(log(A(i+1))-log(A(i)))/dx;
dT_dx=(T(i+1)-T(i))/dx;
% continuity equation
drho_dt_p(i)=-rho(i)*dv_dx-rho(i)*v(i)*dlogA_dx-v(i)*drho_dx;
% momentum equation
dv_dt_p(i)=-v(i)*dv_dx-(1/gamma)*(dT_dx+((T(i)/rho(i))*drho_dx));
% energy equation
dT_dt_p(i)=-v(i)*dT_dx-(gamma-1)*T(i)*(dv_dx+(v(i)*dlogA_dx));
%solution update
rho(i)=rho(i)+drho_dt_p(i)*dt;
v(i)=v(i)+dv_dt_p(i)*dt;
T(i)=T(i)+dT_dt_p(i)*dt;
end
% Corrector step
for j=2:n-1
dv_dx_c=(v(j)-v(j-1))/dx;
drho_dx_c=(rho(j)-rho(j-1))/dx;
dlnA_dx_c=(log(A(j))-log(A(j-1)))/dx;
dT_dx_c=(T(j)-T(j-1))/dx;
% continuity equation
drho_dt_c(j)=-rho(j)*dv_dx_c-rho(j)*v(j)*dlnA_dx_c-v(j)*drho_dx_c;
% momentum equation
dv_dt_c(j)=-v(j)*dv_dx_c-(1/gamma)*(dT_dx_c+((T(j)/rho(j))*drho_dx_c));
% energy equation
dT_dt_c(j)=-v(j)*dT_dx_c-(gamma-1)*T(j)*(dv_dx_c+(v(j)*dlnA_dx_c));
end
% compute the average time derivative
drho_dt=0.5.*(drho_dt_p+drho_dt_c);
dv_dt=0.5.*(dv_dt_p+dv_dt_c);
dT_dt=0.5.*(dT_dt_p+dT_dt_c);
% Final solution update
for l=2:n-1
rho(l)=rho_ini(l)+drho_dt(l)*dt;
v(l)=v_ini(l)+dv_dt(l)*dt;
T(l)=T_ini(l)+dT_dt(l)*dt;
end
% Apply boundary conditions
% Inlet
v(1)=2*v(2)-v(3);
rho(1)=1;
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);
%Computing Pressure(P), Mach number(M) and mass flow rate(m)
P=rho.*T;
M=(v./(T).^0.5);
m = rho.*v.*A;
% Capturing the values at throat for each iterations
rho_throat(k)=rho(throat);
v_throat(k)=v(throat);
T_throat(k)=T(throat);
P_throat(k)=P(throat);
M_throat(k)=M(throat);
m_throat(k)=m(throat);
% Ploting mass flow rate at different times
a=figure(1);
hold on;
if k ==50
plot(x,m,'b','linewidth',2);
else if k == 100
plot(x,m,'g','linewidth',2);
else if k == 150
plot(x,m,'r','linewidth',2);
else if k == 200
plot(x,m,'m','linewidth',2);
else if k == 700
plot(x,m,'-.','linewidth',2);
legend('50Deltat','100Deltat','150Deltat','200Deltat','700Deltat');
title('Variation of Mass flow rate at different time step(Non-Conservative form)','FontWeight', 'bold', 'FontSize',12);
xlabel('Non-Dimensional Length of Nozzle (x/L)','FontWeight', 'bold', 'FontSize', 11 );
ylabel('Non -Dimensional Mass flow rate', 'FontWeight', 'bold', 'FontSize', 11)
grid minor
saveas(a, 'Non-con_mass.png')
end
end
end
end
end
end
endPlease note that the above scripts are uploaded here.