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

MacCormack method

Inputs :- n = 31
x = linspace(0,3,n)
dx = x(2) - x(1)
gamma = 1.4 C = 0.5
R = 8.314
nt = 1400

• Non Conservative form :

continuity equation (non dimensional) :

`(∂ρ)/(∂t) = - ρ(∂v)/(∂x) - ρv(∂(ln a)) / (∂x) - v(∂ρ)/(∂x)`

Momentum equation (non dimensional) :

`(∂v)/(∂t) = - v(∂v)/(∂x) - 1/γ ((∂t)/(∂x) - t/ρ (∂ρ)/(∂x))`

Energy equation (non dimensional):

`(∂t)/(∂t) = - v(∂t)/(∂x) - (γ-1)t ( (∂v)/(∂x) - v (∂(ln a)) / (∂x))`

calculation of time step(dt) :

dt =`C dx/(a + v)`

where, C = courant number;

a = speed of sound in reservoir =`√(γRt)`

v = Dimensionless velocity.

Initial Profiles :-

ρ = 1-0.3146x
t = 1-0.2314x
v = (0.1 + 1.09x)t^0.5

a = 1+2.2(x-1.5)^2

Area is same for both conservative and non conservative form.

Boundary Conditions :-

for inlet -
v(1) = 2*v(2) - v(3);

for outlet -
v(n) = 2v(n-1) - v(n-2);
ρ(n) = 2ρ(n-1) - ρ(n-2);
t(n) = 2t(n-1) - t(n-2);

To slove this we use predictor step and corrector step in which we use forward differencing for predictor step and backward differencing for correcter step. Below is the function code for the non conservative form.

Function code for Non conservative form :-

function [x,rho_non_conservative,v_non_conservative,t_non_conservative,p,m1,mf1,non_conservative_time] = non_conservative_macormack

clear all 
close all
clc

% Inputs
n = 31;
x = linspace(0,3,n);
dx = x(2) - x(1);
gamma = 1.4;

% calculating Initial Profiles
rho = 1-0.3146*x;
t = 1-0.2314*x; 
v = (0.1 + 1.09*x).*t.^0.5;

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

C = 0.5;
R = 8.314;
% time steps
nt = 1400;

tic
% outer time loop 
for k = 1:nt
    
    % time step
    dt = min((C*dx)./(((gamma*R*t).^.5)+v));
    
    rho_old = rho;
    v_old = v;
    t_old = t;
    
    % predictor method
    for j = 2:n-1
        
        % debugging the equations
        dvdx = (v(j+1)-v(j))/dx;
        drhodx = (rho(j+1)-rho(j))/dx;
        dtdx = (t(j+1)-t(j))/dx;
        dlogadx = (log(a(j+1))-log(a(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) = rho(j) + drhodt_p(j)*dt;
        v(j) = v(j) + dvdt_p(j)*dt;
        t(j) = t(j) + dtdt_p(j)*dt;
    end
    
    % Corrector method
    for j = 2:n-1
        
        % debugging the equations
        dvdx = (v(j)-v(j-1))/dx;
        drhodx = (rho(j)-rho(j-1))/dx;
        dtdx = (t(j)-t(j-1))/dx;
        dlogadx = (log(a(j))-log(a(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
    
     % 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;
         t(i) = t_old(i) + dtdt(i)*dt;
     end
     
     % apply boundary conditions
     
     % inlet
     v(1) = 2*v(2) - v(3);
     
     % 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);

     rho_non_conservative = rho;
     v_non_conservative = v;
     t_non_conservative = t;
     p = rho.*t;
     m1 = rho.*a.*v;
     inflow = rho(1)*v(1)*a(1);
     outflow = rho(n)*v(n)*a(n);
     mf1(k) = abs(outflow-inflow);
   
end
non_conservative_time = toc;
end

• Conservative form :

continuity equation (non dimensional):

`(∂(ρa)) / (∂t) + (∂(ρav)) / (∂x)`

momentum equation (non dimensional) :

`(∂(ρav)) / (∂t) + (∂[ρav^2 + (1/γ) ρa]) / (∂x) = (1/γ)p(∂a) / (∂x)`

Energy equation (non dimensional) :

`(∂ [ρ (e/(γ-1) + (γ/2)v^2) a]) / (∂t) + (∂ [ρ (e/(γ-1) + (γ/2)v^2) va + pav]) / (∂x) = 0`

elements definition :

u1 = ρa

u2 = ρav

u3 = ρ (e/(γ-1) + (γ/2)v^2) a

f1 = ρav

f2 = ρav^2 + (1/γ) ρa

f3 = ρ (e/(γ-1) + (γ/2)v^2) va + pav

j2 = (1/γ)p [(∂a) / (∂x)]

After substituting the elements in equations, we get,

`(∂u1) / (∂t) = -(∂f1) / (∂x)`

`(∂u2) / (∂t) = -(∂f2) / (∂x) + j2`

`(∂u3) / (∂t) = -(∂f3) / (∂x)`

now writing primitive variables in form of solution vector U

ρ = u1/a

v = u1/u2

t = e = (γ -1)[ (u3/u1) - (γ/2)v^2 ]

p = ρt

Initial profile :-

if, 0 <= x <= 0.5
then, ρ = 1
t = 1

if, 0.5 <= x <= 1.5
then, ρ = 1-0.366(x-0.5)
t = 1-0.167(x-0.5)

if, 1.5 <= x <= 3.5
then, ρ = 0.634-0.3879(x-1.5)
t = 0.833-0.3507*(x(s)-1.5)

v = 0.59/ρa

a = 1+2.2(x-1.5)^2

Boundary Conditions:-

Inlet
u1(1) = a(1)
u2(1) = 2u2(2) - u2(3)
u3(1) = u1(1)(1/(gamma-1) + (gamma/2)v(1)^2)

Outlet
u1(n) = 2u1(n-1) - u1(n-2)
u2(n) = 2u2(n-1) - u2(n-2)
u3(n) = 2u3(n-1) - u3(n-2)

To slove this we use predictor step and corrector step in which we use forward differencing for predictor step and backward differencing for correcter step. Below is the function code for the conservative form.

Function code for conservative form :-

function [x,rho_conservative,v_conservative,t_conservative,p,m2,mf2,conservative_time] = Conservative_macCormack

clear all 
close all
clc

% Inputs
n = 31;
x = linspace(0,3,n);
dx = x(2)-x(1);
gamma = 1.4;
C = 0.5;
nt = 1400;
R = 8.314;

tic
% Initial Profiles
for s = 1:n
    
    if x(s)>=0 && x(s)<=0.5
        rho(s) = 1;
        t(s) = 1;
        
    elseif x(s)>=0.5 && x(s)<=1.5
        rho(s) = 1-0.366*(x(s)-0.5);
        t(s) = 1-0.167*(x(s)-0.5);
        
    elseif x(s)>=1.5 && x(s)<=3.5
        rho(s) = 0.634-0.3879*(x(s)-1.5);
        t(s) = 0.833-0.3507*(x(s)-1.5);
    end
end

a = 1+(2.2*(x-1.5).^2);
v = 0.59./(rho.*a);
p = rho.*t;
u1 = rho.*a;
u2 = rho.*a.*v;
u3 = rho.*(t/(gamma-1)+(gamma/2)*v.^2).*a;

% time loop
for k = 1:nt
    
   % time step
   dt = min((C*dx)./((t.^.5)+v));
   
   u1_old = u1;
   u2_old = u2;
   u3_old = u3;
   rho_old = u1_old./a;
   v_old = u2_old./u1_old;
   t_old = (gamma-1)*((u3_old./u1_old)-(gamma/2)*v.^2);
   
   f1 = u2;
   f2 = (u2.^2)./u1+((gamma-1)/gamma)*(u3-(gamma/2)*(u2.^2)./u1);
   f3 = (gamma*u2.*u3)./u1 - ((gamma*(gamma-1)/2)*(u2.^3)./u1.^2);
   
   % predictor step
   for j = 2:n-1
       dadx = (a(j+1)-a(j))/dx;
       j2(j) = (1/gamma)*rho(j)*t(j)*dadx;
       
       % 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;
       
       % solution update
       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

   rho = u1./a;
   v = u2./u1;
   t = (gamma-1)*((u3./u1)-(gamma/2)*(u2./u1).^2);
   
   f1 = u2;
   f2 = (u2.^2)./u1+((gamma-1)/gamma)*(u3-(gamma/2)*(u2.^2)./u1);
   f3 = (gamma*u2.*u3)./u1 - ((gamma*(gamma-1)/2)*(u2.^3)./u1.^2);
   
   
   % corrector step
   for j = 2:n-1
       dadx = (a(j)-a(j-1))/dx;
       j2(j) = (1/gamma)*rho(j)*t(j)*dadx;
       
       % 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
   
    % computing the 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
  
   % applying Boundary Conditions
   % Inlet
   u1(1) = a(1);
   u2(1) = 2*u2(2) - u2(3);
   u3(1) = u1(1)*(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);
   
   rho = u1./a;
   v = u2./u1;
   t = (gamma-1)*((u3./u1)-(gamma/2)*v.^2);
   
   rho_conservative = rho;
   v_conservative = v;
   t_conservative = t;
   p = rho.*t;
   m2 = rho.*a.*v;
   inflow = rho(1)*v(1)*a(1);
   outflow = rho(n)*v(n)*a(n);
   mf2(k) = abs(outflow-inflow);

end
conservative_time = toc;
end

Main Program :-

clear all
close all
clc

[x,rho_non_conservative,v_non_conservative,t_non_conservative,p,m1,mf1,non_conservative_time] = non_conservative_macormack;
[x,rho_conservative,v_conservative,t_conservative,p,m2,mf2,conservative_time] = Conservative_macCormack;

figure(1)
plot(x,m1,'color','r',"LineWidth",1.5)
hold on
plot(x,m2,'color','b')
grid on
xlabel('length of nozzle')
legend('non conservative','conservative')
title('mass flow rate')

figure(2)
plot(x,rho_non_conservative,'color','y',"LineWidth",4)
hold on
plot(x,rho_conservative,'color','r',"LineWidth",0.5)
ylim([0 1.02])
grid on
xlabel('length of nozzle')
legend('non conservative','conservative')
title('Density')

figure(3)
plot(x,v_non_conservative,'color','g',"LineWidth",4)
hold on
plot(x,v_conservative,'--b',"LineWidth",0.5)
ylim([0 2.1])
grid on
xlabel('length of nozzle')
legend('non conservative','conservative')
title('Velocity')

figure(4)
plot(x,t_non_conservative,'color','b')
hold on
plot(x,t_conservative,'--r',"LineWidth",2)
ylim([0.3 1.02])
grid on
xlabel('length of nozzle')
legend('non conservative','conservative')
title('Temperature')

figure(5)
plot(x,p,'color','g')
hold on
plot(x,p,'*')
ylim([0 1.02])
grid on
xlabel('length of nozzle')
legend('non conservative','conservative')
title('Pressure')

figure(6)
plot(mf1,'color','b')
hold on
plot(mf2,'r')
grid on
xlabel('no. of Iterations')
legend('non conservative','conservative')
title('Net mass flow')

figure(7)
bar(1,non_conservative_time)
hold on
bar(2,conservative_time)
xticks(1:2)
ylabel('time')
legend('non conservative','conservative')
title('Simulation Time')

Following are the outputs of mass flow rate, density, velocity, temperature, pressure, net flow, simulaton time.

Outputs :-

• at n = 31

1) mass flow rate at n = 31

2) Density at n = 31

3) Velocity at n = 31

4) Temperature at n = 31

5) Pressure at n = 31

6) Net mass Flow at n = 31

7) Simulation Time at n = 31

reference to above figures we can see the following results:

• Mass flow rate of conservative form is more smooth and close to analytic solution compared to non conservative form.
• Conservative's primitive variables profiles are almost same to non conservative profiles, but there are small amount of fluctuation which is almost negligible.
• The small amount of fluctuation in conservative profiles are produced as because conservative form don't use primitive variables directly like non conservative form, as it uses generic form for those equations, in which it uses solution vector U flux vector F and source term J.
• We can see that conservative form takes less no of iterations around 500 iterations to converge and Non conservative form takes more compared to conservative form, which is 1000+ iterations to converge.

Below we have doubled the grid points(n= 61) to check whether these forms are grid dependent.

Outputs :-

• at n = 61

1) mass flow rate at n = 61

2) Density at n = 61

3) Velocity at n = 61

4) Temperature at n = 61

5) Pressure at n = 61

6) Net mass Flow at n = 61

7) Simulation Time at n = 61

From the above figures we can say that,

• mass flow rate and primitive variables profiles are almost same (i.e. they are grid independent)
• only no. of iterations to converge have increased for both method, as we also double the iterations (nt = 2800) to solve the solution.
• Here conservative form takes around 1000 iterations and vnon conservative form takes around 2200 iterations, actually doubles the iterations than in n= 31 and nt =1400.
• Also simulation time for both have also increased as the increase in no of iterations.

Conclusion :-

• Mass flow rate for Conservative is smoothest compared to Non conservative form.
• Looking to above results we can say that the simulation time for both methods is almost same when n= 31 and nt = 1400, but non conservative form is bit faster as it directly uses primitive varibles for calculations.
• And the conservative form is a bit slower as it uses solution and flux vectors to solve solutions and uses primitive variables indirectly.
• Also we can say that, both the conservative and non conservative forms are grid independent as there is no difference in result profiles as show in above outputs.