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

Objective:-

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.

You need to show 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

Procedure:-

1D Super-Sonic Nozel Flow

A nozel is a device which is used to control the characteristics of fluid flow as it exits the closed chamber. Now in the above figure we are seeing a convergent-divergent nozzel in which we have assumed the flow parameters are only changing in the X direction. Such a condition where parameters change only in one direction is know as ' Quasi one dimensional flow'.

Governing equation for non conservative form

Program:-

Main program

clear all
close all
clc

% Initial Conditions
L = 3; % nozzel length
n = 31; % Grid points
x = linspace(0,L,n);
dx = L/(n-1); % grid space
nt = 1400; % time step
throat = (n-1)/2; % Throat
gamma = 1.4; % cp/cv
A = 1 + 2.2*(x-1.5).^2; % Area
C = 0.5; % courant number

% Solution for conservative form

tic
[m_c, M_c, V_c, T_c, P_c, rho_c, T_t_c, rho_t_c, P_t_c, M_t_c] = conservative_form(x,dx,n,C,A,nt,gamma,throat);
time_taken_conser = toc;
fprintf('The computation time by conservative form: %o.4f',time_taken_conser);

% solution for non-conservative form
tic
[m, M, V, T, P, rho, T_t, rho_t, P_t, M_t] = non_conservative_form(x,dx,n,C,A,nt,gamma,throat);
time_taken_non_conser = toc;
fprintf('The computation time by non-conservative form: %o.4f',time_taken_non_conser);

%Plots
% Convergence of parameters at nozzel throat at different time steps in conservative and non conservative form.

figure(1);
subplot(4,1,1);
plot(1:nt,rho_t_c,'b');
hold on;
plot(1:nt,rho_t,'r');
title('Variation of density w.r.t time');
legend('Conservative form','Non-conservative form');
ylabel('rho/rho_{o}');
xlabel('Time steps');
axis([0 1400 0 1.3]);

subplot(4,1,2);
plot(1:nt,T_t_c,'b');
hold on;
plot(1:nt,T_t,'r');
title('Variation of Temperature w.r.t time');
legend('Conservative form','Non-conservative form');
ylabel('T/T_{o}');
xlabel('Time steps');
axis([0 1400 0.6 1]);

subplot(4,1,3);
plot(1:nt,P_t_c,'b');
hold on;
plot(1:nt,P_t,'r');
title('Variation of Pressure w.r.t time');
legend('Conservative form','Non-conservative form');
ylabel('P/P_{o}');
xlabel('Time steps');
axis([0 1400 0.4 1.1]);

subplot(4,1,4);
plot(1:nt,M_t_c,'b');
hold on;
plot(1:nt,M_t,'r');
title('Variation of MAch number w.r.t time');
legend('Conservative form','Non-conservative form');
ylabel('M');
xlabel('Time steps');
axis([0 1400 0.6 1.8]);


% variation of non dimensional flow parameters along x axis in conservative and non-conservative form.

figure(2);
subplot(4,1,1);
plot(x,M_c,'b');
hold on;
plot(x,M,'*r');
title('Variation of non dimensional mach number along x axis');
legend('conservative form','non-conservative form');
ylabel('M');
xlabel('Length(x)');
grid minor;

subplot(4,1,2);
plot(x,P_c,'b');
hold on;
plot(x,P,'*r');
title('Variation of non dimensional pressure ratio along x axis');
legend('conservative form','non-conservative form');
ylabel('P/P_{o}');
xlabel('Length(x)');
grid minor;
axis([0 3 0 1]);

subplot(4,1,3);
plot(x,rho_c,'b');
hold on;
plot(x,rho,'*r');
title('Variation of non dimensional density ratio along x axis');
legend('conservative form','non-conservative form');
ylabel('rho/rho_{o}');
xlabel('Length(x)');
grid minor;
axis([0 3 0 1]);

subplot(4,1,4);
plot(x,T_c,'b');
hold on;
plot(x,T,'*r');
title('Variation of non dimensional Temperature ratio along x axis');
legend('conservative form','non-conservative form');
ylabel('T/T_{o}');
xlabel('Length(x)');
grid minor;
axis([0 3 0 1]);

% Comparison of mass flow rate quantity in conervative and non-conservative form against exact form

figure(3);
hold on;
plot(x,m_c,'r');
plot(x,m,'b');
line([0 3], [0.579 0.579],'linewidth', 2);
legend('Conservative form','Non-conservative form','Exact form');
xlabel('x/L');
ylabel('rhoVA/rho_{o}V_{o}A_{0}');
title('Non dimensional mass flow quantity in conervative and non-conservative form against exact form');

Function program for conservative form:-

function [m_c, M_c, V_c, T_c, P_c, rho_c, T_t_c, rho_t_c, P_t_c, M_t_c] = conservative_form(x,dx,n,C,A,nt,gamma,throat);

% Density and Temperature initialise
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

% Velocity and Pressure initialise
V_c = 0.59./(rho_c.*A);
P_c = rho_c.*T_c;

% solution vector parameter
U1 = rho_c.*A;
U2 = rho_c.*A.*V_c;
U3 = rho_c.*A.*((T_c./(gamma-1)) + (gamma/2)*V_c.^2);

% Time iteration loop
for k = 1:nt
    U1_old = U1;
    U2_old = U2;
    U3_old = U3;
    dt = min(C.*dx./(sqrt(T_c) + V_c));
    
    % flux vector
    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 calculation
    for j = 2:n-1
        
        J2(j) = (1/gamma).*rho_c(j).*T_c(j).*((A(j+1)-A(j))/dx);
        
        % continuity equation
        dU1dt_p(j) = -(F1(j+1) - F1(j))/dx;
        % momentum equation
        dU2dt_p(j) = -((F2(j+1) - F2(j))/dx) + J2(j);
        % energy equation
        dU3dt_p(j) = -(F3(j+1) - F3(j))/dx;
        
        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 values
    
    rho_c = U1./A;
    V_c = U2./U1;
    T_c = (gamma-1)*(U3./U1 -(gamma/2)*V_c.^2);
    P_c = rho_c.*T_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
        
        
        J2(j) = (1/gamma).*rho_c(j).*T_c(j).*((A(j)-A(j-1))/dx);
        
        % continuity equation
        dU1dt_c(j) = -(F1(j) - F1(j-1))/dx;
        % momentum equation
        dU2dt_c(j) = -((F2(j) - F2(j-1))/dx) + J2(j);
        % energy equation
        dU3dt_c(j) = -(F3(j) - F3(j-1))/dx;
    end
    
    % finding average predictor and corrector step
    dU1dt = 0.5*(dU1dt_p + dU1dt_c);
    dU2dt = 0.5*(dU2dt_p + dU2dt_c);
    dU3dt = 0.5*(dU3dt_p + dU3dt_c);
    
    % final solution vector
    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 boundary condition
    % 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));
    
    % 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);
    
    rho_c = U1./A;
    V_c = U2./U1;
    T_c = (gamma-1)*(U3./U1 - (gamma/2)*V_c.^2);
    P_c = rho_c.*T_c;
    M_c = V_c./sqrt(T_c);
    m_c = rho_c.*A.*V_c;
    
    % calculation of vallues at throat
    M_t_c(k) = M_c(throat);
    P_t_c(k) = P_c(throat);
    T_t_c(k) = T_c(throat);
    rho_t_c(k) = rho_c(throat);
    
    % plotting mass flow rate at different time step
    
    figure(4);
    hold on;
    if k == 50
        plot(x,m_c,'b','linewidth',2);
    elseif k == 100
        plot(x,m_c,'g','linewidth',2);
    elseif k == 150
        plot(x,m_c,'r','linewidth',2);
    elseif k == 200
        plot(x,m_c,'c','linewidth',2);
    elseif k == 400
        plot(x,m_c,'m','linewidth',2);
    elseif k == 700
        plot(x,m_c,'y','linewidth',2);
    elseif k == 1400
        plot(x,m_c,'k','linewidth',2);
        
    end
end

title('comparison of non dimensional mass flow rate at different time steps : Conservative form');
xlabel('x/L (non-dimensional distance)');
ylabel('rhoVA/rho_{o}V_{o}A_{o}');
legend('50','100','150','200','400','700','1400');
axis([0 3 0.4 1.2]);
end

Function program for non-conservative form:-

function [m, M, V, T, P, rho, T_t, rho_t, P_t, M_t] = non_conservative_form(x,dx,n,C,A,nt,gamma,throat);

% initial condition
rho = 1-0.3146*x;
T = 1-0.2314*x;
V = (0.1 + 1.09*x).*T.^0.5;
dt = min(C.*dx./(sqrt(T)+V));

% time loop
for k = 1:nt
    rho_old = rho;
    T_old = T;
    V_old = V;
    
    for j = 2:n-1
        dvdx = (V(j+1)-V(j))/dx;
        drhodx = (rho(j+1)-rho(j))/dx;
        dTdx = (T(j+1)-T(j))/dx;
        dlog_a_dx = (log(A(j+1))-log(A(j)))/dx;
        
        % 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)*(dTdx + (T(j)/rho(j))*drhodx);
        % energy equation
        dTdt_p(j) = -V(j)*dTdx - (gamma-1)*T(j)*(dvdx + V(j)*dlog_a_dx);
        
        % 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
        dvdx = (V(j)-V(j-1))/dx;
        drhodx = (rho(j) - rho(j-1))/dx;
        dTdx = (T(j)-T(j-1))/dx;
        dlog_a_dx = (log(A(j)) - log(A(j-1)))/dx;
        drhodt_c(j) = -rho(j)*dvdx - rho(j)*V(j)*dlog_a_dx - V(j)*drhodx;
        dvdt_c(j) = -V(j)*dvdx - (1/gamma)*(dTdx + (T(j)/rho(j))*drhodx);
        dTdt_c(j) =  -V(j)*dTdx - (gamma-1)*T(j)*(dvdx + V(j)*dlog_a_dx);
    end
    
     % calculating average of time derivative values
    drhodt = 0.5*(drhodt_p + drhodt_c);
    dvdt = 0.5*(dvdt_p + dvdt_c);
    dTdt = 0.5*(dTdt_p + dTdt_c);
    
    % final solution
    for i = 2:n-1
        rho(i) = rho_old(i) + drhodt(i)*dt;
        V(i) = V_old(i) + dvdt(i)*dt;
        T(i) = T_old(i) + dTdt(i)*dt;
    end
    
    % Applying boundary conditions
    % inlet
    V(1) = 2*V(2) - V(3);
    rho(1) = 1;
    T(1) = 1;
    % outlet
    V(n) = 2*V(n-1) - V(n-2);
    rho(n) = 2*rho(n-1) - rho(n-2);
    T(n) = 2*T(n-1) - T(n-2);
    
    P = rho.*T;
    m = rho.*A.*V;
    M = V./sqrt(T);
    
    % calculation of throat values
    rho_t(k) = rho(throat);
    T_t(k) = T(throat);
    P_t(k) = P(throat);
    M_t(k) = M(throat);
    
    figure(5);
    hold on;
    if k == 50
        plot(x,m,'b','linewidth',2);
        hold on;
        grid on;
    elseif k == 100
        plot(x,m,'g','linewidth',2);
    elseif k == 150
        plot(x,m,'r','linewidth',2);
    elseif k == 200
        plot(x,m,'c','linewidth',2);
    elseif k == 400
        plot(x,m,'m','linewidth',2);
    elseif k == 700
        plot(x,m,'y','linewidth',2);
    elseif k == 1400
        plot(x,m,'k','linewidth',2);
        
    end
end

title('comparison of non-dimensional mass flow rates at different time steps: Non-conservative form');
xlabel('x/L (non-dimensional distance through the nozzel)');
ylabel('rhoVA/rho_{o}V_{o}A_{o}');
legend('50','100','150','200','400','700','1400');
axis([0 3 0 2]);
end

Result:-

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:-

Conclusion:-

Conservative form produce slightly more accurate results for the flux quantities and the non-conservative form produce slightly more accurate results for primitive quantities. The non-conservative form is faster as compared to the conservative form.