Objective:
To simulagte the isentropic flow through a qasi 1D subsonic-supersonic nozzle using both the conservation and non-conservation forms of the
governing equations and solving them the MacCormack's technique and comparing there solutions. 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.
Problem statement:
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
Description:
Consider the study, isentropic flow through a convergent divergent nozzle. The flow at the inlet to the nozzle comes from reserviour The cross
section area of the reservoir is very large and hence velocity is very small. The flow is locally subsonic the convergent section of thenozzle,
Sonic at the throat and Supersonic at the divergent section.
Consider a control voluem of the nozzle as shown :
The flow properties are uniform across the section. Hence, although the area of the nozzle changes as a function of distance along the nozzle,
X.The cross section area is 'A' ,and therefore in reality the flow field is two dimensional(the flow varies in the two dimensional x-y space), we
make the assumptions that the flow properties vary only with X; this is tantamount to asuming uniform flow properties across any given
cross-section. Such flow is defined as Quasi one dimensional flow.
The Governing Equations continuty, momentum and energy equations in terms of non dimensional variables for this quassi one dimensional
staedy and isentropc flow can be expressed respectively as:
Non-Conservation form:
Continuity Equation:
`(del rho')/(del t') = -rho' (del V')/(delx') -rho' V'((del(lnA'))/(delx'))-V'((del rho')/(del x'))`
Momentum Equation:
`(del V')/(del t') = -V' ((del V')/(delx')) -(1/ gamma) ((delT')/(delx')+(T')/(rho')((del rho')/(del x')))`
Energy Equation :
`(del T')/(del t') = -V (del T')/(del x') - (gamma-1)T'[(delV')/(delx')+V' (del(lnA'))/(delx')]`
Conservation form :
Continuity Equation:
`(del (rho' A'))/(del t') +(del( rho' A'V'))/(delx') =0`
Momentum Equation :
`(del(rho' A' V'))/(del t') + (del[ rho' A' (V')^2 +(1/gamma) rho'A'])/ (delx') = 1/gamma rho'(del A')/(del x')`
Energy Equation :
`(del [rho' ((e')/(gamma-1) + ((gamma)/(2))(V')^2)A'])/(del t') + (del[rho' ((e')/(gamma-1) + ((gamma)/(2))(V')^2)A'+rho'A'V'])/(delx')=0`
Where
gamma is air standard `gamma = 1.4`
These above equations are solved by Macormack method using time marching process
Macormack Method:
Let f(x) is the equation to be solved.Then the following steps are followed
i. inital profile
Before the calculation, the initial profile for all the flow field variables which are velocity, density and temperature has to be determined.
These are guess values but should be specified with more awareness on the problem.
ii. predector step
`((delf)/(delt))_p = (f(i+1)-f(i))/(Deltax)`
iii. solution update
`f_p = f_(Initial)+((delf)/(delt))_p*Deltat`
`f_p =` Flow field variables calculated after the predictor step
`f_(Initial) =`flow field variable in the initial step
`Deltat =`time step size
iv. Corrector step
`((delf_p)/(delt))_c = (f(i)-f(i-1))/(Deltax)`
`((delf)/(delt))_(avg)= 0.5(((delf)/(delt))_(p)-((delf_p)/(delt))_(c))`
v. Final solution update
`f_(sol) = f_(Initial) +((delf)/(delt))_(avg)*Deltat`
`f_(sol) =` Final solution of the flow field variables for the first time step
Matlab Coding:
clear all
close all
clc
%SOLVING QUASI 1D SUPERSONIC NOZZLE FLOW EQUATIONS
%INPUT PARAMETERS
n = 31; %No. of grid points
L = 3; %length of nozzle
gamma = 1.4;
x = linspace(0,3,n);
dx = x(2) - x(1);
nt = 1400; %number of time steps
C = 0.5; %Courant Number
%FUNCTION FOR NON CONSERVATIVE FORM
tic;
[Mass_flow_rate_NC, Pressure_NC, Mach_number_NC, rho_NC, v_NC, T_NC, rho_throat_NC, v_throat_NC, T_throat_NC, Mass_flow_rate_throat_NC, Pressure_throat_NC, Mach_number_throat_NC] = Non_conservative_form(x,dx, n, gamma, nt, C);
Time_elapsed_NC = toc;
fprintf('Time Elasped for Non Conservative form: %0.4f', Time_elapsed_NC);
hold on
%PLOTTING
%Timestep variation of properties at throat area in Non Conservative form
figure(2)
%Density
subplot(4,1,1)
plot(linspace(1,nt,nt),rho_throat_NC,'color','m')
ylabel('Density')
legend({'Density at throat'});
axis([0 1400 0 1.3]);
grid minor;
title('Timestep variation at throat area - NON CONSERVATIVE FORM')
%Temperature
subplot(4,1,2)
plot(linspace(1,nt,nt),T_throat_NC,'color','c')
ylabel('Temperature')
legend({'Temperature at throat'});
axis([0 1400 0.6 1]);
grid minor;
%Pressure
subplot(4,1,3)
plot(linspace(1,nt,nt),Pressure_throat_NC,'color','g')
ylabel('Pressure')
legend({'Pressure at throat'});
axis([0 1400 0.4 1.1]);
grid minor
%Mach Number
subplot(4,1,4)
plot(linspace(1,nt,nt),Mach_number_throat_NC,'color','y')
xlabel('Timesteps')
ylabel('Mach Number')
legend({'Mach Number at throat'});
axis([0 1400 0.6 1.4]);
grid minor;
%Steady state simulation for primitve values for Non Conservative form
figure(3);
%Density
subplot(4,1,1);
plot(x,rho_NC,'color','m')
ylabel('Density')
legend({'Density'});
axis([0 3 0 1])
grid minor;
title('Variation of Flow Rate Distribution - NON CONSERVATIVE')
%Temperature
subplot(4,1,2);
plot(x,T_NC,'color','c')
ylabel('Temperature')
legend({'Temperature'});
axis([0 3 0 1])
grid minor;
%Pressure
subplot(4,1,3);
plot(x,Pressure_NC,'color','g')
ylabel('Pressure')
legend({'Pressure'});
axis([0 3 0 1])
grid minor;
%Mach number
subplot(4,1,4);
plot(x,Mach_number_NC,'color','y')
xlabel('Nozzle Length')
ylabel('Mach number')
legend({'Mach number'});
axis([0 3 0 4])
grid minor;
%FUNCTION FOR CONSERVATIVE FORM
tic;
[Mass_flow_rate_C, Pressure_C, Mach_number_C, rho_C, v_C, T_C, rho_throat_C, v_throat_C, T_throat_C, Mass_flow_rate_throat_C, Pressure_throat_C, Mach_number_throat_C] = Conservative_form(x, dx, n, gamma, nt, C);
Time_elasped_C = toc;
fprintf('Time Elasped for Conservative form: %0.4f', Time_elasped_C);
%PLOTTING
%Timestep variation of properties at throat area in Conservative form
figure(5)
%Density
subplot(4,1,1)
plot(rho_throat_C,'color','m')
ylabel('Density')
legend({'Density at throat'});
axis([0 1400 0 1]);
grid minor;
title('Timestep variation at throat area - CONSERVATIVE FORM')
%Temperature
subplot(4,1,2)
plot(T_throat_C,'color','c')
ylabel('Temperature')
legend({'Temperature at throat'});
axis([0 1400 0.6 1]);
grid minor;
%Pressure
subplot(4,1,3)
plot(Pressure_throat_C,'color','g')
ylabel('Pressure')
legend({'Pressure at throat'});
axis([0 1400 0.2 0.8]);
grid minor
%Mach Number
subplot(4,1,4)
plot(Mach_number_throat_C,'color','y')
xlabel('Timesteps')
ylabel('Mach Number')
legend({'Mach Number at throat'});
axis([0 1400 0.6 1.4]);
grid minor;
%Steady state simulation for primitve values for Conservative form
figure(6);
%Density
subplot(4,1,1);
plot(x,rho_C,'color','m')
ylabel('Density')
legend({'Density'});
axis([0 3 0 1])
grid minor;
title('Variation of Flow Rate Distribution - CONSERVATIVE FORM')
%Temperature
subplot(4,1,2);
plot(x,T_C,'color','c')
ylabel('Temperature')
legend({'Temperature'});
axis([0 3 0 1])
grid minor;
%Pressure
subplot(4,1,3);
plot(x,Pressure_C,'color','g')
ylabel('Pressure')
legend({'Pressure'});
axis([0 3 0 1])
grid minor;
%Mach number
subplot(4,1,4);
plot(x,Mach_number_C,'color','y')
xlabel('Nozzle Length')
ylabel('Mach number')
legend({'Mach number'});
axis([0 3 0 4])
grid minor;
%COMPARISON PLOT of Normalized Mass Flow Rate of both forms
figure(7)
hold on
plot(x,Mass_flow_rate_NC,'r','linewidth',1.25)
hold on;
plot(x,Mass_flow_rate_C,'b','linewidth',1.25)
hold on;
grid on
legend('Non Conservative','Conservative');
xlabel('Length of the Domain')
ylabel('Mass Flow Rate')
title('Comparison of Normalized Mass Flow Rate of both forms')
Non-Conservation form :
%Function for Non-Conservative form
function [Mass_flow_rate_NC, Pressure_NC, Mach_number_NC, rho_NC, v_NC, T_NC, rho_throat_NC, v_throat_NC, T_throat_NC, Mass_flow_rate_throat_NC, Pressure_throat_NC, Mach_number_throat_NC] = Non_conservative_form(~,~, n, ~, nt, C)
% Inputs
n = 31; %No. of grid points
gamma = 1.4;
x = linspace(0,3,n);
dx = x(2) - x(1);
nt = 1400; %number of time steps
C = 0.5; %Courant Number
%Intial boundary conditions
rho_NC = 1 - 0.3146*x;
T_NC = 1 - 0.2341*x;
v_NC = (0.1 + 1.09*x).*T_NC.^0.5;
%Area & Throat profiles
a = 1 + 2.2*(x - 1.5).^2;
throat = find(a==1);
%Outer Time Loop
for k = 1:nt
%Time-step control
dt = min(C.*dx./(sqrt(T_NC) + v_NC));
% copy of Old Values
rho_old = rho_NC;
T_old = T_NC;
v_old = v_NC;
%Predictor Method
for i = 2:n-1
%Simplifying the Equation's spatial terms
dv_dx = (v_NC(i + 1) - v_NC(i))/dx;
drho_dx = (rho_NC(i + 1) - rho_NC(i))/dx;
dT_dx = (T_NC(i + 1) - T_NC(i))/dx;
dlog_a_dx = (log(a(i + 1)) - log(a(i)))/dx;
%Continuity Equation
drhodt_P(i) = -rho_NC(i)*dv_dx - rho_NC(i)*v_NC(i)*dlog_a_dx - v_NC(i)*drho_dx;
%Momentum Equation
dvdt_P(i) = -v_NC(i)*dv_dx - (1/gamma)*(dT_dx + (T_NC(i)*drho_dx/rho_NC(i)));
%Energy Equation
dTdt_P(i) = -v_NC(i)*dT_dx - (gamma - 1)*T_NC(i)*(dv_dx + v_NC(i)*dlog_a_dx);
%Solution Update
rho_NC(i) = rho_NC(i) + drhodt_P(i)*dt;
v_NC(i) = v_NC(i) + dvdt_P(i)*dt;
T_NC(i) = T_NC(i) + dTdt_P(i)*dt;
end
%Corrector Method
for j = 2:n-1
%Simplifying the Equation's spatial terms
dv_dx = (v_NC(j) - v_NC(j - 1))/dx;
drho_dx = (rho_NC(j) - rho_NC(j - 1))/dx;
dT_dx = (T_NC(j) - T_NC(j - 1))/dx;
dlog_a_dx = (log(a(j)) - log(a(j - 1)))/dx;
%Continuity Equation
drhodt_C(j) = - rho_NC(j)*dv_dx - rho_NC(j)*v_NC(j)*dlog_a_dx - v_NC(j)*drho_dx;
%Momentum Equation
dvdt_C(j) = - v_NC(j)*dv_dx - (1/gamma)*(dT_dx + (T_NC(j)*drho_dx/rho_NC(j)));
%Energy Equation
dTdt_C(j) = - v_NC(j)*dT_dx - (gamma - 1)*T_NC(j)*(dv_dx + v_NC(j)*dlog_a_dx);
end
% Average Time derivative
drho_dt = 0.5.*(drhodt_P + drhodt_C);
dv_dt = 0.5.*(dvdt_P + dvdt_C);
dT_dt = 0.5.*(dTdt_P + dTdt_C);
% Solution Update
for m = 2:n-1
rho_NC(m) = rho_old(m) + drho_dt(m)*dt;
v_NC(m) = v_old(m) + dv_dt(m)*dt;
T_NC(m) = T_old(m) + dT_dt(m)*dt;
end
%Applying Boundary Conditions
%At Inlet
rho_NC(1) = 1;
v_NC(1) = 2*v_NC(2) - v_NC(3);
T_NC(1) = 1;
%At Outlet
rho_NC(n) = 2*rho_NC(n-1) - rho_NC(n - 2);
v_NC(n) = 2*v_NC(n - 1) - v_NC(n - 2);
T_NC(n) = 2*T_NC(n - 1) - T_NC(n - 2);
%Defining Mass flow rate, Pressure and Mach number
Mass_flow_rate_NC = rho_NC.*a.*v_NC;
Pressure_NC = rho_NC.*T_NC;
Mach_number_NC = (v_NC./sqrt(T_NC));
%Calculating Variables at throat section
rho_throat_NC(k) = rho_NC(throat);
T_throat_NC(k) = T_NC(throat);
v_throat_NC(k) = v_NC(throat);
Mass_flow_rate_throat_NC(k) = Mass_flow_rate_NC(throat);
Pressure_throat_NC(k) = Pressure_NC(throat);
Mach_number_throat_NC(k) = Mach_number_NC(throat);
%PLOTTING
%Comparison of non-dimensional Mass flow rates at different time steps
figure(1);
if k == 50
plot(x, Mass_flow_rate_NC, 'm', 'linewidth', 1.25);
hold on
elseif k == 100
plot(x, Mass_flow_rate_NC, 'c', 'linewidth', 1.25);
hold on
elseif k == 200
plot(x, Mass_flow_rate_NC, 'g', 'linewidth', 1.25);
hold on
elseif k == 400
plot(x, Mass_flow_rate_NC, 'y', 'linewidth', 1.25);
hold on
elseif k == 800
plot(x, Mass_flow_rate_NC, 'k', 'linewidth', 1.25);
hold on
end
%Labelling
title('Mass flow rates at different Timesteps - NON CONSERVATIVE FORM)')
xlabel('Distance')
ylabel('Mass flow rate')
legend({'50^t^h Timestep','100^t^h Timestep','200^t^h Timestep','400^t^h Timestep','800^t^h Timestep'})
axis([0 3 0 2])
grid on
end
end
Conservation form:
%Function for Conservative form
function [Mass_flow_rate_C, Pressure_C, Mach_number_C, rho_C, v_C, T_C, rho_throat_C, v_throat_C, T_throat_C, Mass_flow_rate_throat_C, Pressure_throat_C, Mach_number_throat_C] = Conservative_form(x, dx, n, gamma, nt, C)
%Determining Initial profile of Density and Temperature
for i = 1:n
if (x(i) >= 0 && x(i) <= 0.5)
rho_C(i) = 1;
T_C(i) = 1;
elseif (x(i) >= 0.5 && x(i) <= 1.5)
rho_C(i) = 1 - 0.366*(x(i) - 0.5);
T_C(i) = 1 - 0.167*(x(i) - 0.5);
elseif (x(i) >= 1.5 && x(i) <= 3)
rho_C(i) = 0.634 - 0.3879*(x(i) - 1.5);
T_C(i) = 0.833 - 0.3507*(x(i) - 1.5);
end
end
%Area & Throat profiles
a = 1 + 2.2*(x - 1.5).^2;
throat = find(a==1);
%Initial profile fot Velocity and Pressure
v_C = 0.59./(rho_C.*a);
Pressure_C = rho_C.*T_C;
%Calculate the Initial Solution Vectors(U)
U1 = rho_C.*a;
U2 = rho_C.*a.*v_C;
U3 = rho_C.*a.*(T_C./(gamma - 1)) + ((gamma/2).*v_C.^2);
%Outer Time Loop
for k =1:nt
%Time-step control
dt = min(C.*dx./(sqrt(T_C) + v_C));
%Saving a copy of Old Values
U1_old = U1;
U2_old = U2;
U3_old = U3;
%Calculating the Flux vector(F)
F1 = U2;
F2 = ((U2.^2)./U1) + ((gamma - 1)/gamma)*(U3 - (gamma/2)*((U2.^2)./U1));
F3 = ((gamma*U2.*U3)./U1) - (gamma*(gamma - 1)*((U2.^3)./(2*U1.^2)));
%Predictor Method
for j = 2:n-1
%Calculating the Source term(J)
J2(j) = (1/gamma)*rho_C(j)*T_C(j)*((a(j + 1) - a(j))/dx);
%Continuity Equation
dU1_dt_P(j) = -((F1(j + 1) - F1(j))/dx);
%Momentum Equation
dU2_dt_P(j) = -((F2(j + 1) - F2(j))/dx) + J2(j);
%Energy Equation
dU3_dt_P(j) = -((F3(j + 1) - F3(j))/dx);
%Updating New values
U1(j) = U1(j) + dU1_dt_P(j)*dt;
U2(j) = U2(j) + dU2_dt_P(j)*dt;
U3(j) = U3(j) + dU3_dt_P(j)*dt;
end
%Updating Flux vectors
F1 = U2;
F2 = ((U2.^2)./U1) + ((gamma - 1)/gamma)*(U3 - ((gamma*0.5*(U2.^2))./(U1)));
F3 = ((gamma*U2.*U3)./U1) - ((gamma*0.5*(gamma - 1)*(U2.^3))./(U1.^2));
%Corrector Method
for m = 2:n-1
%Updating the Source term(J)
J2(m) = (1/gamma)*rho_C(m)*T_C(m)*((a(m) - a(m - 1))/dx);
%Simplifying the Equation's spatial terms
dU1_dt_C(m) = -((F1(m) - F1(m - 1))/dx);
dU2_dt_C(m) = -((F2(m) - F2(m - 1))/dx) + J2(m);
dU3_dt_C(m) = -((F3(m) - F3(m - 1))/dx);
end
%Compute Average Time derivative
dU1_dt = 0.5*(dU1_dt_P + dU1_dt_C);
dU2_dt = 0.5*(dU2_dt_P + dU2_dt_C);
dU3_dt = 0.5*(dU3_dt_P + dU3_dt_C);
%Final Solution Update
for i = 2:n-1
U1(i) = U1_old(i) + dU1_dt(i)*dt;
U2(i) = U2_old(i) + dU2_dt(i)*dt;
U3(i) = U3_old(i) + dU3_dt(i)*dt;
end
% Boundary Conditions
%At Inlet
U1(1) = rho_C(1)*a(1);
U2(1) = 2*U2(2) - U2(3);
v_C(1) = U2(1)./U1(1);
U3(1) = U1(1)*((T_C(1)/(gamma - 1)) + ((gamma/2)*(v_C(1)).^2));
%At 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);
%Final Primitive value update
rho_C = U1./a;
v_C = U2./U1;
T_C = (gamma - 1)*((U3./U1 - (gamma/2)*(v_C).^2));
%Defining Mass flow rate, Pressure and Mach number
Mass_flow_rate_C = rho_C.*a.*v_C;
Pressure_C = rho_C.*T_C;
Mach_number_C = (v_C./sqrt(T_C));
%Calculating Variables at throat section
rho_throat_C(k) = rho_C(throat);
T_throat_C(k) = T_C(throat);
v_throat_C(k) = v_C(throat);
Mass_flow_rate_throat_C(k) = Mass_flow_rate_C(throat);
Pressure_throat_C(k) = Pressure_C(throat);
Mach_number_throat_C(k) = Mach_number_C(throat);
%PLOTTING
%Comparison of non-dimensional Mass flow rates at different time steps
figure(4)
if k == 50
plot(x, Mass_flow_rate_C, 'm', 'linewidth', 1.25);
hold on
elseif k == 100
plot(x, Mass_flow_rate_C, 'c', 'linewidth', 1.25);
hold on
elseif k == 200
plot(x, Mass_flow_rate_C, 'g', 'linewidth', 1.25);
hold on
elseif k == 400
plot(x, Mass_flow_rate_C, 'y', 'linewidth', 1.25);
hold on
elseif k == 800
plot(x, Mass_flow_rate_C, 'k', 'linewidth', 1.25);
hold on
end
%Labelling
title('Mass flow rates at different Timesteps - CONSERVATIVE FORM')
xlabel('Distance')
ylabel('Mass flow rate')
legend({'50^th Timestep','100^th Timestep','200^th Timestep','400^th Timestep','800^th Timestep'})
axis([0 3 0.4 1])
grid on
end
end
Plots :
From the above plot the Non-conservative of mass flow rate is stable from 400th to 800th time step .
The Non-conservative time step variation at throat area.From the graph it is showing that
flow field variables after 400 iterations reached to steady state sloution.
The above graph shows that the variation of flow field varialbles along the length of the nozzle
for the non conservative form. In the above graph it is clear that the flow field variables such as
density, pressure, temperature decreases and mach number increases with a increase
in length of the nozzle.
From the above plot in Conservative of mass flow rate is stable from 400th to 800th timestep.
The Conservative time step variation at throat area.From the graph it is showing that
flow field variables after 400 iterations reached to steady state sloution.
The above graph shows that the variation of flow field varialbles along the length of the nozzle
for the Conservative form. In the above graph it is clear that the flow field variables such as
density, pressure, temperature decreases and mach number increases with a increase
in length of the nozzle.
Mass flow rate at Non-conserevative & conservative form.The above graph shows that the mass flow
rate of conservative forms produces a stable result when compare to that of non conservative form.
Mass flow rate is minimum at the throat then it expands in the non conservative form.
Conclusion :
Simulation of isentropic flow through a Quasi 1D subsonic-supersonic nozzle by deriving both Non-conservative and
conservative forms of governing equation and solve them by using Macormack Method in matlab and results