Simulation of a 1D Supersonic nozzle flow simulation using Macormack Method

Aim -

Using both the conservation and non-conservation forms of the governing equations, create a Matlab programme to simulate the isentropic flow through a quasi-1D subsonic-supersonic nozzle and answer the problems using MacCormack's method.

Objective -

1. Creating a Matlab programme to mimic the isentropic flow via a nozzle that is nearly one-dimensional and subsonic-supersonic.

2. Using MacCormack's method, the problem is solved by taking into account both the conservation and non-conservation forms of the governing equations.

3. To chart the time-varying primitive variables inside the nozzle and the steady state distribution of those variables.

4. To contrast the steady-state mass flow rates in the conservation and non-conservation forms.

5. To chart the changes in the nozzle's mass flow rate distribution at various time intervals while time-marching.

6. To put the solution through a grid independence test.

Theory -

• Non-Conservation Form

1.Continuity Equation

∂ρ/∂t=−ρ∂V/∂x−ρV∂lnA/∂x−V∂ρ/∂x

2. Momentum Equation

∂V/∂t=−V∂V/∂x−1γ(∂T/∂x+Tρ∂ρ/∂x)

3. Energy Equation

∂T/∂t=−V∂T/∂x−(γ−1)T(∂V/∂x+V∂(lnA)/∂x)

4. Initial Conditions(t=0)

ρ=1-0.3146x">ρ=1−0.3146x�=1-0.3146�

T=1-0.2314x">T=1−0.2314x�=1-0.2314�

V=(0.1+1.09x)T12">V=(0.1+1.09x)T12�=(0.1+1.09�)�12

Nozzle shape A=1+2.2(x-1.5)2">A=1+2.2(x-1.5)2">A=1+2.2(x−1.5)2

A=1+2.2(x-1.5)2">5. Boundary Conditions

A=1+2.2(x-1.5)2">Inlet - Two variables are fixed and one variable float that is density and temperature is fixed and velocity is float.

Therefore,

V1=2V2−V3

Outlet - All the three variables are float

VN=2VN-1-VN-2">VN=2VN−1−VN−2��=2��-1-��-2

ρN=2ρN-1-ρN-2">ρN=2ρN−1−ρN−2��=2��-1-��-2

TN=2TN-1-TN-2">TN=2TN−1−TN−2

• Conservation Form

1.Continuity Equation

∂(ρA)/∂t+∂(ρAV)/∂x=0

2. Momentum Equation

∂(ρAV)/∂t+∂(ρAV2+(1γ)pA)/∂x=1γp∂A/∂x

3. Energy Equation

∂(ρ(T/γ−1+γ2V2)A)/∂t+∂(ρ(Tγ−1+γ2V2)AV+pAV)/∂x=0

Lets define the elements of the solutions vector U, the flux vector F and the source term J as follows

U1=ρA">U1=ρA�1=��

U2=ρAV">U2=ρAV�2=���

U3=ρ(Tγ-1+γ2V2)A">U3=ρ(Tγ−1+γ2V2)A�3=�(��-1+�2�2)�

F1=ρAV">F1=ρAV�1=���

F2=ρAV2+1γpA">F2=ρAV2+1γpA�2=���2+1���

F3=ρ(Tγ-1+γ2V2)VA+pAV">F3=ρ(Tγ−1+γ2V2)VA+pAV�3=�(��-1+�2�2)��+���

J2=1γp∂A∂x">J2=1γp∂A∂x

By using these elements the continuity, the momentum and the energy equation is rewritten as,

∂U1∂t=-∂F1∂x">∂U1∂t=−∂F1∂x∂�1∂�=-∂�1∂�

∂U2∂t=-∂F2∂x+J2">∂U2∂t=−∂F2∂x+J2∂�2∂�=-∂�2∂�+�2

∂U3∂t=-∂F3∂x">∂U3∂t=−∂F3∂x

∂U3∂t=-∂F3∂x">4. Initial Conditions (t=0)

[ρ=1T=1]}for 0≤x≤0.5">[ρ=1T=1]}for 0≤x≤0.5[�=1�=1]}for0≤�≤0.5

[ρ=1-0.366(x-0.5)T=1-0.167(x-0.5)]}for 0.5≤x≤1.5">[ρ=1−0.366(x−0.5)T=1−0.167(x−0.5)]}for 0.5≤x≤1.5[�=1-0.366(�-0.5)�=1-0.167(�-0.5)]}for0.5≤�≤1.5

[ρ=0.634-0.3879(x-1.5)T=0.833-0.3507(x-1.5)]}for 1.5≤x≤3.5">[ρ=0.634−0.3879(x−1.5)T=0.833−0.3507(x−1.5)]}for 1.5≤x≤3.5[�=0.634-0.3879(�-1.5)�=0.833-0.3507(�-1.5)]}for1.5≤�≤3.5

V=0.59ρA">V=0.59ρA�=0.59��

Nozzle shape A=1+2.2(x-1.5)2">A=1+2.2(x−1.5)2

A=1+2.2(x-1.5)2">5. Boundary Conditions

A=1+2.2(x-1.5)2">Inlet - Two variables are fixed and one variable float that is density and temperature is fixed and velocity is float.

therefore,

U2(i=1)=2U2(i=2)−U2(i=3)

U3 is also updated by using the value of updated V

Outlet - All the three variables are float

(U1)N=2(U1)N-1-(U1)N-2">(U1)N=2(U1)N−1−(U1)N−2(�1)�=2(�1)�-1-(�1)�-2

(U2)N=2(U2)N-1-(U2)N-2">(U2)N=2(U2)N−1−(U2)N−2(�2)�=2(�2)�-1-(�2)�-2

(U3)N=2(U3)N-1-(U3)N-2">(U3)N=2(U3)N−1−(U3)N−2

(U3)N=2(U3)N-1-(U3)N-2">Procedure(U3)N=2(U3)N-1-(U3)N-2"> -

1. To begin the programme, use the general commands clear all, close all, and clc to end all of Matlab's pre-occupied operations.

2. Deliver the incoming data.

3. To run the simulation using both conservation and non-conservation versions, use the appropriate routines. The problem is resolved by applying the aforementioned equations plus MacCormack's method.

4. Plot the comparison graphs that were discussed during the goal session.

5. Launch and save the file.

Main Code -

% Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
close all
clear all
clc

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

% Function call
[mdot_nc] = non_cons(n, x, dx, nt, gamma);
[mdot_c] = cons(n, x, dx, nt, gamma);

% Plotting Mass flow rate 
figure(7)
plot(x,mdot_nc(end,:), '-r')
hold on
plot(x,mdot_c(end,:), '-k')
xlabel('Nozzle Length');
ylabel('Massflow Rate')
legend('Non conservative form', 'Conservative form')
title('Comparison of steady state mass flow rate distributions of both forms')

Function Script for Non-Conservation Form -

% Function for non-conservative form
function [mdot] = non_cons(n, x, dx, nt, gamma)

% initializing the varibales 
rho = 1-(0.3146*x);
t = 1-(0.2314*x);
v = (0.1+(1.09*x)).*t.^0.5;
a = (1+2.2*(x-1.5).^2);
throat = find(a==1);

dt = min(0.5*dx./(t.^0.5+v));

% looping starts here 
for k = 1:nt

    rho_old = rho;
    t_old = t;
    v_old = v;

    % Predictor
    for j = 2:n-1
    
        dvdx = (v(j+1)-v(j))/dx;
        drhodx = (rho(j+1)-rho(j))/dx;
        dtdx = (t(j+1)-t(j))/dx;
        dloga_dx = (log(a(j+1))-log(a(j)))/dx;

        drhodt_p(j) = -(rho(j)*dvdx)-(rho(j)*v(j)*dloga_dx)-(v(j)*drhodx);
        dvdt_p(j) = -(v(j)*dvdx)-(1/gamma)*((dtdx)+(t(j)/rho(j))*drhodx);
        dtdt_p(j) = -(v(j)*dtdx)-((gamma-1)*t(j)*(dvdx+(v(j)*dloga_dx)));

        % Solution Update
        rho(j) = rho(j)+(drhodt_p(j)*dt);
        v(j) = v(j)+(dvdt_p(j)*dt);
        t(j) = t(j)+(dtdt_p(j)*dt);
    end

    % Corrector
    for j = 2:n-1

        dvdx = (v(j)-v(j-1))/dx;
        drhodx = (rho(j)-rho(j-1))/dx;
        dtdx = (t(j)-t(j-1))/dx;
        dloga_dx = (log(a(j))-log(a(j-1)))/dx;

        drhodt_c(j) = -rho(j)*dvdx-rho(j)*v(j)*dloga_dx-v(j)*drhodx;
        dvdt_c(j) = -v(j)*dvdx-(1/gamma)*((dtdx)+(t(j)/rho(j))*drhodx);
        dtdt_c(j) = -v(j)*dtdx-(gamma-1)*t(j)*(dvdx+(v(j)*dloga_dx));

    end
    
    % Average
    drhodt = 0.5*(drhodt_p+drhodt_c);
    dvdt = 0.5*(dvdt_p+dvdt_c);
    dtdt = 0.5*(dtdt_p+dtdt_c);

    % Final Solution Update
    for i =2:n-1
        rho(i) = rho_old(i)+drhodt(i)*dt;
        v(i) = v_old(i)+dvdt(i)*dt;
        t(i) = t_old(i)+dtdt(i)*dt;
    end

    % Boundary condition 
    % Inlet
    v(1) = 2*v(2)-v(3);

    % Outlet
    v(n) = 2*v(n-1)-v(n-2);
    rho(n) = 2*rho(n-1)-rho(n-2);
    t(n) = 2*t(n-1)-t(n-2);

    p = rho.*t;
    mach = v./(t.^0.5);
    mdot(k,:) = rho.*a.*v;

    % Varibales at throat
    mach_throat(k) = mach(throat);
    p_throat(k) = p(throat);
    t_throat(k) = t(throat);
    rho_throat(k) = rho(throat);
   
    
end

% Plotting
figure(1)
subplot(4,1,1)
plot(x, p, '-r')
title('Steady-state distribution of primitive variables inside the nozzle (Nondimensional) (Non-Conservative form)')
xlabel('Nozzle Length')
ylabel('Pressure')
subplot(4,1,2)
plot(x, rho, '-b')
xlabel('Nozzle Length')
ylabel('Density')
subplot(4,1,3)
plot(x, t, '-k')
xlabel('Nozzle Length')
ylabel('Temperature')
subplot(4,1,4)
plot(x, mach, '-g')
xlabel('Nozzle Length')
ylabel('Mach Number')

figure(2)
subplot(4,1,1)
plot(1:nt, rho_throat, '-r')
title('Time-wise variation of the primitive variables (Nondimensional) at Throat (Non-Conservative form)')
xlabel('Time Step')
ylabel('Density')
subplot(4,1,2)
plot(1:nt, t_throat, '-b')
xlabel('Time Step')
ylabel('Temperature')
subplot(4,1,3)
plot(1:nt, p_throat, '-k')
xlabel('Time Step')
ylabel('Pressure')
subplot(4,1,4)
plot(1:nt, mach_throat, '-g')
xlabel('Time Step')
ylabel('Mach Number')

figure(3)
plot(x,mdot(50,:))
hold on
plot(x,mdot(100,:))
hold on
plot(x,mdot(200,:))
hold on
plot(x,mdot(600,:))
hold on
plot(x,mdot(1400,:))
xlabel('Nozzle Length')
ylabel('Mass-flow rate')
title('Variation of Mass flow rate (Nondimensional) distribution inside the nozzle at different time steps (Non-Conservative form)')
legend('Time Step = 50','Time Step = 100','Time Step = 200','Time Step = 600','Time Step = 1400')

end

Function Script for Conservation Form -

% Function for conservative form
function [mdot] = cons(n, x, dx, nt, gamma)


% initializing the varibales 
for i = 1:n
    if (x(i)>= 0 && x(i)<0.5)
        rho(i) = 1;
        t(i) = 1;
    elseif (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));
    elseif (x(i)>= 1.5 && x(i)<=3.5)
             rho(i) = 0.634-(0.3879*(x(i)-1.5));
             t(i) = 0.833-(0.35*(x(i)-1.5));

    end 

end
a = 1+2.2*(x-1.5).^2;
cfl = 0.5;
v = (0.59)./(rho.*a);
dt = min(cfl*dx./(t.^0.5+v));
throat = find(a==1);
    u1 = rho.*a;
    u2 = rho.*a.*v;
    u3 = rho.*a.*(t./(gamma-1))+(0.5*gamma.*v.^2);

% looping starts here 
for k =1:nt
    u1_old = u1;
    u2_old = u2;
    u3_old = u3;

 % Evaluationg fluxes
    f1 = u2;
    f2 = (u2.^2)./(u1)+((gamma-1)/gamma).*(u3-0.5*gamma.*((u2.^2)./u1));
    f3 = ((gamma.*u2.*u3)./(u1))-((gamma*(gamma-1))*0.5.*(u2.^3)./u1.^2);

    for j =2:1:(n-1) % preictor loop 
        du1dt_p(j) = -(f1(j+1)-f1(j))/(dx);
        du2dt_p(j) = -(f2(j+1)-f2(j))/(dx)+((1/gamma)*rho(j)*t(j)*((a(j+1)-a(j))/dx));
        du3dt_p(j) = -(f3(j+1)-f3(j))/(dx);
 
        % Solution Update
        u1(j) = u1(j)+(du1dt_p(j)*dt);
        u2(j) = u2(j)+(du2dt_p(j)*dt);
        u3(j) = u3(j)+(du3dt_p(j)*dt);

    end 
        rho = u1./a;
        v = u2./u1;
        t = (gamma-1).*(((u3)./(u1))-((gamma*0.5).*v.^2));

 % Updating fluxes again
    f1 = u2;
    f2 = (u2.^2)./(u1)+((gamma-1)/gamma).*(u3-0.5*gamma.*((u2.^2)./u1));
    f3 = ((gamma.*u2.*u3)./(u1))-((gamma*(gamma-1))*0.5.*(u2.^3)./u1.^2);

    for j =2:1:(n-1) % corrector  loop 
        du1dt_c(j) = -(f1(j)-f1(j-1))/(dx);
        du2dt_c(j) = -(f2(j)-f2(j-1))/(dx)+((1/gamma)*rho(j)*t(j)*((a(j)-a(j-1))/dx));
        du3dt_c(j) = -(f3(j)-f3(j-1))/(dx);
    end

    % finding the averagre derivatives 
        du1dt = 0.5*(du1dt_c+du1dt_p);
        du2dt = 0.5*(du2dt_c+du2dt_p);
        du3dt = 0.5*(du3dt_c+du3dt_p);

    for i =2:1:(n-1) % Final update of the solution U 
        u1(i) = u1_old(i)+du1dt(i)*dt;
        u2(i) = u2_old(i)+du2dt(i)*dt;
        u3(i) = u3_old(i)+du3dt(i)*dt;
    end 

    %% update the boundary conditions
    % Inlet
    u1(1) = rho(1)*a(1);
    u2(1) = (2*u2(2)) - u2(3);
    v(1) = u2(1)/u1(1);
    u3(1) = u1(1)*((t(1)/(gamma-1))+((gamma/2)*v(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);
    
    rho = u1./a;
    v = u2./u1;
    t = (gamma-1).*(((u3)./(u1))-((gamma/2).*(u2./u1).^2));
    mach = v./(t.^0.5);
    p = rho.*t;
    mach_throat(k) = mach(throat);
    p_throat(k) = p(throat);
    t_throat(k) = t(throat);
    rho_throat(k) = rho(throat);
    mdot(k,:) = rho.*a.*v;
    mdot_c = mdot(throat);

end 

% Plotting
figure(4)
subplot(4,1,1)
plot(x, p, '-r')
title('Steady-state distribution of primitive variables inside the nozzle (Nondimensional) (Conservative form)')
xlabel('Nozzle Length')
ylabel('Pressure')
subplot(4,1,2)
plot(x, rho, '-b')
xlabel('Nozzle Length')
ylabel('Density')
subplot(4,1,3)
plot(x, t, '-k')
xlabel('Nozzle Length')
ylabel('Temperature')
subplot(4,1,4)
plot(x, mach, '-g')
xlabel('Nozzle Length')
ylabel('Mach Number')

figure(5)
subplot(4,1,1)
plot(1:nt, rho_throat, '-r')
title('Time-wise variation of the primitive variables (Nondimensional) at Throat (Conservative form)')
xlabel('Time Step')
ylabel('Density')
subplot(4,1,2)
plot(1:nt, t_throat, '-b')
xlabel('Time Step')
ylabel('Temperature')
subplot(4,1,3)
plot(1:nt, p_throat, '-k')
xlabel('Time Step')
ylabel('Pressure')
subplot(4,1,4)
plot(1:nt, mach_throat, '-g')
xlabel('Time Step')
ylabel('Mach Number')

figure(6)
plot(x,mdot(50,:))
hold on
plot(x,mdot(100,:))
hold on
plot(x,mdot(200,:))
hold on
plot(x,mdot(600,:))
hold on
plot(x,mdot(1400,:))
xlabel('Nozzle Length')
ylabel('Mass-flow rate')
title('Variation of Mass flow rate (Nondimensional) distribution inside the nozzle at different time steps (Conservative form)')
legend('Time Step = 50','Time Step = 100','Time Step = 200','Time Step = 600','Time Step = 1400')

end

Results -

• Non-Conservation form (Number of nodes = 31)

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.

• Conservation form (Number of nodes = 31)

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

Comparison of steady state mass flow rate distributions of both forms (n = 31)

It is clear from the numbers that the conservative and non-conservative forms produce the same outcomes. The conditions that are accessible, the nature of the problems, and the parameters that are being measured all influence how these methods are used.

It is evident from the plots that neither style is clearly superior in terms of outcome accuracy. For the primitive variables, the non-conservation form appears to offer somewhat more accurate findings; for the flux variables, the conservation form appears to yield slightly more accurate results and a better distribution of mass flow.

Conservation form demands somewhat more work than non-conservation form when comparing the amount of processing time and calculation.

• Grid Independence Test

ρ∗ρ0">ρ∗ρ0	p∗p0">p∗p0	T∗T0">T∗T0	M
Non-Conservation form (n = 31)	0.6387	0.5342	0.8365	0.9994
Non-Conservation form (n = 61)	0.6377	0.5327	0.8354	0.9998
Non-Conservation form (n = 91)	0.6404	0.5359	0.8368	1.0026
Conservation form (n = 31)	0.6520	0.5470	0.8389	0.9786
Conservation form (n = 61)	0.6388	0.5325	0.8335	0.9983
Conservation form (n = 91)	0.6379	0.5316	0.8333	1.0007
Exact Solution	0.634	0.528	0.833	1

• Comparison of steady state mass flow rate distribution of both form (n = 31)

• Comparison of steady state mass flow rate distribution of both form (n = 61)

• Comparison of steady state mass flow rate distribution of both form (n = 91)

It is clear from the grid independence test table that as the number of nodes increases from 31 to 61, the value gets closer to the exact values. However, the value once more begins to veer towards the non-conservative form at 91. In order to obtain a solution that is closer to the exact answer, it may be concluded that the ideal grid points for the non-conservative form are 61. When the number of grid points increases from 31 to 91 in the conservation form, the values of the primitive variables approach the precise solution more closely. Therefore, it can be concluded that the ideal grid points for the conservation form are 91 in order to provide a solution that is closer to the exact answer.

Concluion -

A Matlab program to simulate the isentropic flow through a quasi 1D subsonic-supersonic nozzle by both the conservation and non-conservation forms of the governing equations by using the MacCormack's technique has been created and a study is also carried out by comparing the results obtained from conservation form and non-conservation form. It is then concluded that there is no clear superiority of either form in terms of accuracy of the results and the usage of these method totally depend upon the available conditions, the nature of problems and the parameters being measured.