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

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

OBJECTIVE:- To write code solve the 1D supersonic nozzle flow equations using the Macormack Method.

• Implementing both the conservative and non-conservative forms of the governing equations.
• Performing the grid dependence test.

TO SHOW:-

• To write separate functions for conservative and non-conservative forms.
• To figure out the minimum number of cycles for which the simulation should be run in order for convergence.
• Comparing the normalized mass flow rate between the conservative forms and non-conservative forms - make suitable conclusions
• Which method is faster?

SOLUTIONS:-

In this problem we consider a super-sonic isentropic flow through a convergent-divergent nozzle.

• The initial flow at the inlet to the nozzle comes from a reservoir where the pressure and temperature are represemted as P0 & T0, respectively.
• As the crossectional area of the tank/reservoir is large, hence velocity less. The flow expands isentropically to supersonics speeds at the exit of nozzle.
• Here the flow at inlet is subsonic(M<1) & the flow at outlet is supersonic (M>1). As given in the above figure A is the crossectional area, the flow properties are uniform across the section of nozzle.
• The area of the nozzle changes as a function of distance along the nozzle,x. and hence the flow varies in two dimensional space . The convergence section is connected with reservoir , hence the Mach number is less than or equal to 1.
• The governing continuity, momentum & energy equations for the quasi-one-dimensional, steady, isentropic flow is expressed as;

Where A- Area, rho = density, T-Temperature, V = Velocity.

• For solving a 1D super sonic nozzle flow equations , Macormack method is used to solve the governing flow equations. This method is widely used for discretization scheme for the numerical solution of a hyperbolic partial differential equations.
• In this method , we divide the Time Marching Technique into two phases i.e.

• Predictor step.
• Corrector step.

• In the Predictor step, forward differencing is done for space terms where as for Corrector step we follow backward differencing for the space terms. Further, we take average value of both the terms.
• For this problem we should implement both conservative and non-conservative forms of the governing equations and compare them .

Programme for implementing the non-conservative forms of the governing equations for solving the 1D supersonic nozzle flow equations using the Macormack Method.

% Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
% Non-Conservative form of the governing equations.


clear variables 
close all
clc

%Basic Inputs.
n = 31;
x = linspace(0,3,n);
dx = x(2) - x(1);
gamma = 1.4;
error = 1;
throat = (n + 1)/2;  % initialize throat.

%Calculate Initial Profiles.
rho = 1-0.3146*x; % Initial value of density.
t = 1-0.2314*x; % Initial value of temperature.
v = (0.1 + 1.09*x).*t.^0.5;  % Initial value of velocity. Here '.*t' for element wiseoperator as its an array and '.^0.5' to do element wise exponent.
a = 1+2.2*(x-1.5).^2;  % Area of nozzle.
%dt = min(0.5.*dx./(sqrt(t) + v));

% Time steps.
nt = 1200;   
dt = 0.01;
error_tolerance = 1e-4;


for k = 1:nt
    
    % We need the copy of old values. We require flow data to compare.
    rho_initial = rho;    
    v_initial= v;
    t_initial = t;
    
    % Predictor Step,
    
    for i=2:n-1
        
        dvdx = (v(i+1)-v(i))/dx;
        drhodx = (rho(i+1) - rho(i))/dx;
        dtdx = (t(i+1) - t(i))/dx;
        dlogadx = (log(a(i+1))-log(a(i)))/dx;
        
        % continuity equation:
        drhodt_p(i) = -rho(i)*dvdx -rho(i)*v(i)*dlogadx -v(i)*drhodx;
        % momentum equation:
        dvdt_p(i) = -v(i)*dvdx - (1/gamma)*(dtdx + (t(i)/rho(i))*drhodx);
        % Energy equation:
        dtdt_p(i) = -v(i)*dtdx - (gamma-1)*t(i)*(dvdx +(v(i)*dlogadx));
    
        % Solution Update;
        rho(i) = rho(i) + drhodt_p(i)*dt;
        v(i) = v(i) + dvdt_p(i)*dt;
        t(i) = t(i) + dtdt_p(i)*dt;
        
    end    
         
    
    
     % Corrector Method,  (backward differencing)
     for j=2:n-1
        
        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); 
     
     % solution update,
     for l =2:n-1
         rho(l) = rho_initial(l) + drhodt(l)*dt;
         v(l) = v_initial(l) + dvdt(l)*dt;
         t(l) = t_initial(l) +dtdt(l)*dt;
     end
     
     % Applying boundary conditions, 
     
     %Inlet
     rho(1) = 1;
     v(1) = 2*v(2) - v(3);
     t(1) = 1;
     
     %Outlet
     rho(n) = 2*rho(n-1) - rho(n-2); 
     v(n) = 2*v(n-1) - v(n-2);
     t(n) = 2*t(n-1) - t(n-2);
     
     % Calculating mass flow rate(m), Mach number(M) & Pressure(P)
     
     % Pressure(P)
     P =rho.*t;
     % Mach number
     M = (v./(t).^0.5);
     % mass flow rate(m)
     m = rho.*v.*a;
     
     % To capture values at the throat of nozzle for each iteration.
     rho_throat(k) = rho(throat);
     v_throat(k) = v(throat);
     t_throat(k) = t(throat);
     P_throat(k) = P(throat);
     M_throat(k) = M(throat);
     m_throat(k) = m(throat);
     
     
     % Plotting mass flow rate 
     figure(1)
     hold all
     grid on
     
      if k == 1
         plot(x,m,'b');
     elseif k == 50
            plot(x,m,'g');
     elseif k == 100
            plot(x,m,'r');
     elseif k == 150
            plot(x,m,'m');
     elseif k == 200
            plot(x,m,'y');
     elseif k == 700
            plot(x,m,'p');
             
             xlabel('Non-dimensional distance through nozzle')
             ylabel('Non-dimensional Mass flow rate')
             legend('1-Time step', '50-Time step', '100-Time step', '150-Time step', '200-Time step','700-Time step');
             title(' Mass flow rate at different time step- Non conservative form')
         
     end
     
end

figure(2)
subplot(2,3,1)
plot(x,rho,'r');
xlabel('Length of nozzle')
ylabel('Density')
title('Non-Dimensional Density (Non conservative form)')
grid minor

subplot(2,3,2)
plot(x,t,'b');
xlabel('Length of nozzle')
ylabel('Temperature')
title('Non-Dimensional Temperature (Non conservative form)')
grid minor

subplot(2,3,3)
plot(x,P,'g');
xlabel('Length of nozzle')
ylabel('Pressure')
title('Non-Dimensional Pressure (Non conservative form)')
grid minor

subplot(2,3,4)
plot(x,v,'m');
xlabel('Length of nozzle')
ylabel('velocity')
title('Non-Dimensional velocity (Non conservative form)')
grid minor

subplot(2,3,5)
plot(x,M,'y');
xlabel('Length of nozzle')
ylabel('Mach number')
title('Mach number (Non conservative form)')
grid minor

subplot(2,3,6)
plot(x,m,'k');
xlabel('Length of nozzle')
ylabel('Mass flow rate')
title('Mass flow rate (Non conservative form)')
grid minor

figure (3)
%title('Timewise variations of the Density, Temperature, Pressure, Mach number at Nozzle Throat')
subplot(6,1,1)
plot(1:nt,rho_throat,'r')
xlabel('Time steps')
ylabel('Density')
title('Variation of Density at the nozzle throat (Non conservative form)')
grid minor

subplot(6,1,2)
plot(1:nt,t_throat,'b')
xlabel('Time steps')
ylabel('Temperature')
title('Variation of Temperature at the nozzle throat (Non conservative form)')
grid minor

subplot(6,1,3)
plot(1:nt,P_throat,'g')
xlabel('Time steps')
ylabel('Pressure')
title('Variation of Pressure at the nozzle throat (Non conservative form)')
grid minor

subplot(6,1,4)
plot(1:nt,P_throat,'m')
xlabel('Time steps')
ylabel('Velocity')
title('Variation of Velocity at the nozzle throat (Non conservative form)')
grid minor

subplot(6,1,5)
plot(1:nt,M_throat,'y')
xlabel('Time steps')
ylabel('Mach number')
title('Variation of Mach number at the nozzle throat (Non conservative form)')
grid minor

subplot(6,1,6)
plot(1:nt,m_throat,'k')
xlabel('Time steps')
ylabel('Mass flow rate')
title('Variation of Mass flow rate at the nozzle throat (Non conservative form)')
grid minor

figure(4)
hold all
grid on
plot(x,m,'r')
plot(x,linspace(0.579,0.579,n),'g')
xlabel('Non dimensional Length')
ylabel('Mass flow rate')
legend('Mass flow rate in Non conservative form',' Analytical solution')

Programme for implementing the conservative forms of the governing equations for solving the 1D supersonic nozzle flow equations using the Macormack Method.

% Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
% Conservative form of the governing equations.


clear variables 
close all
clc

% Basic Inputs.
n = 31;
x = linspace(0,3,n);
dx = x(2) - x(1);
gamma = 1.4;
error = 1;
%  Time steps.
nt = 1400;   
dt = 0.01;
error_tolerance = 1e-4;
throat = (n + 1)/2;  % initialize throat.

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.5
          rho(i) = 0.634-(0.3879*(x(i)-015));
          t(i) = 0.833-(0.3507*(x(i)-1.5));
    end
end
    %Calculate Initial Profiles.
    a = 1+2.2*(x-1.5).^2;  % Area of nozzle
    rho = 1-0.3146*x; % Initial value of density.
    t = 1-0.2314*x; % Initial value of temperature.
    v = 0.59./(rho.*a);  % Initial value of velocity. Here '.*t' for element wiseoperator as its an array and '.^0.5' to do element wise exponent.
    P = rho.*t;


    % Defining elements of solution vector U,
    U1 = rho.*a;
    U2 = rho.*a.*v;
    U3 = rho.*a.*((t./(gamma - 1)) + ((gamma/2).*(v.^2)));
 
%Outer time loop

for k = 1:nt
    
    % We need the copy of old values.
    U1_old = U1;    
    U2_old = U2;
    U3_old = U3;
    
    
     % Defining elements of Flux vector F,
     F1 = U2;
     F2 = ((U2.^2)./(U1)) + ((gamma-1)/gamma).*(U3 - (0.5*gamma.*U2.^2)./U1);
     F3 = (gamma.*U2.*U3./U1) - ((0.5*gamma*(gamma-1).*U2.^3)./(U1.^2));
     

     % predictor method, (Forward differencing)
     
     for i = 2:n-1
         
         dF1dx = (F1(i+1)-F1(i))/dx;
         dF2dx = (F2(i+1) -F2(i))/dx; 
         dF3dx = (F3(i+1) -F3(i))/dx; 
         dadx = (a(i+1) - a(i))/dx;
         
         
        
         % continuity equation:
         dU1dt_p(i) = -dF1dx;
         % momentum equation:
         J2(i) = (1/gamma).*rho(i).*t(i).*((a(i+1) - a(i))/dx);
         dU2dt_p(i) =-dF2dx + J2(i);   %Where J is the source term.
         % Energy equation:
         dU3dt_p(i) = -dF3dx;
         
         
         % Solution Update;
         U1(i) = U1(i) + dU1dt_p(i)*dt;
         U2(i) = U2(i) + dU2dt_p(i)*dt;
         U3(i) = U3(i) + dU3dt_p(i)*dt;   
     end
     
     % updating values,
     P = rho.*t;  % Non dimensional pressure
     v = U2./U1;
     rho = U1./a;
     t = (gamma -1).*((U3./U1) - ((gamma./2).*(v.^2)));
     
     % Updating elements of Flux vector F,
     F1 = U2;
     F2 = ((U2.^2)./(U1)) + ((gamma-1)/gamma).*(U3 - (0.5*gamma.*U2.^2)./U1);
     F3 = (gamma.*U2.*U3./U1) - ((0.5.*gamma.*(gamma-1).*U2.^3)./(U1.^2));
     
     % Corrector Method,  (backward differencing)
     
    for j = 2:n-1
         
         dF1dx = (F1(j)-F1(j-1))/dx;
         dF2dx = (F2(j) -F2(j-1))/dx; 
         dF3dx = (F3(j) -F3(j-1))/dx;
        
         % continuity equation:
         dU1dt_c(j) = -dF1dx;
         % momentum equation:
         J2(j) = (1/gamma).*rho(j).*t(j).*((a(j) - a(j-1))/dx);
         dU2dt_c(j) =-dF2dx + J2(j);
         % Energy equation:
         dU3dt_c(j) = -dF3dx;
     end
     
     % Compute the average time derivative
     dU1dt_average = 0.5.*(dU1dt_p + dU1dt_c);
     dU2dt_average = 0.5.*(dU2dt_p + dU2dt_c);
     dU3dt_average = 0.5.*(dU3dt_p + dU3dt_c);
     
     % Updating dependent variables U1,U2&U3.
     for l = 2:n-1
         U1(l) = U1_old(l) + dU1dt_average(l)*dt;
         U2(l) = U2_old(l) + dU2dt_average(l)*dt;
         U3(l) = U3_old(l) + dU3dt_average(l)*dt;
     end
     
     % Applying boundary conditions,
     
     %Inlet
     U2(1) = 2*U2(2) - U2(3);
    
     %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);
    
     % Updating the flow field variables,
     rho = U1./a;
     v = U2./U1;
     T = (gamma -1).*((U3./U1) - ((gamma*v.^2)/2));
     
     
     % Calculating mass flow rate(m), Mach number(M) & Pressure(P)
     
     % Pressure(P)
     P =rho.*t;
     % Mach number
     M = v./((t).^0.5);
     % mass flow rate(m)
     m = rho.*a.*v;
     
     % To capture values at the throat of nozzle for each iteration.
     rho_throat(k) = rho(throat);
     v_throat(k) = v(throat);
     t_throat(k) = t(throat);
     P_throat(k) = P(throat);
     M_throat(k) = M(throat);
     m_throat(k) = m(throat);


   % Plotting mass flow rate 
   figure(1)
   hold all
   grid on
   
     if k == 1
         plot(x,m,'b');
     elseif k == 50
            plot(x,m,'g');
     elseif k == 100
            plot(x,m,'r');
     elseif k == 150
            plot(x,m,'m');
     elseif k == 200
            plot(x,m,'y');
     elseif k == 700
            plot(x,m,'p');
     
             
             xlabel('Non-dimensional distance through nozzle')
             ylabel('Non-dimensional Mass flow rate')
             legend('1-Time step', '50-Time step', '100-Time step', '150-Time step', '200-Time step','700-Time step');
             title(' Mass flow rate at different time step-  conservative form')
         
     end
     
end   

figure(2)
subplot(2,3,1)
plot(x,rho,'b');
xlabel('Length of nozzle')
ylabel('Density')
title('Non-Dimensional Density (conservative form)')
grid minor

subplot(2,3,2)
plot(x,t,'y');
xlabel('Length of nozzle')
ylabel('Temperature')
title('Non-Dimensional Temperature (conservative form)')
grid minor

subplot(2,3,3)
plot(x,P,'r');
xlabel('Length of nozzle')
ylabel('Pressure')
title('Non-Dimensional Pressure (conservative form)')
grid minor

subplot(2,3,4)
plot(x,v,'g');
xlabel('Length of nozzle')
ylabel('velocity')
title('Non-Dimensional velocity (conservative form)')
grid minor

subplot(2,3,5)
plot(x,M,'k');
xlabel('Length of nozzle')
ylabel('Mach number')
title('Mach number (conservative form)')
grid minor

subplot(2,3,6)
plot(x,m,'m');
xlabel('Length of nozzle')
ylabel('Mass flow rate')
title('Mass flow rate (conservative form)')
grid minor



figure (3)
%title('Timewise variations of the Density, Temperature, Pressure, Mach number at Nozzle Throat')
subplot(6,1,1)
plot(1:nt,rho_throat,'b')
xlabel('Time steps')
ylabel('Density')
title('Variation of Density at the nozzle throat (conservative form)')
grid minor

subplot(6,1,2)
plot(1:nt,t_throat,'y')
xlabel('Time steps')
ylabel('Temperature')
title('Variation of Temperature at the nozzle throat (conservative form)')
grid minor

subplot(6,1,3)
plot(1:nt,P_throat,'r')
xlabel('Time steps')
ylabel('Pressure')
title('Variation of Pressure at the nozzle throat (conservative form)')
grid minor

subplot(6,1,4)
plot(1:nt,v_throat,'g')
xlabel('Time steps')
ylabel('Velocity')
title('Variation of Velocity at the nozzle throat (conservative form)')
grid minor

subplot(6,1,5)
plot(1:nt,M_throat,'k')
xlabel('Time steps')
ylabel('Mach number')
title('Variation of Mach number at the nozzle throat (conservative form)')
grid minor

subplot(6,1,6)
plot(1:nt,m_throat,'m')
xlabel('Time steps')
ylabel('Mass flow rate')
title('Variation of Mass flow rate at the nozzle throat (conservative form)')
grid minor

figure(4)
hold all
grid on
plot(x,m,'b')
plot(x,linspace(0.579,0.579,n),'g')
xlabel('Non dimensional Length')
ylabel('Mass flow rate')
legend('Mass flow rate in conservative form',' Analytical solution')

RESULTS:-

For n = 31;

Mass flow rate at different time steps for Non-conservative & conservative forms:

>> The non-dimensional mass flow rate vs non- dimensional distance through nozzle is plotted.

>> The mass flow rate at various time steps for a non-conservative method is observed that, untill 50th time step it is noticed that the mass flow rate is fluctuating extremely whereas after the 100th or 150th time steps suddenly it started dropping down due to the transient variation of the flow field variables.

>> In further stages after 200th time steps , the mass flow distribution started settling down slowly and at 700 time step its clearly observed that it started following a straight horizontal path.

>>Comparing the conservative method to the non conservative method it is noted that, ther wer no much deflections in the low as compared to non conservative method in the initial methods. It too little more time to stablize the flow through the nozzle.

Variation of variables across the length of nozzle for Non-conservative & conservative forms:

Variation of variable w.r.to. time steps at the Nozzle throat for Non-conservative & conservative forms:

>> On comparing both the methods the variations for all the properties are noted. It is observed that the variation is much lower for the conservative form as compared to non-conservative form. Initiall there are very small deviations but immediately they will come to stable position.

>> The variation of mass flow rate along the length of the nozzle is comparitively less for conservative form . Where as for the non-conservative flow during inlet and outlet some disturbances were observed i.e. we can see the oscillations at inlet and out flow flows for non-conservative form.

Comparing the mass flow rate of Non- conservative form with the analytical solution and the mass flow rate of conservative form with the analytical solution:

>> For the mass flow rate of the non conservative form, it is observed that there is too much of disturbances in the flow at inlet and outlet boundaries. When considering the conservative form, initially there is a huge drop and slowly it started following the constant path.

>> Hence compared to non-conservative form , results obtained for conservative form are comparitively closer to the analytical solution.

Grid Dependence Test:

Here by varying the grid points we can check the accuracy of the result .So as we increase the number of grid points , the accuracy of the result will also be increased. On increasing the grid points it is observed that it is reaching close to the actual /analytical solution.

For n = 61;

Variation of variables across the length of nozzle for Non-conservative & conservative forms:

Variation of variable w.r.to. time steps at the Nozzle throat for Non-conservative & conservative forms:

Comparing the mass flow rate of Non- conservative form with the analytical solution and the mass flow rate of conservative form with the analytical solution:

CONCLUSIONS:-

>> It is observed that there is good accuracy on increasing the grid points.

>> The mass flow rate is better in conservative form as it reaches close to the analytical solution when compared to the non- conservative form solution.

>> For mass flow rate in non- conservative form , the solution needs much higher time steps as at lower time steps its still having the disturbances at inflow and outflow.

>> As we know non-conservative form includes primitive variables and conservative form includes flux variables, both tend to give the better outputs for their variables at higher time steps with a very small difference between them.