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

Introduction

Your objective in this project is to write code solve the 1D supersonic nozzle flow equations using the Macormack Method.

• You will be implementing both the conservative and non-conservative forms of the governing equations
• You will perform grid dependence test

Questions

• 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?

The following assumptions are made the flow is steady and isentropic through a convergent-divergent nozzle. First of all its quasi 1d because there is no flow change in the y-direction and flow varies only in x-axis but in reality, it is 2-dimensional flow. The flow at the inlet is subsonic ( M<1 ) and the flow at the outlet is supersonic ( M>1). Moreover, temperature and pressure at the inlet are stagnated. Also, it is assumed the cross-sectional area is very large ( infinity) and the velocity is approximated to zero. For this simulation, we are using two approaches namely conservative form ( control volume stationary) and non-conservative form ( control volume moving). Indeed, the Mach number value at the throat section is one ( M = 1) it pretends to be sonic. For solving this problem we are using the MacCormack explicit technique with the help of the finite difference equation. We divide the x-axis into several discrete points the grid points along the x-axis is uniform because the flow is uniform all over the nozzle. The first grid point is assumed to be a reservoir and the last point namely the nozzle exit is denoted by N. The MacCormack technique is a predictor-corrector method.

Subsonic Supersonic Isentropic Nozzle:

The following variables are used

Temperature

`T\' = frac{T}{T0}`

Density

`rho\' = frac{rho}{rho0}`

Non-Dimensional Distance

`x\' = frac{x}{L}`

Non-Dimensional Velocity

`V\' = frac{V}{a0}`

Non-Dimensional Time

`t\' = frac{t}{ frac{t}{a0}}`

GOVERNING EQUATIONS

Conservative form

Continuity Equation

`frac{del (rho^\' A^\')}{delt^\'} + frac{del (rho^\' A^\' V^\')}{delx^\'} = 0`

Momentum Equation

`frac{del (rho^\' A^\')}{delt^\'} + frac{del (rho^\' A^\' V^\'2 + (frac{1}{gamma} ) (rho^\' A^\'))}{delx^\'} = - (frac{1}{gamma} )rho^\' frac{del A^\'}{delx^\'}`

Energy Equation

`frac{del (rho^\'(frac{e^\'}{gamma-1} + frac{gamma}{2} V^\'2) A^\') } {delt^\'} + frac{del (rho^\'(frac{e^\'}{gamma-1} + frac{gamma}{2} V^\'2) V^\' A^\' + rho^\' V^\' A^\') } {delx^\'} = 0`

Non-Conservative Form

Continuity Equation

`frac{del rho^\' }{delt^\'} = - rho^\' frac{del V^\' }{delx^\'} - rho^\' V^\' frac{del (ln A^\')}{delx^\'} - V^\' frac{del rho^\' }{delx^\'}`

Momentum Equation

`frac{del V^\' }{delt^\'} = - V^\' frac{del V^\' }{delx^\'} - frac{1 }{gamma} ( frac{del T^\' }{delx^\'}+ frac{T^\' }{rho^\'} frac{del rho^\' }{delx^\'})`

Energy Equation

`frac{del T^\' }{delt^\'} = - V^\' frac{del T^\' }{delx^\'} - (gamma -1) T^\'( frac{del V^\' }{delx^\'}+ V^\'frac{del (ln A^\') }{delx^\'})`

Let us define the elements of the solution vectors U, the flux vector F and the source term J as follows.

Solution Vectors

`U 1 =rho^\' A^\'`

`U2 =rho^\' A^\' V^\'`

`U3 = rho^\'(frac{e}{gamma - 1} )+ (frac{gamma}{2} ) (V^\'2) A^\'`

Flux Terms

`F1 =rho^\' A^\' V^\'`

`F2 =rho^\' A^\' V^\' + frac{1}{gamma} rho^\' A^\'`

`F3 = rho^\'(frac{e}{gamma - 1} )+ (frac{gamma}{2} ) (V^\'2) A^\' + rho^\' A^\' V^\'`

Source Terms

`J2 = -frac{1}{gamma} rho^\' frac{∂A^\'}{∂x^\'}`

Gridpoint distribution along with the nozzle

Main program

clc
clear all
close all

%inputs
n=31; %number of nodes
x=linspace(0,3,n);
dx=x(2)-x(1);
gamma=1.4;

throat= ((n-1)/2)+1;% Throat point

C = 0.5; % Courant Number (specified for the problem)

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

%Time step
nt=1400;
dt=0.01;

tic
[rho,T,v,P,M,m,rho_throat,T_throat,v_throat,P_throat,M_throat,m_throat]= non_conservative_form(n,x,dx,gamma,throat,nt,dt,A);

Time_non_conservative=toc;
fprintf(\'Execution time for Non-Conservative form is %0.3g Sec. \\n\' ,Time_non_conservative)

tic
[rho_c,T_c,v_c,P_c,M_c,m_c,rho_throat_c,T_throat_c,v_throat_c,P_throat_c,M_throat_c,m_throat_c]= conservative_form(n,x,dx,gamma,throat,nt,dt,A);

Time_conservative=toc;
fprintf(\'Execution time for Conservative form is %0.3g Sec. \\n \',Time_conservative)

figure(3)
subplot(2,3,1)
plot(x,rho,\'r\',x,rho_c,\'b*\')
xlabel(\'Length of nozzle(m)\')
ylabel(\'Non-Dimensional density\')
title(\'Density Profile\')
grid minor
legend(\'Non-Conservative\',\'Conservative\')
subplot(2,3,2)
plot(x,T,\'r\',x,T_c,\'b*\')
xlabel(\'Length of nozzle(m)\')
ylabel(\'Non-Dimensional Temperature\')
title(\'Temperature Profile\')
grid minor
legend(\'Non-Conservative\',\'Conservative\')
subplot(2,3,3)
plot(x,v,\'r\',x,v_c,\'b*\')
xlabel(\'Length of nozzle(m)\')
ylabel(\'Non-Dimensional Velocity\')
title(\'Velocity Profile\')
grid minor
legend(\'Non-Conservative\',\'Conservative\')
subplot(2,3,4)
plot(x,P,\'r\',x,P_c,\'b*\')
xlabel(\'Length of nozzle(m)\')
ylabel(\'Non-Dimensional Pressure\')
title(\'Pressure Profile\')
grid minor
legend(\'Non-Conservative\',\'Conservative\')
subplot(2,3,5)
plot(x,M,\'r\',x,M_c,\'b*\')
xlabel(\'Length of nozzle(m)\')
ylabel(\'Mach number\')
title(\'Mach number Profile\')
grid minor
legend(\'Non-Conservative\',\'Conservative\')
subplot(2,3,6)
plot(x,m,\'r\',x,m_c,\'b*\')
xlabel(\'Length of nozzle(m)\')
ylabel(\'Non-Dimensional Mass Flow Rate\')
title(\'Mass Flow Rate Profile\')
hold on
line([0 3],[0.579 0.579], \'color\', \'y\')
grid minor
legend(\'Non-Conservative\',\'Conservative\')




figure(4);
subplot(2,3,1)

plot(1:nt,rho_throat,\'r\',1:nt,rho_throat_c,\'b\')
hold on
line([0 nt],[0.634 0.634], \'color\', \'g\')
xlabel(\'Time(s)\')
ylabel(\'Non-Dimensional Density\')
title(\'Variation of Density at the Nozzle Throat\')
grid minor
legend(\'Non-Conservative\',\'Conservative\',\'Exact Solution\')
subplot(2,3,2)
plot(1:nt,T_throat,\'r\',1:nt,T_throat_c,\'b\')
hold on
line([0 nt],[0.833 0.833], \'color\', \'g\')
xlabel(\'Time(s)\')
ylabel(\'Non-Dimensional Temperature\')
title(\'Variation of Temperature at the Nozzle Throat\')
grid minor
legend(\'Non-Conservative\',\'Conservative\',\'Exact Solution\')
subplot(2,3,3)
plot(1:nt,v_throat,\'r\',1:nt,v_throat_c,\'b\')
xlabel(\'Time(s)\')
ylabel(\'Non-Dimensional Velocity\')
title(\'Variation of Velocity at the Nozzle Throat\')
grid minor
legend(\'Non-Conservative\',\'Conservative\')
subplot(2,3,4)
plot(1:nt,P_throat,\'r\',1:nt,P_throat_c,\'b\')
hold on
line([0 nt],[0.528 0.528], \'color\', \'g\')
xlabel(\'Time(s)\')
ylabel(\'Non-Dimensional Pressure\')
title(\'Variation of Pressure at the Nozzle Throat\')
grid minor
legend(\'Non-Conservative\',\'Conservative\',\'Exact Solution\')
subplot(2,3,5)
plot(1:nt,M_throat,\'r\',1:nt,M_throat_c,\'b\')
hold on
line([0 nt],[1 1], \'color\', \'g\')
xlabel(\'Time(s)\')
ylabel(\' Mach number\')
title(\'Variation of Mach number at the Nozzle Throat\')
grid minor
legend(\'Non-Conservative\',\'Conservative\',\'Exact Solution\')
subplot(2,3,6)
plot(1:nt,m_throat,\'r\',1:nt,m_throat_c,\'b\')
hold on 
line([0 nt],[0.579 0.579], \'color\', \'g\')
xlabel(\'Time(s)\')
ylabel(\'Non-Dimensional Mass Flow Rate\')
title(\'Variation of Mass Flow Rate at the Nozzle Throat\')
grid minor
legend(\'Non-Conservative\',\'Conservative\',\'Exact Solution\')

NON CONSERVATIVE FORM

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

%Calculating initial profiles

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

dt=min(0.5.*dx./(sqrt(T) + v));

%Outer time loop
 for k=1:nt
     dt=min(0.5.*dx./(sqrt(T) + v));
     rho_ini=rho;
     T_ini=T;
     v_ini=v;
     
     % Time-step:
     dt = 0.01;
     
     % Predictor step 
     for i=2:n-1
     
     dv_dx=(v(i+1)-v(i))/dx;
     drho_dx=(rho(i+1)-rho(i))/dx;
     dlogA_dx=(log(A(i+1))-log(A(i)))/dx;
     dT_dx=(T(i+1)-T(i))/dx;
     
     % continuity equation
     drho_dt_p(i)=-rho(i)*dv_dx-rho(i)*v(i)*dlogA_dx-v(i)*drho_dx;
     
     % momentum equation
     dv_dt_p(i)=-v(i)*dv_dx-(1/gamma)*(dT_dx+((T(i)/rho(i))*drho_dx));
     
     % energy equation
     dT_dt_p(i)=-v(i)*dT_dx-(gamma-1)*T(i)*(dv_dx+(v(i)*dlogA_dx));
     
     %solution update
     rho(i)=rho(i)+drho_dt_p(i)*dt;
     v(i)=v(i)+dv_dt_p(i)*dt;
     T(i)=T(i)+dT_dt_p(i)*dt;
     
     
     end
      
     % Corrector step      
     for j=2:n-1
     
     dv_dx_c=(v(j)-v(j-1))/dx;
     drho_dx_c=(rho(j)-rho(j-1))/dx;
     dlnA_dx_c=(log(A(j))-log(A(j-1)))/dx;
     dT_dx_c=(T(j)-T(j-1))/dx;
     
     % continuity equation     
     drho_dt_c(j)=-rho(j)*dv_dx_c-rho(j)*v(j)*dlnA_dx_c-v(j)*drho_dx_c;
     
     % momentum equation     
     dv_dt_c(j)=-v(j)*dv_dx_c-(1/gamma)*(dT_dx_c+((T(j)/rho(j))*drho_dx_c));
     
     % energy equation
     dT_dt_c(j)=-v(j)*dT_dx_c-(gamma-1)*T(j)*(dv_dx_c+(v(j)*dlnA_dx_c));
     
     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 l=2:n-1
     rho(l)=rho_ini(l)+drho_dt(l)*dt;
     v(l)=v_ini(l)+dv_dt(l)*dt;
     T(l)=T_ini(l)+dT_dt(l)*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);


     %Computing Pressure(P), Mach number(M) and mass flow rate(m)
     P=rho.*T;
     M=(v./(T).^0.5);
     m = rho.*v.*A;
     
     % Capturing the values at throat for each iterations
     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);
     
     % Ploting mass flow rate at different times
    a=figure(1);    
    hold on;
    if k ==1
        plot(x,m,\'b\',\'linewidth\',2);
     else if k == 50
            plot(x,m,\'g\',\'linewidth\',2);
             else if k == 100
                plot(x,m,\'magenta\',\'linewidth\',2);
        else if k == 250
                plot(x,m,\'r\',\'linewidth\',2);
            else if k == 500
                    plot(x,m,\'m\',\'linewidth\',2);
                else if k == 750
        plot(x,m,\'-.\',\'linewidth\',2);
         else if k == 1000
        plot(x,m,\'y\',\'linewidth\',2);
        
                        legend(\'1th time step\',\'50th time step\',\'100th time step\',\'250th time step\',\'500th time step\', \'750th time step\',\'1000th time step\');
                        title(\'Mass flow rate at different time step(Non-Conservative form)\', \'FontSize\',12);
                        xlabel(\'Non-Dimensional Length of Nozzle (x/L)\', \'FontSize\', 11 );
                        ylabel(\'Non -Dimensional Mass flow rate\',  \'FontSize\', 11) 
                        grid minor
                       
    end
    end
    end
    end
    end
    end 
 end
end

CONSERVATIVE FORM

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


%Calculating initial profiles
for i=1:n
    if x(i)>=0 && x(i)<0.5
    rho(i)=1;
    T(i)=1;
        else if 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));

            else if x(i)>=1.5 && x(i)<3.5
            rho(i)=0.634-(0.3879*(x(i)-1.5));
            T(i)=0.833-(0.3507*(x(i)-1.5));

    end
    end
    end
    end


% Defining velocity and pressure    
v=0.59./(rho.*A); 
P=rho.*T;
 
% Defining solution vectors U1, U2 and U3
U1=rho.*A;
U2=rho.*A.*v;
U3 = rho.*A.*((T./(gamma-1))+((gamma/2).*(v.^2)));

%Outer time loop
 for k=1:nt
     
     % Copies of solution vectors
     U1_old=U1;
     U2_old=U2;
     U3_old=U3;
    
     % Adaptable time step 
     dt=min(0.5.*dx./(sqrt(T) + v));
     
     % Pure form of flux terms
     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 step 
     %forward differences
     for i=2:n-1
     
     dF1_dx=(F1(i+1)-F1(i))/dx;
     dF2_dx=(F2(i+1)-F2(i))/dx;
     dF3_dx=(F3(i+1)-F3(i))/dx;
     
     % continuity equation
     dU1_dt_p(i)=-dF1_dx;
     
     % momentum equation
     
     J2(i) = (1/gamma).*rho(i).*T(i).*((A(i+1)-A(i))/dx);
     dU2_dt_p(i)=-dF2_dx+J2(i);
     
     % energy equation
     dU3_dt_p(i)=-dF3_dx;
     
     %solution update
     U1(i)=U1(i)+dU1_dt_p(i)*dt;
     U2(i)=U2(i)+dU2_dt_p(i)*dt;
     U3(i)=U3(i)+dU3_dt_p(i)*dt;
     
     
     end
     
     %updating the values
% non dimensional pressure
  p = rho.*T ;
%velocity
    v = U2./U1; 
 % density
   rho = U1./A;
   % temperature
   T = (gamma-1).*((U3./U1)-((gamma./2).*(v.^2)));
    
    %Update the values of F1,F2,F3 after predictor step  
     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 step 
     %backward differences
     
     for j=2:n-1
         % source term
     dF1_dx=(F1(j)-F1(j-1))/dx;
     dF2_dx=(F2(j)-F2(j-1))/dx;
     dF3_dx=(F3(j)-F3(j-1))/dx;
     
     % continuity equation
     dU1_dt_c(j)=-dF1_dx;
     
     % momentum equation
     
     J2(j) = (1/gamma).*rho(j).*T(j).*((A(j)-A(j-1))/dx);
     dU2_dt_c(j)=-dF2_dx+J2(j);
     
     % energy equation
     dU3_dt_c(j)=-dF3_dx;
     
     end
     
     % compute the average time derivative
     dU1_dt_avg=0.5.*(dU1_dt_p+dU1_dt_c);
     dU2_dt_avg=0.5.*(dU2_dt_p+dU2_dt_c);
     dU3_dt_avg=0.5.*(dU3_dt_p+dU3_dt_c);
     
     % Final solution update
     
     for l=2:n-1
     U1(l)=U1_old(l)+dU1_dt_avg(l)*dt;
     U2(l)=U2_old(l)+dU2_dt_avg(l)*dt;
     U3(l)=U3_old(l)+dU3_dt_avg(l)*dt;

     
     end

     % Apply boundary conditions 
     % Inlet
     U2(1) = 2*U2(2) - U2(3);
     % Outlet
     U1(n) = 2*U1(n-1)-U1(n-2);
     U2(n) = 2*U1(n-1)-U1(n-2);
     U3(n) = 2*U1(n-1)-U1(n-2);
 
     
     % Updation of flow field variables
     rho= U1./A; 
     v = U2./U1;
     T = (gamma-1).*((U3./U1)-((gamma*v.^2)/2));
     P=rho.*T;
     
     % Calculating Mach  number(M) and mass flow rate(m)
     M= v./((T).^0.5);
     m = rho.*v.*A;
     
     % Capturing the values at throat for each iterations
     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);
     
    % Ploting mass flow rate at different times
    b=figure(2);    
    hold on;
    if k ==1
        plot(x,m,\'b\',\'linewidth\',2);
    else if k == 50
            plot(x,m,\'g\',\'linewidth\',2);
        else if k == 100
                plot(x,m,\'r\',\'linewidth\',2);
            else if k == 250
                    plot(x,m,\'m\',\'linewidth\',2);
                else if k == 500
                        plot(x,m,\'c\',\'linewidth\',2);
                 else if k == 750
                    plot(x,m,\'*\',\'linewidth\',2);
                else if k == 1000
                        plot(x,m,\'y\',\'linewidth\',2);
                        
                        legend(\'1th time step\',\'50th time step \',\'100th time step\',\'250th time step \',\'500th time step\',\'750th time step\',\'1000th time step\');
                        title(\'Mass flow rate at different time step(Conservative form)\', \'FontSize\',12);
                        xlabel(\'Non-Dimensional Length of Nozzle (x/L)\', \'FontSize\', 11 );
                        ylabel(\'Non -Dimensional Mass flow rate\', \'FontSize\', 11) 
                        grid minor
          
    end
    end
    end
    end
    end
    end
    end
 end
end

Graph:

Grid ( n = 31)

Grid ( n = 61)

n = 61

Execution time for Non-Conservative form is 6.57 Sec.
Execution time for Conservative form is 5.61 Sec.

n = 31

Execution time for Non-Conservative form is 5.91 Sec.
Execution time for Conservative form is 4.25 Sec.

RESULTS:

From the above graph, It is observed that variation of mass flow rate at conservative show slight variation compared to the non-conservative from. Also, we can easily find out that conservative form converges quicker ( approximately 600 iterations ) compared to non-conservative form. The minimum number of cycles required for convergence for the non-conservative form is around 1000 cycles.

The accuracy of the solution improves as we increase the grid points. Also, the computational time required for the conservative form is much lesser than the non-conservative form. Moreover, the grid points produce exact value closer to analytical value.

The conservative form preserves more mass compared to non-conservative form.