Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
Introduction:
In this challenge, we consider the flow through a convergent-divergent nozzle. The flow at the inlet to the nozzle comes from a reservoir where the cross-sectional area of the reservoir is large(Theoretically `Trarroo`) and hence the velocity is small (`v rarr0`). The flow expands isentropically to supersonic speeds at the nozzle exit. The flow is locally subsonic at a convergent section of the nozzle, sonic at the throat of the nozzle.
Fig.1 Convergent divergent nozzle
We assume that the flow properties are the same across any given cross-section. the assumption of the quasi-one-dimensional flow dictates that the flow properties are the function of x only. The quasi-one-dimensional flow through the nozzle has been modeled using MacCormack\'s technique for both conservative and nonconservative form of the equation.
The non-conservative form:
The nondimensional continuity equation
`(delrho\')/(delt\')=-rho\'*(delv\')/(delx\')-rho\'*v\'*(delln(A\'))/(delx\')-v\'*(delrho\')/(delx\')`
The nondimensional momentum equation
`(delv\')/(delt\')=-v\'*(delv\')/(delx\')-(1/gamma)*((delT\')/(delx\')+((T\')/(rho\'))*(delrho\')/(delx\'))`
The nondimensional energy eqaution
`(delT\')/(delt\')=-v\'*(delT\')/(delx\')-(gamma-1)*T\'*((delv\')/(delx\')+v\'*(delln(A\'))/(delx\'))`
The conservative form:
The nondimensional continuity equation
`(delU1)/(delt\')=-(delF1)/(delx\')`
The nondimensional momentum equation
`(delU2)/(delt\')=-(delF2)/(delx\')+J`
The nondimensional energy eqaution
`(delU3)/(delt\')=-(delF3)/(delx\')`
Where,
Solution vector, `U1=rho\'*A\'`
`U1=rho\'*A\'*v\'`
`U3=rho\'*A\'*((e\')/(gamma-1)+(gamma/2)*v\'^2)`
Flux vector, `F1=rho\'*A\'*v\'`
`F2=rho\'*A\'*v\'^2+(1/gamma)*rho\'*A\'`
`F3=rho\'*v\'*A\'*((e\')/(gamma-1)+(gamma/2)*v\'^2)+rho\'*v\'*A\'`
Source term `J=(1/gamma)*p\'*(delA\')/(delx\')`
Mesh details:
Total length of the nozzle = 3m
Number of grid points = 31
Step size dx = 0.1m
Profile of the nozzle:
Cross section area `A=1+2.2*(x-1.5)^2`
Time step size calculation:
The governing equations are hyperbolic with respect to time. hence, the stability criteria is based on the Courant number.
`Deltat=C*(Deltax)/(a+V)`
where a = Local speed of sound
v= Flow velocity
For a stable solution(explicit scheme) C<=1
MacCormack\'s technique:
It is an explicit finite-difference technique which is second-order accurate in both space and time. A two-step Predictor-corrector sequence is used with the forward difference on predictor and with a backward difference on a corrector.
Matlab code for solving 1D supersonic nozzle flow simulation using MacCormack\'s technique:
Main code:
clear all
close all
clc
%inputs
n=31;
x=linspace(0,3,n);
dx=x(2)-x(1);
th=(n+1)/2;
a=1+2.2*(x-1.5).^2;
%constants
gamma = 1.4;
%time steps
C = 0.5;
nt=1000;
tic
[v_nc,rho_nc,t_nc,m_nc,M_nc,P_nc]=non_conservative(a,gamma,n,dx,nt,x,C);
t1=toc;
tic
[v_c,rho_c,t_c,U2,M_c,P_c]=conservative(a,gamma,n,dx,nt,x,C);
t2=toc;
%Grid dependence test
%Caluculating primitive values at throat
terms={\'Nondimensional density\';\'Nondimensional temperature\';\'Nondimensional pressure\';\'Mach number\'};
exact_solution=[0.634;0.833;0.528;1];
conservative_terms=[rho_c(th);t_c(th);P_c(th);M_c(th)];
non_conservative_terms=[rho_nc(th);t_nc(th);P_nc(th);M_nc(th)];
table(terms, exact_solution, conservative_terms, non_conservative_terms)
figure(3)
plot(x,m_nc)
hold on
plot(x,U2)
title(\'Comparison of non-dimensional mass flow rates\')
legend(\'Non-Conservative form\',\'Conservative form\')
xlabel(\'x/L\')
ylabel(\'Nondimensional mass flow rate\')
figure(4)
subplot(4,1,1)
plot(x,M_nc)
xlabel(\'x-axis\')
ylabel(\'Mach number\')
subplot(4,1,2)
plot(x,rho_nc)
xlabel(\'x-axis\')
ylabel(\'Density ratio\')
subplot(4,1,3)
plot(x,P_nc)
xlabel(\'x-axis\')
ylabel(\'Pressure ratio\')
subplot(4,1,4)
plot(x,t_nc)
xlabel(\'x-axis\')
ylabel(\'Temperature ratio\')
suptitle(\'Solution for Conservative form at steady state\')
figure(5)
subplot(4,1,1)
plot(x,M_c)
xlabel(\'x-axis\')
ylabel(\'Mach number\')
subplot(4,1,2)
plot(x,rho_c)
xlabel(\'x-axis\')
ylabel(\'Density ratio\')
subplot(4,1,3)
plot(x,P_c)
ylim([0 1])
xlabel(\'x-axis\')
ylabel(\'Pressure ratio\')
subplot(4,1,4)
plot(x,t_c)
ylim([0 1])
xlabel(\'x-axis\')
ylabel(\'Temperature ratio\')
suptitle(\'Solution for Non-Conservative form at steady state\')Matlab function for conservative form:
function [v_c,rho_c,t_c,U2,M_c,P_c]=conservative(a,gamma,n,dx,nt,x,C)
for i=1:n
if (0<=x(i))&&(x(i)<=0.5)
rho_c(i)=1;
t_c(i)=1;
end
if (0.5<=x(i))&&(x(i)<=1.5)
rho_c(i)=1-0.366*(x(i)-0.5);
t_c(i)=1-0.167*(x(i)-0.5);
end
if (1.5<=x(i))&&(x(i)<=3.5)
rho_c(i)=0.634-0.3879*(x(i)-1.5);
t_c(i)=0.833-0.3507*(x(i)-1.5);
end
end
v_c=0.59./(rho_c.*a);
U1=rho_c.*a;
U2=U1.*v_c;
U3=U1.*((t_c/(gamma-1))+(gamma/2).*(v_c.^2));
% Non conservative method
for k=1:nt
%Determinig the time step
A=t_c.^(1/2);
dt=C*(dx./(A+v_c));
dt=min(dt);
%Updating the old values
U1old=U1;
U2old=U2;
U3old=U3;
%Determinig the flux values
F1=U2;
term1=(U2.^2)./U1;
term2=U3-((gamma/2)*term1);
F2=term1+((gamma-1)/gamma)*term2;
term3=(U2.*U3)./U1;
term4=(U2.^3)./(U1.^2);
F3=gamma.*term3-((gamma*(gamma-1))/2).*term4;
%Predictor method
for j=2:n-1
dU1dt_p(j) = -(F1(j+1)-F1(j))/dx;
dU2dt_p(j) = (-(F2(j+1)-F2(j))/dx)+(1/gamma)*(rho_c(j)*t_c(j))*((a(j+1)-a(j))/dx);
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
%Determining the predicted primitive values
rho_c=U1./a;
v_c=U2./U1;
t_c=(gamma-1)*((U3./U1)-(gamma/2)*(v_c.^2));
%Determinig the flux values
F1=U2;
term1=(U2.^2)./U1;
term2=U3-((gamma/2)*term1);
F2=term1+((gamma-1)/gamma)*term2;
term3=(U2.*U3)./U1;
term4=(U2.^3)./(U1.^2);
F3=gamma*term3-((gamma*(gamma-1))/2)*term4;
%Corrector method
for j=2:n-1
dU1dt_c(j) = -(F1(j)-F1(j-1))/dx;
dU2dt_c(j) = (-(F2(j)-F2(j-1))/dx)+(1/gamma)*(rho_c(j)*t_c(j))*((a(j)-a(j-1))/dx);
dU3dt_c(j) = -(F3(j)-F3(j-1))/dx;
end
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 l=2:n-1
U1(l)=U1old(l)+dU1dt(l)*dt;
U2(l)=U2old(l)+dU2dt(l)*dt;
U3(l)=U3old(l)+dU3dt(l)*dt;
end
%Inlet Boundary condition
U1(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 Boundary condition
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);
%Determining the primitive variables
rho_c=U1./a;
v_c=U2./U1;
t_c=(gamma-1)*((U3./U1)-(gamma/2)*(v_c.^2));
M_c=v_c./(t_c.^(1/2)); %Mach number
P_c=rho_c.*t_c; %nondimensional pressure
figure(2)
hold on
if k==100
plot(x,U2,\'b\')
else if k==200
plot(x,U2,\'g\')
else if k==700
plot(x,U2,\'k\')
title(\'Distributions of the dimensionless mass flow rate across the length of the nozzle for conservation form\')
legend(\'100dt\',\'200dt\',\'700dt\')
xlabel(\'x/L\')
ylabel(\'Dimensionless mass flow rate\')
end
end
end
end
return
endMatlab function for non-conservative form:
function [v_nc,rho_nc,t_nc,m_nc,M_nc,P_nc]=non_conservative(a,gamma,n,dx,nt,x,C)
%Conservative
%initial conditions profile
rho_nc=1-0.3146*x;
t_nc=1-0.2314*x;
v_nc=(0.1+1.09*x).*t_nc.^(0.5);
%Outer time loop
for k=1:nt
rho_old=rho_nc;
v_old=v_nc;
t_old=t_nc;
%Determinig the time step
A=t_nc.^(1/2);
dt=C*(dx./(A+v_nc));
dt=min(dt);
%Predictor method
for j=2:n-1
dvdx = (v_nc(j+1)-v_nc(j))/dx;
drhodx = (rho_nc(j+1)-rho_nc(j))/dx;
dtdx = (t_nc(j+1)-t_nc(j))/dx;
dlogadx = (log(a(j+1))-log(a(j)))/dx;
%continuity equation
drhodt_p(j) = -rho_nc(j)*dvdx-rho_nc(j)*v_nc(j)*dlogadx-v_nc(j)*drhodx;
%momentum equation
dvdt_p(j) = -v_nc(j)*dvdx-(1/gamma)*(dtdx+(t_nc(j)/rho_nc(j))*drhodx);
%energy equation
dtdt_p(j) = -v_nc(j)*dtdx-(gamma-1)*t_nc(j)*(dvdx+v_nc(j)*dlogadx);
%solution update
v_nc(j) = v_nc(j)+dvdt_p(j)*dt;
rho_nc(j) = rho_nc(j)+drhodt_p(j)*dt;
t_nc(j) = t_nc(j)+dtdt_p(j)*dt;
end
%Corrector method
for j=2:n-1
dvdx = (v_nc(j)-v_nc(j-1))/dx;
drhodx = (rho_nc(j)-rho_nc(j-1))/dx;
dtdx = (t_nc(j)-t_nc(j-1))/dx;
dlogadx = (log(a(j))-log(a(j-1)))/dx;
%continuity equation
drhodt_c(j) = -rho_nc(j)*dvdx-rho_nc(j)*v_nc(j)*dlogadx-v_nc(j)*drhodx;
%momentum equation
dvdt_c(j) = -v_nc(j)*dvdx-(1/gamma)*(dtdx+(t_nc(j)/rho_nc(j))*drhodx);
%energy equation
dtdt_c(j) = -v_nc(j)*dtdx-(gamma-1)*t_nc(j)*(dvdx+v_nc(j)*dlogadx);
end
% Average time derivative
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
v_nc(i)=v_old(i)+dvdt(i)*dt;
rho_nc(i)=rho_old(i)+drhodt(i)*dt;
t_nc(i)=t_old(i)+dtdt(i)*dt;
end
% Boundary condition
%inlet
v_nc(1)=2*v_nc(2)-v_nc(3);
%outlet
v_nc(n)=2*v_nc(n-1)-v_nc(n-2);
rho_nc(n)=2*rho_nc(n-1)-rho_nc(n-2);
t_nc(n)=2*t_nc(n-1)-t_nc(n-2);
m_nc=rho_nc.*a.*v_nc; %nondimensional mass flow rate
M_nc=v_nc./(t_nc.^(1/2)); %Mach number
P_nc=rho_nc.*t_nc; %nondimensional pressure
figure(1)
hold on
if k==50
plot(x,m_nc,\'r\')
else if k==100
plot(x,m_nc,\'b\')
else if k==150
plot(x,m_nc,\'g\')
else if k==200
plot(x,m_nc,\'y\')
else if k==700
plot(x,m_nc,\'k\')
title(\'Distributions of the dimensionless mass flow rate across the length of the nozzle for non conservation form\')
legend(\'50dt\',\'100dt\',\'150dt\',\'200dt\',\'700dt\')
xlabel(\'x/L\')
ylabel(\'Dimensionless mass flow rate\')
end
end
end
end
end
end
return
endResults:
1)
The above graph was attained after the steady-state reached by solving the conservative form of an equation. It was noticed that the Mach number increases across the nozzle and pressure, temperature, density ratio decreases across the nozzle.
2)
The above graph was attained after the steady-state reached by solving the non-conservative form of an equation. It was noticed that the Mach number increases across the nozzle and pressure, temperature, density ratio decreases across the nozzle.
3)
Here, the nondimensional mass flow rate is plotted as a function of nondimensional distance through the nozzle. Each curve represents the mass flow rate across the length of the nozzle at different time steps(100dt,200dt,700dt). Dimensionless mass flow rate attained convergence after the 700th time step.
4)
Here, the nondimensional mass flow rate is plotted as a function of nondimensional distance through the nozzle. Each curve represents the mass flow rate across the length of the nozzle at different time steps(50dt,100dt,150dt,200dt,700dt). Dimensionless mass flow rate attained convergence after the 700th time step.
5)
The conservation form gives a distribution that is much closer to being a constant. On another hand, the nonconservation results have some oscillations at both inflow and outflow boundaries. But practically, the mass flow rate appears to be constant. The nonconservation results much closer to the exact solution at the throat.
6) grid dependence test:
The grid independence is a serious consideration in CFD. When we solve the problem using CFD, there is a finite number of grid points distributed over the flow field. the solution can be attained at that finite grid points. However, the solution may not be attained at a different number of grid points. So, it is necessary to perform the grid dependence test.
n=31
n=61
n=91
From the above tables, it was noted that the solution varies slightly while increasing the number of grid points for both conservation and nonconservation solution.
7) Time taken for the simulation to run 1000 time steps:
Conservation method =2.63 sec
Nonconcervation method = 3.1268 sec
Hence, the conservation method is faster than nonconservation method