Simulation of a 1D Super-sonic nozzle flow using McCormack's Method

Objective

• Implement both the conservative and non-conservative forms of the governing equations
• Perform grid dependence test
• Write separate functions for conservative and non-conservative forms
• Figure out the minimum number of cycles for which the simulation should be run in order for convergence
• Compare the normalized mass flow rate between the conservative forms and non-conservative forms - make suitable conclusions
• Which method is faster?

Assumption: 

• Flow property varies along x axis only and do not vary along y axis, the flow is quasi 1D flow.

Solution

The main objective of this project is to solve 1D con-di supersonic nozzle by numerical methods. We are using McCormack’s method to solve the continuity, momentum and energy equation for 1D supersonic nozzle. McCormack’s method is an explicit finite difference technique which is second order accurate in both space and time. Here in this project we are calculating the important parameters of nozzle flow by solving the equations numerically. Although, the exact analytical solution of 1D supersonic nozzle is already present but we are using the numerical method for solving the equations. Just to appreciate the accuracy of numerical method.

Reference: All the equations used in the below script are from chapter 7 of Computational fluid dynamics by John D. Anderson Jr. The complete procedure for McCormack’s method is in chapter 6 of the book.

Main matlab code

clear all
close all
clc

%inputs
n = 31;
x = linspace(0,3,n); % mesh
dx = x(2) - x(1);

%number of time steps
nt = 1500;

% defining the node for throat
node_th = ((n-1)/2) + 1;

% calling the function for conservation form
tic
[rho_c,v_c,t_c,mfl_c,M_c,p_c,rho_th_c,v_th_c,t_th_c,mfl_th_c,M_th_c,p_th_c] = conservation(n,nt,node_th);
time_c = toc

% calling the function for non conservation form
tic
[rho_nc,v_nc,t_nc,mfl_nc,M_nc,p_nc,rho_th_nc,v_th_nc,t_th_nc,mfl_th_nc,M_th_nc,p_th_nc] = non_conservation(n,nt,node_th);
time_nc = toc

figure(3)

plot(x,mfl_nc,'r','linewidth',1.5)
hold on
plot(x,mfl_c,'bla','linewidth',1.5)
xlabel('Non-dimensional distance (x)')
ylabel('Mass flow rate')
legend('Non-conservation','Conservation','location','best')
title(sprintf('Comparison of mass flow for \nnon-conservation and conservation analysis \nNon-conservation analysis time = %f \nConservation analysis time =%f',time_nc,time_c))

figure(4)

subplot(4,1,1)
plot(x,M_nc,'r','linewidth',1.5)
hold on
plot(x,M_c,'.','color','bla','markersize',15)
ylabel('Mach number (M)')
legend('Non-conservation','Conservation','location','northwest')
title('Variation of Non-Dimensioanl flow parameters along x axis')
grid on
subplot(4,1,2)
plot(x,p_nc,'r','linewidth',1.5)
hold on
plot(x,p_c,'.','color','bla','markersize',15)
ylabel('Pressure (P)')
legend('Non-conservation','Conservation')
grid on
subplot(4,1,3)
plot(x,rho_nc,'r','linewidth',1.5)
hold on
plot(x,rho_c,'.','color','bla','markersize',15)
ylabel('Desity (rho)')
legend('Non-conservation','Conservation')
grid on
subplot(4,1,4)
plot(x,t_nc,'r','linewidth',1.5)
hold on
plot(x,t_c,'.','color','bla','markersize',15)
xlabel('Non-dimensional distance (x)')
ylabel('Temperature (T)')
legend('Non-conservation','Conservation')
grid on

figure(5)

subplot(4,1,1);
plot(linspace(1,nt,nt),M_th_nc,'r','linewidth',1.5)
hold on
plot(linspace(1,nt,nt),M_th_c,'bla','linewidth',1.5)
ylabel('Mach number (M)')
legend('Non-conservation','Conservation')
title(sprintf('Variation of Non-dimensional flow parameters at throat \naccording to time steps'))
grid on
subplot(4,1,2);
plot(linspace(1,nt,nt),p_th_nc,'r','linewidth',1.5)
hold on
plot(linspace(1,nt,nt),p_th_c,'bla','linewidth',1.5)
ylabel('Pressure (P)')
legend('Non-conservation','Conservation')
grid on
subplot(4,1,3);
plot(linspace(1,nt,nt),rho_th_nc,'r','linewidth',1.5)
hold on
plot(linspace(1,nt,nt),rho_th_c,'bla','linewidth',1.5)
ylabel('Density (rho)')
legend('Non-conservation','Conservation')
grid on
subplot(4,1,4);
plot(linspace(1,nt,nt),t_th_nc,'r','linewidth',1.5)
hold on
plot(linspace(1,nt,nt),t_th_c,'bla','linewidth',1.5)
xlabel('Time steps')
ylabel('Temperature (T)')
legend('Non-conservation','Conservation')
grid on

Function for Non-Conservation analysis

function [rho_nc,v_nc,t_nc,mfl_nc,M_nc,p_nc,rho_th_nc,v_th_nc,t_th_nc,mfl_th_nc,M_th_nc,p_th_nc] = non_conservation(n,nt,node_th)

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

% initial profiles in non dimentional form
rho_nc = 1-(0.3146*x); % density
t_nc = 1-(0.2314*x); % t = temperature
v_nc = (0.1 + (1.09*x)).*t_nc.^0.5; % velocity
a = 1 + (2.2*(x-1.5).^2); % area
c = 0.5; % courant number
dt = min(c.*(dx./((t_nc.^0.5) + v_nc))); % time step

    % outer time loop
    
    for k = 1:nt
    
        
        
        rho_old = rho_nc;
        v_old = v_nc;
        t_old = t_nc;
        
        % predictor method
        for j = 2:n-1
            dvdx = (v_nc(j+1) - v_nc(j))/dx;
            drhodx = (rho_nc(j+1) - rho_nc(j))/dx;
            dtdx = (t_nc(j+1) - t_nc(j))/dx;
            dlogadx = (log(a(j+1)) - log(a(j)))/dx;
            
            % continuity equation
            drhodt_p(j) = -rho_nc(j)*dvdx -rho_nc(j)*v_nc(j)*dlogadx -v_nc(j)*drhodx;
            
            % momentum equation
            dvdt_p(j) = -v_nc(j)*dvdx -(1/gamma)*(dtdx +(t_nc(j)/rho_nc(j))*drhodx);
            
            % energy equation
            dtdt_p(j) = -v_nc(j)*dtdx -(gamma-1)*t_nc(j)*(dvdx + v_nc(j)*dlogadx);
            
            % solution update
            rho_nc(j) = rho_nc(j) + drhodt_p(j)*dt;
            v_nc(j) = v_nc(j) + dvdt_p(j)*dt;
            t_nc(j) = t_nc(j) + dtdt_p(j)*dt;
            
        end
        
        % corrector method
        for j = 2:n-1
            dvdx = (v_nc(j) - v_nc(j-1))/dx;
            drhodx = (rho_nc(j) - rho_nc(j-1))/dx;
            dtdx = (t_nc(j) - t_nc(j-1))/dx;
            dlogadx = (log(a(j)) - log(a(j-1)))/dx;
            
            % continuity equation
            drhodt_c(j) = -rho_nc(j)*dvdx -rho_nc(j)*v_nc(j)*dlogadx -v_nc(j)*drhodx;
            
            % momentum equation
            dvdt_c(j) = -v_nc(j)*dvdx -(1/gamma)*(dtdx +(t_nc(j)/rho_nc(j))*drhodx);
            
            % energy equation
            dtdt_c(j) = -v_nc(j)*dtdx -(gamma-1)*t_nc(j)*(dvdx + v_nc(j)*dlogadx);
            
        end
        % computing 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 j = 2:n-1
            
            rho_nc(j) = rho_old(j) + drhodt(j)*dt;
            v_nc(j) = v_old(j) + dvdt(j)*dt;
            t_nc(j) = t_old(j) + dtdt(j)*dt;
        end
        % boundary condition
        
        %inlet
        v_nc(1) = 2*v_nc(2) - v_nc(3);
        
        %outlet
        rho_nc(n) = 2*rho_nc(n-1) - rho_nc(n-2);
        v_nc(n) = 2*v_nc(n-1) - v_nc(n-2);
        t_nc(n) = 2*t_nc(n-1) - t_nc(n-2);
        
        % calculating solution
        mfl_nc = rho_nc.*v_nc.*a; % mass flow rate
        M_nc = v_nc./((t_nc).^0.5); % Mach number
        p_nc = rho_nc.*t_nc; % pressure
        
        % calculating solution at throat
        rho_th_nc(k) = rho_nc(node_th);
        v_th_nc(k) = v_nc(node_th);
        t_th_nc(k) = t_nc(node_th);
        mfl_th_nc(k) = mfl_nc(node_th);
        M_th_nc(k) = M_nc(node_th);
        p_th_nc(k) = p_nc(node_th);
              
            % plotting mass flow rate at different time steps
            figure(1)
            if k == 1
                plot(x, mfl_nc,'--','color','bla','linewidth',1.5)
                hold on
            elseif k == 50
                plot(x, mfl_nc,'b','linewidth',1.5)
                hold on
            elseif k == 100
                plot(x, mfl_nc,'g','linewidth',1.5)
                hold on
            elseif k == 150
                plot(x, mfl_nc,'y','linewidth',1.5)
                hold on
            elseif k == 200
                plot(x, mfl_nc,'c','linewidth',1.5)
                hold on
            elseif k == 700
                plot(x, mfl_nc,'r','linewidth',1.5)
                hold on
            elseif k == 1000
                plot(x, mfl_nc,'.','color','bla','markersize',10)
   
               xlabel('Non-dimansional distance through nozzle (x)')
               ylabel('Non-dimensional mass flow')
               legend('1 time step','50 time steps','100 time steps','150 time steps','200 time steps','700 time steps','1000 time steps','location','north')
               title(sprintf('(Non-Conservation analysis) \nInstantaneous distribution of the non-dimensional mass flow \nas a function of distance through the nozzle \nat seven different times'))
               grid on
            end
    
        
    end

end

Function for Conservation analysis

function [rho_c,v_c,t_c,mfl_c,M_c,p_c,rho_th_c,v_th_c,t_th_c,mfl_th_c,M_th_c,p_th_c] = conservation(n,nt,node_th)

x = linspace(0,3,n); % mesh
dx = x(2)-x(1);
gamma = 1.4;
a = 1 + 2.2*(x-1.5).^2; % area
c = 0.5; % courant number

% initial profiles in non dimentional form
% defining density and temperature according to the condition given in Chapter 7 John
% D Anderson
    for i = 1:n
        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.366*(x(i) - 0.5);
            t_c(i) = 1 - 0.167*(x(i) - 0.5);
        elseif (x(i)>=1.5 && x(i)<=3.5)
            rho_c(i) = 0.634 - 0.3879*(x(i) - 1.5);
            t_c(i) = 0.833 - 0.3507*(x(i) - 1.5);
        end
        
    end

v_c = (0.59)./(rho_c.*a); 
dt = min(c.*(dx./((t_c.^0.5) + v_c))); % time step

% Defining the solution vector u, to represent quasi 1D flow in generic
% form

p_c = rho_c.*t_c;
u1 = rho_c.*a;
u2 = rho_c.*a.*v_c;
u3 = u1.*((t_c./(gamma - 1)) + ((gamma/2).*(v_c.^2)));

for k = 1:nt
   
    u1_old = u1;
    u2_old = u2;
    u3_old = u3;
   
    % flux terms defined in terms of u1,u2 and u3
    f1 = u2;
    f2 = ((u2.^2)./u1) + (((gamma - 1)/gamma)*(u3 - ((gamma/2)*((u2.^2)./u1))));
    f3 = ((gamma*u2.*u3)./u1) - (gamma*(gamma-1)*(u2.^3))./(2*u1.^2);
    
   % predictor method
   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;
       
       dadx(j) = (a(j+1) - a(j))/dx;
       j3(j) = (1/gamma)*(rho_c(j)*t_c(j))*dadx(j);
       
       % continuity equation
       du1dt_p(j) = -df1dx(j);
       
       % momentum equation
       du2dt_p(j) = -df2dx(j) + j3(j);
       
       % energy equation
       du3dt_p(j) = -df3dx(j);
       
       % 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
   
   % calculating the primitive terms from updated values of u1,u2,u3. Also
   % calculating f1,f2,f3.
   
   rho_c = u1./a;
   v_c = u2./u1;
   t_c = (gamma - 1)*((u3./u1) - ((gamma/2)*(v_c.^2)));
   f1 = u2;
   f2 = ((u2.^2)./u1) + (((gamma - 1)/gamma)*(u3 - (gamma*u2.^2)./(2*u1)));
   f3 = ((gamma*u2.*u3)./u1) - (gamma*(gamma-1)*(u2.^3))./(2*u1.^2);
   
   % corrector method
   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;
       
       dadx_1(j) = (a(j) - a(j-1))/dx;
       j3(j) = (1/gamma)*(rho_c(j)*t_c(j))*dadx_1(j);
       
       % continuity equation
       du1dt_c(j) = -df1dx(j);
       
       % momentum equation
       du2dt_c(j) = -df2dx(j) + j3(j);
       
       % energy equation
       du3dt_c(j) = -df3dx(j);
       
   end
   
   % coomputing 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 p_c = 2:n-1
       u1(p_c) = u1_old(p_c) + du1dt(p_c)*dt;
       u2(p_c) = u2_old(p_c) + du2dt(p_c)*dt;
       u3(p_c) = u3_old(p_c) + du3dt(p_c)*dt;
   end
      
   % boundary condition
   
   % inlet
   u1(1) = rho_c(1)*a(1);
   u2(1) = 2*u2(2) - u2(3);
   v_c(1) = u2(1)/u1(1);
   u3(1) = u1(1)*((t_c(1)/(gamma-1)) + ((gamma/2)*v_c(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);
   
   % calculating the primitive terms from the final updated solution
   rho_c = u1./a;
   v_c = u2./u1;
   t_c = (gamma - 1)*((u3./u1) - ((gamma*v_c.^2)/2));

   
   % calculating mass flow rate, mach number and pressure
   mfl_c = rho_c.*v_c.*a;
   M_c = v_c./((t_c).^0.5);
   p_c = rho_c.*t_c;
   
   % calculating solution at throat
   rho_th_c(k) = rho_c(node_th);
   v_th_c(k) = v_c(node_th);
   t_th_c(k) = t_c(node_th);
   mfl_th_c(k) = mfl_c(node_th);
   M_th_c(k) = M_c(node_th);
   p_th_c(k) = p_c(node_th);

   % plotting mass flow rate at different time steps
            figure(2)
            if k == 1
                plot(x, mfl_c,'--','color','bla','linewidth',1.5)
                hold on
            elseif k == 50
                plot(x, mfl_c,'b','linewidth',1.5)
                hold on
            elseif k == 100
                plot(x, mfl_c,'g','linewidth',1.5)
                hold on
            elseif k == 150
                plot(x, mfl_c,'y','linewidth',1.5)
                hold on
            elseif k == 200
                plot(x, mfl_c,'c','linewidth',1.5)
                hold on
            elseif k == 700
                plot(x, mfl_c,'r','linewidth',1.5)
                hold on
            elseif k == 1000
                plot(x, mfl_c,'.','color','bla','markersize',10)
   
               xlabel('Non-dimansional distance through nozzle (x)')
               ylabel('Non-dimensional mass flow')
               legend('1 time step','50 time steps','100 time steps','150 time steps','200 time steps','700 time steps','1000 time steps','location','north')
               title(sprintf('(Conservation analysis) \nInstantaneous distribution of the non-dimensional mass flow \nas a function of distance through the nozzle \nat seven different times'))
               grid on
            end
       
 end

end

Report

• Minimum number of cycles required for convergence is 700 cycles, this can be seen from the graph of mass flow at different time steps. The mass flow at 700 time steps and 1000 time steps coincides each other, hence it can be said that 700 cycles are required for convergence.
• below is the graph for comparison of mass flow rate between the conservative form and non-conservative form.
• Analysis time for Non-Conservation form is 0.520425 seconds and analysis time for Conservation form is 1.247168 seconds. Non-Conservation analysis is faster than conservation analysis.
• This is because non-conservation form deals with the primitive terms directly and solve them directly. While, the conservation form deals with some extra variables which are defined in terms of primitive terms. These extra variables are solution terms and flux terms. As a result of that more work and time is needed to decode those extra variables to calculate the primitive terms (primitive terms are: density, velocity, temperature, mach no. pressure)

First, define the the input data. In the main code we define number of grid points (n), mesh (x) and grit size (dx). Also define the node for throat area in the nozzle, to calculate the solution at throat. After defining the input data we call the functions for conservative and non-conservative analysis. These two functions calculate all the solutions.

Inside the function for non-conservation analysis, we define the initial conditions. first we define density, velocity, temperature of the fluid flow and area. All these parameters are in non dimensional form. Then calculate the appropriate value of "dt". We can calculate this value with the help of courant number "c". The formula for calculating "dt" is given in the above code. Now we have to start the time loop which runs from 1 to nt (total time steps). inside this time loop we have predictor step defined in a loop which runs from 2 to n-1 node in the mesh. After predictor step we update the values of density, velocity and temperature. these updated values are used in the corrector step. After corrector step we finally have the final update values of density, velocity and temperature. After corrector step we define boundary condition at inlet and outlet by extrapolation from the neighbouring points. After defining the boundary condition we calculate the other variables of interest like mach number, mass flow and pressure. We also calculate these value at the throat area and plot there timewise variation on the graph.

The script for conservation analysis is similar to the non-conservation analysis, except here, we define the solution terms u1, u2, u3 in terms of the primitive variables like density, velocity and temperature. And we also define the pure flux terms in terms of u1, u2 and u3. these solution terms and flux terms are use in the predictor step and corrector step to find the solution. After getting the solution we decode density, velocity and temperature from these solution terms u1,u2,u3.

Graphs

In the above graph we can see that the first line is way off from convergence, this is because it is at first time step and the solution has not yet converged. the cyan colored line at 200 time steps is close to convergence. The red line at 700 time steps and the dotted line at 1000 time steps coincides each other. these two line indicates that the solution has converged and minimum 700 cycles are required to achieve convergence.

Here, we can see that the lines are pretty close to convergence in the starting of the time loop. The lines are not scattered as we saw in non-conservation analysis of mass flow through the nozzle. here, we can also see that the solution converges at 700 time steps. And the line at 700 time steps and 1000 time steps coincides each other.

In the graph above we can see that plot for non-conservation is varying a lot it decreases at throat area and then increases after that. But, for conservation analysis the mass flow in almost constant, the plot s almost straight. This is because the the control volume is moving with the flow and hence there are a lot of variation in the parameters as the flow advances. But in conservation analysis the controlvolume is fixed, that is why it is easy to maintain the flow variable at a fix point in the fluid flow. And hence we see an almost staright line for mass flow in conservation analysis.

In the title of the graph we can see that analysis time is also shown for both non-conservation and conservation form of analysis. here, we can see that non-censervation analysis is faster than the conservation analysis. This is because in conservation analysis there are some extra solution terrms which are needed to be decoded to get the primitive terms and hence require some extra work and time.

Above is the garph for variation of mach number. pressure, density and temperature along the non dimenstional distance "x" through the nozzle. Here, we can see that the plot for non-conservation and conservation analysis is nearly identical. In the graph we can see that all the parameters follows the conditions for con-di supersonic nozzle. The mach number is increasing and becomes greater than two at the divergent side of the nozzle as expected for a supersonic nozzle. Pressure, temperature and density is very high at first but keeps on decreasing at the end of nozzle.

Above is the graph for timewise variation of flow parameters at the throat. The plot for the non-conservation is varying at first but then soon it becomes stable, but in case of conservation form the plot is fairly stable right from the starting of the cycle. Both the plots almost coincides each other. The mach number is one at the throat as it was supposed to be. The rest of the parameters also achived stable condition at around 200 cycles and keeps stable for rest of the time cycles.

GRID DEPENDENCE TEST

The grid dependence test is necessary for chacking the solution is dependent on number of grid points or it is independent. In CFD it is always desireble to get a grid independent result. For this purpose we ran our code at two different number of grid points, n=31 and n=61. We just double the number of grid points to check the results. In our analysis the solution is independent of number of grid points, the results are fairly similar according to our level of accuracy. Here in our analysis we can say that non-conservational analysis is more accurate than conservation analysis. As there is very less difference in the solution for 31 grid points and 61 grid points. 

Conditions at nozzle throat:-