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

INTRODUCTION:

We consider the steady, isentropic flow through a convergent-divergent nozzle. The flow at the inlet to the nozzle comes from a reservoir where pressure and temperatures are at stagnation. Also it is assumed that the crosss sectional area of the nozzle is very large so that the velocity is almost V ~ 0. The flow expands adiabatically to supersonic speeds at the nozzle exit. the flow is locally subsonic in the convergent section (M < 1), sonic at the throat section (M = 1) and supersonic at the divergent section ( M> 1).

We assume that at a given section, where the cross sectional area is A, the flow properties are uniform across the section. Hence although the area of the nozzle changes as the function of distance along the nozzle x, and therefore in reality the flow field is two-dimensional(the flow varies in two dimensional xy space), we make the assumption that the flow properties vary only with x; this is tantamount to assuming uniform flow properties across any given cross section.Such flow is defined as Quasi one dimensional flow.

MacCormack technique:

The MacCormack method is an explicit finite-difference technique which is second order accurate in both space and time. This method is a popular technique as it provides satisfactory results for many fluid flow applications.MacCormack method is a variant of the Lax- Wendroff but is much simpler in its application because there is no need to evaluate the second order time derivatives with was the case for Lax-Wendroff method..MacCormack method is a two-step predictor- corrector sequence is used with forward differences on the predictor and with rearward differences on the corrector is a second order accurate method.
As we assume that the flow field at each grid point in time t is known, the flow field variables at each grid point at time `t+trianglet` can be calculated using maccormack method in a predictor-corrector sequence.

Main Program

clear all
close all
clc

n = 31;
x = linspace(0,3,n);
dx = x(2)-x(1);
gamma = 1.4;
nt = 1400;
C = 0.5;
throat = (n-1)/2;
a = (1+2.2*(x-1.5).^2);

tic;
[rho,v,T,P,m,M,T_throat,rho_throat,P_throat,M_throat] = non_conserv(n,x,dx,gamma,nt,C,a,throat);
time_elapsed = toc;
fprintf('computation time = %f',time_elapsed);

time_variation = figure(1);

subplot(4,1,1)
plot(rho_throat,'b')
title('Time-wise variations in quantity in non conservative form');
legend('density')
ylabel('rho/rho_{o}');
grid minor

subplot(4,1,2)
plot(T_throat,'b')
legend('temp')
ylabel('T/T_{o}');
grid minor

subplot(4,1,3)
plot(P_throat,'b')
legend('pressure')
ylabel('p/p_{o}');
grid minor

subplot(4,1,4)
plot(M_throat,'b')
legend('mach no')
grid minor
ylabel('M');
xlabel('no.of iter')

qualitative_non_cons = figure(2);
subplot(4,1,1)
plot(x,M,'b')
title('qualitative aspects of 1D supersonic nozzle:non-conservative form');
ylabel('M')
legend('mach no.')
grid minor

subplot(4,1,2)
plot(x,T,'b')
ylabel('T')
legend('temp')
grid minor

subplot(4,1,3)
plot(x,P,'b')
ylabel('P')
legend('press')
grid minor

subplot(4,1,4)
plot(x,rho,'b')
ylabel('rho')
legend('density')
grid minor

% conservation form
tic
[m_c,v_c,T_c,P_c,M_c,rho_c,T_throat_c,rho_throat_c,P_throat_c,M_throat_c] = conserv_form(x,dx,n,C,a,nt,gamma,throat);
time_elapsed = toc;
time_variation_conser = figure(3);

subplot(4,1,1)
plot(rho_throat_c,'b')
title('Time-wise variations in conservative form');
axis([0 1500 0.6 1.3])
legend('density')
ylabel('rho/rho_{o}');
grid minor

subplot(4,1,2)
plot(T_throat_c,'b')
axis([0 1500 0.6 1])
legend('temp')
ylabel('T/T_{o}');
grid minor

subplot(4,1,3)
plot(P_throat_c,'b')
axis([0 1500 0.4 1.1])
legend('pressure')
ylabel('p/p_{o}');
grid minor

subplot(4,1,4)
plot(M_throat_c,'b')
axis([0 1500 0.6 1.4])
legend('mach no')
ylabel('M');
grid minor
xlabel('no.of iter')

qualitative_cons = figure(4);
subplot(4,1,1)
plot(x,M_c,'b')
title('qualitative aspects of 1D supersonic nozzle flow: Conservative form')
ylabel('M')
legend('mach no.')

subplot(4,1,2)
plot(x,T_c,'b')
ylabel('T/T_{o}');
legend('temp')

subplot(4,1,3)
plot(x,P_c,'b')
ylabel('p/p_{o}');
legend('press')

subplot(4,1,4)
plot(x,rho_c,'b')
ylabel('rho/rho_{o}');
legend('density')

%comparison of mass flow artes by conservative and non-conservative
%methods and comaprison with exact form
mass_comparison = figure(5);
hold on;
plot(x,m_c,'g','LineWidth',1.75);
plot(x,m,'m','LineWidth',1.75);
line([0 3],[0.579 0.579],'LineWidth',1.75);
legend('Conservative Form','Non-Conservative Form','Exact Solution');
grid minor;
xlabel('x/L');
ylabel('rhoVA/rho_{o}V_{o}A_{0}');
title('non-dimensional mass flow comparison with exact solution');

Conservative form:

function [m_c,M_c,v_c,T_c,P_c,rho_c,T_throat_c,rho_throat_c,P_throat_c,M_throat_c] = conserv_form(x,dx,n,C,a,nt,gamma,throat)

for i = 1:length(x)
    if (x(i) >=0 && x(i) < 0.5)
        rho_c(i) = 1;
        T_c(i) = 1;
    else if (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);
        else if (x(i) >= 1.5 && x(i) <= 3)
                rho_c(i) = 0.634 - 0.3879*(x(i) - 1.5);
                T_c(i) = 0.833 - 0.3507*(x(i) - 1.5);
            end
        end
    end
end

v_c = 0.59./(rho_c.*a);
P_c = rho_c.*T_c;

U1 = rho_c.*a;
U2 = rho_c.*a.*v_c;
U3 = rho_c.*a.*((T_c./(gamma-1))+(gamma/2)*(v_c.^2));

for k = 1:1400
    
    U1_old = U1;
    U2_old = U2;
    U3_old = U3;
    
     dt = min(C*(dx./(sqrt(T_c)+v_c)));
     
     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
         
         J2(j) = (1/gamma)*rho_c(j)*T_c(j)*((a(j+1)-a(j))/dx);
         
         %predictor step
         dU1_dt_p(j) = -(F1(j+1)-F1(j))/dx;
         dU2_dt_p(j) = -((F2(j+1) - F2(j))/dx) + J2(j);
         dU3_dt_p(j) = -(F3(j+1)-F3(j))/dx;
         
         U1(j) = U1(j) + dU1_dt_p(j)*dt;
         U2(j) = U2(j) + dU2_dt_p(j)*dt;
         U3(j) = U3(j) + dU3_dt_p(j)*dt;
     end
     
     % updating vriables and flux
      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)/2)*((U2.^3)./(U1.^2));
     
      for j = 2:n-1
         
         J2(j) = (1/gamma)*rho_c(j)*T_c(j)*((a(j)-a(j-1))/dx);
         
         %corrector step
         dU1_dt_c(j) = -(F1(j)-F1(j-1))/dx;
         dU2_dt_c(j) = -((F2(j) - F2(j-1))/dx) + J2(j);
         dU3_dt_c(j) = -(F3(j)-F3(j-1))/dx;
      end
      
      %updating the sol
      %avg val
      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);
       
       % final sol 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
       
       % flow field variables update
         rho_c = U1./a;
         v_c = U2./U1;
         T_c = (gamma-1)*((U3./U1) - ((gamma/2)*(v_c.^2)));
         
       %boundary cond
       %inlet
       U1(1) = rho_c(1)*a(1);
       U2(1) = 2*U2(2) - U2(3);
       U3(1) =  U1(1)*((T_c(1)/(gamma-1))+(0.5*gamma*(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);
         

         
         % p,mach,mass
         P_c = rho_c.*T_c;
         M_c = v_c./(sqrt(T_c));
         m_c = rho_c.*a.*v_c;
         
         % @throat
         M_throat_c(k) = M_c(throat);
         P_throat_c(k) = P_c(throat);
         T_throat_c(k) = T_c(throat);
         rho_throat_c(k) = rho_c(throat);
         
         %comparison of non-dimensional mass flow rates at diffferent
         %times:
         mass2 = figure(7);
         hold on;
         if k ==100
             plot(x,m_c,'b','linewidth',1.75);
         else if k ==200
                 plot(x,m_c,'g','linewidth',1.75);
             else if k ==700
                     plot(x,m_c,'r','linewidth',1.75);
                 end
             end
         end
end
title('Comparison of non-dimensional mass flow rates at different time-instants(conservative form');
xlabel('x/L (Non-dimensional distance through the nozzle)');
ylabel('rhoVA/rho_{o}V_{o}A_{o}'); 
legend('1000Deltat','200Deltat','700Deltat'); 
axis([0 3 0.4 1.6]); 
grid on; 
grid minor; 
end
                 
       
                              

Non-conservative form:

function [rho,v,T,P,m,M,T_throat,rho_throat,P_throat,M_throat] = non_conserv(n,x,dx,gamma,nt,C,a,throat)


rho = 1-0.3146*x;
T = 1-0.2314*x;
v = (0.1+1.09*x).*T.^0.5;

for k = 1:nt
    
    rho_old = rho;
    v_old = v;
    T_old = T;
    
    dt = min(C.*dx./(sqrt(T)+v));
    
    % predictor method

    for j = 2:n-1
        
        dv_dx = (v(j+1)-v(j))/dx;
        drho_dx = (rho(j+1)-rho(j))/dx;
        dT_dx = (T(j+1)-T(j))/dx;
        dloga_dx = (log(a(j+1))-log(a(j)))/dx;
        
        % continuity eq
        drho_dt_p(j) = -rho(j)*dv_dx-rho(j)*v(j)*dloga_dx-v(j)*drho_dx;
        
        %momentum eq
        dv_dt_p(j) = -v(j)*dv_dx-(1/gamma)*(dT_dx+(T(j)/rho(j))*drho_dx);
        
        % energy eq
        dT_dt_p(j) = -v(j)*dT_dx-(gamma-1)*T(j)*(dv_dx+v(j)*dloga_dx);
        
        % solution update
    rho(j) = rho(j) + drho_dt_p(j)*dt;
    v(j) = v(j) + dv_dt_p(j)*dt;
     T(j) = T(j) + dT_dt_p(j)*dt;
    
    end


% corrector method
    
    % continuity eq
    
    for j = 2:n-1
        
        dv_dx = (v(j)-v(j-1))/dx;
        drho_dx = (rho(j)-rho(j-1))/dx;
        dT_dx = (T(j)-T(j-1))/dx;
        dloga_dx = (log(a(j))-log(a(j-1)))/dx;
        
        % continuity eq
        drho_dt_c(j) = -rho(j)*dv_dx-rho(j)*v(j)*dloga_dx-v(j)*drho_dx;
        
        %momentum eq
        dv_dt_c(j) = -v(j)*dv_dx-(1/gamma)*(dT_dx+(T(j)/rho(j))*drho_dx);
        
        % energy eq
        dT_dt_c(j) = -v(j)*dT_dx-(gamma-1)*T(j)*(dv_dx+v(j)*dloga_dx);
        
    end
   
    
    % compute the average time derivative
    
    drho_dt = 0.5*(drho_dt_p+drho_dt_c);
    dv_dt = 0.5*(dv_dt_p+dv_dt_c);
    dT_dt = 0.5*(dT_dt_p+dT_dt_c);
    
    % final solution update
    
    for i = 2:n-1
    rho(i) = rho_old(i) + drho_dt(i)*dt;
    v(i) = v_old(i) + dv_dt(i)*dt;    
    T(i) = T_old(i) + dT_dt(i)*dt;   
    end
    
    % apply boundary conditions
    %inlet
    v(1) = 2*v(2)-v(3);
    rho(1) = 1;
    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);
    
  % defining the values for p,m,
  P = rho.*T;
  m = rho.*a.*v;
  M = v./sqrt(T);
  
  % @throat
  T_throat(k) = T(throat);
  P_throat(k) = P(throat);
  rho_throat(k) = rho(throat);
  M_throat(k) = M(throat);
    
% comparison of non-dimensional mass flow rates at different times:
mass = figure(6);
hold on;
if k == 50
    plot(x,m,'b','linewidth',1.75);
else if k == 100
        plot(x,m,'g','linewidth',1.75);
    else if k == 150
            plot(x,m,'r','linewidth',1.75);
        else if k == 200
                plot(x,m,'m','linewidth',1.75);
            else if k == 700
                    plot(x,m,'-','linewidth',1.75);
                    plot(x(end),m(end),'b*');
                end
            end
        end
    end
end
end

title('Comparison of non-dimensional mass flow rates at different time-instants');
xlabel('x/L (Non-dimensional distance through the nozzle)');
ylabel('rhoVA/rho_{o}V_{o}A_{o}'); 
legend('50Deltat','100Deltat','150Deltat','200Deltat','700Deltat'); 
axis([0 3 0 2.5]); 
grid on; 
grid minor; 
end



                
        
    
    
    

Results:

1.Time wise variation of pressure,temperature,velocity and mach number at nozzle throat.

a) Conservative form

b) Non-conservative form

2. Qualitative aspects of quasi 1D nozzle flow at steady state.

a) conservative form

b) non-conservative form

3. Non dimensional mass flow rate distribution at various time steps

a) Conservative form

b)non-conservative form

4. Comparison of  non dimensional mass flow rates of both forms

As evident from the plot the conservative form gives a mass flowrate variation that is close to a constant value and simultaneously close to the exact value of 0.579 at the steady state. The non-conservative form shows some oscillations around the inflow and outflow boundaries and is less satisfactory compared to the conservation form.

5. Grid independency test

From the above results we can see that for a grid-size of 30 points the Non-Conservative form gives a more accurate result compared to the Conservative form. However when we change the grid-size to 60 points the values of the flow quantites are quite close when calculated by either of the methods. Additionally both the methods give results which are closer to the exact value in this case.