Simulation of Super Sonic Nozzle Flow using Mac-Cormack Method

 Calculating Time Step 

   For the solution to be stable, the concept of CFL number is used. For linear Hyperbolic equation, the value of CFL < 1 will be stable solution. The time step for each iteration is calculated so that the result must be stable.

                      dt = c * dx / v+a

        for non dimensional form a = t^0.5 is substitued.

The minimum value of the time step is used to find the time derivative at each iteration, so the value of that variables are correct for all the iterations.

The value of dt must be less than or equal to the time it takes a sound wave to move from one grid to next the dt value must statisfy the condition for stability of equation.

Main Program

  1. The program grid point is set @ 31 for the both the methods for the Macormack method.

  2. The tic & toc command is used to find the simulation time of the program & bar chart to compare          the results.

  3. The comparision of the different profiles at the throat section for conservative & non conservative          form.

  4. The plot on the mass flow rate in the nozzle at the steady state with the exact solution.

% Program to comparision of both the methods

clear all
close all
clc

%inputs
l = 3; % length of nozzle
n = 31;% Grid space
x = linspace(0,l,n); % create x 
dx = x(2)-x(1);  % difference in x axis
gamma = 1.4;     
m_exact = 0.579;  

% calculate initial profiles 

a = 1+2.2*(x-1.5).^2;
nt = 1400;   % no of iterations
throat=find(a-1<0.001);

% courant number
c = 0.5;     


tic
[k,rho,t,v,M,m,p,rho_throat,p_throat,t_throat,M_throat] = non_conservative(x,n,dx,gamma,a,nt,throat,c)
solution_nc = toc;

tic
[k,rho_c,t_c,v_c,M_c,m_c,p_c,rho_throat_c,p_throat_c,t_throat_c,M_throat_c] = conservative(x,n,dx,gamma,a,nt,throat,c)
solution_c = toc;

g = ([solution_nc;solution_c])
figure(10)
bar (g)   % bar graph for the simulation time 
legend('Different forms')
xlabel('Methods')
ylabel('Simulation time')

% Comparing the Conservative & Non Conservative at the Throat section 
plot_throat(rho_throat,p_throat,t_throat,M_throat,rho_throat_c,p_throat_c,t_throat_c,M_throat_c)

% Steady state plot for mass flow rate 
figure(9)
plot(x,m,'b')
hold on
plot(x,m_c,'r')
hold on
plot(x,m_exact,'*')
xlabel('Nozzle length')
ylabel('Mass flow rate')
title('Distribution of Mass flow rate in the nozzle at steady state')
legend('Numerical solution Non conservative','Numerical solution conservative','Exact solution')
grid on

Plot comparision with the Conservative & Non Conservative with Exact solution

% Function for 1D supersonic simulation nozzle flow using Macormack method
% using non conservative form

function [k,rho,t,v,M,m,p,rho_throat,p_throat,t_throat,M_throat] = non_conservative(x,n,dx,gamma,a,nt,throat,c)


% Non conservative solution
% Calculate initial profiles 
rho = 1 - 0.3146*x;
t = 1-0.2314*x;
v = (0.1 + 1.09*x).*t.^0.5;

% Time-step control:
for i = 1:31
dt_value(i) = min(c.*dx./(sqrt(t(i)) + v(i)));
dt = min(dt_value);
end

% outer loop 

for k = 1:nt 
   % while (error>tol)
             
      % property copy
      rho_old =rho;
      t_old =t;
      v_old =v;
      
      % predictor method 

      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;
          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
            dvdx_c =( v(j)-v(j-1))/dx;
            drhodx_c =( rho(j)-rho(j-1))/dx;
            dtdx_c =( t(j)-t(j-1))/dx;
            dlogadx_c = (log(a(j))-log(a(j-1)))/dx;
           
            % continuity equation
            drhodt_c(j) = - rho(j)*dvdx_c - rho(j)*v(j)*dlogadx_c -v(j)*drhodx_c;
           
            % momentum equation
            dvdt_c(j) = - v(j)*dvdx_c - (1/gamma )*(dtdx_c +(t(j)/rho(j))*drhodx_c);
            
            % energy  equation
             dtdt_c(j) = -v(j)*dtdx_c - (gamma-1)*t(j) *(dvdx_c +v(j)*dlogadx_c);
        end
        % compute 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
        % Boudary conditions -Inlet
        v(1)=2*v(2) -v(3);
        rho(1) = 1;
        t(1) = 1;
        % Boudary conditions -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);
        
    % calculation of other variables of interest
    M=v./sqrt(t); %Mach Number
    m=rho.*v.*a; %Non-dimensional Mass Flow Rate
    p=rho.*t; %Non-dimensional pressure
    
    
   % Transient Throat Values
     rho_throat(k)=rho(throat);
     t_throat(k)=t(throat);
     p_throat(k)=p(throat);
     M_throat(k)=M(throat);
     figure (1)
     Mass_nc(x,k,m)
end 
   % Plot of different variables
     steady_state_nc (x,M,p,rho,t)

Conservative Form 

1. Function for 1D supersonic simulation nozzle flow using Macormack method for conservative form.

2. Formula for the conservative alone changes and all the procedures remains the same.

% Function for 1D supersonic simulation nozzle flow using Macormack method
% using conservative form

function [k,rho_c,t_c,v_c,M_c,m_c,p_c,rho_throat_c,p_throat_c,t_throat_c,M_throat_c] = conservative(x,n,dx,gamma,a,nt,throat,c)

% conservative solution
% Initial profiles for  density and temperature:
for i = 1:length(x)
    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 -0.366*(x(i) - 0.5);
            t_c(i) = 1.0 - 0.167*(x(i) - 0.5);
    elseif x(i) >= 1.5 && x(i) <= 3.0
                rho_c(i) = 0.634 - 0.3879*(x(i) - 1.5);
                t_c(i) = 0.833 - 0.3507*(x(i) - 1.5);
    end
end
    

% Initial profiels for velocity :
v_c = 0.59./(rho_c.*a);

% Time-step control:
for i = 1:31
dt_value(i) = min(c.*dx./(sqrt(t_c(i)) + v_c(i)));
dt = min(dt_value);
end

% Preallocation the different variable to iteration simple

dU1dt_p = zeros(1,30);
dU2dt_p = zeros(1,30);
dU3dt_p = zeros(1,30);
dU1dt_c = zeros(1,30);
dU2dt_c = zeros(1,30);
dU3dt_c = zeros(1,30);
J2 = zeros(1,30);
dt_value = zeros(1,31);
p_throat_c = zeros(1,1400);
rho_throat_c = zeros(1,1400);
t_throat_c = zeros(1,1400);
M_throat_c = zeros(1,1400);

% Initializing the solution vectors:
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
    
    % Preserving the old solution vectors:
    U1_old = U1;
    U2_old = U2;
    U3_old = U3;
    
    % Defining the flux vectors:
    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);
    
    % Solution by Mc-Cormack Method:1
    for j = 2:n-1
        
        % Defining the source term:
        J2(j) = (1/gamma).*rho_c(j).*t_c(j).*((a(j+1)-a(j))/dx);
        % Predictor step:
        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 the new 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;
    end
    
    % Updating the primitive variables and flux vectors:
    
    rho_c = U1./a;
    t_c = (gamma -1)*(U3./U1 - (gamma/2)*((U2./U1).^2));
        
    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
        % Updating the source term:
        J2(j) = (1/gamma)*rho_c(j).*t_c(j)*((a(j)-a(j-1))/dx);
        % Corrector step:
        dU1dt_c(j) = -(F1(j) - F1(j-1))/dx;
        dU2dt_c(j) = -((F2(j) - F2(j-1))/dx) + J2(j);
        dU3dt_c(j) = -(F3(j) - F3(j-1))/dx;
    end
    % Updating the solution vectors after corrector step:
    % Averaging the predictor and corrector steps:
    dU1dt = 0.5*(dU1dt_p + dU1dt_c);
    dU2dt = 0.5*(dU2dt_p + dU2dt_c);
    dU3dt = 0.5*(dU3dt_p + dU3dt_c);
    
    % Calculating the new solution vectors:
    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 the boundary conditions:
    % inlet conditions
    U1(1) = rho_c(1).*a(1);
    U2(1) = 2*U2(2) - U2(3);
    U3(1) = U1(1)*(t_c(1)/(gamma-1)+(gamma/2)*v_c(1).^2);
    
    % outlet conditions:
    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);
    
    % Final updation of the primitive variables:
    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;
    
    % Defining the Values of mach Number and mass flow rate:
    M_c = v_c./sqrt(t_c);
    m_c = rho_c.*a.*v_c;
    
    % Calculating values of specific quantites at the throat:
    M_throat_c(k) = M_c(throat);  % Mach Number
    p_throat_c(k) = p_c(throat);  % Pressure Ratio
    t_throat_c(k) = t_c(throat);  % Non-dimensional Temperature
    rho_throat_c(k) = rho_c(throat);  % Density Ratio
    
    % Comparison of non-dimensional mass-flow rates at different times:
      Mass_c(x,k,m_c)
end
     % Steady state plot
     steady_state_c (x,M_c,p_c,rho_c,t_c)

Throat 

function[]=plot_throat(rho_throat,p_throat,t_throat,M_throat,rho_throat_c,p_throat_c,t_throat_c,M_throat_c)

% Exact analytical solution at throat at steady state
rho_exact=0.634*ones(1,1400);
t_exact=0.833*ones(1,1400);
p_exact=0.528*ones(1,1400);
M_exact=1*ones(1,1400);

% Plots for convergence (at throat section)
figure(3)
plot(rho_throat)
hold on 
plot(rho_exact,'r')
hold on
plot(rho_throat_c,'g')
xlabel('Number of iterations')
ylabel('Density')
legend('Numerical solution Non conservative','Exact solution','Numerical solution conservative')
title('Density at Throat section')
    
figure(4)
plot(t_throat)
hold on 
plot(t_exact,'r')
hold on
plot(t_throat_c,'g')
xlabel('Number of iterations')
ylabel('Temperature')
legend('Numerical solution Non conservative','Exact solution','Numerical solution conservative')
title('Temperature at Throat section')
    
figure(5)
plot(p_throat)
hold on 
plot(p_exact,'r')
hold on
plot(p_throat_c,'g')
xlabel('Number of iterations')
ylabel('Pressure')
legend('Numerical solution Non conservative','Exact solution','Numerical solution conservative')
title('Pressure at Throat section')

figure(6)
plot(M_throat)
hold on 
plot(M_exact,'r')
hold on
plot(M_throat_c,'g')
xlabel('Number of iterations')
ylabel('Mach number')
legend('Numerical solution Non conservative','Exact solution','Numerical solution conservative')
title('Mach number at Throat section')

end

Steady State Sub-plot

Non conservative form

function[]=steady_state_nc(x,M,p,rho,T)

figure(7)
% Plots at Steady-state
subplot(4,1,1)
plot(x,M)
xlabel('Nozzle length')
ylabel('Mach number')
grid minor
axis ([0 3 0 4])
subplot(4,1,2)
plot(x,p)
xlabel('Nozzle length')
ylabel('P ratio')
grid minor
axis ([0 3 0 1])
subplot(4,1,3)
plot(x,rho)
xlabel('Nozzle length')
ylabel('rho ratio')
grid minor
axis ([0 3 0 1])
subplot(4,1,4)
plot(x,T)
xlabel('Nozzle length')
ylabel('Temp ratio')
grid minor
axis ([0 3 0 1])
suptitle('Qualitative aspects of Quasi 1-D Supersonic nozzle flow for NC')

end

Conservative Form