BlogFor BusinessTalk to a counsellor
CoursesIIT ProgramsPlacementsLog in
BlogFor Business
Talk to a counsellor
Challenge

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

AAkshai A22 Jun 2022
MATLAB

Aim:-

To simulate Quasi 1D supersonic nozzle flow simulation using MacCormack method.

Objective:-

1. Steady-state distribution of primitive variables inside the nozzle
2. Time-wise variation of the primitive variables
3. Variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process
4. Comparison of Normalized mass flow rate distributions of both forms.

To simulate the isentropic flow through a Quasi 1D nozzle by deriving both the conservation and non conservation forms of the governing equations and solve them by using the MacCormacks method using MATLAB.

One dimensional, steady, isentropic flow through the Convergent-Divergent nozzle, it is assumed that the nozzle is connected to a large reservoir, such that its area is infinite. Velocity at the inlet is zero or negligible. Pressure and temperature at the inlet of the nozzle are P0 and T0 .It is assumed that the flow properties are varying only along the X-direction and constant along Y- direction.

Nozzle diagram with sections:-

  • Convergent Section &
  • Divergent Section.

Governing Equations:-

Soving the above equations analytically we get the following results shown below;

Governing Equations in Non- Conservative form is shown below;

Governing Equations in Conservative form is shown below;

Non Dimensional quantities are;

  • Temperature:-

  • Density:-

  • Pressure:-

  • Non dimensional distance:-

  • Non dimensional Velocity:-

(a0 is the local sonic velocity).

  • Non dimensional time:-

(L/a0 is the Characteristic time scale).

  • Non dimensional Area:-

To Discretize the numerical solution the technique that iam used is MacCormacks method. The MacCormack method is a variation of the two-step Lax-Wendroff scheme and simpler in application. The method proceeds in two steps called Predictor step and which is followed by Corrector step.

Predictor step:-

In predictor step the spatial derivatives are replaced with forward differences..Discritizaton and solution of governing equations by first order forward differences.

ie; To find

then to obtain the predicted values of temperature, density and velocity;

Corrector Step:-

Discritizaton and solution of governing equations by first order forward differences.

Compute the average time derivative of temperature, densiy and velocity. The average time derivatives are given by;

ie;

To obtain the corrected values of temperature, density and velocity;

Time step- calculation:-

The time step is calculated for all time steps using the following expressions,

(C = CFL number)

Nozzle Shape:-

Nozzle shape is varying parabolicall area distribution is given as w.r.to x as A = 1 + 2.2(x-1.5)2 , 0 ≤ x ≥ 3

Initial conditions:-

Initial conditions for density, velocity and temperature is calculated from the below showing expressions such as;

Boundary Conditions:-

  • The subsonic inflow boundary condition (first node);

  • The Supersonic outflow boundary condition(last node);

Now let us define the elements of solution vector U, the flux vector F and the source vecot J .

while substituting the above equations in the governing differential equation;

we get;

From the definitions of U1,U2 &U3;

Where;

Boundary Conditions;

  • For subsonic inflow boundary condition(First Node)

  • For Supersonic outflow boundary condition(Last node);

Matlab program for Non-Conservative form

clear 
close all
clc

% Input Values
% No of grid points
n = 31;
% Length of domain
l = 3; 
x = linspace(0,l,n);
dx = x(2) - x(1);
gamma = 1.4;

% Initial conditions

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


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

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

dt = min(abs(dt_1));

% Theoretical 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 
           
          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 value 


         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 derivatives

     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,'g','linewidth',1.5)
         hold on
     else if j==100
             plot(x,m,'r','linewidth',1.5)
             hold on
     else if j==200
             plot(x,m,'black','linewidth',1.5)
             hold on
     else if j==300
             plot(x,m,'b','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,'g')
xlabel('Length of Nozzle');
ylabel('Fluid Velocity');
legend('VELOCITY');
title(sprintf('Time Step: %d',j));

figure(2)
plot(x,rho,'r')
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,'y')
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')

Outputs:-

Matlab program for Conservative form

clear 
close all
clc

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

% Initial Conditions

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

% Primitive variables of the Nozzle

% The primitive variable Temperature and Densiy 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 values 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; % Theoretical 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 interms 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 variable 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,'g')
xlabel('Length of Nozzle');
ylabel('Fluid Velocity');
legend('VELOCITY');
title(sprintf('Time Step: %d',j));

figure(2)
plot(x,rho,'r')
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,'y')
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')

Outputs:-

GRID INDEPENDENCE TEST:-

Conclusion:-

  • Here while anlysing I can understand that the non conservative form results are distinctly closer to the exact analytical solution results and while doubling the number of grid points, the conservative form of results are slightly closer to exact analytical results.
  • Mass flow distribution in conservative form is much closer to the analytical value and mostly some what constant, but in non conservative form there is a variation in mass flow distribution.
  • Non conservative forms leads to small residuals, however the computation time is very less.
  • The results in both cases are certainly satisfactory, ie , both the governing equations are satisfactory while comparing it with analytical form.
Your portfolio, next

Ship work like this
in your own course.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.