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

Aim:

To simulate the isentropic flow through a quasi 1D subsonic-supersonic nozzle by deriving both the conservation and non-conservation forms of the governing equation and solve them by implementing the MacCormack’s technique using MATLAB.

Description:

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

Non-Conservative Form:

Continuity Equation:

(∂ρ^')/(∂t^' )=(ρ^' ∂V^')/(∂x^' )-ρ^' V^' (∂(lnA)/(∂x^' ))-V^' ((∂p^')/(∂x^' ))

Momentum Equation:

(∂V^')/(∂t^' )=V^' ((∂V^')/(∂x^' ))-1/γ ((∂T^')/(∂x^' )+(T^'/p^'  ∂p^')/(∂x^' ))

Energy Equation:

(∂T^')/(∂t^' )=V^'((∂T^')/(∂x^' ))-(γ-1) T^' [(∂V^')/(∂x^' )+(V^' ∂(lnA^' ))/(dx^' )]

Conservative Form:

Continuity Equation:

∂(ρ^' A^' V^' )/(∂t^' )+∂(ρ^' A^' V^' )/dx=0

Momentum Equation:

∂(ρ^' A^' V^' )/(∂t^' )+(ρ^' A^' V^'2+1/γ p^' A^')/(∂x^' )=1/γ p^' ((∂A^')/(∂x^' ))

Energy Equation:

∂[ρ^' (T/(γ-1)+γ/2 V^'2) A^']/(∂t^')+∂[p^'(T/(γ-1)+γ/2 V^'2) V^' A^'+p^' A^'V^' ]/(∂x^')=0

Non Dimensional Quantities:

Temperature:

T’ =T/T0

Density:

ρ^'=ρ/ρ_(0) 

Non Dimensional Distance:

x’ = x/L

Non Dimensional Velocity:

V’ = V/a0 

where ‘a₀’ is the local sonic velocity

a0 = √(γRT_0 )

Non Dimensional Time:

t’ = T/(L/a0) 

where L/a₀ is characteristic time scale

Non Dimensional Area:

A’ = A/A*

Nozzle Diagram with sections:

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

1. Convergent Section
2. Divergent Section

MacCarmack Method:

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

Predictor step:

In the predictor step, a “provisional” value of `u`at time level `n+1`(denoted by `u_i^(n+1)`) is estimated obtained by replacing the spatial and temporal derivatives in the previous first order equation using Forward Difference Scheme.

`∆t`is used for Timestep.

Corrector Step:

In the corrector step, the predicted value `u_i^(n+1)` is corrected according to the equation. Note that the corrector step uses Backward Difference Scheme for spatial derivative. The time-step used in the corrector step is `∆t/2` in contrast to the used in the predictor step.

Main program:

clear all 
close all
clc

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

gamma = 1.4;
c = 0.5;
nt = 1400;

% Function for Non Conservative Form
tic;
[mass_flow_rate,pressure,mach_number,rho,v,T,rho_throat,v_throat,T_throat,mass_flow_rate_throat,pressure_throat,mach_number_throat] = non_conserative(x,dx,n,gamma,nt,c);
simulation_time = toc;
fprintf('Time taken for Non Conservative Form: %f',simulation_time);
hold on

% Timestep varition of properties at throat area in Non Conservative form 
figure(2)
%Density
subplot(4,1,1)
plot(linspace(1,nt,nt),rho_throat,'r')
xlabel('No. of grid points')
ylabel('density')
legend({'Density at throat'})
axis([0 1400 0 1.3])
grid minor
title('Timestep variation at Throat Area- Non Conservative Form')
%Temperature
subplot(4,1,2)
plot(linspace(1,nt,nt),T_throat,'g')
xlabel('No. of grid points')
ylabel('Temperature')
legend({'Temperature at throat'})
axis([0 1400 0.6 1])
grid minor
%pressure
subplot(4,1,3)
plot(linspace(1,nt,nt),pressure_throat,'b')
xlabel('No. of grid points')
ylabel('Pressure')
legend({'Pressure at throat'})
axis([0 1400 0.4 1.1])
grid minor
%Mach Number
subplot(4,1,4)
plot(linspace(1,nt,nt),mach_number_throat,'g')
xlabel('No. of grid points')
ylabel('Mach Number')
legend({'Mach Number at throat'})
axis([0 1400 0.6 1.4])
grid minor
%steady State Simulation for Primitive Values for Non Conservative Form
figure(3)
%Density
subplot(4,1,1)
plot(x,rho,'r')
xlabel('No. of grid points')
ylabel('Density')
legend({'Density'})
axis([0 3 0 1])
grid minor
title('Variation of Flow Rate Distribution-Non Conservative')
%Temperature
subplot(4,1,2)
plot(x,T,'g')
xlabel('No. of grid points')
ylabel('Temperature')
legend({'Temperature'})
axis([0 3 0 1])
grid minor
%pressure
subplot(4,1,3)
plot(x,pressure,'b')
xlabel('No. of grid points')
ylabel('pressure')
legend({'pressure'})
axis([0 3 0 1])
grid minor
%Mach Number
subplot(4,1,4)
plot(x,mach_number,'m')
xlabel('No. of grid points')
ylabel('Mach Number')
legend({'Mach Number'})
axis([0 3 0 4])
grid minor
% 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);
simulation_time_c = toc;
fprintf('Time taken for conervative form: %f',simulation_time_c)

% Timestep varition of properties at throat area in Conservative form 
figure(5)
%Density
subplot(4,1,1)
plot(rho_throat_c,'r')
xlabel('No.of grid points')
ylabel('density')
legend({'Density at throat'})
axis([0 1400 0 1.3])
grid minor
title('Timestep variation at Throat Area- Conservative Form')
%Temperature
subplot(4,1,2)
plot(T_throat_c,'g')
xlabel('No. of grid points')
ylabel('Temperature')
legend({'Temperature at Throat'})
axis([0 1400 0.6 1])
grid minor
%pressure
subplot(4,1,3)
plot(pressure_throat_c,'b')
xlabel('No. of grid points')
ylabel('pressure')
legend({'pressure at Throat'})
axis([0 1400 0.2 0.8])
grid minor
%Mach Number
subplot(4,1,4)
plot(mach_number_throat_c,'m')
xlabel('No. of grid points')
ylabel('Mach Number')
legend({'Mach Number at Throat'})
axis([0 1400 0.6 1.4])
grid minor
%steady State Simulation for Primitive Values for Non Conservative Form
figure(6)
%Density
subplot(4,1,1)
plot(x,rho_c,'r')
xlabel('No. of grid points')
ylabel('Density')
legend({'Density'})
axis([0 3 0 1])
grid minor
title('Variation of Flow Rate Distribution-Conservative Form')
%Temperature
subplot(4,1,2)
plot(x,T_c,'g')
xlabel('No. of grid points')
ylabel('Temperature')
legend({'Temperature'})
axis([0 3 0 1])
grid minor
%pressure
subplot(4,1,3)
plot(x,pressure_c,'b')
xlabel('No. of grid points')
ylabel('pressure')
legend({'pressure'})
axis([0 3 0 1])
grid minor
%Mach Number
subplot(4,1,4)
plot(x,mach_number_c,'m')
xlabel('No. of grid points')
ylabel('Mach Number')
legend({'Mach Number'})
axis([0 3 0 4])
grid minor
%comparison of Nomalized Mass Flow Rate of both Forms
figure(7)
hold on
plot(x,mass_flow_rate,'r','LineWidth',1.25)
hold on
plot(x,mass_flow_rate_c,'b','LineWidth',.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')

Function Code for Non-Conservative Form:

function [mass_flow_rate,pressure,mach_number,rho,v,T,rho_throat,v_throat,T_throat,mass_flow_rate_throat,pressure_throat,mach_number_throat] = non_conserative(x,dx,n,gamma,nt,c);

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

A = 1+2.2*(x-1.5).^2;
throat = find(A==1);

for k = 1:nt
    dt = min(c.*dx./(sqrt(T)+v));
    
    rho_old = rho;
    T_old = T;
    v_old = v;
    
    %Predictor Method
    for i = 2:n-1
        
        dv_dx = (v(i+1)-v(i))/dx;
        drho_dx = (rho(i+1)-rho(i))/dx;
        dT_dx = (T(i+1)-T(i))/dx;
        dlog_A_dx = (log(A(i+1))-log(A(i)))/dx;
        
        %Continuity Equation
        drho_dt_P(i) = -rho(i)*dv_dx-rho(i)*v(i)*dlog_A_dx-v(i)*drho_dx;
        
        %momentum Equation
        dv_dt_P(i) = -v(i)*dv_dx-(1/gamma)*(dT_dx+(T(i)*drho_dx/rho(i)));
        
        %Energy Equation
        dT_dt_P(i) = -v(i)*dT_dx-(gamma-1)*T(i)*(dv_dx+v(i)*dlog_A_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 Method
    
    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;
        dlog_A_dx = (log(A(j))-log(A(j-1)))/dx;
        
        %Continuity Equation
        drho_dt_C(j) = -rho(j)*dv_dx-rho(j)*v(j)*dlog_A_dx-v(j)*drho_dx;
        
        %momentum Equation
        dv_dt_C(j) = -v(j)*dv_dx-(1/gamma)*(dT_dx+(T(j)*drho_dx/rho(j)));
        
        %Energy Equation
        dT_dt_C(j) = -v(j)*dT_dx-(gamma-1)*T(j)*(dv_dx+v(j)*dlog_A_dx);
        
    end
    
    %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
    for g = 2:n-1
        rho(g) = rho_old(g)+drho_dt(g)*dt;
        v(g) = v_old(g)+dv_dt(g)*dt;
        T(g) = T_old(g)+dT_dt(g)*dt;
        
    end
    %Boundary Conditions
    % at Inlet
    rho(1) = 1;
    v(1) = 2*v(2)-v(3);
    T(1)= 1;
    
    %at 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);
    
    mass_flow_rate = rho.*A.*v;
    pressure = rho.*T;
    mach_number = (v./sqrt(T));
    
    rho_throat(k) = rho(throat);
    T_throat(k) = T(throat);
    v_throat(k) = v(throat);
    mass_flow_rate_throat(k) = mass_flow_rate(throat);
    pressure_throat(k) = pressure(throat);
    mach_number_throat(k) = mach_number(throat);
    
    % Comparison of non-dimensional Mass Flow Rate at Different timesteps
    
    figure(1)
   if k==50
        plot(x,mass_flow_rate,'r','LineWidth',1.25)
         hold on
     elseif k==100
         plot(x,mass_flow_rate,'g','LineWidth',1.25)
         hold on
     elseif k==200
         plot(x,mass_flow_rate,'b','LineWidth',1.25)
         hold on
     elseif k==400
         plot(x,mass_flow_rate,'y','LineWidth',1.25)
         hold on
     elseif k==800
         plot(x,mass_flow_rate,'m','LineWidth',1.25)
         hold on
     end
     
     title('Mass Flow Rate at Different Timesteps - Non 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 2])
     grid on
end
end

Function Code 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);

for g = 1:n
    if (x(g) >= 0 && x(g) <= 0.5)
        rho_c(g) = 1;
        T_c(g) = 1;
    elseif (x(g) >= 0.5 && x(g) <= 1.5)
        rho_c(g) = 1-0.366*(x(g)-0.5);
        T_c(g) = 1-0.167*(x(g)-0.5);
    elseif (x(g) >= 1.5 && x(g) <= 3.5)
        rho_c(g) = 0.634-0.3879*(x(g)-1.5);
        T_c(g) = 0.833-0.3507*(x(g)-1.5);
    end
    
end
A = 1+2.2*(x-1.5).^2;
throat = find(A==1);

v_c = 0.59./(rho_c.*A);
pressure_c = rho_c.*T_c;

%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);
dt = min(c.*dx./(sqrt(T_c)+v_c));
for k = 1:nt
    
    U1_old = U1;
    U2_old = U2;
    U3_old = U3;
    
    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
        %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);
        
        %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
    
    %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 q = 2:n-1
        %source term(J)
        J2(q) = (1/gamma)*rho_c(q)*T_c(q)*((A(q)-A(q-1))/dx);
        
        dU1_dt_C(q) = -((F1(q)-F1(q-1))/dx);
        dU2_dt_C(q) = -((F2(q)-F2(q-1))/dx)+J2(q);
        dU3_dt_C(q) = -((F3(q)-F3(q-1))/dx);
    end
    
    %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
    %Apply Boundary Condition 
    %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 Perimitive Value Update
    rho_c = U1./A;
    v_c = U2./U1;
    T_c = (gamma-1)*((U3./U1-(gamma/2)*(v_c).^2));
    
    mass_flow_rate_c = rho_c.*v_c.*A;
    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,'r','LineWidth',1.25);
         hold on
     elseif k==100
         plot(x,mass_flow_rate_c,'g','LineWidth',1.25);
         hold on
     elseif k==200
         plot(x,mass_flow_rate_c,'b','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,'m','LineWidth',1.25);
         hold on
     end
     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'});
     grid on
     
end

end

Code Explanation:

• Starting few commands like ‘clear all’, ‘close all’, ‘clc’ are entered to clear any existing data from previous program or command section.
• Enter the input variables x,dx,n,gamma,nt and c.
• Create separate function codes for Non Conservative and Conservative forms with the appropriate inputs and equations.
• In those function code the values of Mass flow rate, Pressure, Mach Number, Density, Velocity, Temperature,Density at throat, Velocity at throat, Temperature at throat, Mass flow rate at throat, Pressure at throat, Mach number at throat were calculated by the equation for all grid points.
• The initial boundary conditions are established for Density, Temperature and Velocity.
• MacCarmack’s technique is used to solve the numerical solution of governing equations such as Continuity, Momentum and Energy Equation.
• Predictor and corrector method are assigned within the for loop individually for calculating values.
• The calculated values are averaged in the nest step and updating the solution with old values.
• Boundary conditions are applied within for loop and subsequently used for calculating the Mass flow rate, Mach number,Pressure nf also calculating the values at throat section.
• Use tic toc command in main code for calculating time taken between functions.
• Enter the plot command for the followings,

1. Variation of Flow Rate Distribution
2. Time-step variation of primitive variable at Throat area
3. Mass Flow Rates at Different Timesteps
4. Comparison of Normalized Mass Flow Rate of both Non Conservative & Conservative Form

• Save and run the program

Outputs:

Variation of Flow Rate Distribution:

Non-Conservative Form:

Conservative Form:

The above plots shows that the both forms are in same phase giving the exact same solution for all the values.

Time-step variation of primitive variable at Throat area:

Non-Conservative Form:

Conservative Form:

The above plots shows that the time taken for the solution became stable in non conservative form is higher when compare to conservative form.

Mass Flow Rates at Different Timesteps:

Non-Conservative Form:

Conservative Form:

The above plots shows that the solution became stable after the 400­^th timestep for both forms attaining a constant value.

Comparison Plot of Normalized Mass Flow Rate of both Non Conservative & Conservative Form:

The above plot shows the solution in non conservative form oscillates much when compare to conservative form.