Week 7 - Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
- In this challenge you are going to simulate the isentropic flow through a quasi 1D subsonic-supersonic nozzle. You will derive both the conservation and non-conservation forms of the governing equations and sovle them using the MacCormack's technique. You 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.
- Objective- To simulate the isentropic flow through a quasi 1D subsonic-supersonic nozzle.
- Derive both the conservation and non-conservation forms of the governing equations and sovle them using the MacCormack's technique.
- 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.
-
A nozzle is a device designed to control the direction or characteristics of fluid flow (especially to increase velocity) as it exits (or enters) a closed chamber or pipe.
-
A nozzle is often a pipe or tube with a variable cross-sectional area and can be used to direct or modify the flow of a fluid (liquid or gas). Nozzles are often used to control the flow rate, speed, direction, weight, shape, and/or pressure of the stream emerging from them. In the nozzle, the velocity of the fluid increases at the expense of its pressure energy.
- Often the goal of the nozzle is to increase the kinetic energy of the flowing medium at the expense of its pressure and internal energy.
- Isentropic flow occurs when the change in flow variables is small and gradual, such as ideal jet flow.
- In computational fluid dynamics, the MacCormack method is a widely used discretization scheme for the numerical solution of hyperbolic partial differential equations
-
MacCormack's method is a variation of the two-step Lax-Wendroff scheme, but is much simpler to apply.
The application of the MacCormack method takes place in two steps; a predictor step followed by a corrector step.
- It is an explicit finite difference technique.
- The solution can be obtained through three steps:
- Predictor step- In this step an itermediate solution with initial conditions and solution gradient can be canculated.
- Corrector step- In this step , corrected solution with predicted values and solution gradient can be obtained.
- Final solution- In this step, averaged solution gradient can be obtained.
- For convergent divergent nozzle ,the inlet condition would be in the subsonic region.
- Outlet condition would be in the supersonic region.
-
Non Conservative Form:
Conservative Form:
-
The main program of the subsonic-supersonic jet Quasi 1D
Specifically, two functions
1.Non-conservative form
2.Conservative form.
Since the problem setting is 1D, we only consider the x-axis distribution. Here the nozzle can be divided into two partsConvergent section
Divergent section
- Grid Independence test-
- To check whether the solution is dependent on the mesh size, we perform a mesh independence test.
- It shows us the limit to which we need to refine our mesh to get accurate results, as further refinement of the mesh would only increase the computation time.
- When grids are changed,we need to ensure grid independence is applied.
- Code-
clear all
close all
clc
%Inputs
n = 41
C = 0.5
x = linspace(0,3,n)
dx = x(2)- x(1)
gamma = 1.4;
nt = 1400;
t = 20
%fUNCTION FOR NON CONSERVATIVE FORM
tic;
% [Mass_flow_rate_C, 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);
[Mass_flow_rate_throat_NC, Pressure_throat_NC, Mach_number_throat_NC,Mass_flow_rate_NC, Pressure_NC, Mach_number_NC, rho_NC, v_NC, T_NC, rho_throat_NC, v_throat_NC, T_throat_NC ] = Non_conservative_form(n, gamma, nt,x,dx,C)
Time_elapsed_NC = toc;
fprintf('Time Elasped for Non Conservative form: %0.4f', Time_elapsed_NC);
hold on
%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);
[Mass_flow_rate_throat_C, Pressure_throat_C, Mach_number_throat_C,Mass_flow_rate_C, Pressure_C, Mach_number_C, rho_C, v_C, T_C, rho_throat_C, v_throat_C, T_throat_C] = Conservative_form(n, gamma, nt,x,dx,C)
Time_elasped_C = toc;
fprintf('Time Elasped for Conservative form: %0.4f', Time_elasped_C);
%Plotting
%Density
figure(1)
subplot(4,1,1)
plot(rho_throat_NC,'color','r')
ylabel('Density')
legend({'Density at throat'});
grid minor;
title('Timestep variation at throat area - NON CONSERVATIVE FORM')
%Temperature
subplot(4,1,2);
plot(T_throat_NC,'color','g')
ylabel('Temperature')
legend({'Temperature'});
grid minor;
%Pressure
subplot(4,1,3)
plot(Pressure_throat_NC,'color','b')
ylabel('Pressure')
legend({'Pressure at throat'});
grid minor
%Mach Number
subplot(4,1,4)
plot(Mach_number_throat_NC,'color','y')
xlabel('Timesteps')
ylabel('Mach Number')
legend({'Mach Number at throat'});
grid minor;
figure(2)
%Density
subplot(4,1,1)
plot(rho_throat_C,'color','r')
ylabel('Density')
legend({'Density at throat'});
grid minor;
title('Timestep variation at throat area - CONSERVATIVE FORM')
%Temperature
subplot(4,1,2)
plot(T_throat_C,'color','g')
ylabel('Temperature')
legend({'Temperature at throat'});
grid minor;
%Pressure
subplot(4,1,3)
plot(Pressure_throat_C,'color','b')
ylabel('Pressure')
legend({'Pressure at throat'});
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'});
grid minor;
figure(3)
plot(x,Mass_flow_rate_NC,'b',x,Mass_flow_rate_C,'g','linewidth',2)
xlabel('Length of nozzle')
ylabel('Mass flow rate')
legend('Conservative form','Non-Conservative form')
title('Mass flow rate using "Non-Conservative form vs Conservative form" in the nozzle')
figure(4)
subplot (2,1,1)
plot(Mass_flow_rate_throat_C,'linewidth',2)
hold on
xlabel('Time steps (nt)')
ylabel('Mass flow rate')
% legend('@31 nodes','@52 nodes','@63 nodes')
title('Conservative form: Grid independency test at throat point i=16')
subplot (2,1,2)
plot(Mass_flow_rate_throat_NC,'linewidth',2)
hold on
xlabel('Time steps (nt)')
ylabel('Mass flow rate')
% legend('@31 nodes','@52 nodes','@63 nodes')
title('Non-Conservative form: Grid independency test at throat point i=16')
- Function for conservative Form
function[Mass_flow_rate_throat_C, Pressure_throat_C, Mach_number_throat_C,Mass_flow_rate_C, Pressure_C, Mach_number_C, rho_C, v_C, T_C, rho_throat_C, v_throat_C, T_throat_C ] = Conservative_form(n, gamma, nt,x,dx,C)
%Determining Initial profile of Density and Temperature
for i = 1:length(x)
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 j = 2:n-1
%Updating the Source term(J)
J2(j) = (1/gamma)*rho_C(j)*T_C(j)*((a(j) - a(j - 1))/dx);
%Simplifying the Equation's spatial terms
dU1_dt_C(j) = -((F1(j) - F1(j - 1))/dx);
dU2_dt_C(j) = -((F2(j) - F2(j - 1))/dx) + J2(j);
dU3_dt_C(j) = -((F3(j) - F3(j - 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
%Applying 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);
hold on
elseif k == 100
plot(x, Mass_flow_rate_C, 'c', 'linewidth', 1);
hold on
elseif k == 200
plot(x, Mass_flow_rate_C, 'g', 'linewidth', 1);
hold on
elseif k == 400
plot(x, Mass_flow_rate_C, 'y', 'linewidth', 1);
hold on
elseif k == 800
plot(x, Mass_flow_rate_C, 'k', 'linewidth', 1);
hold on
end
end
%Labelling
title('Mass flow rates at different Timesteps - 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.4 1])
grid on
figure(1000)
title('Conservative form - Steady-state distribution of primitive variables inside the nozzle')
% plot(x,v,'b',x,p,'r',x,T,'k',x,rho,'m')
ylim([-0.1 3.7])
title (sprintf('Time step = %d',i))
subplot(4,1,1)
plot(rho_C)
xlabel('Length of nozzle')
ylabel('Non dimensional parameters')
legend('rho_{C}')
subplot(4,1,2)
plot(v_C)
xlabel('Length of nozzle')
ylabel('Non dimensional parameters')
legend('v_{C}')
subplot(4,1,3)
plot(T_C)
xlabel('Length of nozzle')
ylabel('Non dimensional parameters')
legend('T_{C}')
subplot(4,1,4)
plot(Mach_number_C)
xlabel('Length of nozzle')
ylabel('Non dimensional parameters')
legend('Mach number_{C}')
xlabel('Length of nozzle')
ylabel('Non dimensional parameters')
legend('Density','Velocity','Temperature','Mach number','location','northeast')
end- Function for Non- conservative form.
%function[Mass_flow_rate_throat_NC, Pressure_throat_NC, Mach_number_throat_NC,Mass_flow_rate_NC, Pressure_NC, Mach_number_NC, rho_NC, v_NC, T_NC, rho_throat_NC, v_throat_NC, T_throat_NC] = Non_conservative_form(n, gamma, nt,x,dx,C)
function[Mass_flow_rate_throat_NC, Pressure_throat_NC, Mach_number_throat_NC,Mass_flow_rate_NC, Pressure_NC, Mach_number_NC, rho_NC, v_NC, T_NC, rho_throat_NC, v_throat_NC, T_throat_NC ] = Non_conservative_form(n, gamma, nt,x,dx,C)
%Calculate initial profiles
rho_NC = 1-0.3146*x
T_NC = 1-0.2314*x
v_NC = (0.1+1.09.*x).*T_NC.^0.5
%Area & Throat profiles
a = 1 + 2.2*(x - 1.5).^2;
throat = (n+1)/2;
%outer time loop
for k=1:nt
%Time-step control
dt = min(C.*dx./(sqrt(T_NC) + v_NC));
%Saving a copy of Old Values
rho_old = rho_NC;
T_old = T_NC;
v_old = v_NC;
%predetor method
for j= 2:n-1
dv_dx = (v_NC(j + 1) - v_NC(j))/dx;
drho_dx = (rho_NC(j + 1) - rho_NC(j))/dx;
dT_dx = (T_NC(j + 1) - T_NC(j))/dx;
dlog_a_dx = (log(a(j + 1)) - log(a(j)))/dx;
%Continuity Equation
drho_dt_P(j) = -rho_NC(j)*dv_dx - rho_NC(j)*v_NC(j)*dlog_a_dx - v_NC(j)*drho_dx;
%Momentum Equation
dv_dt_P(j) = -v_NC(j)*dv_dx - (1/gamma)*(dT_dx + (T_NC(j)*drho_dx/rho_NC(j)));
%Energy Equation
dT_dt_P(j) = -v_NC(j)*dT_dx - (gamma - 1)*T_NC(j)*(dv_dx + v_NC(j)*dlog_a_dx);
%Solution Update
rho_NC(j) = rho_NC(j) + drho_dt_P(j)*dt;
v_NC(j) = v_NC(j) + dv_dt_P(j)*dt;
T_NC(j) = T_NC(j) + dT_dt_P(j)*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
drho_dt_C(j) = - rho_NC(j)*dv_dx - rho_NC(j)*v_NC(j)*dlog_a_dx - v_NC(j)*drho_dx;
%Momentum Equation
dv_dt_C(j) = - v_NC(j)*dv_dx - (1/gamma)*(dT_dx + (T_NC(j)*drho_dx/rho_NC(j)));
%Energy Equation
dT_dt_C(j) = - v_NC(j)*dT_dx - (gamma - 1)*T_NC(j)*(dv_dx + v_NC(j)*dlog_a_dx);
end
%Compute 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 i = 2:n-1
rho_NC(i) = rho_old(i) + drho_dt(i)*dt;
v_NC(i) = v_old(i) + dv_dt(i)*dt;
T_NC(i) = T_old(i) + dT_dt(i)*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);
hold on
elseif k == 100
plot(x, Mass_flow_rate_NC, 'c', 'linewidth', 1);
hold on
elseif k == 200
plot(x, Mass_flow_rate_NC, 'g', 'linewidth', 1);
hold on
elseif k == 400
plot(x, Mass_flow_rate_NC, 'y', 'linewidth', 1);
hold on
elseif k == 800
plot(x, Mass_flow_rate_NC, 'k', 'linewidth', 1);
hold on
end
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
figure(100)
title('Non Conservative form - Steady-state distribution of primitive variables inside the nozzle')
%plot(x,v,'b',x,p,'r',x,T,'k',x,rho,'m')
%ylim([-0.1 3.7])
%title (sprintf('Time step = %d',i))
subplot(4,1,1)
plot(v_NC)
xlabel('Length of nozzle')
ylabel('Non dimensional parameters')
legend('v_{NC}')
subplot(4,1,2)
plot(rho_NC)
xlabel('Length of nozzle')
ylabel('Non dimensional parameters')
legend('rho_{NC}')
subplot(4,1,3)
plot(T_NC)
xlabel('Length of nozzle')
ylabel('Non dimensional parameters')
legend('T_{NC}')
subplot(4,1,4)
plot(Mach_number_NC)
xlabel('Length of nozzle')
ylabel('Non dimensional parameters')
legend('Mach number_{NC}')
- Stepwise Explaination of the code-
- Code started with clear all ,close all , clc .
- Inupts parameters are assigned for the variables x, dx, n, gamma, nt and C.
- Array are created for the function x .
- Functions are created for the conservative form and non conservative form with necessary input parameters.
- To calculate the initial profiles ,flow parameters values are assigned i.e, density , temperature and velocity.
- For the throat profile,
throat = (n+1)/2; formula are used.
- tic toc command is used for calculating the time elasped between the functions.
- For outer time loop for loop command is used is used.
- MacCarmack technique is used to solve the numerical solution of governing equations such as Continuity, Momentum & Energy Equation.
- The technique implies two methods for solving namely Predictor method followed by Corrector method.
- The calculated values are averaged in the next step and updating the solution with old values.
- Predator method is used for continuity ,momentum and energy equations.
- Corrector method are used to simplyfying the governing equations.
- Then average time derivatives are computed.
- Final solution are updated for the flow parameters and boundary conditions are assigned at inlet and outlet of nozzle.
- Defined the mass flow rate , mach number and pressure values.
- Flow parameters are calculated for throat section of the nozzle.
- Plots are done for various parameters.
- 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 - Each plot is calculated as sub-plots of four primitive variables such as density, temperature, pressure and throat Mach number, as well as convergent/divergent sections.
- Comparison of non-dimensional Mass flow rates at different time steps is done.
- Outputs-
- Time steps variation of primitive variable at Throat area
- Conservative form
- Non conservative form
- From the above graphs, we can conclude that the initial deviations observed in the non-conservative form over a period of time require a higher amount of time to reach a stable solution compared to the conservative form.
- Grid independency test
- The total time required to solve the conservative form is longer compared to the non-conservative form.
- Although the final solutions of the two forms are almost similar, the non-conservative form oscillates considerably compared to the conservative form. The latter also provides an almost stable solution along the entire length. From this we can conclude that the conservative form reaches stability in fewer iterations with the earliest possible time.
- Here by observing the above graphs we could see that both forms are in the same phase and give exactly the same solution for all values.
- The above graphs show that the mass flow rate at different time steps along the distance initially oscillates, but the solution settles to a constant value after the 400th time step for both forms.
- Conclusion-
- we can obviously conclude that the conservative form converges slightly faster than the non-conservative form which offers a stable solution along the simulation.
- The complexity in solving the conservative form is more compared to the non-conservative form because the conservative part has more mathematical calculations compared to other forms