Project

1-D Supersonic Nozzle Flow Simulation

JJay Gite18 Jul 2019

INTRODUCTION

The aim of this project is to determine the numerical solution of a 1-D Supersonic Nozzle Flow using the Macormack method. Both Conservative, as well as Non-Conservative Approach, are discussed in this project. Although analytically identical, the numerical solution of both these approaches tends to vary a little from each other. Grid dependence test and CFL based adaptive time control were programmed as a part of this project to ascertain the accuracy of obtained results.

Maccormack Method has an advantage over First Order Forward and Backward Differencing schemes that it does not introduce the diffusive errors into the numerical solution. Thus, this method is very useful to dump out diffusive errors in applications where no diffusion is desired to be a part of the numerical solution. Supersonic Nozzle flow simulation is one of those applications where the introduction of diffusive errors could spoil the entire solution and make it unrealistic. To avoid this, Maccormack method is used in this project.

The governing equation for this project comes from the finite control volume represented by the nozzle. A steady, isentropic flow through a convergent-divergent (C-D) Nozzle. The nozzle inlet comes from a reservoir where the pressure and temperature are at a stagnation. Also, the cross-section area of the nozzle is large enough to assume that inlet flow velocity (V) ~ 0. The flow expands isentropically to supersonic speeds at the nozzle exit, i.e. Mach Number (M)>1. The flow is defined as a quasi-1-D flow meaning that the flow properties vary only in x-direction. The governing equations are written such that the non-dimensionalized form of flow parameters are used.

As a general rule of thumb, when the control volume is fixed in space with the fluid particles moving through the space defined by the governing equations, a conservative approach is used. Likewise, when the control volume is moving with the fluid such that same fluid particle is being traced throughout its path, a non-conservative approach is used.

Initial Conditions, Boundary Conditions, Nozzle Geometry are in accordance with the data provided in Chapter 7 of the book: John D. Anderson Jr. Computational Fluid Dynamics: The Basics with Applications.

PROGRAMMING

Main Script

%%% Macormack Method to solve 1-D Supersonic Nozzle Flow %%%
% Main Program %

clear all
close all
clc

% Input Parameters
% Mesh Generation
n = 31;
x = linspace(0,3,n);
dx = x(2)-x(1);
gamma = 1.4;

% Area Profile

A = 1 + (2.2*(x-1.5).^2);
n_throat = find((A-1)<0.001);

Type_of_Solver = input(\'Select The Type of Solver: n=1 for Non-Conservative and n=2 for Conservative\');

% Non-Conservative Form Solver

if Type_of_Solver == 1
    
tic

[rho,v,T,p,M,m,p_throat,rho_throat,v_throat,T_throat,M_throat,k] = NC(n,x,dx,A,gamma,n_throat);

time_nc = toc;
fprintf(\'Computational Time: Non-Conservative Approach = %f sec\', toc);
fprintf(\'\\n Number of Iterations Required = %f\', k);

% Transient Properties at throat

figure(2);
% Density Ratio
subplot(4,1,1);
plot(rho_throat,\'b\');
title(\'Transient Analysis of Properties at Throat\');
legend(\'Density Ratio\');
ylabel(\'\\rho/\\rho_{o}\');
axis([0 700 0.5 1]);
grid minor;
% Temperature
subplot(4,1,2);
plot(T_throat,\'m\');
legend(\'Temperature\');
ylabel(\'T/T_{o}\');
axis([0 700 0.5 1]);
grid minor;
% Pressure
subplot(4,1,3);
plot(p_throat,\'g\');
legend(\'Pressure\');
ylabel(\'P/P_{o}\');
axis([0 700 0.2 1]);
grid minor;
% Mach-Number
subplot(4,1,4);
plot(M_throat,\'r\');
legend(\'Mach-Number\');
ylabel(\'M/M_{o}\');
axis([0 700 0.8 1.5]);
grid minor;
xlabel(\'Number of Iterations\');

% Steady Profiles

figure(3);
% Density Ratio
subplot(4,1,1);
plot(x,rho,\'b\');
title(\'Steady State Profiles of Nozzle\');
legend(\'Density Ratio\');
ylabel(\'\\rho/\\rho_{o}\');
grid minor;
% Temperature
subplot(4,1,2);
plot(x,T,\'m\');
legend(\'Temperature\');
ylabel(\'T/T_{o}\');
grid minor;
% Pressure
subplot(4,1,3);
plot(x,p,\'g\');
legend(\'Pressure\');
ylabel(\'P/P_{o}\');
grid minor;
% Mach-Number
subplot(4,1,4);
plot(x,M,\'r\');
legend(\'Mach-Number\');
ylabel(\'M/M_{o}\');
grid minor;
xlabel(\'Space Domain\');

% Conservative Form Solver

elseif Type_of_Solver == 2

tic

[rho_con,v_con,T_con,p_con,M_con,m_con,rho_throat_con,v_throat_con,T_throat_con,M_throat_con,p_throat_con] = C(n,x,dx,A,gamma,n_throat);

time_nc = toc;
fprintf(\'Computational Time: Conservative Approach = %f sec\', toc);

% Transient Properties at throat

figure(2);
% Density Ratio
subplot(4,1,1);
plot(rho_throat_con,\'b\');
title(\'Transient Analysis of Properties at Throat\');
legend(\'Density Ratio\');
ylabel(\'\\rho/\\rho_{o}\');
axis([0 700 0.5 1]);
grid minor;
% Temperature
subplot(4,1,2);
plot(T_throat_con,\'m\');
legend(\'Temperature\');
ylabel(\'T/T_{o}\');
axis([0 700 0.5 1]);
grid minor;
% Pressure
subplot(4,1,3);
plot(p_throat_con,\'g\');
legend(\'Pressure\');
ylabel(\'P/P_{o}\');
axis([0 700 0.2 1]);
grid minor;
% Mach-Number
subplot(4,1,4);
plot(M_throat_con,\'r\');
legend(\'Mach-Number\');
ylabel(\'M/M_{o}\');
axis([0 700 0.8 1.5]);
grid minor;
xlabel(\'Number of Iterations\');

% Steady Profiles

figure(3);
% Density Ratio
subplot(4,1,1);
plot(x,rho_con,\'b\');
title(\'Steady State Profiles of Nozzle\');
legend(\'Density Ratio\');
ylabel(\'\\rho/\\rho_{o}\');
grid minor;
% Temperature
subplot(4,1,2);
plot(x,T_con,\'m\');
legend(\'Temperature\');
ylabel(\'T/T_{o}\');
grid minor;
% Pressure
subplot(4,1,3);
plot(x,p_con,\'g\');
legend(\'Pressure\');
ylabel(\'P/P_{o}\');
grid minor;
% Mach-Number
subplot(4,1,4);
plot(x,M_con,\'r\');
legend(\'Mach-Number\');
ylabel(\'M/M_{o}\');
grid minor;
xlabel(\'Space Domain\');

end

Non-Conservative Approach Solver

%%% Macormack Method to solve 1-D Supersonic Nozzle Flow %%%
%%% Non-Conservative Approach %%%

function [rho,v,T,p,M,m,p_throat,rho_throat,v_throat,T_throat,M_throat,k] = NC(n,x,dx,A,gamma,n_throat)

% Initial Profiles

rho = 1-(0.3146*x);
T = 1-(0.2314*x);
v = (0.1+(1.09*x)).*T.^0.5;
nt = 700;
C = 0.5;
error = 1e9;
tol = 1e-5;


% Initialization of Time Loop
  k = 0;

% Time Marching Approach

  while error > tol

    k = k+1;
      
    rho_old = rho;
    v_old = v;
    T_old = T;
    m_old = rho_old.*A.*v_old;
    a_old = sqrt(T_old);
    
    
    % Time Step
    TS = C*dx./(a_old+v_old); % Courant Number time-step criteria
    dt = min(TS);
    
    % Predictor Step
    
    for j = 2 : n-1
        
        drhodx_p = (rho(j+1)-rho(j))/dx;
        dvdx_p = (v(j+1)-v(j))/dx;
        dlnA_dx_p = (log(A(j+1))-log(A(j)))/dx;
        dTdx_p = (T(j+1)-T(j))/dx;
        
        
        % COntinuity Equation
        
        drhodt_p(j) = (-rho(j)*dvdx_p) - ((rho(j)*v(j))*dlnA_dx_p) - (v(j)*drhodx_p);
        
        % Momentum Equation
        
        dvdt_p(j) = (-v(j)*dvdx_p) - ((1/gamma)*((dTdx_p)+((T(j)/rho(j))*drhodx_p)));
        
        % Energy Equation
        
        dTdt_p(j) = (-v(j)*dTdx_p) - ((gamma-1)*T(j))*(dvdx_p + (v(j)*dlnA_dx_p));    
        
        % Solution Update
        rho(j) = rho(j) + (drhodt_p(j)*dt);
        v(j) = v(j) + (dvdt_p(j)*dt);
        T(j) = T(j) + (dTdt_p(j)*dt);
        
    end
    
    % Corrector Step
    for j = 2 : n-1
        
        drhodx_c = (rho(j)-rho(j-1))/dx;
        dvdx_c = (v(j)-v(j-1))/dx;
        dTdx_c = (T(j)-T(j-1))/dx;
        dlnA_dx_c = (log(A(j))-log(A(j-1)))/dx;
        
        % COntinuity Equation
        
        drhodt_c(j) = (-rho(j)*dvdx_c) - ((rho(j)*v(j))*dlnA_dx_c) - (v(j)*drhodx_c);
        
        % Momentum Equation
        
        dvdt_c(j) = (-v(j)*dvdx_c) - ((1/gamma)*((dTdx_c)+((T(j)/rho(j))*drhodx_c)));
        
        % Energy Equation
        
        dTdt_c(j) = (-v(j)*dTdx_c) - ((gamma-1)*T(j))*(dvdx_c + (v(j)*dlnA_dx_c)); 
        
    end
    
    % Average Time Derivative
    
    drhodt_avg = 0.5*(drhodt_p+drhodt_c);
    dvdt_avg = 0.5*(dvdt_p+dvdt_c);
    dTdt_avg = 0.5*(dTdt_p+dTdt_c);
    
    % Final Update
    
    for i = 2:n-1
        
        rho(i) = rho_old(i) + (drhodt_avg(i)*dt);
        v(i) = v_old(i) + (dvdt_avg(i)*dt);
        T(i) = T_old(i) + (dTdt_avg(i)*dt);
        
    end
    
    % Boundary Conditions
    
    % Inlet
    rho(1) = 1;
    v(1) = (2*v(2))-v(3);
    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);
    
    % Calculating Important Parameters
    
    % Mach Number
    M = v./sqrt(T);
    % Mass-Flow
    m = rho.*A.*v;
    % Pressure
    p = rho.*T;
    
    % Throat Values
    
    rho_throat(k) = rho(n_throat);
    T_throat(k) = T(n_throat);
    M_throat(k) = M(n_throat);
    v_throat(k) = v(n_throat);
    p_throat(k) = p(n_throat);
    
    % Time-Loop error criteria for convergence
    % Mass-Flow is used as the comparing criteria, when mass-flow becomes
    % constant, it can be said that solution has converged
    
    error = max(abs(m-m_old));
    
    % Capturing change in mass-flow before steady state is reached
    
    figure(1);
    hold on;
   
    
    if k == 1
        plot(x,m,\'*\');
    elseif k == 50
        plot(x,m,\'--\');
    elseif k == 100
        plot(x,m,\'o\');
    elseif k == 150
        plot(x,m,\'r\');
    elseif k == 200
        plot(x,m,\'b\');
    elseif k ==400
        plot(x,m,\':\');
    elseif k == 700
        plot(x,m,\'g\');
    end
    
       
end
title(\'Transient Non-Dimensional Mass Flow Rate : Non-Conservative Approach\')
xlabel(\'x/L\');
ylabel(\'m/m{o}\')
legend(\'t=0\',\'t=50\',\'t=100\',\'t=150\',\'t=200\',\'t=400\',\'t=700\')
axis([0 x(n) 0.2 2]);
grid minor;

end

Conservative Approach Solver

%%% Macormack Method to solve 1-D Supersonic Nozzle Flow %%%
%%% Conservative Approach %%%

function [rho_con,v_con,T_con,p_con,M_con,m_con,rho_throat_con,v_throat_con,T_throat_con,M_throat_con,p_throat_con] = C(n,x,dx,A,gamma,n_throat)

nt = 700;
C = 0.5;

% Initial Profiles

for i = 1:length(x)
    if x(i)>=0 && x(i)<0.5
        rho_con(i) = 1;
        T_con(i) = 1;
    elseif x(i)>=0.5 && x(i)<1.5
        rho_con(i) = 1 - 0.366*(x(i)-0.5);
        T_con(i) = 1 - 0.167*(x(i)-0.5);
    elseif x(i)>=1.5 && x(i)<=3
        rho_con(i) = 0.634 - 0.3879*(x(i)-1.5);
        T_con(i) = 0.833 - 0.3507*(x(i)-1.5);
    end
end


v_con = (0.59)./(rho_con.*A);
rho_old_con = rho_con;
    v_old_con = v_con;
    T_old_con = T_con;

% Time Marching Solution

U1 = rho_con.*A;
U2 = rho_con.*A.*v_con;
U3 = rho_con.*A.*((T_con./gamma-1)) + ((gamma/2).*v_con.^2);

% Time Loop

for k = 1:nt
    
   
    U1_old = U1;
    U2_old = U2;
    U3_old = U3;
    
    
    % Time Step
    a_con = sqrt(T_con);
    TS = C*dx./(a_con+v_con); % Courant Number time-step criteria
    dt = 0.001;
    
    % Flux Vectors
    F1 = U2;
    F2 = (U2.^2./U1) + ((gamma-1)/gamma)*(U3-(gamma*U2.^2)./(2.*U1));
    F3 = (gamma*U2.*U3./U1) - ((gamma*(gamma-1)/2)*(U2.^3./U1.^2));
    
    % Predictor Step
    
    for j = 2:n-1
        
        J(j) = (1/gamma).*rho_con(j).*T_con(j).*((A(j+1)-A(j))/dx);
        
        % Continuity Equation
        dU1dt_con_p(j) = -(F1(j+1)-F1(j))/dx;
        
        % Momentum Equation
        dU2dt_con_p(j) = (-(F2(j+1)-F2(j))/dx) + J(j);
        
        % Energy Equation
        dU3dt_con_p(j) = -(F3(j+1)-F3(j))/dx;
        
        % Predictor Step Solution Update

        U1(j) = U1(j) + (dU1dt_con_p(j)*dt);
        U2(j) = U2(j) + (dU2dt_con_p(j)*dt);
        U3(j) = U3(j) + (dU3dt_con_p(j)*dt);
        
    end 
    
    % Predicted Value of Properties (Variables) of Interest and Flux
    
    rho_con = U1./A;
    v_con = U2./U1;
    T_con = (gamma-1)*((U3./U1)-((gamma/2)*(U2./U1).^2));
    p_con = rho_con.*T_con;
    F1 = U2;
    F2 = (U2.^2./U1) + (((gamma-1)/gamma)*(U3-(((gamma/2)*U2.^2)./U1)));
    F3 = (gamma*U2.*U3./U1) - ((gamma*(gamma-1)/2)*(U2.^3./U1.^2));
    
    
    % Corrector Step
    
    for j = 2:n-1
        
        J(j) = (1/gamma).*rho_con(j).*T_con(j).*((A(j)-A(j-1))/dx);
        
        % Continuity Equation
        dU1dt_con_c = -(F1(j)-F1(j-1))/dx;
        
        % Momentum Equation
        dU2dt_con_c = -(F2(j)-F2(j-1))/dx + J(j);
        
        % Energy Equation
        dU3dt_con_c = -(F3(j)-F3(j-1))/dx;
    
    end
    
    % Average Time Derivative
    
    dU1dt_con_avg = 0.5*(dU1dt_con_p + dU1dt_con_c);
    dU2dt_con_avg = 0.5*(dU2dt_con_p + dU2dt_con_c);
    dU3dt_con_avg = 0.5*(dU3dt_con_p + dU3dt_con_c);
    
    % Final Solution Update
    for i = 2:n-1
    
    U1(i) = U1_old(i) + (dU1dt_con_avg(i)*dt);
    U2(i) = U2_old(i) + (dU2dt_con_avg(i)*dt);
    U3(i) = U3_old(i) + (dU3dt_con_avg(i)*dt);
    
    end
    
    % Boundary Conditions 
    % Inlet
    rho_con(1) = 1;
    T_con(1) = 1;
    U1(1) = rho_con(1)*A(1);
    U2(1) = 2*U2(2)-U2(3);
    U3(1) = U1(1)*(T_con(1)/(gamma-1)+(gamma/2)*v_con(1).^2);
    
    % 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);
    
    % Calculation of Properties (Variables) of Interest
    
    rho_con = U1./A;
    v_con = U2./U1;
    T_con = (gamma-1)*(U3./U1-(gamma/2)*v_con.^2);
    p_con = rho_con.*T_con;
    
    % Calculation of Remaining Important Properties
    
    % Mach Number
    M_con = v_con./sqrt(T_con);
    %Mass-Flow
    m_con = rho_con.*A.*v_con;
    
    % Throat Values
    
    rho_throat_con(k) = rho_con(n_throat);
    T_throat_con(k) = T_con(n_throat);
    M_throat_con(k) = M_con(n_throat);
    v_throat_con(k) = v_con(n_throat);
    p_throat_con(k) = p_con(n_throat);
    
    % Capturing change in mass-flow before steady state is reached
    
    figure(1);
    hold on;
   
    
    if k == 1
        plot(x,m_con,\'*\');
    elseif k == 50
        plot(x,m_con,\'--\');
    elseif k == 100
        plot(x,m_con,\'o\');
    elseif k == 150
        plot(x,m_con,\'r\');
    elseif k == 200
        plot(x,m_con,\'b\');
    elseif k ==400
        plot(x,m_con,\':\');
    elseif k == 700
        plot(x,m_con,\'g\');
    end
    
       
end
title(\'Transient Non-Dimensional Mass Flow Rate : Conservative Approach\')
xlabel(\'x/L\');
ylabel(\'m/m{o}\')
legend(\'t=0\',\'t=50\',\'t=100\',\'t=150\',\'t=200\',\'t=400\',\'t=700\')
axis([0 x(n) 0.2 2]);
grid minor;

                 
end

Non-Conservative Approach Analysis

TRANSIENT VARIATION

a) Throat Section

Here, we analyze the variations of flow properties at the throat section of the nozzle as the iterations move ahead. The expected non-dimensional flow properties at the nozzle throat section are known in advance and the accuracy of their prediction is comparatively higher than other sections. This is why this section is chosen to check the stability of the numerical scheme. Any strange deviations either high or low of the numerical results from the expected results can be studied.

For a C-D nozzle, the flow attains M=1 (sonic level) at nozzle throat section. It can be seen from the graph that there are initial disturbances which get dumped out and the Mach Number (M) value gets stabilized at unity. Thus, it is evident that the numerical solution obtained is stable. Just like Mach Number, Non-dimensionalized Pressure, Temperature and Density become stable after initial disturbances are dumped out.

b) Normalized Mass Flow Rate

Mass Flow Rate is a useful parameter to check the consistency of the numerical scheme. The governing equation is solved for steady-state conditions. Thus, the mass flow rate should remain constant across the entire space domain (all the nodes) of the nozzle.

As can be seen from the graph, the initial distorted mass flow profile can be attributed to the multiplication of density, velocity and parabolic profile of the nozzle cross-section. The profile undergoes radical changes for every 50-time steps until time step of 200. which have been plotted. The mass flow rate profile attains fairly consistent value by 400-time steps and after 700 time-steps, it can be clearly seen that the mass-flow rate values have converged to a final solution. The final converged solution is extremely close to the analytical result of 0.579.

STEADY STATE ANALYSIS

a) USEFUL FLOW PROPERTIES

The Steady-State Analysis of crucial flow properties helps us to determine the behaviour of fluid flow through the C-D nozzle. This results obtained from the steady-state analysis can be trusted as the numerical scheme is stable and consistent. As can be seen from the graphs, the Mach Number (M) increases along the course of flow through the nozzle. This goes to show that the fluid accelerates through the nozzle and reaches super-sonic levels (M>1) as it exits the nozzle. This acceleration is at the expense of a simultaneous drop in Pressure, Temperature and Density of the fluid flow along the course of its flow through the nozzle.

Conservative Approach Analysis

TRANSIENT VARIATION

a) Throat Section

Once again, the throat section was chosen to check the stability of the numerical method for the same reasons as before. Graphs similar to those obtained from non-conservative approach were obtained. The peaks obtained during the start of the numerical solution were much flatter than those obtained from a non-conservative approach. Just like before, the numerical solution quickly stabilizes and converges to the final solution.

b) Normalized Mass Flow Rate

Mass Flow Rate, just as before, proves to be an important parameter to check the consistency of the numerical scheme. A much faster and smoother convergence to the stable and consistent mass flow rate was observed for the Conservative Approach.

This might lead us to think that the Conservative approach is best suited for this kind of problems. However, there are other criteria which need to be evaluated before we jump to any conclusion. One of those criteria is the grid-dependence test.

STEADY STATE ANALYSIS

The results obtained from steady-state analysis once again proved to stable, consistent and very close to the expected results from the analytical solution. The plots obtained were identical to those obtained from a non-conservative approach.

GRID DEPENDENCE TEST

The basic fundamental of this test is to alter the number of grid points and check if there are any changes in the results obtained. A stable numerical scheme will give very little fluctuation in the final result when subjected to changing grid points. The objective of any numerical approach is to have a scheme which continues to be stable and consistent irrespective of the size of the grid. In this project, we have changed the number of grid points from 31 to 61 and observed the result at the throat section. The results obtained are shown in the table below. It can be seen that the numerical scheme is independent of the grid size and continue to remain stable. For evaluating fundamental flow parameters, Non-Conservative approach is a better choice as it gives results close to those predicted from the analytical analysis.

Analytical Result:

Density Ratio	Temperature Ratio	Pressure Ratio	Mach Number (M)
0.634	0.833	0.528	1.0

For n = 31,

Non-Conservative Approach:

Density Ratio	Temperature Ratio	Pressure Ratio	Mach Number (M)
0.6386	0.8364	0.5341	0.9993

Conservative Approach:

Density Ratio	Temperature Ratio	Pressure Ratio	Mach Number (M)
0.6478	0.8387	0.5437	0.9872

For n = 61,

Non-Conservative Approach:

Density Ratio	Temperature Ratio	Pressure Ratio	Mach Number (M)
0.6377	0.8355	0.5328	0.9998

Conservative Approach:

Density Ratio	Temperature Ratio	Pressure Ratio	Mach Number (M)
0.6437	0.8378	0.5402	0.9887

1-D Supersonic Nozzle Flow Simulation

Ship work like this
in your own course.

Talk to a career counsellor

Ship work like thisin your own course.

Talk to a career counsellor

Ship work like this
in your own course.