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

OBJECTIVE:

To simulate the isentropic flow through a Quasi 1D subsonic-supersonic nozzle by deriving both the Conservation & Non-Conservation forms of the governing equations and solve them by applying MacCormack's technique using MATLAB.

Problem statement

Determine the steady-state temperature distribution for the flow-field variables and investigate the difference between the two forms of governing equations by comparing their solutions along with the plots stated below.

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 form

THEORY:

Any numerical solution of isentropic flow through a quasi 1D subsonic-supersonic nozzle is unreasonable, while an analytic solution exists for these kinds of problems. However, we are trying here to establish the broadened capability of CFD in various streams rather in representing the approximated form of the analytic solution. Here we implement the numerical solution with Conservative and Non-Conservative forms.

Governing Equations are marched in terms of PDE for both forms in order to use Time marching solution for the problem. These forms represent in numerical aspect only, if we take analytic formation all the methods are the same.

Non Conservative Form:

Conservative Form:

Non-Dimensional Quantities:

Also as part of the problem setup, we represent quantities in terms of non-dimensional quantities.

Nozzle diagram with sections:

Since the problem setup in 1D, we consider only distribution in the x-axis. Here Nozzle can be split into two sections

• Convergent Section
• Divergent Section

MacCormack method:

To discretize the numerical solution we use the technique called the MacCormack method. The MacCormack method is a variation of the two-step Lax-Wendroff Scheme but is much simpler in the application.

The application of the MacCormack method proceeds in two steps; a predictor step which is followed by a corrector step.

Predictor step: In the predictor step, a "provisional" value of at time level n+1 (denoted by `u_i^{\bar{n+1}}` ) is estimated obtained by replacing the spatial and temporal derivatives in the previous first-order equation using Forward Difference Scheme. `\mathrm{\Delta t}` is used for Timestep.

Corrector step: In the corrector step, the predicted value `u_i^{\bar{n+1}}` is corrected according to the equation. Note that the corrector step uses the Backward Difference Scheme for spatial derivative. The time-step used in the corrector step is `\mathrm{\Delta t}/2` in contrast to the used in the predictor step.

MAIN CODE

clear all
close all
clc


%INPUT PARAMETERS
%No. of grid points    
n = 31;             

gamma = 1.4;
x = linspace(0,3,n);
dx = x(2) - x(1);

%No. of time step
nt = 1400;

%Courant Number
C = 0.5;            


%FUNCTION FOR NON CONSERVATIVE FORM
tic;
[Mass_flow_rate_NC, Pressure_NC, Mach_number_NC, rho_NC, v_NC, T_NC, rho_throat_NC, v_throat_NC, T_throat_NC, Mass_flow_rate_throat_NC, Pressure_throat_NC, Mach_number_throat_NC] = Non_conservative_form(x, dx, n, gamma, nt, C);
Time_elapsed_NC = toc;
fprintf('Time Elasped for Non Conservative form: %0.4f', Time_elapsed_NC);
hold on

%PLOTTING 
%Timestep variation of properties at throat area in Non Conservative form
figure(2)

%Density
subplot(4,1,1)
plot(linspace(1,nt,nt),rho_throat_NC,'color','m')
ylabel('Density')
legend({'Density at throat'});
axis([0 1400 0 1.3]);
grid on;
title('Timestep variation at throat area - NON CONSERVATIVE FORM')

%Temperature
subplot(4,1,2)
plot(linspace(1,nt,nt),T_throat_NC,'color','c')
ylabel('Temperature')
legend({'Temperature at throat'});
axis([0 1400 0.6 1]);
grid on;

%Pressure
subplot(4,1,3)
plot(linspace(1,nt,nt),Pressure_throat_NC,'color','g')
ylabel('Pressure')
legend({'Pressure at throat'});
axis([0 1400 0.4 1.1]);
grid on;

%Mach Number
subplot(4,1,4)
plot(linspace(1,nt,nt),Mach_number_throat_NC,'color','y')
xlabel('Timesteps')
ylabel('Mach Number')
legend({'Mach Number at throat'});
axis([0 1400 0.6 1.4]);
grid on;


%Steady state simulation for primitve values for Non Conservative form
figure(3);
%Density

subplot(4,1,1);
plot(x,rho_NC,'color','m')
ylabel('Density')
legend({'Density'});
axis([0 3 0 1])
grid on;
title('Variation of Flow Rate Distribution - NON CONSERVATIVE')

%Temperature
subplot(4,1,2);
plot(x,T_NC,'color','c')
ylabel('Temperature')
legend({'Temperature'});
axis([0 3 0 1])
grid on;

%Pressure
subplot(4,1,3);
plot(x,Pressure_NC,'color','g')
ylabel('Pressure')
legend({'Pressure'});
axis([0 3 0 1])
grid on;

%Mach number
subplot(4,1,4);
plot(x,Mach_number_NC,'color','y')
xlabel('Nozzle Length')
ylabel('Mach number')
legend({'Mach number'});
axis([0 3 0 4])
grid on;


%FUNCTION FOR CONSERVATIVE FORM
tic;
[Mass_flow_rate_C, Pressure_C, Mach_number_C, rho_C, v_C, T_C, rho_throat_C, v_throat_C, T_throat_C, Mass_flow_rate_throat_C, Pressure_throat_C, Mach_number_throat_C] = Conservative_form(x, dx, n, gamma, nt, C);
Time_elasped_C = toc;
fprintf('Time Elasped for Conservative form: %0.4f', Time_elasped_C);

%PLOTTING 
%Timestep variation of properties at throat area in Conservative form
figure(5)

%Density
subplot(4,1,1)
plot(rho_throat_C,'color','m')
ylabel('Density')
legend({'Density at throat'});
axis([0 1400 0 1]);
grid on;
title('Timestep variation at throat area - CONSERVATIVE FORM')

%Temperature
subplot(4,1,2)
plot(T_throat_C,'color','c')
ylabel('Temperature')
legend({'Temperature at throat'});
axis([0 1400 0.6 1]);
grid on;

%Pressure
subplot(4,1,3)
plot(Pressure_throat_C,'color','g')
ylabel('Pressure')
legend({'Pressure at throat'});
axis([0 1400 0.2 0.8]);
grid on;

%Mach Number
subplot(4,1,4)
plot(Mach_number_throat_C,'color','y')
xlabel('Timesteps')
ylabel('Mach Number')
legend({'Mach Number at throat'});
axis([0 1400 0.6 1.4]);
grid on;


%Steady state simulation for primitve values for Non Conservative form
figure(6);

%Density 
subplot(4,1,1);
plot(x,rho_C,'color','m')
ylabel('Density')
legend({'Density'});
axis([0 3 0 1])
grid on;
title('Variation of Flow Rate Distribution - CONSERVATIVE FORM')

%Temperature 
subplot(4,1,2);
plot(x,T_C,'color','c')
ylabel('Temperature')
legend({'Temperature'});
axis([0 3 0 1])
grid on;

%Pressure
subplot(4,1,3);
plot(x,Pressure_C,'color','g')
ylabel('Pressure')
legend({'Pressure'});
axis([0 3 0 1])
grid on;

%Mach number
subplot(4,1,4);
plot(x,Mach_number_C,'color','y')
xlabel('Nozzle Length')
ylabel('Mach number')
legend({'Mach number'});
axis([0 3 0 4])
grid on;


%COMPARISON PLOT of Normalized Mass Flow Rate of both forms
figure(7)
hold on
plot(x,Mass_flow_rate_NC,'r','linewidth',1.25)
hold on;
plot(x,Mass_flow_rate_C,'b','linewidth',1.25)
hold on;
grid on

legend('Non Conservative','Conservative');
xlabel('Length of the Domain')
ylabel('Mass Flow Rate')
title('Comparison of Normalized Mass Flow Rate of both forms')

NON CONSERVATIVE FORM

%Function for non conservative form
function[Mass_flow_rate_NC,Pressure_NC,Mach_number_NC,rho_NC,v_NC,T_NC,rho_throat_NC,v_throat_NC,T_throat_NC,Mass_flow_rate_throat_NC,Pressure_throat_NC,Mach_number_throat_NC]=Non_conservative_form(x,dx,n,gamma,nt,C)
%Initial Boundary conditions
rho_NC=1-0.3146*x;
T_NC=1-0.2341*x;
v_NC=(0.1+1.09*x).*T_NC.^0.5;
%Area & throat profiles
a=1+2.2*(x-1.5).^2;
throat=find(a==1);
%outer time loop
for k=1:nt
    %time step control
    dt=min(C.*dx./(sqrt(T_NC)+v_NC));
     rho_old=rho_NC;
     T_old=T_NC;
     v_old=v_NC;
     
     % Predictor Method
     for j=2:n-1
         %simpliying the equation;s spatial terms
         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)*drho_dx/rho_NC(j)));
          % 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 Method
     for j=2:n-1
         %simpliying the equation;s spatial terms
         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)*drho_dx/rho_NC(j)));
          % 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 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_NC(i)=rho_old(i)+drho_dt(i)*dt;
               v_NC(i)=v_old(i)+dv_dt(i)*dt;
               T_NC(i)=T_old(i)+dT_dt(i)*dt;
           end
           %% Applying Boundary conditions
           % At Inlet
           rho_NC(1)=1;
           v_NC(1)=2*v_NC(2)-v_NC(3);
           T_NC(1)=1;
           % At Outlet
           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);
           %Defining mass flow rate , Mach no and pressure
           Mass_flow_rate_NC=rho_NC.*a.*v_NC;
           Pressure_NC=rho_NC.*T_NC;
           Mach_number_NC=(v_NC./sqrt(T_NC));
           % Calculating variables at the throat
           rho_throat_NC(k)=rho_NC(throat);
           T_throat_NC(k)=T_NC(throat);
           v_throat_NC(k)=v_NC(throat);
           Mass_flow_rate_throat_NC(k)=Mass_flow_rate_NC(throat);
           Pressure_throat_NC(k)=Pressure_NC(throat);
           Mach_number_throat_NC(k)=Mach_number_NC(throat);
           %%Plotting
           %Comparison of mass flow rate at different time steps
           figure(1);
           if k==50
               plot(x,Mass_flow_rate_NC,'m','linewidth',1.25);
               hold on
           elseif k==100
               plot(x,Mass_flow_rate_NC,'c','linewidth',1.25);
               hold on
           elseif k==200
               plot(x,Mass_flow_rate_NC,'g','linewidth',1.25);
               hold on
           elseif k==400
               plot(x,Mass_flow_rate_NC,'y','linewidth',1.25);
               hold on
           elseif k==800
               plot(x,Mass_flow_rate_NC,'k','linewidth',1.25);
               hold on
           end
           %labelling
           title('Mass flow rate at different time steps-(NON CONSERVATIVE FORM)')
           xlabel('Distance')
           ylabel('Mass flow rate')
           legend({'50^t^h Timesteps','100^t^h Timesteps','200^t^h Timesteps','400^t^h Timesteps','800^t^h Timesteps'})
           axis([0 3 0 2])
           grid on
end
end

CONSERVATIVE FORM

%Function for Conservative form

function [Mass_flow_rate_C, Pressure_C, Mach_number_C,  rho_C, v_C, T_C, rho_throat_C, v_throat_C, T_throat_C, Mass_flow_rate_throat_C, Pressure_throat_C, Mach_number_throat_C] = Conservative_form(x, dx, n, gamma, nt, C)

%Determining Initial profile of Density and Temperature
  for i = 1:n
    
    if (x(i) >= 0 && x(i) <= 0.5)
        rho_C(i) = 1;
        T_C(i) = 1;
    elseif (x(i) >= 0.5 && x(i) <= 1.5)
        rho_C(i) = 1 - 0.366*(x(i) - 0.5);
        T_C(i) = 1 - 0.167*(x(i) - 0.5);
    elseif (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

%Area & Throat profiles
a = 1 + 2.2*(x - 1.5).^2;
throat = find(a==1);

%Initial profile fot Velocity and Pressure
v_C = 0.59./(rho_C.*a);
Pressure_C = rho_C.*T_C;

%Calculate the Initial Solution Vectors(U)
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-step control
    dt = min(C.*dx./(sqrt(T_C) + v_C));
    
    %Saving a copy of Old Values
    U1_old = U1;
    U2_old = U2;
    U3_old = U3;

    %Calculating the Flux vector(F)
    F1 = U2;
    F2 = ((U2.^2)./U1) + ((gamma - 1)/gamma)*(U3 - (gamma/2)*((U2.^2)./U1));
    F3 = ((gamma*U2.*U3)./U1) - (gamma*(gamma - 1)*((U2.^3)./(2*U1.^2)));
    
    
    %Predictor Method
    for j = 2:n-1
    
        %Calculating the Source term(J)
        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);
        
        
        %Updating New values
        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 Flux vectors
    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 j = 2:n-1
        
        %Updating the Source term(J)
        J2(j) = (1/gamma)*rho_C(j)*T_C(j)*((a(j) - a(j - 1))/dx);

        %Simplifying the Equation's spatial terms
        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

    %Compute 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);

    %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
    
    
    %Applying Boundary Conditions
    %At Inlet
    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));
    
    %At 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);
    
    %Final Primitive value update
    rho_C = U1./a;
    v_C = U2./U1;
    T_C = (gamma - 1)*((U3./U1 - (gamma/2)*(v_C).^2));
    
    %Defining Mass flow rate, Pressure and Mach number
    Mass_flow_rate_C = rho_C.*a.*v_C;
    Pressure_C = rho_C.*T_C;
    Mach_number_C = (v_C./sqrt(T_C));
    
    %Calculating Variables at throat section
    rho_throat_C(k) = rho_C(throat);
    T_throat_C(k) = T_C(throat);
    v_throat_C(k) = v_C(throat);
    Mass_flow_rate_throat_C(k) = Mass_flow_rate_C(throat);
    Pressure_throat_C(k) = Pressure_C(throat);
    Mach_number_throat_C(k) = Mach_number_C(throat);

    %PLOTTING
    %Comparison of non-dimensional Mass flow rates at different time steps
    figure(4)
    if k == 50
       plot(x, Mass_flow_rate_C, 'm', 'linewidth', 1.25);
       hold on
    elseif k == 100
       plot(x, Mass_flow_rate_C, 'c', 'linewidth', 1.25);
       hold on
    elseif k == 200
       plot(x, Mass_flow_rate_C, 'g', 'linewidth', 1.25);
       hold on
    elseif k == 400
       plot(x, Mass_flow_rate_C, 'y', 'linewidth', 1.25);
       hold on
    elseif k == 800
       plot(x, Mass_flow_rate_C, 'k', 'linewidth', 1.25);
       hold on
    end
           
    
    %Labelling
    title('Mass flow rates at different Timesteps - CONSERVATIVE FORM')
    xlabel('Distance')
    ylabel('Mass flow rate')
    legend({'50^t^h Timestep','100^t^h Timestep','200^t^h Timestep','400^t^h Timestep','800^t^h Timestep'})
    axis([0 3 0.4 1])
    grid on
    
end

end

RESULTS

Time step variation at throat for Non-Conservative and Conservative forms

On comparing the two graphs it can be visualized that the Non-conservative form consumes a higher amount of time than the Conservative form to reach a stable solution.

Flow rate distribution at throat for Non-Conservative and Conservative forms

It can be predicted by both the Mass flow rate graph that initially it oscillates but becomes stable after 400-time steps for both the norms.

COMPARISON PLOT

It can be depicted from the above figure that the Non-conservative form oscillates more when compared to the conservative form.And the Conservative form attains stability in lesser no of iterations.

GRID INDEPENDENCE TEST

Types of Form	Grid Number	Variables
		Density	Pressure	Temperature	Mach Number
Non-Conservative Form	N=31	0.6387	0.5342	0.8365	0.9994
	N=61	0.6377	0.5329	0.8355	0.9998
Conservative Form	N=31	0.6516	0.5478	0.8406	0.9794
	N=61	0.6383	0.5331	0.8352	0.9993

From the table, it can be depicted that increasing the grid points leads to more accurate results at expense of time.

CONCLUSION

On cumulative comparing the results, we can say that Conservative form converges faster than the Non-conservative form offering the stable solution.