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

Objective

In the project, we simulate isentropic flow through a quasi 1D supersonic converging-diverging nozzle using the MacCormack method and adaptive time steps. The simulation is done and analysed for both the conservation and the non-conservation forms of the governing equations. Specifically, the following analyses were the focus of the project:

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

Theory

A supersonic nozzle is a converging-diverging nozzle. In reality, the flow through the nozzle is 2-D because of the change in the area along the `x` direction. However, the assumption that the flow properties is uniform across any given cross-section of the nozzle makes the analyses quasi 1-D flow.

In Compuational Fluid Dynamics, the governing equations of flow can be cast in two forms namely,

• Conservative form
• Non-Conservative form.

Theoretically, there's no difference between these two forms. But in CFD, the distinction is in how easy it is to pose the numerical problem with one form or the other. Also, numerical solutions are more accurate for some problems when cast in one form compared to the other.

The derivation of the quasi 1-D equations in conservation and non-conservation forms are well explained in the Chapters 2 and 7 of the book Computational Fluid Dynamics: the basics with applications by J. D. Anderson, Jr. Readers are encouraged to peruse these chapters. However, the final forms of the equations are displayed below for the two forms.

Conservation form

`(del(rho'A'))/(delt') + (del(rho'A'u'))/(delx') = 0`

`(del(rho'A'u'))/(delt') + (del[rho'A'u'^2 +1/gammaP'A'])/(delx') = 1/gammaP'(delA')/(delx')`

`(del[rho'((T')/(gamma-1) +gamma/2u'^2)A'])/(delt') + (del[rho'((T')/(gamma-1) +gamma/2u'^2)u'A' + P'A'u'])/(delx') = 0`

These equations can be re-cast in the form

`(delU)/(delt) = -(delF)/(delx) +J`

where

`U = {[rho'A'],[rho'A'u'],[rho'((T')/(gamma-1) +gamma/2u'^2)A']}` the solution terms

`F = {[rho'A'u'],[rhoA'u'^2 + P'A'], [rho'((T')/(gamma-1) +gamma/2u'^2)u'A' + P'A'u']}` the convective terms

`J = {[0],[1/gammaP'(delA')/(delx')], [0]}` the source terms

`P' = rho'T'`

and in the Pure flux form, the convective terms and the source terms can be re-cast as as functions of the solution terms. i.e.

`F = {[U_2], [(U_2^2)/U_1 + (gamma-1)/gamma(U_3 - gamma/2(U_2^2)/U_1)],[gamma(U_2U_3)/U_1 - (gamma(gamma-1))/2(U_2^3)/U_1^2]}`

`J = {[0],[(gamma-1)/gamma(U_3 - gamma/2(U_2^2)/U_1)(dellnA')/(delx')],[0]}`

These equations are used in the discretization.

Non-Conservation Form

`(delrho')/(delt') = -rho'(delu')/(delx') - rho'u'(dellnA')/(delx') - u'(delrho')/(delx')`

`(delu')/(delt') = -u'(delu')/(delx') - 1/gamma((delT')/(delx') + (T')/(rho')(delrho')/(delx'))`

`(delT')/(delt') = -u'(delT')/(delx') -(gamma-1)T'((delu')/(delx') + u'(dellnA')/(delx'))`

Boundary Conditions

At the inlet, the velocity is the variable that is chosen to allow to float because, physically, the mass flow rate should be allowed to adjust to proper steady state. Thus allowing the velocity to vary seems appropriate for this adjustment. The velocity is therefore calculated from linear extrapolation of the nearest two internal points. Other parameters are kept fixed. i.e.

`u'_1 = 2u'_2 - u'_3`

`rho'_1 = 1`

`T'_1 = 1`

At the supersonic outlet, all flow field variables are alloed to float. And using linear extrapolation as well,

`u'_n = 2u'_(n-1) - u'_(n-2)`

`rho'_n = 2rho'_(n-1) - rho'_(n-2)`

`T'_n = 2T'_(n-1) - T'_(n-2)`

Initial Conditions

For the non-conservation form, the initial distribution of the flow field inside the nozzle parameters are given as

`rho' = 1 - 0.3246x'`

`T' = 1 - 0.2314x'`

`u' = (0.1 + 1.09x')T'^(1/2)`

For the conservation form, the initial distribution of the flow field inside the nozzle parameters are given as

`0 le x' le 0.5`

`rho' = 1`

`T' = 1`

`0.5 le x' le 1.5`

`rho' = 1 - 0.366(x' -0.5)`

`T' = 1 - 0.167(x' - 0.5)`

`1.5 le x' le 3`

`rho' = 0.634 - 0.3879(x' - 1.5)`

`T' = 0.833 - 0.3507(x' - 1.5)`

An intial constant mass flow of `0.59` was assumed throughout the nozzle so that

`u' = 0.59/(rho'A')`

where

`A' = 1 + 2.2(x' - 1.5)^2`

The figure below is from the reference book Computational Fluid Dynamics: the basics with applications by J. D. Anderson, Jr. It represents the plot of the Area profile equation above.

Adaptive Time Step

The time step is chosen so that the CFL condition is true. The courant number for 1-D simulation is given as

`C = Deltat(a +u)/(Deltax) le 1`

so that

`Delta t' = C(Deltax')/(a'+u')`

This equation means that the time for a sound wave to tavel from one grid point to the other must be less than or equal to 1. However, since `Deltax'` is fixed and `a' & u'` varies, it means that `Deltat_i` varies from one grid point to the next. As such a constant timestep is might not make the time -marching solution converge faster. The adaptive technique used here was to estimate the `Deltat` at all the grid points for every time step. The minimum is then chosen as the value of `Deltat` for that time step. i.e.

`Delta t^t = min(Deltat_1^t, Deltat_2^t, Deltat_3^t, ..., Deltat_(n-1)^t, Deltat_n^t)`

With this technique, the time-marching solutions follow the actual unsteady flow variations to give a time-accurate solution of the actual flow fields and avoids artificial time-warping consistent with the other technique that makes use of the local `Deltat^t`.

The MacCormack Method

In CFD, the MacCormack method is a widely used discretization scheme for the numerical solution of hyperbolic partial differential equations. This second-order finite difference method was introduced by Robert W. MacCormack in 1969.

The MacCormack method is elegant and easy to understand and program. It is a two-step variation of the Lax-Wendroff method but much simpler. The first step is the predictor step where a "provisional" value of `u` at time level `n+1` is estimated. The next step is the corrector step where the predicted value from the previous step is corrected.

The MacCormack method is well suited for nonlinear equations. For linear equations, the MacCormack scheme is equivalent to the Lax–Wendroff method. It is known to introduce dispersive errors (Gibbs phenomenon) in the region where the gradient is high but not diffusive errors in the solution.

MacCormack method is the discretization scheme used. It is a time-marching solution that starts with the initial conditions and continues until there is no change in the solutions any more. At the state where there is no change in solutions, the simulation is said to have reached steady state.

Three octave scripts were written, two of which simulated the conservation and the non-conservation forms using the time-marching technique. Each of these scripts could be run separately. But there is a third interactive script that is able to run both scripts for the sake of automation and comparison of results. Furthermore, the conservation form code was written in a modular form in order to avoid boilerplate codes and facilitate code reuse.

Several plots were generated for primitive variables in both transient and steady states for both forms. The transient variations were monitored at the nozzle throat while the steady state results were observed as a distribution within the nozzle.

The Codes and the plots are displayed below.

The General Code

fprintf("There are two scripts for the two forms of the supersonic nozzle simulation:");
fprintf("\n");
fprintf("conservative or non-conservative");
fprintf("\n");
fprintf("****************************************************************************");
fprintf("\n");
scriptname = input('what script do you want to run first? ', 's');
if strcmp(scriptname, 'conservative')
    clear all; close all; clc;
    fprintf("running conservative form first...");
    fprintf("\n");
    run('supersonic_nozzle_1D_C.m')
    rho_c = rho;
    v_c = v;
    T_c = T;
    P_c = P;
    Ma_c = Ma;
    mdot_c = mdot;
    
    clear rho v T P Ma mdot A C;
    clear i k n nt gamma;
    clear d* *throat;
    clear U* F* 
    %close all;
    
    fprintf("running non-conservative form now...");
    fprintf("\n");
    run('supersonic_nozzle_1D_NC.m')
    rho_nc = rho;
    v_nc = v;
    T_nc = T;
    P_nc = P;
    Ma_nc = Ma;
    mdot_nc = mdot;
    
    clear rho v T P Ma mdot A C;
    clear i k n nt  gamma;
    clear d* *throat;
    clear U* F* 

elseif strcmp(scriptname, 'non-conservative')

    clear all; close all; clc;
    fprintf("running non-conservative form first...");
    fprintf("\n");
    run('supersonic_nozzle_1D_NC.m')
    rho_nc = rho;
    v_nc = v;
    T_nc = T;
    P_nc = P;
    Ma_nc = Ma;
    mdot_nc = mdot;
    
    clear rho v T P Ma mdot A C;
    clear i k n nt gamma;
    clear d* *throat;
    clear U* F* 
    %close all;
    
    fprintf("running conservative form now...");
    fprintf("\n");
    run('supersonic_nozzle_1D_C.m')
    rho_c = rho;
    v_c = v;
    T_c = T;
    P_c = P;
    Ma_c = Ma;
    mdot_c = mdot;
    
    clear rho v T P Ma mdot A C;
    clear i k n nt gamma;
    clear d* *throat;
    clear U* F* 
endif
figure(7, 'position', [100, 100, 800, 400])
plot(x, mdot_c, x, mdot_nc, "linewidth", 2)
title("Comparing Mass Flow Variations for Conservative and Non-Conservative Form")
legend("conservation", "non-conservation")
xlabel("Nozzle Length")
ylabel("N-D massflow")
saveas(7, "NC-C_mdot_plots.png");

Non-Conservation Form Code

%clear all; close all; clc;

%inputs
n=31;
x=linspace(0,3,n);
dx = x(2) - x(1);
C = 0.5;
dt_array = [1, 50, 100, 150, 200, 700];

%calculate initial profiles
rho = 1-0.3146*x;
T = 1-0.2314*x;
v=(0.1+1.09*x).*(T.^0.5);
Ma = v./sqrt(T);

%%calculate initial profiles
%for i = 1:length(x)
 % if x(i) <= 0.5
 %   rho(i) = 1;
  %  T(i) = 1;
  %elseif x(i) > 0.5 && x(i) <=1.5
  %  rho(i) = 1-0.366*(x(i)-0.5);
  %  T(i) = 1-0.167*(x(i)-0.5);
  %else
  %  rho(i) = 0.634-0.3879*(x(i)-1.5);
   % T(i) = 0.833-0.3507*(x(i)-1.5);
  %endif
%end
A=1+2.2*(x-1.5).^2; 
ithroat = find(A==min(A));

%v=0.59./(rho.*A);
%Ma = v./sqrt(T);

gamma = 1.4;


%timesteps
nt = 1400;

for k = 1:nt
  rho_old =rho;
  v_old = v;
  T_old = T;
  
  % calculate dt
  dt = min(C*dx./(sqrt(T) + v));
  
  %predictor
  for j = 2:n-1
    dvdx = (v(j+1) - v(j))/dx;
    drhodx = (rho(j+1) - rho(j))/dx;
    dlogAdx = (log(A(j+1))-log(A(j)))/dx;
    dTdx = (T(j+1)-T(j))/dx;
    

    
    %continuity equation
    drhodt_p(j) = -rho(j)*dvdx - rho(j)*v(j)*dlogAdx - 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)*dlogAdx);
    
   
    
    % solution update
    rho(j) += drhodt_p(j)*dt;
    v(j) += dvdt_p(j)*dt;
    T(j) += dTdt_p(j)*dt;
    
   end
   
   %corrector
   for j = 2:n-1
     dvdx = (v(j) - v(j-1))/dx;
     drhodx = (rho(j) - rho(j-1))/dx;
     dlogAdx = (log(A(j))-log(A(j-1)))/dx;
     dTdx = (T(j)-T(j-1))/dx;
    
     %continuity equation
     drhodt_c(j) = -rho(j)*dvdx - rho(j)*v(j)*dlogAdx - v(j)*drhodx;
    
     %momentum equation
     dvdt_c(j) = -v(j)*dvdx - (1/gamma)*(dTdx + (T(j)/rho(j))*drhodx);
    
     %energy equation
     dTdt_c(j) = -v(j)*dTdx - (gamma-1)*T(j)*(dvdx + v(j)*dlogAdx);
    
    end
   
    %comput 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;
      T(i) = T_old(i) + dTdt(i)*dt;
    end

    % Apply the 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);
    
    P = rho.*T;
    Ma = v./sqrt(T);
    mdot = rho.*v.*A;
    
    %plot_variables(rho, v, T, P, mdot, Ma, x, k, ithroat)
    
  rhothroat(k) = rho(ithroat);
  vthroat(k) = v(ithroat);
  Tthroat(k) = T(ithroat);
  Pthroat(k) = P(ithroat);
  Mathroat(k) = Ma(ithroat);
  mdotthroat(k) = mdot(ithroat);
  
  figure(4, 'position', [100, 100, 1000, 1000])
  cla
  subplot(6,1,1)
  plot(rhothroat, "linewidth", 2, "color", "g")
  title("Transient Variations of Primitive Parameters at the Nozzle Throat for the Non-Conservative Form")
  xlabel("Time steps")
  ylabel("Density")
  subplot(6,1,2)
  plot(vthroat, "linewidth", 2, "color", "y")
  xlabel("Time steps")
  ylabel("Velocity")
  subplot(6,1,3)
  plot(Tthroat, "linewidth", 2, "color", "r")
  xlabel("Time steps")
  ylabel("Temperature")
  subplot(6,1,4)
  plot(Pthroat, "linewidth", 2, "color", "k")
  xlabel("Time steps")
  ylabel("Pressure")
  subplot(6,1,5)
  plot(Mathroat, "linewidth", 2, "color", "m")
  xlabel("Time steps")
  ylabel("Mach Number")
  subplot(6,1,6)
  plot(mdotthroat, "linewidth", 2, "color", "c")
  xlabel("Time steps")
  ylabel("mass flow rate")
  if k==1400
    hold on
    saveas(4, "transient_plots_NC.png");
    hold off
  endif
  
  figure(5, 'position', [100, 100, 1000, 1000])
  subplot(6,1,1)
  plot(x, rho, "linewidth", 2, "color", "g")
  title("Steady State Distribution of Primitive Parameters within the Nozzle for the Non-Conservative Form")
  xlabel("nozzle Length")
  ylabel("Density")
  subplot(6,1,2)
  plot(x, v, "linewidth", 2, "color", "y")
  xlabel("nozzle Length")
  ylabel("Velocity")
  subplot(6,1,3)
  plot(x, T, "linewidth", 2, "color", "r")
  xlabel("nozzle Length")
  ylabel("Temperature")
  subplot(6,1,4)
  plot(x, P, "linewidth", 2, "color", "k")
  xlabel("nozzle Length")
  ylabel("Pressure")
  subplot(6,1,5)
  plot(x, Ma, "linewidth", 2, "color", "m")
  xlabel("nozzle Length")
  ylabel("Mach Number")
  subplot(6,1,6)
  plot(x, mdot, "linewidth", 2, "color", "c")
  xlabel("nozzle Length")
  ylabel("mass flow rate")
  if k==1400
    hold on
    saveas(5, "steady_state_plots_NC.png");
    hold off
  endif
  %legend("v","rho","T","Ma")
  
  figure(6, 'position', [100, 100, 800, 400])
  if ismember(k, dt_array)==1
    plot(x, mdot, 'linewidth', 2, 'DisplayName', ['\Deltat = ' num2str(k)])
    title("Mass Flow Distribution for Different Timesteps within the Nozzle for the Non-Conservative Form")
    xlabel("Nozzle Length")
    ylabel("N-D massflow")
    hold on 
  endif
  pause(0.003)
end
legend('show')
saveas(6, "mdot_plots_NC.png");
hold off;

The plots above are steady state distributions of the primitive variables inside the supersonic nozzle for the time-marching solution of the non-conservative form of governing equations. This state was achieved after 1400 time steps. They are consistent with the those derived from analytic results in the reference book.

The plots above are the transient variations of the primitive parameters observed at the nozzle throat also for the time-marching solution of the non-conservation form of the governing equations. The x-axis is the time axis and it started from `0` which represents the initial conditions. It is observed that the largest changes took place at the early times after which the steady state is approached asymptotically. This numerical time-marching procedure converged to the analytical steady state answer.

The plot above is another visualization into the transient variations and approach to steady state. It is the plot of non-dimensional mass flow distribution at different time steps. Radical changes are aobserved in the distributions from `Deltat = 0` up until the time `Deltat = 700` where the distribution became uniform indicating the mass flow distribution has converged to a constant, steady state value throughout the nozzle.

Conservation Form Codes

function [F1 F2 F3] = flux_terms_F(U1, U2, U3, gamma)
  
  F1 = U2;

  F2a = (U2.^2)./U1;
  F2b = ((gamma-1)/gamma).*(U3 - (gamma/2).*F2a);
  F2 =  F2a + F2b;

  F3a = gamma.*(U2.*U3)./U1;
  F3b = ((gamma^2-gamma)/2).*(U2.^3)./(U1.^2);
  F3 = F3a - F3b;
end

function J2 = source_term(U1, U2, U3, gamma, dlogAdx)
  J2a = (U3 - (gamma/2)*(U2.^2)/U1);
  J2 = (gamma-1)/gamma*J2a.*dlogAdx;
end

function [U1 U2 U3 dU1dt_p dU2dt_p dU3dt_p] = predictor(F1, F2, F3, A, U1, U2, U3, gamma, n, dt, dx)
  
  for j = 2:n-1
    
    dF1dx(j) = (F1(j+1) - F1(j))/dx;
    dF2dx(j) = (F2(j+1) - F2(j))/dx;
    dF3dx(j) = (F3(j+1) - F3(j))/dx;
    dlogAdx(j) = (log(A(j+1)/A(j)))/dx;
    
    J2(j) = source_term(U1(j), U2(j), U3(j), gamma, dlogAdx(j));
    
    %continuity equation
    dU1dt_p(j) = -dF1dx(j);
    
    %momentum equation
    dU2dt_p(j) = -dF2dx(j) + J2(j);
    
    %energy equation
    dU3dt_p(j) = -dF3dx(j);
    
    
    % Predicted values for U
    U1(j) += dU1dt_p(j)*dt;
    
    U2(j) += dU2dt_p(j)*dt;
    
    U3(j) += dU3dt_p(j)*dt;
    
   end
end

function [dU1dt_c dU2dt_c dU3dt_c] = corrector(F1, F2, F3, A, U1, U2, U3, gamma, n, dt, dx)
  
  for j = 2:n-1
    
    dF1dx(j) = (F1(j) - F1(j-1))/dx;
    dF2dx(j) = (F2(j) - F2(j-1))/dx;
    dF3dx(j) = (F3(j) - F3(j-1))/dx;
    dlogAdx(j) = (log(A(j)/A(j-1)))/dx;
    
    J2(j) = source_term(U1(j), U2(j), U3(j), gamma, dlogAdx(j));
    
    %continuity equation
    dU1dt_c(j) = -dF1dx(j);
    
    %momentum equation
    dU2dt_c(j) = -dF2dx(j) + J2(j);
    
    %energy equation
    dU3dt_c(j) = -dF3dx(j);
       
   end
end

function [rho v T P Ma mdot] = primitive_parameters(U1, U2, U3, A, gamma)
  rho = U1./A;
  v = U2./U1;
  T = (gamma-1)*((U3./U1) - (gamma/2)*(U2./U1).^2);
  % Mach Number
  Ma = v./sqrt(T);
  P = rho.*T;
  mdot = rho.*A.*v;
end

%%inputs
n=31;
x=linspace(0,3,n);
dx = x(2) - x(1);
gamma = 1.4;
C = 0.5;
dt_array = [1, 50, 100, 150, 200, 700];

%timesteps
nt = 1400;

%%calculate initial profiles
for i = 1:length(x)
  if x(i) <= 0.5
    rho(i) = 1;
    T(i) = 1;
  elseif x(i) > 0.5 && x(i) <=1.5
    rho(i) = 1-0.366*(x(i)-0.5);
    T(i) = 1-0.167*(x(i)-0.5);
  else
    rho(i) = 0.634-0.3879*(x(i)-1.5);
    T(i) = 0.833-0.3507*(x(i)-1.5);
  endif
end
A=1+2.2*(x-1.5).^2; 
ithroat = find(A==min(A));
 
v=0.59./(rho.*A);
Ma = v./sqrt(T);


U1 = rho.*A;
U2 = rho.*A.*v;
U3 = rho.*A.*((T./(gamma-1)) + ((gamma/2).*v.^2));


for k = 1: nt
  U1_old =U1;
  U2_old = U2;
  U3_old = U3;
  
  %Corrected values for F
  [F1 F2 F3] = flux_terms_F(U1, U2, U3, gamma);
  
  % calculate dt
  dt = min(C*dx./(sqrt(T) + v));
  
  %predictor
  [U1 U2 U3 dU1dt_p dU2dt_p dU3dt_p] = predictor(F1, F2, F3, A, U1, U2, U3, gamma, n, dt, dx);
    
  %Predicted values for F
  [F1 F2 F3] = flux_terms_F(U1, U2, U3, gamma);
    
  %corrector
  [dU1dt_c dU2dt_c dU3dt_c] = corrector(F1, F2, F3, A, U1, U2, U3, gamma, n, dt, dx);
   
  %compute average time derivative
  dU1dt = 0.5*(dU1dt_p + dU1dt_c);
  dU2dt = 0.5*(dU2dt_p + dU2dt_c);
  dU3dt = 0.5*(dU3dt_p + dU3dt_c);
   
  %Final Solution update
  for i = 2:n-1
    U1(i) = U1_old(i) + dU1dt(i)*dt;
    U2(i) = U2_old(i) + dU2dt(i)*dt;
    U3(i) = U3_old(i) + dU3dt(i)*dt;
      
  end
    
  % Apply the Boundary Conditions
  %inlet
  U1(1) = rho(1)*A(1);
  U2(1) = 2*U2(2) - U2(3);
  v(1) = U2(1)/U1(1);
  U3(1) = U1(1)*(T(1)/(gamma-1) + (gamma/2)*v(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);
    
  %primitive parameters
  [rho v T P Ma mdot] = primitive_parameters(U1, U2, U3, A, gamma);
  rhothroat(k) = rho(ithroat);
  vthroat(k) = v(ithroat);
  Tthroat(k) = T(ithroat);
  Pthroat(k) = P(ithroat);
  Mathroat(k) = Ma(ithroat);
  mdotthroat(k) = mdot(ithroat);
  
  figure(1, 'position', [100, 100, 1000, 1000])
  cla
  subplot(6,1,1)
  plot(rhothroat, "linewidth", 2, "color", "g")
  title("Transient Variations of Primitive Parameters at the Nozzle Throat for the Conservative Form")
  xlabel("Time steps")
  ylabel("Density")
  subplot(6,1,2)
  plot(vthroat, "linewidth", 2, "color", "y")
  xlabel("Time steps")
  ylabel("Velocity")
  subplot(6,1,3)
  plot(Tthroat, "linewidth", 2, "color", "r")
  xlabel("Time steps")
  ylabel("Temperature")
  subplot(6,1,4)
  plot(Pthroat, "linewidth", 2, "color", "k")
  xlabel("Time steps")
  ylabel("Pressure")
  subplot(6,1,5)
  plot(Mathroat, "linewidth", 2, "color", "m")
  xlabel("Time steps")
  ylabel("Mach Number")
  subplot(6,1,6)
  plot(mdotthroat, "linewidth", 2, "color", "c")
  xlabel("Time steps")
  ylabel("mass flow rate")
  if k==1400
    hold on
    saveas(1, "transient_plots_C.png");
    hold off
  endif
  
  figure(2, 'position', [100,100, 1000, 1000])
  subplot(6,1,1)
  plot(x, rho, "linewidth", 2, "color", "g")
  title("Steady State Distribution of Primitive Parameters within the Nozzle for the Conservative Form")
  xlabel("nozzle Length")
  ylabel("Density")
  subplot(6,1,2)
  plot(x, v, "linewidth", 2, "color", "y")
  xlabel("nozzle Length")
  ylabel("Velocity")
  subplot(6,1,3)
  plot(x, T, "linewidth", 2, "color", "r")
  xlabel("nozzle Length")
  ylabel("Temperature")
  subplot(6,1,4)
  plot(x, P, "linewidth", 2, "color", "k")
  xlabel("nozzle Length")
  ylabel("Pressure")
  subplot(6,1,5)
  plot(x, Ma, "linewidth", 2, "color", "m")
  xlabel("nozzle Length")
  ylabel("Mach Number")
  subplot(6,1,6)
  plot(x, mdot, "linewidth", 2, "color", "c")
  xlabel("nozzle Length")
  ylabel("mass flow rate")
  %legend("v","rho","T","Ma")
  if k==1400
    hold on
    saveas(2, "steady_state_plots_C.png");
    hold off
  endif
  
  figure(3, 'position', [100, 100, 800, 400])
  if ismember(k, dt_array)==1
    plot(x, mdot, 'linewidth', 2, 'DisplayName', ['dt = ' num2str(k)])
    title("Mass Flow Distribution for Different Timesteps within the Nozzle for the Conservative Form")
    xlabel("Nozzle Length")
    ylabel("N-D massflow")
    hold on 
    
  endif
  pause(0.003)
end
legend('show')
saveas(3, "mdot_plots_C.png");
hold off;

The plots above are steady state distributions of the primitive variables inside the supersonic nozzle for the time-marching solution of the conservative form of governing equations. This state was achieved after 1400 time steps. They are consistent with the those derived from analytic results in the reference book.

The plots above are the transient variations of the primitive parameters observed at the nozzle throat also for the time-marching solution of the conservation form of the governing equations. The x-axis is the time axis and it started from `0` which represents the initial conditions. It is observed that the largest changes took place at the early times after which the steady state is approached asymptotically. This numerical time-marching procedure converged to the analytical steady state answer.

The plot above is another visualization into the transient variations and approach to steady state. It is the plot of non-dimensional mass flow distribution at different time steps. Apart from the high inlet and outlet flows at `Deltat = 50`, there is no much radical changes are observed in the distributions at different times. the mass flow distributions are near close to the converged steady state solutions at different times. This is partly due to the constant initial mass flow assumption.

Comparing the Conservation and Non-Conservation Form

From the plot for both non-conservation and conservation forms, one observes that there plots are similar for the primitive parameter for both steady state and transient. This should be expected since the present problem is an isentropic subsonic-supersonic flow through the nozzle which does not involve any shock capturing.

Comparing the steady state mass flow distribution, however, one can draw certain conclusions, all be it slight.

• There is no clear superiority in terms of accuracy of results. The non-conservation form produces more accurate primitive variables while the conservation form produces more accurate flux.
• The mass flow distribution from the conservation form is better than from the non-conservation.
• The non-conservation form does better in terms of quality of numerical algorithm as it yields smaller residue.
• conservation form is more computationally intensive due to the fact that primitive variables have to be decoded at every timestep. The non-conservation form is not like that.

Establishing Grid Independence

Generally in CFD analysis, it is essential to ensure that solutions obtain are independent of the number of grids used in domain discretization. The way this is done is to double the number of grids and run another simulation on the new grid. If the solutions do not change, then one has established grid independence.

For the solutions above, 31 grid points have been used. To establish grid independence, the number of grid is doubled to `61` grid points and the simulation run again. The non-conservation form of governing equation is used in this study. The table below presents the results tand compare The last row calculated the percentage absolute difference as

`"%difference" = (V_(31pts) - V_(61pts))/V_(31pts) times 100`

	`rho'^star`	`T'^star`	`P'^star`	`Ma'^star`
Exact Analytical Solution	`0.634`	`0.833`	`0.528`	`1.00`
Conservation Form (31 points)	`0.69206`	`0.85842`	`0.59408`	`0.91524`
Non-conservation Form (31 points)	`0.63869`	`0.83646`	`0.53424`	`0.99937`
Non-conservation Form (61 points)	`0.63756`	`0.83539`	`0.53261`	`0.99970`
%Absolute Difference	`0.177%`	`0.128%`	`0.305%`	`0.033%`

From the table above, the steady state values for the primitive parameters on `31` grid points were within `le 1%` of those on `61` grid points. Thus, it is safe to conclude that the simulation reached grid independence on `31` grid points.