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

AIM:- To determine the steady-state temperature distribution for the flow-field variables of 1-D super-sonic nozzle flow and investigate the difference between the two forms of governing equations by comparing their solutions by using the MacCormack's method.

Governing Equations:- To complete this simulation, the basic concepts of the conservative and non-conservative forms of the governing equations to be used must be acknowledged. The governing equations can be expressed either in conservative or non-conservative forms which also may be expressed in integral or differential forms.

Though all the equations may look very different from each other, they are actually the same. Just by doing some mathematical manipulations, the equations can be made to look the same. That is, the conservative integral form can be converted to conservative differential form using some mathematical manipulations and which can also be converted to non-conservative differential form. This means that all the equations are the same. The below figure shows the different forms of the continuity equation. But, it is to be noted that the conservative and non-conservative form will differ when considered under numerical calculations.

For the 1-D nozzle flow, the continuity, momentum and energy equations can be simplified as,

1. Conservative Form:

Continuity: `del/(delt)(rhoA) + del/(delx)(rhoAV) = 0`

Momentum: `del/(delt)(rhoAV) + del/(delx)(rhoAV^2) = -(delP)/(delx)`

Energy: `del/(delt)(rhoAV^2/2) + del/(delx)(rhoAV^3/2) = -AV(delP)/(delx)`

2. Non-Conservative Form:

Continuity: `del/(delt)(rhoA) + Vdel/(delx)(rhoA) = -rhoAdel/(delx)V`

Momentum: `rho(delV)/(delt) + PV(delV)/(delx) = -(delP)/(delx)`

Energy: `rho(dele)/(delt) + rhoV(dele)/(delx) = -P(delV)/(delx) -PVdel/(delx) (lnA)`

These above equations are to be converted into non-dimensional form to normalize these equations so that variables do not have units. The non-dimensional form of the governing equations can be written as,

1. Conservative Form:

Continuity: `del/(delt')(rho'A') + del/(delx')(rho'A'V') = 0`

Momentum: `del/(delt')(rho'A'V') + del/(delx')(rho'A'V'^2+1/gammaP'A') = -1/gammaP'(delA')/(delx)`

Energy: `del/(delt')(rho'((e')/(gamma-1)+gamma/2V'^2)A') + del/(delx')(rho'((e')/(gamma-1)+gamma/2V'^2)V'A'+P'A'V') =0`

The governing equations in the conservative can be expressed in simple form using the solution vector U, flux vector F and source term J as,

With the help of the above equations, the governing equations are simplified as,

Continuity:`(delU1)/(delt') = -(delF1)/(delx')`

Momentum:`(delU2)/(delt') = -(delF2)/(delx')+J2`

Energy:`(delU3)/(delt') = -(delF3)/(delx')`

2. Non-Conservative Form:

Continuity: `(delrho')/(delt')= -rho'(delV')/(delx')-rho'V'(del(lnA'))/(delx')-V'(delrho')/(delx')`

Momentum: `(delV')/(delt')= -V'(delV')/(delx')-1/gamma((delT')/(delx')+(T')/(rho')(delrho')/(delx'))`

Energy: `(delT')/(delt')= -V'(delT')/(delx')-(gamma-1)T'((delV')/(delx')+V'(dellnA')/(delx'))`

Where,

Dimensionless temperature`T'=T/(To)` , T is Flow temperature and To is reservoir temperature

Dimensionless density`rho'=rho/(rhoo)`, `rho`is flow density and `rhoo`is reservoir density

Dimensionless length`x'=x/(L)`, L is length of nozzle

Dimensionless velocity `V'=V/(ao)`, V is flow velocity and ao is speed of sound

Dimensionless Area `A'=A/(Ao)`, A is area of nozzle and Ao is sonic throat area

Dimensionless internal energy `e'=e/(eo)`, `e = (Cv)/T` and `eo = (Cv)/(To)`

Nozzle Shape:-

Timestep:- `Deltat = min(C*(Deltax)/(sqrt(T')+V'))`

Boundary Conditions and Initial Conditions:-

1. Conservative:

Subsonic Inlet:

Supersonic Outlet:

2. Non-Conservative:

Subsonic Inlet:

Supersonic Outlet:

The above equations are used in the calulcations variables in the nozzle flow using the Mac-Cormack's method by setting up the finite-difference expressions.

Mac-Cormack Method:-

MacCormack's method is an explicit finite difference technique of second-order-accurate in both space and time. It is a variant of the Lax-Wendroff approach but is much simpler compared to it. Briefly, it uses a predictor and corrector step. In the predictor step, the space terms are replaced with forward differencing and in the corrector step, the space terms are replaced with backward differencing. The solution is updated and average time derivatives of predictor and corrector step are calculated for every grid point and timestep.

E.g;

Predictor Step:

Continuity:

Momentum:

Energy:

Main Code:

clear all
close all
clc

%inputs
n = 31;
x = linspace(0,3,n);
dx = x(2) - x(1);
gamma = 1.4;
%area of nozzle
a = 1+2.2*(x-1.5).^2;
%area at throat
throat = find(a==1);
%Time steps
nt = 1400;
%courant number
C = 0.5;
tic
%Non-Conservative form
[rho_nc,T_nc,v_nc,p_nc,mfr_nc,M_nc,dt,rho_nct,T_nct,p_nct,M_nct,k] = non_conservative(n,x,dx,gamma,nt,C,a,throat);
time_nc = toc;

tic
%Conservative form
[rho_c,T_c,v_c,p_c,mfr_c,M_c,rho_ct,T_ct,p_ct,M_ct] = conservative(n,x,dx,gamma,nt,C,a,throat);
time_c = toc;

%Steady State    
    figure(1)
    subplot(4,1,1)
    plot(x,rho_nc)
    hold on
    plot(x,rho_c,'.')
    grid minor
    title('Steady-state distribution of primitive variables inside the nozzle')
    xlabel('Non-Dimensional Length')
    ylabel('Density')
    legend('Non-conservative','Conservative')
    subplot(4,1,2)
    plot(x,M_nc)
    hold on
    plot(x,M_c,'.')
    grid minor
    xlabel('Non-Dimensional Length')
    ylabel('Mach Number')
    legend('Non-conservative','Conservative')
    subplot(4,1,3)
    plot(x,T_nc)
    hold on
    plot(x,T_c,'.')
    grid minor
    xlabel('Non-Dimensional Length')
    ylabel('Temperature')
    legend('Non-conservative','Conservative')
    subplot(4,1,4)
    plot(x,p_nc)
    hold on
    plot(x,p_c,'.')
    grid minor
    xlabel('Non-Dimensional Length')
    ylabel('Pressure')
    legend('Non-conservative','Conservative')

%Time wise variations in variables at throat section of nozzle
figure(2)
    subplot(4,1,1)
    plot(linspace(1,nt,nt),rho_nct)
    hold on
    plot(linspace(1,nt,nt),rho_ct)
    grid minor
    title('Time-wise variation of the primitive variables at Throat')
    xlabel('Timestep')
    ylabel('Density')
    legend('Non-conservative','Conservative')
    subplot(4,1,2)
    plot(linspace(1,nt,nt),p_nct)
    hold on
    plot(linspace(1,nt,nt),p_ct)
    grid minor
    xlabel('Timestep')
    ylabel('Pressure')
    legend('Non-conservative','Conservative')
    subplot(4,1,3)
    plot(linspace(1,nt,nt),T_nct)
    hold on
    plot(linspace(1,nt,nt),T_ct)
    grid minor
    xlabel('Timestep')
    ylabel('Temperature')
    legend('Non-conservative','Conservative')
    subplot(4,1,4)
    plot(linspace(1,nt,nt),M_nct)
    hold on
    plot(linspace(1,nt,nt),M_ct)
    grid minor
    xlabel('Timestep')
    ylabel('Mach Number')
    legend('Non-conservative','Conservative')
%Comparison of Normalized mass flow rate distributions
figure(3)
    plot(x,mfr_c)
    hold on
    plot(x,mfr_nc)
    grid minor
    xlabel('Non-Dimensional Length')
    ylabel('Normalized Mass Flow')
    legend('Non-Conservative','Conservative')
    title(' Comparison of Normalized mass flow rate distributions of both forms')

Function code for the non-conservative form:

function[rho_nc,T_nc,v_nc,p_nc,mfr_nc,M_nc,dt,rho_nct,T_nct,p_nct,M_nct] = non_conservative(n,x,dx,gamma,nt,C,a,throat)

%calculate initial profiles
rho_nc = 1-0.3146*x; %density
T_nc = 1-0.2314*x; %temperature
v_nc = (0.1 + 1.09*x).*T_nc.^0.5; %velocity

%outer time loop
for k = 1:nt
    
    rho_old_nc = rho_nc;
    v_old_nc = v_nc;
    T_old_nc = T_nc;
   
    %Time Step
    dt = min((C*dx)./(sqrt(T_nc) + v_nc));
 
   
    %PREDICTOR STEP
    for j = 2:n-1
        %drho_dt_p = -rho(j)*( v(j+1) - v(j) )/dx - rho(j)*v(j)*( log(a(j+1)) - log(a(j)) )/dx - v(j)*( rho(j+1) - rho(j) )/dx;
        dv_dx = (v_nc(j+1)-v_nc(j))/dx;
        drho_dx = (rho_nc(j+1) - rho_nc(j))/dx;
        dT_dx = (T_nc(j+1) - T_nc(j))/dx;
        dlog_a_dx = (log(a(j+1)) - log(a(j)))/dx;
        
        %Continuity equation
        drho_dt_p(j) = -rho_nc(j)*dv_dx - rho_nc(j)*v_nc(j)*dlog_a_dx - v_nc(j)*drho_dx;
        
        %Momentum Equation
        dv_dt_p(j) = -v_nc(j)*dv_dx - (1/gamma)*(dT_dx + (T_nc(j)/rho_nc(j))*drho_dx);
        
        %Energy Equation
        dT_dt_p(j) = -v_nc(j)*dT_dx - ((gamma-1)*T_nc(j)*(dv_dx + v_nc(j)*dlog_a_dx));
        
        %Solution Update
        rho_nc(j) = rho_nc(j) + drho_dt_p(j)*dt;
        v_nc(j) = v_nc(j) + dv_dt_p(j)*dt;
        T_nc(j) = T_nc(j) + dT_dt_p(j)*dt;    
    end
    
    %Corrector STEP
    for j = 2:n-1
        
        %drho_dt_p = -rho(j)*( v(j+1) - v(j) )/dx - rho(j)*v(j)*( log(a(j+1)) - log(a(j)) )/dx - v(j)*( rho(j+1) - rho(j) )/dx;
        dv_dx = (v_nc(j)-v_nc(j-1))/dx;
        drho_dx = (rho_nc(j) - rho_nc(j-1))/dx;
        dT_dx = (T_nc(j) - T_nc(j-1))/dx;
        dlog_a_dx = (log(a(j)) - log(a(j-1)))/dx;
        
        %Continuity equation
        drho_dt_c(j) = -rho_nc(j)*dv_dx - rho_nc(j)*v_nc(j)*dlog_a_dx - v_nc(j)*drho_dx;
        
        %Momentum Equation
        dv_dt_c(j) = -v_nc(j)*dv_dx - (1/gamma)*(dT_dx + (T_nc(j)/rho_nc(j))*drho_dx);
        
        %Energy Equation
        dT_dt_c(j) = -v_nc(j)*dT_dx - ((gamma-1)*T_nc(j)*(dv_dx + v_nc(j)*dlog_a_dx));
        
    end
    
    %Compute the avereage 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_nc(i) = rho_old_nc(i) + drho_dt(i)*dt;
        v_nc(i) = v_old_nc(i) + dv_dt(i)*dt;
        T_nc(i) = T_old_nc(i) + dT_dt(i)*dt; 
        M_nc(i) = v_nc(i)/(T_nc(i)^0.5);
    end
    
    %Apply Boundary Conditions
    %Subsonic Inflow
    rho_nc(1) = 1;
    v_nc(1) = 2*v_nc(2) - v_nc(3);
    T_nc(1) = 1;
    
    %Supersonic Outflow 
    rho_nc(n) = 2*rho_nc(n-1) - rho_nc(n-2);
    v_nc(n) = 2*v_nc(n-1) - v_nc(n-2);
    T_nc(n) = 2*T_nc(n-1) - T_nc(n-2);
    
    %Mach number
    M_nc = v_nc./sqrt(T_nc);
    %Pressure
    p_nc = rho_nc.*T_nc;
    %Mass flow rate
    mfr_nc = rho_nc.*v_nc.*a;
        
    %Variables at throat
    rho_nct(k) = rho_nc(throat);
    T_nct(k) = T_nc(throat);
    p_nct(k) = p_nc(throat);
    M_nct(k) = M_nc(throat);
    %mfr_nct(k) = mfr_nc(throat);
    
 %Mass Flow rate variation
 figure(5)
    if k==50
    plot(x,mfr_nc,'b')
    hold on
    elseif k==100
    plot(x,mfr_nc,'r')
    hold on
    elseif k==200
    plot(x,mfr_nc,'k')
    hold on
    elseif k==400
    plot(x,mfr_nc,'-o')
    hold on
    elseif k==800
    plot(x,mfr_nc,'m')
    hold on
    elseif k==1000
    plot(x,mfr_nc,'-*')
    hold on
    elseif k==1400
    plot(x,mfr_nc,'y')
    grid minor
    xlabel('Non-Dimensional length')
    ylabel('Non-Dimensional Mass Flow Rate')
    title('Variation of mass flow rate at different time steps in nonconservative form')
    legend('50th Timestep','100th Timestep','200th Timestep','400th Timestep','800th Timestep','1000th Timestep','1400th Timestep')                        
    end   
end
end

Function code for the conservative form:

function[rho_c,T_c,v_c,p_c,mfr_c,M_c,rho_ct,T_ct,p_ct,M_ct] = conservative(n,x,dx,gamma,nt,C,a,throat)

%calculate initial conditions
for m = 1:n
    if x(m)>=0 && x(m)<=0.5   
        rho_c(m) = 1;
        T_c(m) = 1;
    end
    if x(m)>=0.5 && x(m)<=1.5
        rho_c(m) = 1.0 - 0.366*(x(m) - 0.5);
        T_c(m) = 1.0 - 0.167*(x(m) - 0.5);
    end
    if x(m)>=1.5 && x(m)<=3.5
        rho_c(m) = 0.634 - 0.3879*(x(m) - 1.5);
        T_c(m) = 0.833 - 0.3507*(x(m) - 1.5);
    end
end
 v_c = (0.59)./(rho_c.*a);
 
%Solutions Vectors
U1 = rho_c.*a;
U2 = rho_c.*a.*v_c;
U3 = rho_c.*a.*((T_c/(gamma-1) + (gamma/2)*v_c.^2));


%outer time loop
for k = 1:nt
    
    %Time Steps
    dt = min((C*dx)./((T_c.^0.5) + v_c));
    
    %Savings old values
    
    U1_old = U1;
    U2_old = U2;
    U3_old = U3;
    
    %Flux vectors
    %Pure form
    
    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));    
    
    %PREDICTOR STEP
   
    for j = 2:n-1
        
       J2(j) = (1/gamma)*rho_c(j).*T_c(j).*((a(j+1) - a(j))/dx);
        
        %Continuity equation
        dU1_dt_p(j) = -(F1(j+1) - F1(j))/dx;
        
        %Momentum Equation
        dU2_dt_p(j) = -((F2(j+1) - F2(j))/dx) + J2(j);
        
        %Energy Equation
        dU3_dt_p(j) = -(F3(j+1) - F3(j))/dx;
        
        %Solution Update
        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

    %Obtaining primitive variables and flux values
   
    rho_c = U1./a;
    v_c = U2./U1;
    T_c = (gamma-1)*((U3./U1)-((gamma/2)*(v_c.^2)));
    p_c = rho_c.*T_c;
    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 STEP
    
    for j = 2:n-1
        
        J2(j) = (1/gamma)*rho_c(j).*T_c(j).*((a(j) - a(j-1))/dx);
        
        %Continuity equation
        dU1_dt_c(j) = -(F1(j) - F1(j-1))/dx;
        
        %Momentum Equation
        dU2_dt_c(j) = -((F2(j) - F2(j-1))/dx) + J2(j);
        
        %Energy Equation
        dU3_dt_c(j) = -(F3(j) - F3(j-1))/dx;
        
    end
    
    %Compute the avereage 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);
    
    %Final Solution Update
    
    for m = 2:n-1
        
        U1(m) = U1_old(m) +(dU1_dt(m)*dt);
        U2(m) = U2_old(m) +(dU2_dt(m)*dt);
        U3(m) = U3_old(m) +(dU3_dt(m)*dt);

    end
    
    %Apply Boundary Conditions
    %Subsonic Inflow
  
    U1(1) = rho_c(1)*a(1);
    U2(1) = 2*U2(2) - U2(3);
    v_c(1) = U2(1)/U1(1);
    U3(1) = U1(1)*(T_c(1)/(gamma-1) + (gamma/2)*v_c(1).^2);
    
    %Supersonic Outflow 
    
    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);
    
   %Calculation of primitive variables
   
   rho_c = U1./a;
   v_c = U2./U1;
   T_c = (gamma-1)*(U3./U1 - (gamma/2)*v_c.^2);
   p_c = rho_c.*T_c;
   
   %Mach number
   M_c = v_c./sqrt(T_c);
   %Mass flow rate
   mfr_c = rho_c.*v_c.*a;
   
   %Variables at throat
    rho_ct(k) = rho_c(throat);
    T_ct(k) = T_c(throat);
    p_ct(k) = p_c(throat);
    M_ct(k) = M_c(throat);
    

%Mass Flow rate variation
 figure(6)
    if k==50
    plot(x,mfr_c,'b')
    hold on
    elseif k==100
    plot(x,mfr_c,'r')
    hold on
    elseif k==200
    plot(x,mfr_c,'k')
    hold on
    elseif k==400
    plot(x,mfr_c,'-o')
    hold on
    elseif k==800
    plot(x,mfr_c,'m')
    hold on
    elseif k==1000
    plot(x,mfr_c,'-*')
    hold on
    elseif k==1400
    plot(x,mfr_c,'y')
    grid minor
    xlabel('Non-Dimensional length')
    ylabel('Non-Dimensional Mass Flow Rate')
    title('Variation of mass flow rate at different time steps in nonconservative form')
    legend('50th Timestep','100th Timestep','200th Timestep','400th Timestep','800th Timestep','1000th Timestep','1400th Timestep')                        
    end
         
end
end

OUTPUT:-

1. Steady-state distribution of primitive variables inside the nozzle

The steady-state distribution of primitive variables inside the nozzle are resulted as expected. As observed in the plot, the variables are almost constant in the initial state but when approached towards the throat of the nozzle, a disturbing change is observed in the variables resulting in a gradual decrease or increase. Irrespective of the form, the distribution of variables inside the nozzle is very much equal. Also, the Mach number as it should be is equal to 1 which can be observed from the plot.

2. Time-wise variation of the primitive variables

The distribution of variables inside the nozzle for both the forms may be very much equal, but a difference in the variation of the variables with respect to the timestep is noticed. From the plot, it can be observed that both the forms reach almost a steady-state from the 400th timestep and proceed linearly along with the simulation.

3. Variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process

Non-Conservative:

The plot above shows the variation of mass flow at different time steps in a non-conservative form of the governing equation. After the 400th timestep, the solution reaches steady state and negligible changes are found after it.Conservative:

The plot above shows the variation of mass flow at different time steps in a conservative form of the governing equation. After the 400th timestep, the solution reaches steady state and negligible changes are found after it.

4. Comparison of Normalized mass flow rate distributions of both forms

From the plot, it is easily identifiable that the mass flow rate in the conservative form varies unlikely when compared to non-conservative form along the length of the nozzle. But the mass flow rate in the conservative form at the throat of the nozzle is less compared to non-conservative form, even though having higher mass flow rates at the starting and end of the nozzle. Also, in the non-conservative form, the mass flow rate does not change much.

5.Comparison of computational time of both forms

Although both the forms mean the same theoretically, there is an identifiable difference when it comes to numerical calculations. Therefore, the computational time of both the forms is compared to know which one is faster.

From the bar graph, it can be said that a conservative form takes less time to compute compared to non-conservative form. Even though it is faster, it cannot be concluded that it is the best one as there are some differences in the variables along the timesteps.

6. Grid dependency test

When the number of grid points is increased from 31 to 91, that is refining the mesh, it is observed that there are negligible changes in the primitive values of the flow, which are shown in the table below.

Non-dimensional Flow variables	Non-Conservative		Conservative
	n=31	n=91	n=31	n=91
Density	0.6387	0.6406	0.6504	0.6380
Pressure	0.5342	0.5361	0.5463	0.5329
Mach number	0.9994	1.0019	0.9818	1.0003
Temperature	0.8365	0.8369	0.8400	0.8353

Hence, it can be concluded that the problem is independent of the grid points. No matter how refined it be, there will no changes in primitive values of the flow inside the nozzle.