Week 7 - Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
Abstract
In this challenge, it is going to simulate the isentropic flow through a quasi 1D subsonic-supersonic nozzle.it will derive both the conservation and non-conservation forms of the governing equations and solve them using the MacCormack's technique. form this need to determine the steady-state temperature distribution for the flow-field variables and investigate the difference between the two forms of governing equations by comparing their solutions.
from this experiment, we will get the following plots inside your report
1. Steady-state distribution of primitive variables inside the nozzle
2. Time-wise variation of the primitive variables
3. Variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process
4. Comparison of Normalized mass flow rate distributions of both forms.
Introduction-:
The nozzle is used to convert the chemical thermal energy generated in the combustion chamber into kinetic energy. The nozzle converts the low velocity, high velocity, high pressure, high-temperature gas in the combustion chamber into high-velocity gas of lower pressure and temperature. The one dimensional in isentropic flow in a convergent-divergent nozzle.
Fig 1:Basic diagram of 1-D-Supersonic Nozzle Flow.
The convergent-divergent nozzle is designed for attaining speeds that are greater than the speed of sound.
The design of this nozzle came from the area-velocity relation
(dA/dV)= - (A/V)(1-M2)
Where,
M = Mach number (which means the ratio of the local speed of flow to the local speed of sound)
A = area of the nozzle.
V= velocity
P= Pressure
MACORMACK METHOD-:
This method is employed to solve the governing equations of quasi 1-D-Supersonic nozzle flow rather than the differencing schemes as we have studied earlier.Macormack method is an explicit finite difference technique
Which is second-order accurate in both time and space. In computational fluid
Dynamics, the macro mask method is a widely used discretization schema for the numerical solution of the hyperbolic partial differential equation. This method essentially uses a predictor step and corrector steps which will be depicted in the code.
According to the macromack method, the values of a variable at the grid point i,j at the time level
Where the average time derivatives are the mean value between time ‘t‘ and ’t+dt‘. The average time derivatives are obtained from the predictor-corrector philosophy.
Predictor Step:
In the predictor step, the spatial terms on the right-hand side of the equation are replaced by forwarding differences. The RHS of the equation is known as all values of variables that are for time steps ’t‘.Next, the values of the variables on the LHS are calculated by employing the time marching solution and for the time step,’ the t+dt‘ value of the variables at all grid points is obtained.
Corrector steps :
In the predictor step, the spatial terms on the right-hand side of the equation are replaced by backward differences. The values of variables on RHS are taken are obtained in the preditor step for the time's steps t+dt. Then the corrected time derivative is calculated.
Now the predicted and corrected value of time derivative is used to calculate the average time derivative which in turn gives the value of the variable at a given point for a time steps. The same steps are repeated until the values of variables are obtained for all grid points for a given number of time steps. In predictor steps, we are using forward differencing method for space terms.
GOVERNING EQUATIONS:
There are 3 equations that govern the flow of the fluid in the nozzle continuity equation momentum equation and energy equation.
Non-Conservation Equation-:
In this type, a control volume is taken into consideration, where the control volume moves with the fluid, and properties of the fluid are taken into consideration.
Conservative Equation-:
In this type, a control volume is taken into consideration where the control volume is fixed and the flow of the fluid in that control volume is studied.
Programming for Non-conservative form and conservative form-
Non-conservative form-:
clear all
close all
clc
n=31;
x=linspace(0,3,n);
gamma=1.4;
dx=x(2)-x(1);
error=1;
%INITIALIZE THROAT
throat=(n+1)/2;
%INITIAL CONDITION
%DENSITY
rho=1-0.3146*x;
%Temprature
t=1-0.23148*x;
%velocity
v=(0.1+1.09*x).*t.^0.5;
%area
a=1+2.2*(x-1.5).^2;
nt=900;
dt=0.01;
error_tol=1e-4;
for k=1:nt
rho_initial=rho;
v_initial=v;
t_initial=t;
%predictor step
for i=2:n-1
dvdx=(v(i+1)-v(i))/dx;
drhodx=(rho(i+1)-rho(i))/dx;
dtdx=(t(i+1)-t(i))/dx;
dlogadx=(log(a(i+1))-log(a(i)))/dx;
% continuity equation
drhodt_p(i)= -rho(i)*dvdx- rho(i)*v(i)*dlogadx - v(i)*drhodx;
% momentum equation
dvdt_p(i)= -v(i)*dvdx - (1/gamma)*(dtdx + (t(i)/rho(i))*drhodx);
% energy equation
dtdt_p(i)=-v(i)*dtdx - (gamma-1)*t(i)*(dvdx+(v(i)*dlogadx));
% solution update
rho(i)=rho(i)+drhodt_p(i)*dt;
v(i)=v(i)+dvdt_p(i)*dt;
t(i)=t(i)+dtdt_p(i)*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));
end
drhodt=0.5*(drhodt_p+drhodt_c);
dvdt=0.5*(dvdt_p+dvdt_c);
dtdt=0.5*(dtdt_p+dtdt_c);
% solution upadte
for i=2:n-1
rho(i)=rho_initial(i)+drhodt(i)*dt;
v(i)=v_initial(i)+dvdt(i)*dt;
t(i)=t_initial(i)+dtdt(i)*dt;
end
%boundary condition
%inlet
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);
%pressure
P=rho.*t;
% Mach number
M=(v./(t).^0.5);
% mass flow rate
m=rho.*v.*a;
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);
figure(1)
hold all
grid on
if k==1
plot(x,m,'c');
elseif k==50
plot(x,m,'r');
elseif k==100
plot(x,m,'b');
elseif k==150
plot(x,m,'y');
elseif k==200
plot(x,m,'m');
elseif k==700
plot(x,m,'g');
xlabel ('distance through nozzle')
ylabel('mass flow rate')
legend('1','50','100','150','200','700');
title('mass flow rate at different time step non conservative form')
end
end
figure(2)
subplot(2,3,1)
plot(x,rho,'r');
xlabel('length of nozzle')
ylabel('density')
title('density non conservative form')
grid minor
subplot(2,3,2)
plot(x,t,'b');
xlabel('length of nozzle')
ylabel('temperature')
title('temperature non conservative form')
grid minor
subplot(2,3,3)
plot(x,P,'g');
xlabel('length of nozzle')
ylabel('pressure')
title('pressure non conservative form')
grid minor
subplot(2,3,4)
plot(x,v,'m');
xlabel('length of nozzle')
ylabel('velocity')
title('velocity non conservative form')
grid minor
subplot(2,3,5)
plot(x,M,'y');
xlabel('length of nozzle')
ylabel('velocity')
title('velocity non conservative form')
grid minor
subplot(2,3,6)
plot(x,m,'k');
xlabel('length of nozzle')
ylabel('mass flow rate')
title('mass flow rate non conservative form')
grid minor
figure(3)
subplot(6,1,1)
plot(1:nt,rho_throat,'r')
xlabel('time steps')
ylabel('density')
title('variation of density at the nozzle throat non conservative form')
subplot(6,1,1)
plot(1:nt,rho_throat,'g')
xlabel('time steps')
ylabel('temperature')
title('variation of density at the nozzle throat non conservative form')
grid minor
subplot(6,1,2)
plot(1:nt,rho_throat,'m')
xlabel('time steps')
ylabel('temperature')
title('variation of temperature at the nozzle throat non conservative form')
grid minor
subplot(6,1,3)
plot(1:nt,P_throat,'y')
xlabel('time steps')
ylabel('Pressure')
title('variation of pressure at the nozzle throat non conservative form')
grid minor
subplot(6,1,4)
plot(1:nt,v_throat,'k')
xlabel('time steps')
ylabel('velocity')
title('variation of velocity at the nozzle throat non conservative form')
grid minor
subplot(6,1,5)
plot(1:nt,M_throat,'b')
xlabel('time steps')
ylabel('Mach number')
title('variation of mach number at the nozzle throat non conservative form')
grid minor
subplot(6,1,6)
plot(1:nt,m_throat,'c')
xlabel('time steps')
ylabel('mass_flow_rate')
title('variation of mass flow rate at the nozzle throat non conservative form')
grid minor
figure(4)
hold all
grid on
plot(x,m,'c')
plot(x,linspace (0.579,0.579,n),'r')
xlabel('Non dimensional length')
ylabel('Mass flow rate')
legend('mass flow rate in Non conservative form','analytical solution')
Conservative form-:
clear all
close all
clc
n=31;
x=linspace(0,3,n);
dx=x(2)-x(1);
gamma=1.4;
error=1;
%INITIALIZE THROAT
throat=(n+1)/2;
nt=900;
dt=0.01;
error_tol=1e-4;
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.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
%initial conditions
% density
rho=1-0.3146*x;
% temperature
t=1-0.2314*x;
% area
a=1+2.2*(x-1.5).^2;
%velocity
v=0.59./(rho.*a);
P=rho.*t;
% defining elements of solution vector U
U1=rho.*a;
U2=rho.*a.*v;
U3=rho.*a.*((t./(gamma-1))+((gamma/2).*(v.^2)));
%outer time loop
for k=1:nt
U1_old=U1;
U2_old=U2;
U3_old=U3;
%flux vector
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 methods
for i=2:n-1
dF1dx=((F1(i+1)-F1(i))/dx);
dF2dx=(F2(i+1)-F2(i))/dx;
dF3dx=((F3(i+1)-F3(i))/dx);
dadx=(a(i+1)-a(i)/dx);
% Continuity euuation
dU1dt_p(i)= -(dF1dx);
%momentum equation
J2(i)=(1/gamma).*rho(i).*t(i).*((a(i+1)-a(i))/dx);
dU2dt_p(i)=-dF2dx+J2(i);
%energy equation
dU3dt_p(i)=-dF3dx;
%solution upadte
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
% update values
P=rho.*t;
v=U2./U1;
rho=U1./a;
t=(gamma-1).*((U3./U1)-((gamma./2).*(v.^2)));
%flux vector
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
for j=2:n-1
dF1dx=((F1(j)-F1(j-1))/dx);
dF2dx=((F2(j)-F2(j-1))/dx);
dF3dx=((F3(j)-F3(j-1))/dx);
% continuity equ
% Continuity euuation
dU1dt_c(j)=-dF1dx;
%momentum equation
J2(j)=(1/gamma).*rho(j).*t(j).*((a(j)-a(j-1))/dx);
dU2dt_c(j)=-dF2dx+J2(j);
%energy equation
dU3dt_c(j)=-dF3dx;
end
%time derivative
dU1dt_average=0.5*(dU1dt_p + dU1dt_c);
dU2dt_average=0.5*(dU2dt_p + dU2dt_c);
dU3dt_average=0.5*(dU3dt_p + dU3dt_c);
% variable U1,U2,U3.
for l=2:n-1
U1(l)= U1_old(l)+dU1dt_average(l)*dt;
U2(l)= U2_old(l)+dU2dt_average(l)*dt;
U3(l)= U3_old(l)+dU3dt_average(l)*dt;
end
%BOUNDRY CONDITION
%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);
% flow filed variable
rho=U1./a;
v=U2./U1;
T=(gamma-1).*((U3./U1)-((gamma*v.^2)/2));
P=rho.*t;
%mach number
M=(v./(t).^0.5);
% mass flow rate
m=rho.*v.*a;
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);
figure(1)
hold all
grid on
if k==1
plot(x,m,'g');
elseif k==50
plot(x,m,'m');
elseif k==100
plot(x,m,'y');
elseif k==150
plot(x,m,'b');
elseif k==200
plot(x,m,'r');
elseif k==700
plot(x,m,'c');
xlabel ('distance through nozzle')
ylabel('mass flow rate')
legend('1','50','100','150','200','700');
title('mass flow rate at different time step conservative form')
end
end
figure(2)
subplot(2,3,1)
plot(x,rho,'r');
xlabel('length of nozzle')
ylabel('density')
title('density conservative form')
grid minor
subplot(2,3,2)
plot(x,t,'b');
xlabel('length of nozzle')
ylabel('temperature')
title('temperature conservative form')
grid minor
subplot(2,3,3)
plot(x,P,'g');
xlabel('length of nozzle')
ylabel('pressure')
title('pressure conservative form')
grid minor
subplot(2,3,4)
plot(x,v,'m');
xlabel('length of nozzle')
ylabel('velocity')
title('velocity conservative form')
grid minor
subplot(2,3,5)
plot(x,M,'y');
xlabel('length of nozzle')
ylabel('velocity')
title('velocity conservative form')
grid minor
subplot(2,3,6)
plot(x,m,'k');
xlabel('length of nozzle')
ylabel('mass flow rate')
title('mass flow rate conservative form')
grid minor
figure(3)
subplot(6,1,1)
plot(1:nt,rho_throat,'r')
xlabel('time steps')
ylabel('density')
title('variation of density at the nozzle throat conservative form')
subplot(6,1,1)
plot(1:nt,rho_throat,'g')
xlabel('time steps')
ylabel('temperature')
title('variation of density at the nozzle throat conservative form')
grid minor
subplot(6,1,2)
plot(1:nt,rho_throat,'g')
xlabel('time steps')
ylabel('temperature')
title('variation of temperature at the nozzle throat conservative form')
grid minor
subplot(6,1,3)
plot(1:nt,P_throat,'m')
xlabel('time steps')
ylabel('Pressure')
title('variation of pressure at the nozzle throat conservative form')
grid minor
subplot(6,1,4)
plot(1:nt,v_throat,'b')
xlabel('time steps')
ylabel('velocity')
title('variation of velocity at the nozzle throat conservative form')
grid minor
subplot(6,1,5)
plot(1:nt,M_throat,'k')
xlabel('time steps')
ylabel('Mach number')
title('variation of mach number at the nozzle throat conservative form')
grid minor
subplot(6,1,6)
plot(1:nt,m_throat,'y')
xlabel('time steps')
ylabel('mass flow rate')
title('variation of mass flow rate at the nozzle throat conservative form')
grid minor
figure(4)
hold all
grid on
plot(x,m,'c')
plot(x,linspace (0.579,0.579,n),'g')
xlabel('Non dimensional length')
ylabel('Mass flow rate')
legend('mass flow rate in conservative form','analytical solution')
RESULTS AND GRAPH PLOT-:
Conservative form-:
The mass flow rate of the graph for both conservative and observed that for the conservative form mass flow rate is maximum at 100th-time steps and it seems that the mass flow rate is fluctuating whereas at the 700th time steps it does not show much more deflection. other then all graphs show a wavy appearance.
Non-conservative form-:
The mass flow rate of the graph for non-conservative and it observed that for the non-conservative form mass flow rate is maximum at 50th and 100th-time steps and it seems that the mass flow rate is fluctuating whereas at the 700th time steps it does not show much more deflection. other then all graphs show a wavy appearance.
Variation of the length of the nozzle with respect to various functions for both forms-:
Non-conservative form-:
Conservative form-:
Variation of flow properties with respect to time steps at throat-:
Conservative form-:
Non-conservative form-:
it seems that the variation is much less for the conservative form as compared to the non-conservative form. at first steps, there is very little depreciation but they transitorily come in a stable position.
Distinguishing between the flow rate of non-conservative form and conservative form with respect to the analytical solution.
Non-Conservative form-:
Conservative form-:
For the mass flow rate of the non-conservative form, it seems that at first half its show that stipulated at lower in the graph but slowly increasing in the amount and reach maximum position.it shows the less disturbance between the inlet and outlet boundaries. In conservative form, it shows wavy scenarios between the inlet and outlet boundaries.
Conclusion-:
- The 1D supersonic nozzle is simulated with the maccormacks technique with conservative form and non-conservative form.
- Grid independence is established with n=31 and n=61.
- accuracy of obtaining a steady-state mass flow rate in the nozzle after iteration is checked for both forms.
- accuracy is better with the conservative form.
- The conservation form yields a better mass flow distribution. The conservation form simply does a better job of conserving mass.