Study of quasi 1-D supersonic flow through a nozzle using Macormack technique

• The objective of this challenge is to study the Macormack technique and set up a simulation case for the quasi 1-D supersonic flow through a nozzle. 
• Both, the conservative and non-conservative forms of the governing equations were set up. The characteristics of both these forms were studied.
• The case set-up and initial conditions were taken from the book \'Computational Fluid Dynamics\' by John. D. Anderson, Jr.  [Chapter 7].

A. Non-Conservative form:

Code Snippet:

Note: To avoid overlapping of results, the code section of obtaining variation of flow-field variables w.r.t time was run separately. (The section is commented) 

clc
clear all
close all

%Initializing the grid
n=31;                       %Number of grid points
x=linspace(0,3,n);
dx=x(2)-x(1);               %Cell size

%Calculating the initial profiles
rho=1-0.3146*x;             %Density value
T=1-0.2314*x;               %Temperature value
v=(0.1+0.9*x).*T.^0.5;      %Velocity value
p= rho.*T;                  %Pressure value

a=1+2.2*(x-1.5).^2;         %Area profile of the nozzle

%Calculating the time step
C=0.5;                      %Courant Number = 0.5
dt=min((C*dx)./(a+v));      %Time step as a function of CFL Number, area and velocity

gamma=1.4;                  %Gamma value

%Outer time loop
for i=1:900
    
    rho_old = rho;
    T_old = T;
    v_old = v;
    
    %Predictor Step
    for j=2:n-1
        
        %Standardizing the values
        dvdx = (v(j+1)-v(j))/dx;
        dlogadx = (log(a(j+1))-log(a(j)))/dx;
        dTdx = (T(j+1)-T(j))/dx;
        drhodx = (rho(j+1)-rho(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));
        
        %Updating the values
        rho(j) = rho(j)+drhodt_p(j)*dt;
        T(j) = T(j)+dTdt_p(j)*dt;
        v(j) = v(j)+dvdt_p(j)*dt;
        p(j) = rho(j).*T(j);
    end
    
    %Corrector Step
    for j=2:n-1
        
        %Standardizing the values
        dvdx = (v(j)-v(j-1))/dx;
        dlogadx = (log(a(j))-log(a(j-1)))/dx;
        dTdx = (T(j)-T(j-1))/dx;
        drhodx = (rho(j)-rho(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
    
    %Average values
    drhodt = 0.5*(drhodt_p + drhodt_c);
    dvdt = 0.5*(dvdt_p + dvdt_c);
    dTdt = 0.5*(dTdt_p + dTdt_c);
    
    %Final value update
    for k=2:n-1
    rho(k) = rho_old(k) + drhodt(k)*dt;
    v(k) = v_old(k) + dvdt(k)*dt;
    T(k) = T_old(k) + dTdt(k)*dt;
    end
    
    %Inlet BC
    v(1) = 2*v(2) - v(3);
    
    %Outlet BCs
    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);
    
    m = rho.*a.*v;
    M = v./(sqrt(T));
    p = rho.*T;
%%%%%Executed separately to prevent overlapping of results%%%%%     
%     v_throat(i) = v(1,16);
%     rho_throat(i) = rho(1,16);
%     T_throat(i) = T(1,16);
%     p_throat(i) = p(1,16);
%     
%     figure(1);
%     subplot(4,1,1);
%     plot(v_throat,\'b\',\'linewidth\',1.5);
%     axis([0 900 0 2]);
%     hold on
%     xlabel(\'Iteration Number\');
%     ylabel(\'Velocity\');
%     title(sprintf(\'Iteration Number = %d\',i));
%     pause(0.03);
% 
%     subplot(4,1,2);
%     plot(rho_throat,\'g\',\'linewidth\',1.5);
%     axis([0 900 0 2]);
%     hold on
%     xlabel(\'Iteration Number\');
%     ylabel(\'Density\');
%     title(sprintf(\'Iteration Number = %d\',i));
%     pause(0.03);
% 
%     subplot(4,1,3);
%     plot(T_throat,\'r\',\'linewidth\',1.5);
%     axis([0 900 0 2]);
%     hold on
%     xlabel(\'Iteration Number\');
%     ylabel(\'Temperature\');
%     title(sprintf(\'Iteration Number = %d\',i));
%     pause(0.03);
%     
%     subplot(4,1,4);
%     plot(p_throat,\'y\',\'linewidth\',1.5);
%     axis([0 900 0 2]);
%     hold on
%     xlabel(\'Iteration Number\');
%     ylabel(\'Pressure\');
%     title(sprintf(\'Iteration Number = %d\',i));
%     pause(0.03);
    
    figure(2);
    subplot(4,1,1);
    plot(x,v);
    xlabel(\'Displacement\');
    ylabel(\'Velocity\');
    title(sprintf(\'Iteration Number = %d\',i));
    pause(0.03);
    grid minor;
    
    subplot(4,1,2);
    plot(x,rho);
    xlabel(\'Displacement\');
    ylabel(\'Density\');
    title(sprintf(\'Iteration Number = %d\',i));
    pause(0.03);
    grid minor;
    
    subplot(4,1,3);
    plot(x,T);
    xlabel(\'Displacement\');
    ylabel(\'Temperature\');
    title(sprintf(\'Iteration Number = %d\',i));
    pause(0.03);
    grid minor;
    
    subplot(4,1,4);
    plot(x,p);
    xlabel(\'Displacement\');
    ylabel(\'Pressure\');
    title(sprintf(\'Iteration Number = %d\',i));
    pause(0.03);
    grid minor;
    
        
figure(3)
if i==100
plot(x,m,\'g\',\'linewidth\',1);
hold on
else if i==250
plot(x,m,\'y\',\'linewidth\',1);
hold on
else if i==400
plot(x,m,\'b\',\'linewidth\',1);
hold on
else if i==500
plot(x,m,\'r\',\'linewidth\',1);
hold on
else if i==600
plot(x,m,\'magenta\',\'linewidth\',1);
hold on
else if i==700
plot(x,m,\'black\',\'linewidth\',1);
hold on
else if i==850
plot(x,m,\'cyan\',\'linewidth\',1);
legend(\'100dt\',\'200dt\',\'700dt\');
title(\'Mass flow rate in a nozzle wrt time step for non-conservative form\');
xlabel(\'Non Dimensional length of the nozzle\')
ylabel(\'mass flow rate\')
hold off
    end
    end
    end
    end
    end
    end
end


end

Results:

1. Profile of various parameters along the nozzle length:

                                                            (a)

2. Mass flow rate w.r.t time:

• As can be seen from the plot, the mass flow rate variation stabilizes after 700 time steps and the steady state has been reached at 850 time steps.
• Theoretically, the solution will take infinite time to reach the steady state but for practical purposes, it is achieved much earlier, albeit with minor deviations from the exact solution.

3. Residuals:

• As can be seen from the above plot, the values of all the parameters show a wavy nature in the beginning till 700 time steps. This means that with the non-conservative approach, the solution is initially unstable but attains steady state after a certain time (here, 700 time steps).

B. Conservative Form:

Code Snippet:

Note: To avoid overlapping of results, the code section of obtaining variation of flow-field variables w.r.t time was run separately. (The section is commented) 

clc
close all
clear all

%Initializing the grid
n=31;                       %Number of grid points
x=linspace(0,3,n);      
dx=x(2)-x(1);               %Cell size

%Initializing the parameter profiles
for i=1:length(x)               %Initializing values of density and temperature based on nozzle profile
if x(i)>=0 && x(i)<=0.5
rho(i) = 1;
T(i) = 1;
end
if 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);
end
if x(i)> 1.5 && x(i)<=3
rho(i) = 0.634-0.3879*(x(i)-1.5);
T(i) = 0.833-0.3507*(x(i)-1.5);
end
end

A = 1+2.2*((x-1.5).^2);                 %Area profile of the nozzle
v = 0.59./(rho.*A);                     %velocity value
C = 0.5;                                %Courant Number
y = 1.4;                                %Gamma value

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

%Outer Time Loop
for i=1:1400

%Time Step Calculation
dt = min(C.*(dx./(sqrt(T)+v)));

U1_old = U1;
U2_old = U2;
U3_old = U3;

F1 = U2;
F2 = ((U2.^2)./U1) + (((y-1)/y)*(U3 - ((y*0.5*(U2.^2))./U1)));
F3 = ((y*U2.*U3)./U1) - ((y*0.5*(y-1)*(U2.^3))./(U1.^2));

%Space marching
for j=2:n-1

J2(j) = (1/y).*rho(j).*T(j).*((A(j+1)-A(j))/dx);

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


% Upadated value of density, velocity and temperature after predicted values.

rho(j) = U1(j)/A(j);
v(j) = U2(j)/U1(j);
T(j) = (y-1)*((U3(j)/U1(j)) - ((y/2)*(v(j)^2)));
end
% Update the values of F1,F2,F3 after predictor step
F1 = U2;
F2 = ((U2.^2)./U1) + (((y-1)/y)*(U3 - ((y*0.5*(U2.^2))./U1)));
F3 = ((y*U2.*U3)./U1) - ((y*0.5*(y-1)*(U2.^3))./(U1.^2));

% Corrector step
for k=2:n-1
J2(k) = (1/y)*rho(k)*T(k)*((A(k)-A(k-1))/dx); 
dU1dt_c(k) = -((F1(k) - F1(k-1))/dx);
dU2dt_c(k) = -((F2(k) - F2(k-1))/dx) + J2(k);
dU3dt_c(k) = -((F3(k) - F3(k-1))/dx);
end

%Average Time Derivatives
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
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

%Updating primitive variables
rho = U1./A;
T = (y-1)*((U3./U1) - ((y/2)*(v.^2)));
v = U2./U1;

%Inlet BC
U1(1) = rho(1)*A(1); 
U2(1) = 2*U2(2) - U2(3);
U3(1) = U1(1)*((T(1)/(y-1))+((0.5*y)*(v(1)^2))); 

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


m = rho.*A.*v;
M = v./(sqrt(T));
p = rho.*T;

%%%%%Executed separately to prevent overlapping of results%%%%%   
%     v_throat(i) = v(1,16);
%     rho_throat(i) = rho(1,16);
%     T_throat(i) = T(1,16);
%     p_throat(i) = p(1,16);
%     
%     figure(1);
%     subplot(4,1,1);
%     plot(v_throat,\'b\',\'linewidth\',1.5);
%     axis([0 1400 0 2]);
%     hold on
%     xlabel(\'Iteration Number\');
%     ylabel(\'Velocity\');
%     title(sprintf(\'Iteration Number = %d\',i));
%     pause(0.03);
% 
%     subplot(4,1,2);
%     plot(rho_throat,\'g\',\'linewidth\',1.5);
%     axis([0 1400 0 2]);
%     hold on
%     xlabel(\'Iteration Number\');
%     ylabel(\'Density\');
%     title(sprintf(\'Iteration Number = %d\',i));
%     pause(0.03);
% 
%     subplot(4,1,3);
%     plot(T_throat,\'r\',\'linewidth\',1.5);
%     axis([0 1400 0 2]);
%     hold on
%     xlabel(\'Iteration Number\');
%     ylabel(\'Temperature\');
%     title(sprintf(\'Iteration Number = %d\',i));
%     pause(0.03);
%     
%     subplot(4,1,4);
%     plot(p_throat,\'y\',\'linewidth\',1.5);
%     axis([0 1400 0 2]);
%     hold on
%     xlabel(\'Iteration Number\');
%     ylabel(\'Pressure\');
%     title(sprintf(\'Iteration Number = %d\',i));
%     pause(0.03);


figure(2)
subplot(4,1,1);
plot(x,v);
xlabel(\'Displacement\');
ylabel(\'Velocity\');
title(sprintf(\'Iteration Number = %d\',i));
pause(0.03);
grid minor;
    
subplot(4,1,2);
plot(x,rho);
xlabel(\'Displacement\');
ylabel(\'Density\');
title(sprintf(\'Iteration Number = %d\',i));
pause(0.03);
grid minor;
  
subplot(4,1,3);
plot(x,T);
xlabel(\'Displacement\');
ylabel(\'Temperature\');
title(sprintf(\'Iteration Number = %d\',i));
pause(0.03);
grid minor;
   
subplot(4,1,4);
plot(x,p);
xlabel(\'Displacement\');
ylabel(\'Pressure\');
title(sprintf(\'Iteration Number = %d\',i));
pause(0.03);
grid minor;
    
figure(3)
if i==100
plot(x,m,\'g\',\'linewidth\',1);
hold on
else if i==250
plot(x,m,\'y\',\'linewidth\',1);
else if i==400
plot(x,m,\'b\',\'linewidth\',1);
else if i==500
plot(x,m,\'r\',\'linewidth\',1);
else if i==600
plot(x,m,\'magenta\',\'linewidth\',1);
else if i==200
plot(x,m,\'cyan\',\'linewidth\',1);
else if i==700
plot(x,m,\'black\',\'linewidth\',1);
else if i==900
plot(x,m,\'-rs\',\'linewidth\',1);
else if i==1000
plot(x,m,\'-bo\',\'linewidth\',1);
else if i==1200
plot(x,m,\'-k^\',\'linewidth\',1);
plot(x,m,\'r\',\'linewidth\',1)
legend(\'100dt\',\'200dt\',\'700dt\');
title(\'Mass flow rate in a nozzle wrt time step for conservative form\');
xlabel(\'Non Dimensional length of the nozzle\')
ylabel(\'mass flow rate\')
hold off
end
end
end
end
end
end
end
end
end
end

end

1. Profile of various parameters along the nozzle length:

                                                               (b)

• As can be seen from figures (a) and (b), alll the parameter profiles are identical in both, the non-conservative as well as the conservative approach.
• Hence, if the objective is to get a profile of any parameter, both the approaches will yield similar results. However, as can be observed, the non-conservative approach is quicker than the conservative approach.

2. Residuals:

• From the above plot, it is observed that the parameters attain a stable value at roughly the 250th time step, much lower as compared to the 700th+ time step in the non-conservation form.
• Hence, the conservative approach is able to reach the steady state value in a much lower time.

3. Mass flow rate w.r.t time:

                                                  (c) 31 grid points

                                                    (d) 61 grid points

• Here, there are less fluctuations when compared with the non-conservation plot.
• This is partly because of the initial consitions, which were selected with constant mass flow due to the higher sensitivity of the conservative approach.
• The solution starts attaining stability at around 650th time step and then sticks with the constant value which is very close to the exact value.

• The figures c,d and the above table can be used to study the concept of grid independence. A case set-up is said to have grid independence if the solution changes only marginally by increasing the number of cells in the grid(thus reducing the individual cell size).
• In this case, the solution values for all the parameters are very close to those of the exact analytical values. The deviations are very small and can be neglected since the accuracy required is not high.
• Hence, this set-up can be said to be grid independent at n=31 and having more number of grid points will not make any considerable difference in the solutions.