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

AIM: Simulation of Quasi 1 D Supersonic nozzle flow simulation using MacCormack Method

OBJECTIVE:

1. Steady-state distribution of primitive variable inside the nozzle.
2. Time-Wise variation of a primitive variable.
3. Variation of mass flow rate distribution inside the nozzle at different time steps during the time marching process.
4. Comparison of the normalized mass flow rate distribution of both forms.

INTRODUCTION OF PROBLEM/PHYSICAL ASPECT:

1D, steady, isentropic flow considered through the Convergent – Divergent nozzle. It is assumed that the nozzle is connected to the very large reservoir whose area is infinite. Hence, the velocity at the inlet is zero or negligible. Pressure and temperature at the inlet of the nozzle are P₀ and T₀.

It is assumed that the flow properties are varying along the X – Direction only and at section flow properties are constant along the Y – Direction.

ANALYTICAL SOLUTION OF PROBEM:

 Governing equations is as follows:

Solving the above equation analytically, we get the following results:

Solution by using CFD technique:

The solution is obtained using both conservative and non-conservative forms of governing equation.

Following are the governing equations in non-conservative form (Dimensional):

Following are the governing equations in non-conservative form (Non – Dimensional):

Where:

Please note that all prime quantities are non-Dimensional in nature.

Now, to solve the above equations using MacCormack Method (Predictor and Corrector method)

1. Predictor step:                                
2. Corrector Step:

Calculation of Time Step:

Time step is calculated for all time steps using the following expression:

Where C = CFL number

Time step is calculated for all nodes and minimum of it is taken as time step for further calculation.

Nozzle Shape:

Nozzle shape is varying parabolically w.r.t.x as A = 1 + 2.2(x – 1.5)² f OR 0 ≤ x ≤ 3

Initial Conditions:

The initial condition for density, velocity, and temperature is calculated from the following expression:

Boundary Conditions:

Solution using a non-conservative form of governing equation:

Now, let us define elements if solution vector U, the solution vector F, and the source term J as follows:

Substituting the above equations in governing differential equation, we get

The above equations are governing differential equations for quasi 1D flow in conservative form.

Now from definitions U1, U2, and U3

Initial Conditions:

Nozzle Shape:

Boundary Conditions:

1. Subsonic Inflow boundary condition (First Node):

     2. Supersonic outflow boundary condition (Last node):

Main Code:

Non - Conservative Part:

clear all
close all
clc

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

%% Initial Values

rho = 1-0.3146.*x;              % Density
T = 1-0.2314.*x;                % Temperature
V = (0.1+1.09.*x).*T.^0.5;      % Velocity

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

%% Time Step calculation by Courant no.
nt = 1400; % Total no. of time steps
%C = linspace(0.1,1,n); % Courant No. defined (0<=C<=1)
C = 0.5;
a = sqrt(T);
for i = 1:n
    dt_1 = C.*(dx./(a+V));
end

dt = min(abs(dt_1));

%% Therotical mass flow rate
for i = 1:n
    m_initial(i) = rho(i)*A(i)*V(i);
end

%% Outer time loop (Solving using McCormack Method)

for j = 1:nt
    
    % Saving the old values
        rho_old = rho;
        V_old = V;
        T_old = T;
    
    % Predictor method    
    for k = 2:(n-1)
        
        % Splitting the values of easyfication
        
        dvdx = (V(k+1)-V(k))/dx;
        drhodx = (rho(k+1)-rho(k))/dx;
        dlogadx = (log(A(k+1))-log(A(k)))/dx;
        dTdx = (T(k+1)-T(k))/dx;
        
        % Continuity Equation       
        drho_dt_p(k) = -rho(k)*(dvdx)-rho(k)*V(k)*(dlogadx)-V(k)*(drhodx);
        
        % Momentum Equation
        dv_dt_p(k) = -V(k)*(dvdx)-(1/gamma)*(dTdx + (T(k)/rho(k))*drhodx);
        
        % Energy Equation
        dT_dt_p(k) = -V(k)*(dTdx)-(gamma-1)*T(k)*(dvdx + V(k)*dlogadx);
        
        
        % Solution Update
        rho(k) = rho(k) + drho_dt_p(k)*dt;
        V(k) = V(k) + dv_dt_p(k)*dt;
        T(k) = T(k) + dT_dt_p(k)*dt;
        
    end
    
    
    % Corrector method    
    for k = 2:(n-1)
        
        % Splitting the values of easyfication
        
        dvdx = (V(k)-V(k-1))/dx;
        drhodx = (rho(k)-rho(k-1))/dx;
        dlogadx = (log(A(k))-log(A(k-1)))/dx;
        dTdx = (T(k)-T(k-1))/dx;
        
        % Continuity Equation       
        drho_dt_c(k) = -rho(k)*(dvdx)-rho(k)*V(k)*(dlogadx)-V(k)*(drhodx);
        
        % Momentum Equation
        dv_dt_c(k) = -V(k)*(dvdx)-(1/gamma)*(dTdx + (T(k)/rho(k))*drhodx);
        
        % Energy Equation
        dT_dt_c(k) = -V(k)*(dTdx)-((gamma-1)*T(k))*(dvdx + V(k)*dlogadx);             
        
    end
    
    % Compute the average time derivative
    
    dvdt = 0.5*(dv_dt_p + dv_dt_c);
    drhodt = 0.5*(drho_dt_p + drho_dt_c);
    dTdt = 0.5*(dT_dt_p + dT_dt_c);
    
    % Final Solution 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
    
    % Defining the Boundary Conditions
    % Inlet
    V(1) = 2*V(2)-V(3);
    rho(1) = 1;
    T(1) = 1;
    
    %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);
    
    % Calculating the Mach No. and pressure
    a = sqrt(T);
    M = V./a;
    p = rho.*T;
    
    % Calculating the mass flow rate along the length of the nozzle
    m = rho.*A.*V;
    
    % Calculation of primitive variables at throat
    rho_t(j) = rho(16);
    V_t(j) = V(16);
    T_t(j) = T(16);
    p_t(j) = p(16);
    M_t(j) = M(16);
    m_t(j) = m(16);
    
    % Plotting
    
    figure(6)
    if j==50    
        plot(x,m,'r','linewidth',1.5)
        hold on
    else if j==100
        plot(x,m,'g','linewidth',1.5)
        hold on
        else if j==200
            plot(x,m,'b','linewidth',1.5)
            hold on
            else if j==300
                plot(x,m,'black','linewidth',1.5)
                hold on
                else if j==400
                    plot(x,m,'y','linewidth',1.5)
                    hold on
                    else if j==500
                        plot(x,m,'linewidth',1.5)
                        hold on
                        else if j==600
                            plot(x,m,'linewidth',1.5)
                            hold on
                            else if j==700
                                plot(x,m,'linewidth',1.5)
                                hold on
                                
                                plot(x,m_initial,'--','linewidth',1.5)
        
                                xlabel('Length of nozzle');
                                ylabel('Mass Flow Rate');
                                legend('50dt','100dt','200dt','300dt','400dt','500dt','600dt','700dt','Initial Mass flow rate')
                                title('Mass Flow rate from the Nozzle for Non-Conservative  part')
                                end
                            end
                        end
                    end
                end
            end
        end
    end
end

figure(1)
plot(x,V,'r')
xlabel('Length of Nozzle');
ylabel('Fluid Velocity');
legend('VELOCITY');
title(sprintf('Time Step: %d',j));
   
figure(2)
plot(x,rho,'g')
xlabel('Length of Nozzle');
ylabel('Fluid Density');
legend('DENSITY');
title(sprintf('Time Step: %d',j));
    
figure(3)
plot(x,T,'b')
xlabel('Length of Nozzle');
ylabel('Fluid Temperature');
legend('TEMPERATURE');
title(sprintf('Time Step: %d',j));
    
figure(4)
plot(x,M,'black')
xlabel('Length of Nozzle');
ylabel('Mach Number');
legend('MACH NUMBER');
title(sprintf('Time Step: %d',j));
    
figure(5)
plot(x,m_initial)
hold on
plot(x,m,'black')
xlabel('Length of Nozzle');
ylabel('Fluid Mass Flow Rate');
legend('Initial Mass Flow Rate','Final Mass flow rate');
title(sprintf('Time Step: %d',j));
    
pause(0.01)
hold off

figure(7)
subplot(511)
plot(1:j,V_t,'linewidth',1.5)
xlabel('No. of Time Step')
ylabel('Velocity at Throat')
txt1 = (sprintf('Time Step: %d',j));
txt2 = ('Time-wise variation of the primitive variables');
title({txt1,txt2})
    
subplot(512)
plot(1:j,rho_t,'linewidth',1.5)
xlabel('No. of Time Step')
ylabel('Density at Throat')
    
subplot(513)
plot(1:j,T_t,'linewidth',1.5)
xlabel('No. of Time Step')
ylabel('Temperature at Throat')

subplot(514)
plot(1:j,p_t,'linewidth',1.5)
xlabel('No. of Time Step')
ylabel('Pressure at Throat')

subplot(515)
plot(1:j,M_t,'linewidth',1.5)
xlabel('No. of Time Step')
ylabel('Mach Number at Throat')

Conservative Part:

clear all
close all
clc

%% Inputs (Defining the mesh size)
n = 31;
l = 3; % Length of domain
x = linspace(0,l,n);
dx = x(2)-x(1);
gamma = 1.4;

%% Initial Values

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

% Primitive variable of the nozzle

% The primitive variable Temperature and Density were kept constant and the
% third primitive variable Velocity was kept floating with the grid point

for i = 1:n
    
    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)
        rho(i) = 0.634-0.3879*(x(i)-1.5);
        T(i) = 0.833-0.3507*(x(i)-1.5);
    end
    
    % For finding the floating variable of velocity, we use the mass flux
    % variable U2 which is nothing but 'm_dot'. Where m_dot is nothing but
    % the initial value of mass flow rate.
    
    % m_initial(i) = rho(i)*A(i)*V(i)
    % V(i) = m_initial(i)/(rho(i)*A(i))
    
    m_initial(i) = 0.59; % Therotical value of mass flow rate
    
    V(i) = m_initial(i)/(rho(i).*A(i));
    
end

%% Speed of sound and time step
a = sqrt(T);
nt = 1400;          % Total no. of time steps

%% Derived flux variable of the nozzle

U1 = rho.*A;
U2 = rho.*A.*V;
U3 = rho.*A.*((T/(gamma-1))+(gamma/2)*V.^2);

%% Outer time Loop (Solving using McCormack Method)

for j = 1:nt
    
    % The governing equation of Continuity, Momentum, and Energy has been
    % expressed in terms of flux variables
    
    % Continuity Equation   =   dU1/dt+dF1/dx=0
    % Momentum Equation     =   dU2/dt+dF2/dx=J2
    % Energy Equation       =   dU3/dt+dF3/dx=0
    
    % Saving the old values
        U1_old = U1;
        U2_old = U2;
        U3_old = U3;
        
        F1 = U2;
        F2 = (U2.^2./U1)+(((gamma-1)/gamma)*(U3-((gamma*0.5*(U2.^2))./U1)));
        F3 = ((gamma*U2.*U3)./U1)-((gamma*0.5*(gamma-1)*(U2.^3))./(U1.^2));
        
        %% Time Step calculation by Courant no.
        %nt = 1400; % Total no. of time steps
        C = 0.5;
        %for i = 1:n
            dt_1 = C*(dx./(a+V));
        %end

        dt = min(abs(dt_1));

                
        % Predictor Method
        for k = 2:(n-1)         
            
            J2(k) = (1/gamma)*rho(k)*T(k)*((A(k+1)-A(k))/dx);
            
            % Continuity Equation            
            dU1_dt_p(k) = -(F1(k+1)-F1(k))/dx;
            
            % Momentum Equation
            dU2_dt_p(k) = -(F2(k+1)-F2(k))/dx+J2(k);
            
            % Energy Equation
            dU3_dt_p(k) = -(F3(k+1)-F3(k))/dx;
            
            
            % Solution Update
            U1(k) = U1(k) + dU1_dt_p(k)*dt;
            U2(k) = U2(k) + dU2_dt_p(k)*dt;
            U3(k) = U3(k) + dU3_dt_p(k)*dt;
            
            rho(k) = U1(k)/A(k);
            V(k) = U2(k)/U1(k);
            T(k) = (gamma-1)*((U3(k)/U1(k)) - ((gamma/2)*(V(k)^2)));
        end
        
        F1 = U2;
        F2 = (U2.^2./U1)+(((gamma-1)/gamma)*(U3-((gamma*0.5*(U2.^2))./U1)));
        F3 = ((gamma*U2.*U3)./U1)-((gamma*0.5*(gamma-1)*(U2.^3))./(U1.^2));
        
        
        % Corrector Method
                
        for k = 2:(n-1)      
                 
            J2(k) = (1/gamma)*rho(k)*T(k)*((A(k)-A(k-1))/dx);
            
            % Continuity Equation            
            dU1_dt_c(k) = -(F1(k)-F1(k-1))/dx;
            
            % Momentum Equation
            dU2_dt_c(k) = -(F2(k)-F2(k-1))/dx+J2(k);
            
            % Energy Equation
            dU3_dt_c(k) = -(F3(k)-F3(k-1))/dx;
        end
        
        % Computing the average time derivative
        
        dU1_dt = 0.5*(dU1_dt_p+dU1_dt_c);
        dU2_dt = 0.5*(dU2_dt_p+dU2_dt_c);
        dU3_dt = 0.5*(dU3_dt_p+dU3_dt_c);
        
        % Performing the final solution update
        
        for i = 2:(n-1)
            U1(i) = U1_old(i) + dU1_dt(i)*dt;
            U2(i) = U2_old(i) + dU2_dt(i)*dt;
            U3(i) = U3_old(i) + dU3_dt(i)*dt;
        end
        
        % Defining the boundary conditions
        % Inlet
        U1(1) = rho(1)*A(1);
        U2(1) = 2*U2(2)-U2(3);
        U3(1) = U1(1)*((T(1)/(gamma-1))+(0.5*gamma*(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);
        
        % Density Update
        rho = U1./A;
        
        % Velocity Update
        V = U2./U1;
        
        % Temperature Update
        T = (gamma-1)*((U3./U1) - ((gamma/2)*(V.^2)));
        
        % Calculating the Mach No.
        a = sqrt(T);
        M = V./a;
        p = rho.*T;

        % Calculating the mass flow rate along the length of the nozzle
        m = rho.*A.*V;
        
        % Calculation of primitive variables at throat
        rho_t(j) = rho(16);
        V_t(j) = V(16);
        T_t(j) = T(16);
        p_t(j) = p(16);
        M_t(j) = M(16);
        m_t(j) = m(16);
        
        % Plotting
    
    figure(6)
    if j==50    
        plot(x,m,'r','linewidth',1.5)
        hold on
    else if j==100
            plot(x,m,'g','linewidth',1.5)
            hold on
        else if j==200
                plot(x,m,'b','linewidth',1.5)
                hold on
            else if j==300
                    plot(x,m,'black','linewidth',1.5)
                    hold on
                else if j==400
                        plot(x,m,'y','linewidth',1.5)
                        hold on
                    else if j==500
                            plot(x,m,'linewidth',1.5)
                            hold on
                        else if j==600
                                plot(x,m,'linewidth',1.5)
                                hold on
                            else if j==700
                                    plot(x,m,'linewidth',1.5)
                                    hold on
                                
                                    plot(x,m_initial,'--','linewidth',1.5)
                                    
                                xlabel('Length of nozzle');
                                ylabel('Mass Flow Rate');
                                legend('50dt','100dt','200dt','300dt','400dt','500dt','600dt','700dt','Initial Mass flow rate')
                                title('Mass Flow rate from the Nozzle for Conservative  part')
                                end
                            end
                        end
                    end
                end
            end
        end
    end
end
                        
figure(1)
plot(x,V,'r')
xlabel('Length of Nozzle');
ylabel('Fluid Velocity');
legend('VELOCITY');
title(sprintf('Time Step: %d',j));
    
figure(2)
plot(x,rho,'g')
xlabel('Length of Nozzle');
ylabel('Fluid Density');
legend('DENSITY');
title(sprintf('Time Step: %d',j));
    
figure(3)
plot(x,T,'b')
xlabel('Length of Nozzle');
ylabel('Fluid Temperature');
legend('TEMPERATURE');
title(sprintf('Time Step: %d',j));
    
figure(4)
plot(x,M,'black')
xlabel('Length of Nozzle');
ylabel('Mach Number');
legend('MACH NUMBER');
title(sprintf('Time Step: %d',j));
    
figure(5)
plot(x,m_initial,'linewidth',1.5)
hold on
plot(x,m,'black','linewidth',1.5)
xlabel('Length of Nozzle');
ylabel('Fluid Mass Flow Rate');
axis([0 1 0 1]);
legend('Initial Mass Flow Rate','Final Mass flow rate');
title(sprintf('Time Step: %d',j));
    
pause(0.01)
hold off

figure(7)
subplot(511)
plot(1:j,V_t,'linewidth',1.5)
xlabel('No. of Time Step')
ylabel('Velocity at Throat')
txt1 = (sprintf('Time Step: %d',j));
txt2 = ('Time-wise variation of the primitive variables');
title({txt1,txt2})
    
subplot(512)
plot(1:j,rho_t,'linewidth',1.5)
xlabel('No. of Time Step')
ylabel('Density at Throat')
    
subplot(513)
plot(1:j,T_t,'linewidth',1.5)
xlabel('No. of Time Step')
ylabel('Temperature at Throat')

subplot(514)
plot(1:j,p_t,'linewidth',1.5)
xlabel('No. of Time Step')
ylabel('Pressure at Throat')

subplot(515)
plot(1:j,M_t,'linewidth',1.5)
xlabel('No. of Time Step')
ylabel('Mach Number at Throat')

Results:

Non-Conservative form of Governing Equations at n = 31:

Mass flow rate from the nozzle for the non-conservative part at various time steps:

Comparisons of initial and Final mass flow rate:

DENSITY:

MACH NUMBER:

TEMPERATURE:

VELOCITY:

Conservative form of Governing Equations at n = 31:

Mass flow rate from the nozzle for the conservative part at various time steps:

DENSITY:

MACH NUMBER:

TEMPERATURE:

VELOCITY:

Discussion and Conclussion:

1. Mass flow distribution in conservative form is much closer to the analytical value and it is somewhat constant. But in the case of non - conservative form, there is a sizeable variation in mass flow distribution.
2. The non - conservative form leads to a smaller residual. Hence, in terms of the quality of the algorithm, non - conservative form does a better job.
3. Both the form of governing differential equations are satisfactory if we compare them with analytical form.
4. It almost takes 600 iterations to have a steady-state solution.
5. The solution is grid-independent.