FINAL INDEPENDENT PROJECT

Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method

Given Problem:-

Our objective in this project is to write code solve the 1D supersonic nozzle flow equations using the Macormack Method.

• You will be implementing both the conservative and non-conservative forms of the governing equations
• You will perform grid dependence test
• Write separate functions for conservative and non-conservative forms
• Figure out the minimum number of cycles for which the simulation should be run in order for convergence
• Compare the normalized mass flow rate between the conservative forms and non-conservative forms - make suitable conclusions
• Which method is faster?

Solution:-

Our Solution will consist of following parts in given series:-

1. What is 1D super sonic nozzle Flow?

2. What is macoromac method?

3. Equations we are using

4. Result Graph for Time-wise variation of pressure, temperature, velocity and Mach number at the nozzle throat(for n=61 and n=31)

5. Result graph for Qualitative aspects of Quasi 1D nozzle flow at steady state.(for n=61 and n=31)

6. Result graph for Non-Dimensional mass flowrate distribution at various time-steps.(for n=61 and n=31)

7. Result graph on Comparison of the non-dimensional mass flowrates of both the forms.(for n=61 and n=31)

8. Grid-dependancy test

9. Analysis of Computation time and number of timesteps required for convergence

10. MATLAB code for MAIN SCRIPT

11. MATLAB code for NON-CONSERVATIVE FORM OF THE GOVERNING EQUATION: (USER-DEFINED FUNCTION)

12. MATLAB code for CONSERVATIVE FORM OF THE GOVERNING EQUATIONS: (USER-DEFINED FUNCTION)

1. 1d Super sonic Nozzle flow

When mach number crosses value 1 then this condition is known as super sonic condition i.e velocity faster than sound. Mach number generally represent speed of aircraft and missile with respect to sound.

This given nozzle has 3 parts mainly i.e. Convergent part(Subsonic), Throat(Sonic) and Divergent Part(Super Sonic). The fluid have subsonic speed in convergent part, Sonic speed in throat and super sonic speed in divergent part.

In 1D case of nozzle, the air or fluid can only change or vary with respect to single axis.

2.Macormac Method:-

MacCormack's technique is used, which provides 2nd order accuracy without the involvement of 2nd derivative, thus reducing the complexity of the mathematical calculation. We will be setting up finite difference solution, using time marching approach.

3. Equations we are using:-

We consider the stady, isentropic flow through a convergent-divergent nozzle. The flow at the inlet to the nozzle comes from a reservoir where pressure and temperatures are at stagnation. Also it is assumed that the crosss sectional area of the nozzle is very large so that the velocity is almost V≅0. The flow expands adiabatically to supersonic speeds at the nozzle exit. the flow is locally subsonic in the convergent section (M<1), sonic at the throat section (M=1) and supersonic at the divergent section (M>1).

We define the necessary non-dimensional variables T',ρ',x',V',t'and A' as:

In the problem specification we have taken the number of timesteps as 1400 as mentioned in Chapter 7 of the book: John D. Anderson Jr. Computational Fluid Dynamics: The Basics with Applications. Other problem specifications were taken directly from the above book and the results obtained was used to compare and analyze the solutions obtained form the programs.

4. Result Graph for Time-wise variation of pressure, temperature, velocity and Mach number at the nozzle throat.

(a) Non-Conservative Form: (Grid-size n=61)

(b) Conservative Form: (Grid-size n = 61)

After analysing the above two cases, we can say that timewise variation of quantities at the nozzle throat (n=15) for conservative form is less then Non conservative form. As we can see clearly in above conservative graph that Conservative form starting was unstable somehow but soon it got stable and become good.

The result is Nearly same as compare to reference book.

5. Result graph for Qualitative aspects of Quasi 1D nozzle flow at steady state

(a). Non-Conservative form: (Grid-size n = 61)

(b). Conservative Form: (Grid-size n = 61)

Its same as compare to refrence given.

6. Result graph for Non-Dimensional mass flowrate distribution at various time-steps

a). Non-Conservative form: (Grid-size n = 61)

(b). Conservative Form: (Grid-size n = 61)

Mass flow rate variation in conservative form is slightly less less than Non conservative form.

7. Result graph on Comparison of the non-dimensional mass flowrates of both the forms

As seen from the plot the conservative form gives a mass flow rate variation that is close to a constant value and simultaneously close to the exact value of 0.579 at the steady state. The non-conservative form shows some oscillations around the inflow and outflow boundaries and is less satisfactory compared to the conservation form.

Now Lets show the graph for Grid size= 31.

(a). Non-Conservative Form:(n=31)

(b). Conservative Form:(n=31)

8. Grid-dependancy test

We would be examining the value of flow quantities by changing the value of grid-size. The results are summarized below:

Exact Analytical Solution:-

ρ/ρo=0.634

T/To=0.833

p/po=0.528

M=1

Index(i)=---

NON CONSERVATIVE FORM-

For N=31

ρ/ρo=0.6387

T/To=0.8365

p/po=0.5342

M=0.9994

Index(i)=16

For N=61

ρ/ρo=0.6376

T/To=0.8354

p/po=0.5326

M=0.9997

Index(i)=31

Conservative Form-

For N=31

ρ/ρo=0.6504

T/To=0.8400

p/po=0.5464

M=0.9817

Index(i)=16

For N=61

ρ/ρo=0.6382

T/To=0.8352

p/po=0.5330

M=0.9993

Index(i)=31

From the above results we can see that for a grid-size of 30 points the Non-Conservative form gives a more accurate result compared to the Conservative form. However when we change the grid-size to 60 points the values of the flow quantites are quite close when calculated by either of the methods. Additionally both the methods give results which are clsoer to the exact value in this case.

9. Analysis of Computation time and number of timesteps required for convergence

NON CONSERVATIVE FORM-

For N=31, Computational Time= 4.5740

Minimum number of Time-steps required for steady state= 680(Approx)

For N=61, Computational Time= 6.2957

Minimum number of Time-steps required for steady state= 925(Approx)

Conservative Form-

For N=31, Computational Time= 5.2501

Minimum number of Time-steps required for steady state= 400(approx)

For N=61, Computational Time= 7.3008

Minimum number of Time-steps required for steady state= 600(approx)

10. MATLAB code for MAIN SCRIPT

%% MATLAB Program to simulate the 1D Supersonic Nozzle flow using Mc-Cormack Method:
%% Analysis of Conservative and Non-conservative forms of governing equations and grid dependancy test:
clear all;
close all;
clc;

%% Defining the initial conditions:
L = 3; %Length of the nozzle
n = 31; %Manually entering the number of grid-points
x = linspace(0,L,n);  % Initializing the 1D grid
dx = L/(n-1); % Grid spacings
gamma = 1.4; % Ratio of Cp and Cv for air assuming it to be ideal
nt = 1400; % number of timesteps
throat = (n-1)/2; % Throat point
a = 1 + 2.2*(x - 1.5).^2; % Area calculation
C = 0.5; % Courant Number (specified for the problem)

%% Solution by Non-Conservative form:
tic
[m, M, v, Temp, p, rho, Temp_throat, rho_throat, p_throat, M_throat] = non_conservative_form(x,dx,n,C,a,nt,gamma,throat);
time_elapsed = toc;
fprintf('The Computation time by Non-Conservative form: %0.4f',time_elapsed);

%% Time-wise variation of density, temperature, pressure and Mach number for the Non-Conservative form:
time_variation = figure(2);

subplot(4,1,1);
plot(rho_throat,'b');
title('Time-wise variations in quantity in Non-Conservative form');
legend('Density');
ylabel('rho/rho_{o}');
axis([0 1400 0 1.3]);
grid minor;

subplot(4,1,2);
plot(Temp_throat,'r');
legend('Temperature');
ylabel('T/T_{o}');
axis([0 1400 0.6 1]);
grid minor;

subplot(4,1,3);
plot(p_throat,'g');
legend('Pressure');
ylabel('p/p_{o}');
axis([0 1400 0.4 1.1]);
grid minor;

subplot(4,1,4);
plot(M_throat,'m');
legend('Mach Number');
ylabel('M');
axis([0 1400 0.6 1.4]);
grid minor;
xlabel('Number of iterations');

%% Qualitative Aspects of 1D flow in Non-Conservative Form:
qualitative_non_cons = figure(3);
subplot(4,1,1);
plot(x,M,'-b');
title('Qualitative Aspects of 1D Supersonic Nozzle flow: Non-Conservative form');
ylabel('M');
legend('Mach Number','location','southeast');
grid minor;

subplot(4,1,2);
plot(x,p,'-r');
ylabel('p/p_{o}');
legend('Pressure ratio');
grid minor;

subplot(4,1,3);
plot(x,rho,'-g');
ylabel('rho/rho_{o}');
legend('Density ratio','outside');
grid minor;

subplot(4,1,4);
plot(x,Temp,'-m');
ylabel('T/T_{o}');
legend('Temperature ratio','outside');
grid minor;

%% Solution by Conservative form:
tic
[m_c, M_c, v_c, Temp_c, p_c, rho_c, Temp_throat_c, rho_throat_c, p_throat_c, M_throat_c] = conservative_form(x,dx,n,C,a,nt,gamma,throat);
time_elapsed_cons = toc;
fprintf('The computation time by Conservative form: %0.4f',time_elapsed_cons);

%% Time-wise variation of density, temperature, pressure and Mach number for the Non-Conservative form:
time_variation_cons = figure(5);

subplot(4,1,1);
plot(rho_throat_c,'b');
title('Time-wise variations in quantity in Conservative Form');
legend('Density');
ylabel('rho/rho_{o}');
axis([0 1400 0 1.3]);
grid minor;

subplot(4,1,2);
plot(Temp_throat_c,'r');
legend('Temperature');
ylabel('T/T_{o}');
axis([0 1400 0.6 1]);
grid minor;

subplot(4,1,3);
plot(p_throat_c,'g');
legend('Pressure');
ylabel('p/p_{o}');
axis([0 1400 0.4 1.1]);
grid minor;

subplot(4,1,4);
plot(M_throat_c,'m');
legend('Mach Number');
ylabel('M');
axis([0 1400 0.6 1.4]);
grid minor;
xlabel('Number of iterations');

%% Qualitative Aspects of 1D flow in Conservative Form:
qualitative_cons = figure(6);
subplot(4,1,1);
plot(x,M_c,'-b');
title('Qualitative Aspects of 1D Supersonic Nozzle flow: Conservative form');
ylabel('M');
legend('Mach Number','location','southeast');
grid minor;

subplot(4,1,2);
plot(x,p_c,'-r');
ylabel('p/p_{o}');
legend('Pressure ratio');
grid minor;

subplot(4,1,3);
plot(x,rho_c,'-g');
ylabel('rho/rho_{o}');
legend('Density ratio','outside');
grid minor;

subplot(4,1,4);
plot(x,Temp_c,'-m');
ylabel('T/T_{o}');
legend('Temperature ratio','outside');
grid minor;

%% Comparison of Mass-Flow rates by Conservative and Non-Conservative methods and comparison with Exact Form:
mass_comparison = figure(7);
hold on;
plot(x, m, 'g','LineWidth', 1.75);
plot(x, m_c, 'm','LineWidth', 1.75);
line([0 3], [0.579 0.579],'LineWidth', 1.75);
legend('Conservative Form', 'Non-Conservative Form', 'Exact Solution');
grid minor;
xlabel('x/L');
ylabel('rhoVA/rho_{o}V_{o}A_{0}');
title('Non-Dimensional Mass Flow comparison with Exact Solution');

11. MATLAB code for NON-CONSERVATIVE FORM OF THE GOVERNING EQUATION: (USER-DEFINED FUNCTION)

function [m, M, v, Temp, p, rho, Temp_throat, rho_throat, p_throat, M_throat] = non_conservative_form(x,dx,n,C,a,nt,gamma,throat)

% Initializing the variables:
rho = 1 - 0.3146*x;
Temp = 1- 0.2314*x;
v = (0.1 + 1.09*x).*Temp.^0.5;

%Outer Time Loop
for k = 1:nt
    rho_old = rho;
    Temp_old = Temp;
    v_old = v;
    %Adaptive time-step control
    dt = min(C.*dx./(sqrt(Temp)+v));
    
    %Predictor Part followed by Solution Update:
    for j = 2: n-1
        dvdx = (v(j+1)-v(j))/dx;
        drhodx = (rho(j+1) - rho(j))/dx;
        dTempdx = (Temp(j+1) - Temp(j))/dx;
        dlog_a_dx = (log(a(j+1))-log(a(j)))/dx;
        % Predictor Step:
        %Continuity Equation:
        drhodt_p(j) = -rho(j)*dvdx - rho(j)*v(j)*dlog_a_dx - v(j)*drhodx;
        % Momentum Equation:
        dvdt_p(j) = -v(j)*dvdx - (1/gamma)*(dTempdx + (Temp(j)/rho(j))*drhodx);
        % Energy Equation:
        dTdt_p(j) = -v(j)*dTempdx - (gamma-1)*Temp(j)*(dvdx + v(j)*dlog_a_dx);
        
        % Solution Update step:
        rho(j) = rho(j) + drhodt_p(j)*dt; %Density Update
        v(j) = v(j) + dvdt_p(j)*dt; %Velocity Update
        Temp(j) = Temp(j) + dTdt_p(j)*dt; %Temperature Update
        
    end
    % Corrector step
    for j = 2: n-1
        dvdx = (v(j)-v(j-1))/dx;
        drhodx = (rho(j) - rho(j-1))/dx;
        dTempdx = (Temp(j) - Temp(j-1))/dx;
        dlog_a_dx = (log(a(j))-log(a(j-1)))/dx;
        %Continuity Equation:
        drhodt_c(j) = -rho(j)*dvdx - rho(j)*v(j)*dlog_a_dx - v(j)*drhodx;
        % Momentum Equation:
        dvdt_c(j) = -v(j)*dvdx - (1/gamma)*(dTempdx + (Temp(j)/rho(j))*drhodx);
        % Energy Equation:
        dTdt_c(j) = -v(j)*dTempdx - (gamma-1)*Temp(j)*(dvdx + v(j)*dlog_a_dx);
        
    end
    
    % Compute the 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
        rho(i) = rho_old(i) + drhodt(i)*dt;
        v(i) = v_old(i) + dvdt(i)*dt;
        Temp(i) = Temp_old(i) + dTdt(i)*dt;
    end
    % Applying the boundary conditions
    % Inlet
    v(1) = 2*v(2) - v(3);
    rho(1) = 1;
    Temp(1) = 1;
    % Outlet
    v(n) = 2*v(n-1) - v(n-2);
    rho(n) = 2*rho(n-1) - rho(n-2);
    Temp(n) = 2*Temp(n-1) - Temp(n-2);
    
    % Defining the values for pressure, mass flow rate and mach number:
    p = rho.*Temp; %Non-Dimensional Pressure
    m = rho.*a.*v; %Non-Dimensional Mass flow rate
    M = v./sqrt(Temp); %Mach Number
    
    %Calculation of Quantities at the Throat point of the nozzle:
    Temp_throat(k) = Temp(throat);
    p_throat(k) = p(throat);
    rho_throat(k) = rho(throat);
    M_throat(k) = M(throat);
    
    % Comparison of non-dimensional mass-flow rates at different times:
    mass = figure(1);
    hold on;
    if k == 50
        plot(x,m,'b','linewidth',1.75);
    else if k == 100
            plot(x,m,'g','linewidth',1.75);
        else if k == 150
                plot(x,m,'r','linewidth',1.75);
            else if k == 200
           

12. MATLAB code for CONSERVATIVE FORM OF THE GOVERNING EQUATIONS: (USER-DEFINED FUNCTION)

function [m_c, M_c, v_c, Temp_c, p_c, rho_c, Temp_throat_c, rho_throat_c, p_throat_c, M_throat_c] = conservative_form(x,dx,n,C,a,nt,gamma,throat)

%% Initialization step:
% Initializing density and temperature:
for i = 1:length(x)
    if x(i) >= 0 && x(i) < 0.5
        rho_c(i) = 1;
        Temp_c(i) = 1;
    else if x(i) >= 0.5 && x(i) < 1.5
            rho_c(i) = 1.0 -0.366*(x(i) - 0.5);
            Temp_c(i) = 1.0 - 0.167*(x(i) - 0.5);
        else if x(i) >= 1.5 && x(i) <= 3.0
                rho_c(i) = 0.634 - 0.3879*(x(i) - 1.5);
                Temp_c(i) = 0.833 - 0.3507*(x(i) - 1.5);
            end
        end
    end
end
% Initializing velocity and pressure:
v_c = 0.59./(rho_c.*a);
p_c = rho_c.*Temp_c;

% Initializing the solution vectors:
U1 = rho_c.*a;
U2 = rho_c.*a.*v_c;
U3 = rho_c.*a.*((Temp_c./(gamma - 1)) + (gamma/2)*v_c.^2);

% Time Iteration Loop
for k = 1:nt
    
    % Preserving the old solution vectors:
    U1_old = U1;
    U2_old = U2;
    U3_old = U3;
    
    % Adaptive Time-step control:
    dt = min(C.*dx./(sqrt(Temp_c) + v_c));
    
    % Defining the 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);
    
    %% Solution by Mc-Cormack Method:
    for j = 2:n-1
        
        % Defining the source term:
        J2(j) = (1/gamma).*rho_c(j).*Temp_c(j).*((a(j+1)-a(j))/dx);
        % Predictor step:
        dU1dt_p(j) = -(F1(j+1) - F1(j))/dx;
        dU2dt_p(j) = -((F2(j+1) - F2(j))/dx) + J2(j);
        dU3dt_p(j) = -(F3(j+1) - F3(j))/dx;
        
        % Updating the new values of solution vectors:
        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
    
    % Updating the primitive variables and flux vectors:
    
    rho_c = U1./a;
    v_c = U2./U1;
    Temp_c = (gamma -1)*(U3./U1 - (gamma/2)*v_c.^2);
    p_c = rho_c.*Temp_c;
    
    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);
    
    for j = 2:n-1
        % Updating the source term:
        J2(j) = (1/gamma)*rho_c(j).*Temp_c(j)*((a(j)-a(j-1))/dx);
        % Corrector step:
        dU1dt_c(j) = -(F1(j) - F1(j-1))/dx;
        dU2dt_c(j) = -((F2(j) - F2(j-1))/dx) + J2(j);
        dU3dt_c(j) = -(F3(j) - F3(j-1))/dx;
    end
    % Updating the solution vectors after corrector step:
    % Averaging the predictor and corrector steps:
    dU1dt = 0.5*(dU1dt_p + dU1dt_c);
    dU2dt = 0.5*(dU2dt_p + dU2dt_c);
    dU3dt = 0.5*(dU3dt_p + dU3dt_c);
    
    % Calculating the new solution vectors:
    for j = 2:n-1
        U1(j) = U1_old(j) + dU1dt(j)*dt;
        U2(j) = U2_old(j) + dU2dt(j)*dt;
        U3(j) = U3_old(j) + dU3dt(j)*dt;
    end
    % Applying the boundary conditions:
    % left side/ inlet conditions
    % U1(1) = rho_c(1).*Temp_c(1);
    U2(1) = 2*U2(2) - U2(3);
    % U3(1) = U1(1)*(Temp_c(1)/(gamma-1)+(gamma/2)*v_c(1).^2);
    
    % right side/outlet conditions:
    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 updation of the primitive variables:
    rho_c = U1./a;
    v_c = U2./U1;
    Temp_c = (gamma -1)*(U3./U1 - (gamma/2)*v_c.^2);
    p_c = rho_c.*Temp_c;
    %%
    % Defining the Values of mach Number and mass flow rate:
    M_c = v_c./sqrt(Temp_c);
    m_c = rho_c.*a.*v_c;
    
    % Calculating values of specific quantites at the throat:
    M_throat_c(k) = M_c(throat);  % Mach Number
    p_throat_c(k) = p_c(throat);  % Pressure Ratio
    Temp_throat_c(k) = Temp_c(throat);  % Non-dimensional Temperature
    rho_throat_c(k) = rho_c(throat);  % Density Ratio
    
    % Comparison of non-dimensional mass-flow rates at different times:
    mass2 = figure(4);
    hold on;
    if k == 100
        plot(x,m_c,'b','linewidth',1.75);
    else if k == 200
            plot(x,m_c,'g','linewidth',1.75);
        else if k == 700
                plot(x,m_c,'r','linewidth',1.75);
            end
        end
    end
end
title('Comparison of non-dimensional mass flow rates at different time-instants (Conservative Form)');
xlabel('x/L (Non-dimensional distance through the nozzle)');
ylabel('rhoVA/rho_{o}V_{o}A_{o}');
legend('100Deltat','200Deltat','700Deltat');
axis([0 3 0.4 1.6]);
grid on;
grid minor;

end

13. Conclusion

1. Conservation form produces slightly more accurate results for the flux quantities (Ex. m') and the Non-conservation form produces slightly more accurate results for the primitive variables (ρ',T',V') but the non-conservative form is more accurate overall.

2. The non-conservative form leads to smaller residuals in computation.

3. The non-conservative form is faster in comparison to the conservative form when it comes to the computation times. Here conservative form requires marginally more time due to increased number of calculations.