Numerical Solution of Quasi 1D Nozzle Flows: MacCormack's Technique

  • I am going to write a program in MATLAB for solving Quasi one dimensional nozzle flows using MacCormack's Technique implementing both of the conservative and non-conservative forms of governing equations and perform the grid dependence test. Note, I will use the set of initial conditions provided in Anderson's textbook (chapter `7`) along with the nozzle's shape and length details.

 

Problem Setup:
1. Form finite-difference expressions applicable to MacCormack's technique.
2. Chosen Courant number is `0.5`.
3. Grid points used for project purpose is `30`.
4. For grid dependence test another grid point used is `60`.


Governing flow equations

  1. Non-conservative form:

    1. Continuity equation: `(delrho)/(delt) = -rho(delV)/(delx) - rhoV(del(lnA))/(delx) - V(delrho)/(delx)`

    2. Momentum equation: `(delV)/(delt) = -V(delV)/(delx) - 1/gamma((delT)/(delx) + T/rho (delrho)/(delx))`

    3. Energy equation: `(delT)/(delt) = -V(delT)/(delx) - T(gamma-1)[(delV)/(delx) + V(del(lnA))/(delx)]`

  2. Conservative form:

    1.
    Continuity equation: `(delU_1)/(delt) = - (delF_1)/(delx)`

    2. Momentum equation: `(delU_2)/(delt) = - (delF_2)/(delx) + J_2`

    3.
    Energy equation: `(delU_3)/(delt) = - (delF_3)/(delx)`

    In the above equations, 

    Solution vectors are defined as:
    `U_1 = rhoA`, `U_2 = rhoAV`, `U_3 = rho(T/(gamma-1)+gamma/2V^2)A`

    Flux vectors are defined as:
    `F_1 = U_2`, `F_2 = (U_2^2)/(U_1) + 1/gammaU_1(gamma-1)[U_3/U_1-gamma/2(U_2/U_1)^2]`, `F_3 = gamma(U_2U_3)/(U_1)-(gamma(gamma-1))/2(U_2^3)/U_1^2`

    Source term is defined as:
    `J_2 = 1/gammarhoT(delA)/(delx)`

 
NOTE: All of the above variables are dimensionless.


MATLAB PROGRAM:

Function files

1. 'fd_dx.m':

% fd_dx - calculates forward difference of a input function

function f1 = fd_dx(f,dx,i)

f1 = (f(i+1) - f(i))/dx;

end

2. 'rd_dx.m':

% rd_dx - calculates reward differences of the input function

function f1 = rd_dx(f,dx,i)

f1 = (f(i) - f(i-1))/dx;

end

 3. 'pred_corr.m':

% pred_corr - calculates predictor/corrector step 

function [f1, f2, f3] = pred_corr(v, t, a, r, y, rho, V, T, dx, i)

term1 = -rho(i) * v;
term2 = -rho(i) * V(i) * a;
term3 = -V(i) * r;
term4 = -V(i) * v;
term5 = -y^-1 * (t + T(i)/rho(i)  * r);
term6 = -V(i) * t;
term7 = -(y-1) * T(i) * (v + V(i) * a);

f1 = term1 + term2  + term3;
f2 = term4 + term5;
f3 = term6 + term7;

end

4. 'conservative.m':

%conservative - Solves conservative form of govering equations 

function [m, M, V, T, p, rho, temperature, density, pressure, mach_n] = conservative(x, dx, n, C, A, nt, y, throat)

%% Initial conditions
for i = 1 : length(x)
     if x(i) >= 0 && 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);
    elseif x(i) >= 1.5 && x(i) <= 3.5 
        rho(i) = 0.634 - 0.3879 * (x(i) - 1.5);
        T(i) = 0.833 - 0.3507 * (x(i) - 1.5);
    end
end 

% Defining velocity and pressure
V = 0.59 ./ (rho .* A); 
p = rho .* T;

% Defining solution vectors
U1 = rho .* A;
U2 = rho .* A .* V;
U3 = rho .* A .* (T / (y - 1) + y / 2 * V.^2); 

% Iteration loop
for k = 1 : nt
     
     % Copies of solution vectors
     U1_old = U1; 
     U2_old = U2; 
     U3_old = U3;
     
     % dt calculations
     dt = min(C * dx ./ (sqrt(T) + V));
     
     % Flux vectors     
     F1 = U2;
     F2=(U2.^2 ./ U1) + ((y - 1) / y) * (U3 - (y * U2.^2) ./ (2 * U1));
     F3 =(y * U2 .* U3 ./ U1) - (y * (y - 1) / 2) * (U2.^3 ./ U1.^2); 
     
     % Predictor step
     for i = 2 : n - 1
          J2(i) = 1 / y * rho(i) .* T(i) * fd_dx(A, dx, i);
          dU1dt_p(i) = -fd_dx(F1, dx, i);
          dU2dt_p(i) = -fd_dx(F2, dx, i) + J2(i);
          dU3dt_p(i) = -fd_dx(F3, dx, i);
          
          % Updating solution vectors
          U1(i)= U1(i) + dU1dt_p(i) * dt;
          U2(i)= U2(i) + dU2dt_p(i) * dt;
          U3(i)= U3(i) + dU3dt_p(i) * dt;
     end
     
     % Updation of primitive variables and flux vectors
      rho = U1 ./ A;
      T = (y - 1) * (U3 ./ U1 - (y / 2) * (U2 ./ U1).^2);
      p = rho .* T;
      F1 = U2;
      F2=(U2.^2 ./ U1) + ((y - 1) / y) * (U3 - (y * U2.^2) ./ (2 * U1));
      F3 =(y * U2 .* U3 ./ U1) - (y * (y - 1) / 2) * (U2.^3 ./ U1.^2);
      
      % Corrector step
      for j = 2 : n - 1
           J2(j) = 1 / y * rho(j) .* T(j) * rd_dx(A, dx, j);       % Updated source term
           dU1dt_c(j) = -rd_dx(F1, dx, j);
           dU2dt_c(j) = -rd_dx(F2, dx, j) + J2(j);
           dU3dt_c(j) = -rd_dx(F3, dx, j);
      end
     
     % Averaging corrector step and predictor step
      dU1dt_a = 0.5 * (dU1dt_p + dU1dt_c);
      dU2dt_a = 0.5 * (dU2dt_p + dU2dt_c);
      dU3dt_a = 0.5 * (dU3dt_p + dU3dt_c);
      
      % Updating solution vectors after averaging
      for m = 2 : n - 1 
           U1(m) = U1_old(m) + dU1dt_a(m) * dt;
           U2(m) = U2_old(m) + dU2dt_a(m) * dt;
           U3(m) = U3_old(m) + dU3dt_a(m) * dt;
      end   

      % Boundary conditions 
      % Inlet
      U2(1) = 2*U2(2) - U2(3);

      % 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);

      % Updation of flow field variables
      rho = U1 ./ A;
      V = U2 ./ U1;
      T = (y - 1) * (U3 ./ U1 - (y / 2) * (U2 ./ U1).^2);
      p = rho .* T;
      
      % Defining mass flow and mach number
      m = rho .* A .* V;
      M = V ./ sqrt(T);
      
      % Capturing the values at throat for each iterations
      mach_n(k) = M(throat);
      pressure(k) = p(throat);
      density(k) = rho(throat);
      temperature(k) = T(throat);
      
      % Ploting mass flow rate at different times
      if n == 31
      
          v = figure(1);
          hold on
          if k == 100
             plot(x,m, 'b', 'linewidth', 2)
          elseif k  == 200
            plot(x, m, 'r', 'linewidth', 2)
          elseif k  == 300
            plot(x, m, 'm', 'linewidth', 2)
          elseif k  == 550
            plot(x, m, 'k', 'linewidth', 2)  
            title('Non - Dimensional Mass Flow Distribution at different Time Steps', 'FontAngle', 'italic', 'FontSize', 13, 'FontWeight', 'bold')
            legend('100Deltat', '200Deltat', '300Deltat', '550Deltat')
            ylim([0.5 0.7])
            xlabel('x/L', 'FontWeight', 'bold', 'FontSize', 11)
            ylabel('rhoVA/rho_{o}V_{o}A_{0}', 'FontWeight', 'bold', 'FontSize', 11) 
            grid minor
            saveas(v, 'cons_mass_flow.png')     
          end
 
     end
      
end     

end

5. 'non_conservative.m':

%nonconservative - Solves non-conservative form of govering equations 

function [m, M, V, T, p, rho, temperature, density, pressure, mach_n] = non_conservative(x, dx, n, C, A, nt, y, throat)

%% Initial conditions for the flow field variables
rho = 1 - 0.3146 * x;                                                   % Density 
T = 1 - 0.2314 * x;                                                      % Temprature
V = (0.1 + 1.09 * x) .* T.^0.5;                                     % Velocity

for k = 1 : nt  

     % Copy of Field Variables
     rho_old = rho;
     V_old = V;
     T_old = T;     
     
     % Time step
     dt = min(C * dx ./ (T.^0.5 + V));
     
     % Predictor step
     for i = 2:n-1
          vf = fd_dx(V,dx,i); tf = fd_dx(T,dx,i); af = fd_dx(log(A),dx,i); rf = fd_dx(rho,dx,i);
          [drdt_p(i), dVdt_p(i), dTdt_p(i)] = pred_corr(vf, tf, af, rf, y, rho, V, T, dx, i);
          
          %solution update
          rho(i) = rho(i) + drdt_p(i)*dt;
          V(i) = V(i) + dVdt_p(i)*dt;
          T(i) = T(i) + dTdt_p(i)*dt;
     end
     
     % Corrector step
     for j = 2:n-1
          vr = rd_dx(V,dx,j); tr = rd_dx(T,dx,j); ar = rd_dx(log(A),dx,j); rr = rd_dx(rho,dx,j); 
          [drdt_c(j), dVdt_c(j), dTdt_c(j)] = pred_corr(vr, tr, ar, rr, y, rho, V, T, dx, j);
     end
     
     % Taking average of corrector and predictor step
     drdt_a = 0.5 * (drdt_c + drdt_p);
     dVdt_a = 0.5 * (dVdt_c + dVdt_p);
     dTdt_a = 0.5 * (dTdt_c + dTdt_p);
     
     % Flow field variables at next time step 
     for s = 2:n-1
          rho(s) = rho_old(s) + drdt_a(s) * dt;
          V(s) = V_old(s) + dVdt_a(s) * dt;
          T(s) = T_old(s) + dTdt_a(s) * dt;
     end  
     
     % Boundary Conditions    
     % subsonic inflow
     V(1) = 2*V(2) - V(3);
     % supersonic outflow
     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);
     
      % Defining pressure, mass flow and mach number
      p = rho .* T;
      m = rho .* A .* V;
      M = V ./ sqrt(T);

      % Capturing values at throat for each iterations
      mach_n(k) = M(throat);
      pressure(k) = p(throat);
      density(k) = rho(throat);
      temperature(k) = T(throat);
      
      % Mass flow at different time steps
      if n == 31 
      
          v = figure(2);
          hold on
          if k  == 50
            plot(x, m, 'r', 'linewidth', 2)
          elseif k  == 100
            plot(x, m, 'm', 'linewidth', 2)
          elseif k  == 150
            plot(x, m, 'c', 'linewidth', 2)
          elseif k  == 700
            plot(x, m, 'k', 'linewidth', 2) 
            title('Non - Dimensional Mass Flow Distribution at different Time Steps', 'FontAngle', 'italic', 'FontSize', 13, 'FontWeight', 'bold')
            legend('50Deltat', '100Deltat', '150Deltat', '700Deltat')       
            xlabel('x/L', 'FontWeight', 'bold', 'FontSize', 11)
            ylabel('rhoVA/rho_{o}V_{o}A_{0}', 'FontWeight', 'bold', 'FontSize', 11) 
            grid minor
            saveas(v, 'non_cons_mass_flow.png')
          end  
      
      end     
      
end
    
end


The driver script 'quasi_one_D_nozzle_flows.m' is as follows

% Quasi 1D Nozzle Flows

clear all
close all
clc

%% Setup
L = 3;                                               %Length of the nozzle
n = input('Choose number of grid points: ') + 1;     % Grid points - choose it manually
C =0.5;                                              % Courant number 
x = linspace(0,L,n);                                 % Creating 1D grid
dx = L/(n-1);                                        % Adjecent grid spacing
y = 1.4;                                             % Gamma                                                  
nt = 1400;                                           % Iterations / Total time steps
throat = (n-1)/2 + 1;                                % Throat point
A = 1 + 2.2 * (x - 1.5).^2;                          % Area

%% Conservative form
tic
[Cm, CM, CV, CT, Cp, Crho, Ctemperature, Cdensity, Cpressure, Cmach_n] = conservative(x, dx, n, C, A, nt, y, throat);
cons = toc;
fprintf('nSolution time for conservative form is %0.3g secondsn', cons)

%% Non-Conservative form
tic
[Nm, NM, NV, NT, Np, Nrho, Ntemperature, Ndensity, Npressure, Nmach_n] = non_conservative(x, dx, n, C, A, nt, y, throat);
non_cons = toc;
fprintf('Solution time for non - conservative form is %0.3g secondsn', non_cons)

%% Ploting the above data
u = figure(3);
subplot(411)
hold on 
plot(Cmach_n, 'b', 'linewidth', 2)
plot(Nmach_n, '--k', 'linewidth', 2)
s = legend('Conservative Form', 'Non-Conservative Form');
set(s, 'Location', 'northeastoutside')
title('Timewise Variation at the Nozzle Throat', 'FontAngle', 'italic', 'FontSize', 13, 'FontWeight', 'bold')
ylabel('M', 'FontSize', 11, 'FontWeight', 'bold')
grid minor

subplot(412)
hold on 
plot(Cpressure, 'm', 'linewidth', 2)
plot(Npressure, '--k', 'linewidth', 2)
e = legend('Conservative Form', 'Non-Conservative Form');
set(e, 'Location', 'northeastoutside')
ylabel('p/p_{0}', 'FontSize', 11, 'FontWeight', 'bold')
grid minor

subplot(413)
hold on
plot(Cdensity, 'r', 'linewidth', 2)
plot(Ndensity, '--k', 'linewidth', 2)
r = legend('Conservative Form', 'Non-Conservative Form');
set(r, 'Location', 'northeastoutside')
ylabel('rho/rho_{0}', 'FontSize', 11, 'FontWeight', 'bold')
grid minor

subplot(414)
hold on
plot(Ctemperature, 'g', 'linewidth', 2)
plot(Ntemperature, '--k', 'linewidth', 2)
t = legend('Conservative Form', 'Non-Conservative Form');
set(t, 'Location', 'northeastoutside')
ylabel('T/T_{o}', 'FontSize', 11, 'FontWeight', 'bold')
xlabel('Number of Iterations', 'FontSize', 11, 'FontWeight', 'bold')
grid minor
saveas(u,'nozzle_throat.png')


b = figure(4);
subplot(411)
hold on
plot(x, CM, 'b', 'linewidth', 2)
plot(x, NM, '-k+', 'linewidth', 1, 'markersize', 5, 'markerfacecolor', 'k')
o = legend('Conservative Form', 'Non-Conservative Form');
set(o, 'Location', 'northeastoutside')
grid minor
hold on
title('Qualitative Aspects of Quasi - 1D Nozzle Flow', 'FontSize', 13, 'FontWeight', 'bold', 'FontAngle', 'italic')
ylabel('Mach number', 'FontWeight', 'bold', 'FontSize', 11)

subplot(412)
hold on
plot(x, CT, 'g', 'linewidth', 2)
plot(x, NT, '-k+', 'linewidth', 1, 'markersize', 5, 'markerfacecolor', 'k')
N = legend('Conservative Form', 'Non-Conservative Form');
set(N, 'Location', 'northeastoutside')
ylabel('Temperature ratio', 'FontWeight', 'bold', 'FontSize', 11)
grid minor

subplot(413)
hold on
plot(x, Crho, 'r', 'linewidth', 2)
plot(x, Nrho, '-k+', 'linewidth', 1, 'markersize', 5, 'markerfacecolor', 'k')
D = legend('Conservative Form', 'Non-Conservative Form');
set(D, 'Location', 'northeastoutside')
ylabel('Density ratio', 'FontWeight', 'bold', 'FontSize', 11)
grid minor

subplot(414)
hold on
plot(x, Cp, 'm', 'linewidth', 2)
plot(x, Np, '-k+', 'linewidth', 1, 'markersize', 5, 'markerfacecolor', 'k')
F = legend('Conservative Form', 'Non-Conservative Form');
set(F, 'Location', 'northeastoutside')
ylabel('Pressure ratio', 'FontWeight', 'bold', 'FontSize', 11)
xlabel('Non-Dimensional Length of the Nozzle', 'FontWeight', 'bold', 'FontSize', 11)
grid minor
saveas(b, 'qualitative_aspects.png')

q = figure(5);
hold on
plot(x, Cm, 'b', 'linewidth', 2)
plot(x, Nm, 'r', 'linewidth', 2)
line([0 3], [0.579 0.579], 'color', 'm', 'linewidth', 2.5)
Z = legend('Conservative Form', 'Non-Conservative Form', 'Exact Solution');
set(Z, 'Location', 'north')
grid minor
xlabel('x/L', 'FontWeight', 'bold', 'FontSize', 11)
ylabel('rhoVA/rho_{o}V_{o}A_{0}', 'FontWeight', 'bold', 'FontSize', 11)
title('Non-Dimensional Mass Flow Distribution', 'FontWeight', 'bold', 'FontSize', 13, 'FontAngle', 'italic')
saveas(q, 'mass_flow_comparission.png')

 

Results

Comparision of steady-state results; conservative vs. non-conservative form - 


                                               `rho/(rho_o)`            `T/(T_o)`          `p/(p_o)`           `M`

Exact analytical solution:
     `0.634`        `0.833`       `0.528`      `1.000`
                                                                                                                

Non-conservative form:       `0.639`        `0.836`       `0.534`       `0.999`
       (30 points)                                                                                        
                                                                                                               
Non-conservative form:       `0.638`        `0.835`       `0.533`      `1.000`
      (60 points) 

Conservative form:               `0.650`        `0.840`       `0.546`      `0.982`
      (30 points)           
                   
                                                                                          
Conservative form:               `0.638`        `0.835`       `0.533`      `0.999`      
     
(60 points)

 

Command window display -

For `30` points:

Solution time for conservative form is 22.9 seconds
Solution time for non - conservative form is 34.7 seconds


For `60` points:

Solution time for conservative form is 53.7 seconds
Solution time for non - conservative form is 87.9 seconds



Plots - 

  1. For `30` points:

    1. Qualitative aspects:


    2. Timewise variation at the throat:


    3. Mass flow distribution at different times: Conservative form:


    4. Mass flow distribution at different times: Conservative form:


    5. Mass flow rate comparison:


  2. For `60` points:

    6.
     Mass flow rate comparison plot - 

 

Conclusions

  1. For `30` points, as seen from the above plots `4` and `5` the minimum number of cycles for convergence in approximately `550` for conservative form and approximately `700` for non-conservative form.

  2. It totally depends on what circumstances you are performing the simulation. If a very minor error is not going to create an issue to your task then `30` grid points are enough else you can increase the number of grid points for better yield of the solution.

  3. In terms of accuracy, there is no clear supremacy in either of the forms. Though the results for each of the forms are certainly satisfactory.


  4. As from the above plot `5` and `6` its seen that conservative form results overall are pretty much closer to the exact analytical solution (`0.579`). Whereas in non-conservative form there are spurious oscillations observed at inflow and outflow boundaries. 

  5. The conservative form fetches a better mass flow distribution. It simply does a better of conserving mass which is obviously guessed from the name of the form itself.

  6. Although it takes more efforts in computing the conservative form, it eventually produces results which are faster than the non-conservative form.

Projects by Rohan Shah

a new title
Rohan Shah · 2018-11-25 19:35:42

Comments from the graderGreat job!But unfortunately, the values for power and specific fuel consumption is wrong Link to challenge:https://projects.skill-lync.com/projects/Data-Visualizer-Tool-69264 Read more

a new title
Rohan Shah · 2018-11-25 19:35:40

Comments from the graderGreat job!But unfortunately, the values for power and specific fuel consumption is wrong Link to challenge:https://projects.skill-lync.com/projects/Data-Visualizer-Tool-69264 Read more

a new title
Rohan Shah · 2018-11-25 19:35:40

Comments from the graderGreat job!But unfortunately, the values for power and specific fuel consumption is wrong Link to challenge:https://projects.skill-lync.com/projects/Data-Visualizer-Tool-69264 Read more

a new title
Rohan Shah · 2018-11-25 19:35:40

Comments from the graderGreat job!But unfortunately, the values for power and specific fuel consumption is wrong Link to challenge:https://projects.skill-lync.com/projects/Data-Visualizer-Tool-69264 Read more

a new title
Rohan Shah · 2018-11-25 19:35:39

Comments from the graderGreat job!But unfortunately, the values for power and specific fuel consumption is wrong Link to challenge:https://projects.skill-lync.com/projects/Data-Visualizer-Tool-69264 Read more

a new title
Rohan Shah · 2018-11-25 19:35:37

Comments from the graderGreat job!But unfortunately, the values for power and specific fuel consumption is wrong Link to challenge:https://projects.skill-lync.com/projects/Data-Visualizer-Tool-69264 Read more

a new title
Rohan Shah · 2018-11-25 19:35:37

Comments from the graderGreat job!But unfortunately, the values for power and specific fuel consumption is wrong Link to challenge:https://projects.skill-lync.com/projects/Data-Visualizer-Tool-69264 Read more

a new title
Rohan Shah · 2018-11-25 19:35:36

Comments from the graderGreat job!But unfortunately, the values for power and specific fuel consumption is wrong Link to challenge:https://projects.skill-lync.com/projects/Data-Visualizer-Tool-69264 Read more

a new title
Rohan Shah · 2018-11-25 19:35:36

Comments from the graderGreat job!But unfortunately, the values for power and specific fuel consumption is wrong Link to challenge:https://projects.skill-lync.com/projects/Data-Visualizer-Tool-69264 Read more

a new title
Rohan Shah · 2018-11-25 19:35:35

Comments from the graderGreat job!But unfortunately, the values for power and specific fuel consumption is wrong Link to challenge:https://projects.skill-lync.com/projects/Data-Visualizer-Tool-69264 Read more


Loading...

The End