Simulation of a Quasi 1D Subsonic-Supersonic Nozzle Flow using the Maccormack s Technique

Objective: Simulation of 1D subsonic-supersonic nozzle flow simulation using the Maccormack’s technique for Non-conservative and conservative forms of governing equations.

Introduction:

Maccormack’s technique is an explicit finite differencing method and a variant of the Lax-Wendroff approach. Maccormack’s method is second order accurate in space as well as in time. It is particularly suited for time marching solutions. This technique will be used to determine the numerical solution of governing equations (conservative and non-conservative form) for a 1D subsonic-supersonic nozzle.

In this study, a convergent-divergent nozzle as depicted in the Figure 1. will be analyzed under steady isentropic flow conditions. The nozzle inlet flow comes from a reservoir with large cross sectional area (A tends to infinity), hence the velocity at the inlet is very small. The inlet pressure and temperature at nozzle is denoted as Po and To. The flow expands isentropically to supersonic speeds at the nozzle exits, where the exit preesure, temperature, velocity and Mach number is denoted by Pe, Te, Ve, Me. The flow is locally subsonic at the inlet (M<1), sonic at the throat (M=1) and supersonic at the exit (M>1). The sonic flow at the exit means that the local velocity at the throat is equal to the local speed of sound. In this study, it is assumed that the flow properties varies only in the x-direction, hence quasi 1D flow.

Figure 1. Schematic for subsonic-supersonic isentropic nozzle flow

Governing equations:

• Non-dimensional non-conservative form of governing equations

Continuity equation

Momentum equation

Energy equation

Non-dimensional length:

Non-dimensional velocity:

`underset(o)(a)`= speed of sound

Non-dimensional time:

Non-dimensional area:

Initial conditions and boundary profiles are expressed as below:

`Density,rho = 1-0.3146*x;`

`Temperature, T = 1 - 0.2314*x;`

`Velocity, v = (0.1 + 1.09*x)*t^0.5;`

`Area, A = 1 + 2.2*(x - 1.5)^2;`

Figure 2. Grid point distribution along the nozzle

Inlet:

Outlet:

• Non-dimensional conservative form of governing equations

Continuity equation

Momentum equation

Energy equation

Defining the elements of the solutions vector U, the flux vector F and the source term J as,

With the above elements, the non dimensional conservative form of governing equations can be expressed as,

Initial conditions and Boundary profiles are as expressed below:

Inlet:

Outlet:

Grid independence test:

The numerical solution achieved thorugh CFD is dependent on the grid system adopted. Moreover, the grid spacing also impacts the accuracy and computational time taken to generate an appreciable numerical solution. Grid independence test is done to check whether the solution alters from change in the number of grid points (N). If the solution doesn\'t change with increase in grid points, then the solution is said to be grid independent. In this study, N=31 and N=61 is used to analyze the flow variable property as well as for conducting the grid independence.

Matlab code for implementation of function:

clear all
close all
clc
% Inputs
L = 3;  % Length of the domain
n = 31; % Number of the grid points
x = linspace(0,L,n);
dx = x(2) - x(1);   % Grid spacing
a = 1+2.2*(x-1.5).^2;   % Area
gamma = 1.4
% node at the throat
throat = ((n-1)/2)+1;
% courant number
c = 0.5
% No. of time steps
nt = 1400;

% non conservative form
[v1,rho1,t1,m1,p1,mach1,v_t1,rho_t1,t_t1,m_t1,p_t1,mach_t1,k1,SimulationTime_NC]=non_conservative(n,x,dx,c,gamma,a,nt,throat);
% conservative form
[v2,rho2,t2,m2,p2,mach2,v_t2,rho_t2,t_t2,m_t2,p_t2,mach_t2,k2,Simulationtime_C]=conservative(n,x,dx,c,gamma,a,nt,throat);

% plots comparing nonconservative and conservative
figure(3)
s1=subplot(4,1,1);
plot(linspace(1,nt,nt), rho_t1);
hold on
plot(linspace(1,nt,nt), rho_t2);
hold on
plot(linspace(nt-400,nt,800),linspace(0.634,0.634,800),\'color\',\'g\');
grid on
ylabel(\"rho\'\",\'fontsize\',12,\'fontweight\',\'bold\')
h=legend(s1,{\'Non-conservative Form\', \'Conservative Form\',\'Analytical solution\'})
set(h, \'Location\',\'northeastoutside\')
title(\'Timewise variation of flow parameters at throat\',\'fontsize\',12,\'fontweight\',\'bold\')

s2=subplot(4,1,2)
plot(linspace(1,nt,nt), t_t1);
hold on
plot(linspace(1,nt,nt), t_t2);
hold on
plot(linspace(nt-400,nt,800),linspace(0.833,0.833,800),\'color\',\'g\');
grid on
ylabel(\"T\'\",\'fontsize\',12,\'fontweight\',\'bold\')
h=legend(s2,{\'Non-conservative Form\', \'Conservative Form\',\'Analytical solution\'})
set(h, \'Location\',\'northeastoutside\')

s3=subplot(4,1,3)
plot(linspace(1,nt,nt), p_t1);
hold on
plot(linspace(1,nt,nt), p_t2);
hold on
plot(linspace(nt-400,nt,800),linspace(0.528,0.528,800),\'color\',\'g\');
grid on
ylabel(\"pr\'\",\'fontsize\',12,\'fontweight\',\'bold\')
h=legend(s3,{\'Non-conservative Form\', \'Conservative Form\',\'Analytical solution\'})
set(h, \'Location\',\'northeastoutside\')
legend(\'Non-conservative Form\', \'Conservative Form\',\'Analytical solution\')

s4=subplot(4,1,4)
plot(linspace(1,nt,nt), mach_t1);
hold on
plot(linspace(1,nt,nt), mach_t2);
hold on
plot(linspace(nt-400,nt,800),linspace(1,1,800),\'color\',\'g\');
grid on
xlabel(\"Number of time steps\",\'fontsize\',12,\'fontweight\',\'bold\')
ylabel(\"Ma\'\",\'fontsize\',12,\'fontweight\',\'bold\')
h=legend(s4,{\'Non-conservative Form\', \'Conservative Form\',\'Analytical solution\'})
set(h, \'Location\',\'northeastoutside\')

% variation in flow parameters along x axis
figure(4)
q1=subplot(4,1,1)
plot(x,rho1,\'r\')
hold on
plot(x,rho2,\'b\')
grid on
ylabel(\"rho\'\",\'fontsize\',12,\'fontweight\',\'bold\')
h = legend(q1,{\'Non conservative Form\',\'Conservative Form\'});
set(h,\'location\',\'northeastoutside\')
title(\'Variation of flow parameters along x axis\',\'fontsize\',12,\'fontweight\',\'bold\')

q2=subplot(4,1,2)
plot(x,t1,\'r\')
hold on
plot(x,t2,\'b\')
grid on
ylabel(\"T\'\",\'fontsize\',12,\'fontweight\',\'bold\')
h = legend(q2,{\'Non conservative Form\',\'Conservative Form\'});
set(h,\'location\',\'northeastoutside\')

q3=subplot(4,1,3)
plot(x,p1,\'r\')
hold on
plot(x,p2,\'b\')
grid on
ylabel(\"pr\'\",\'fontsize\',12,\'fontweight\',\'bold\')
h = legend(q3,{\'Non conservative Form\',\'Conservative Form\'});
set(h,\'location\',\'northeastoutside\')

q4=subplot(4,1,4)
plot(x,mach1,\'r\')
hold on
plot(x,mach2,\'b\')
grid on
xlabel(\'Non-Dimensional Length\',\'fontsize\',12,\'fontweight\',\'bold\')
ylabel(\"Ma\'\",\'fontsize\',12,\'fontweight\',\'bold\')
h = legend(q4,{\'Non conservative Form\',\'Conservative Form\'});
set(h,\'location\',\'northeastoutside\')

Matlab code for non-conservative form of governing equation:

function [v,rho,t,m,p,mach,v_t,rho_t,t_t,m_t,p_t,mach_t,k,Simulationtime_NC] = non_conservative(n,x,dx,c,gamma,a,nt,throat)

tic
% Inputs

L = 3   % Length of the domain
n = 31  % no of the grid points
x =linspace(0,L,n)
dx = x(2) - x(1)    % Grid Spacing
gamma = 1.4

% Initial profiles
rho = 1 - 0.3146*x; % Density profile
t = 1 - 0.2314*x;   % Temperature profile
v = (0.1 + 1.09*x).*t.^0.5;   % Velocity profile
a  = 1 + 2.2*(x-1.5).^2;    % Area profile

% Time step definition

nt = 1400; % Number of time step

% Determination of dt using the CFL stability condition
for i = 1:n
    % delta t = c * ((dx)/(a+v)), where a = sqrt(t)
    del_t(i)= c*(dx/(t(i)^0.5+v(i)));  
end
    % Choose the minimum value of delta t obtained across grid points
    dt = min(del_t(i));

% Outer time loop
for k = 1:nt
    
    rho_old = rho;
    v_old = v;
    t_old = t;
    a_old = a;
    
    % predictor method
    for j = 2:n-1   % Solving only for the inner nodes
        dvdx = (v(j+1)-v(j))/dx;
        drhodx = (rho(j+1)-rho(j))/dx;
        dtdx = (t(j+1)-t(j))/dx;
        dlogadx = (log(a(j+1)) - log(a(j)))/dx;
        % continuity equation
        drhodt_p(j) = -rho(j)*dvdx - rho(j)*v(j)*dlogadx - v(j)*drhodx;
        % Momentum equation
        dvdt_p(j) = -v(j)*dvdx - (1/gamma)*(dtdx + (t(j)/rho(j))*drhodx);
        % Energy equation
        dtdt_p(j) = -v(j)*dtdx - (gamma - 1)*t(j)*(dvdx + v(j)*dlogadx);
        
        % Updating solution variables
        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 method
    for j = 2:n-1   % Solving only for the inner nodes
        dvdx = (v(j)-v(j-1))/dx;
        drhodx = (rho(j)-rho(j-1))/dx;
        dtdx = (t(j)-t(j-1))/dx;
        dlogadx = (log(a(j)) - log(a(j-1)))/dx;
        % continuity equation
        drhodt_c(j) = -rho(j)*dvdx - rho(j)*v(j)*dlogadx - v(j)*drhodx;
        % Momentum equation
        dvdt_c(j) = -v(j)*dvdx - (1/gamma)*(dtdx + (t(j)/rho(j))*drhodx);
        % Energy equation
        dtdt_c(j) = -v(j)*dtdx - (gamma - 1)*t(j)*(dvdx + v(j)*dlogadx);
    end
    % compute the average time derivative
    
    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
    % Applying Boundary conditions
    
    % 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);
    
    % Calculating solutions
    % Mass flow rate
    m = rho.*a.*v;
    % pressure
    p = rho.*t;
    % Mach number
    mach = v./((t).^0.5);
    
    % Throat values
    rho_t(k) = rho(throat);
    v_t(k)=v(throat);
    t_t(k)= t(throat);
    p_t(k)=p(throat);
    m_t(k)=m(throat);
    mach_t(k)= mach(throat);
    
    % Plotting mass flow rate at different time steps
    figure(1);
    if k == 50
    plot(x,m,\'r\',\'linewidth\',2)
    hold on
    else
    if k == 100
    plot(x,m,\'b\',\'linewidth\',2)
    else
    if k == 200
    plot(x,m,\'y\',\'linewidth\',2)
    else
    if k == 400
    plot(x,m,\'g\',\'linewidth\',2)
    else
    if k == 800
    plot(x,m,\'m\',\'linewidth\',2)
    else
    if k == 1400
    plot(x,m,\'c\',\'linewidth\',2)
    end
    end
    end
    end
    end
    end
    grid on
    
    xlabel(\'Non-Dimensional distance through nozzle(x)\',\'fontweight\',\'bold\',\'fontsize\',12)
    ylabel(\'Non-Dimensional mass flow rate through nozzle(m)\',\'fontweight\',\'bold\',\'fontsize\',12)
    title(\'Non-Dimensional mass flow rate(Non-conservative form) along x at different time steps\',\'fontweight\',\'bold\',\'fontsize\',12,\'color\',\'k\')
    legend({\'50 time steps\',\'100 time steps\',\'200 time steps\',\'400 time steps\',\'800 timesteps\',\'1400 time steps\'})
    Simulationtime_NC = toc
end

Matlab code for conservative form of governing equation:

function [v,rho,t,m,p,mach,v_t,rho_t,t_t,m_t,p_t,mach_t,k,Simulationtime_C] = conservative(n,x,dx,c,gamma,a,nt,throat)

tic
% Inputs
L = 3   % Length of the domain
n = 31  % no of the grid points
x =linspace(0,L,n)
dx = x(2) - x(1)    % Grid spacing
gamma = 1.4

% Initial profile
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

% Assuming steady state flow for the initial profile
v = 0.59./(rho.*a);
% Time steps
nt = 1400
% calculating the value of time step using courant criteria
for i = 1:n
    del_t(i)= c*(dx/(t(i)^0.5+v(i)));
end
dt = min(del_t(i));

% Variables
p = rho.*t;
U1 = rho.*a;
U2 = rho.*a.*v;
U3 = rho.*((t./(gamma-1))+((gamma/2).*(v.^2))).*a;

% outer Time loop
for k = 1:nt;
    U1_old = U1;
    U2_old = U2;
    U3_old = U3;
    
    % Flux variables
    F1 = U2;
    F2 = ((U2.^2)./U1) + ((gamma-1)/gamma)*(U3-(((gamma/2)*U2.^2)./U1));
    F3 = ((gamma*U2.*U3)./U1) - (gamma*(gamma-1)/2)*((U2.^3)./(U1.^2));
    
    % predictor step
    for j=2:n-1
        dF1dx = (F1(j+1)-F1(j))/dx;
        dF2dx = (F2(j+1)-F2(j))/dx;
        dF3dx = (F3(j+1)-F3(j))/dx;
        dlogadx = (log(a(j+1))-log(a(j)))/dx
        dadx = (a(j+1)-a(j))/dx
        J2 = (1/gamma).*rho(j).*t(j).*dadx
        
        % continuity equation
        dU1dt_p(j) = -dF1dx
        % Momentum equation
        dU2dt_p(j) = -dF2dx + J2
        % Energy equation
        dU3dt_p(j) = -dF3dx
        
        % 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
    
        % Primitive terms for corrector step
        rho = U1./a;
        v = U2./U1;
        t = (gamma-1)*((U3./U1)-((gamma/2)*((U2./U1).^2)));
        p = rho.*t;
        F1 = U2;
        F2 = ((U2.^2)./U1) + ((gamma-1)/gamma)*(U3-(((gamma/2)*U2.^2)./U1));
        F3 = ((gamma*U2.*U3)./U1)-(gamma*(gamma-1)/2)*((U2.^3)./(U1.^2));
        
        % corrector step
    for j=2:n-1
        dF1dx = (F1(j)-F1(j-1))/dx;
        dF2dx = (F2(j)-F2(j-1))/dx;
        dF3dx = (F3(j)-F3(j-1))/dx;
        dlogadx = (log(a(j))-log(a(j-1)))/dx
        dadx = (a(j)-a(j-1))/dx
        J2 = (1/gamma).*rho(j).*t(j).*dadx
        
        % continuity equation
        dU1dt_c(j) = -dF1dx
        % Momentum equation
        dU2dt_c(j) = -dF2dx + J2
        % Energy equation
        dU3dt_c(j) = -dF3dx
    end
        % Average derivatives
        dU1dt = 0.5*(dU1dt_p + dU1dt_c);
        dU2dt = 0.5*(dU2dt_p + dU2dt_c);
        dU3dt = 0.5*(dU3dt_p + dU3dt_c);
        
        % Final solution update
    for i = 2:n-1
        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
    
        % Applying the Boundary conditions
        % Inlet
        U1(1)=rho(1)*a(1);
        U2(1)= 2*U2(2)-U2(3);
        v(1)= U2(1)/U1(1);
        t(1) = 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)));
        
        % caculating solutions
        % Mass flow rate
        m = rho.*v.*a;
        % pressure
        p = rho.*t;
        % Mach number
        mach = v./((t).^0.5);
        
        % Throat values
        rho_t(k) = rho(throat);
        v_t(k)=v(throat);
        t_t(k)= t(throat);
        p_t(k)=p(throat);
        m_t(k)=m(throat);
        mach_t(k)= mach(throat);
        
        % Plotting mass flow rate at different time steps
        figure(2);
        if k == 50
        plot(x,m,\'r\',\'linewidth\',2)
        hold on
        else
        if k == 100
        plot(x,m,\'b\',\'linewidth\',2)
        else
        if k == 200
        plot(x,m,\'y\',\'linewidth\',2)
        else
        if k == 400
        plot(x,m,\'g\',\'linewidth\',2)
        else
        if k == 800
        plot(x,m,\'m\',\'linewidth\',2)
        else
        if k == 1400
        plot(x,m,\'c\',\'linewidth\',2)
        end
        end
        end
        end
        end
        end
        grid on
        
        xlabel(\'Non-Dimensional distance through nozzle(x)\',\'fontweight\',\'bold\',\'fontsize\',12)
        ylabel(\'Non-Dimensional mass flow rate through nozzle(m)\',\'fontweight\',\'bold\',\'fontsize\',12)
        title(\'Non-Dimensional mass flow rate(conservative form) along x at different time steps\',\'fontweight\',\'bold\',\'fontsize\',12,\'color\',\'k\')
        legend({\'50 time steps\',\'100 time steps\',\'200 time steps\',\'400 time steps\',\'800 timesteps\',\'1400 time steps\'})
        Simulationtime_C = toc
end
        

Results and Discussions:

The numerical solution for quasi 1D nozzle flow was estimated in Non-conservative as well as conservative form using the Maccormack\'s techinque.

• N = 31

Figure 3. Timewise variation of nondimensional mass flow as a function of distance along the nozzle (Non-conservative form)

Figure 4. Timewise variation of nondimensional mass flow as a function of distance along the nozzle (Conservative form)

An insight into the timewise variation of non-dimensional mass flow rate along the length of the nozzle is shown in Figure 3 and Figure 4. Six different curves are plotted in each figure and each one is for a different time during the course of the time marching procedure. It is noticed that there is a considerable change in flow after 50 time steps. Similarily the flow keeps on changing with respect to time. By 400 times steps, the mass flow distribution starts to settle and by 800 time steps the flow distribution is a staright horizontl line. The lines within 800 and 1400 time steps ovelaps. This suggests that the mass flow is converged to a constant, steady state value. throughout the nozzle.

Figure 5. Timewise variation of density, temperature, pressure and Mach number at the nozzle throat

Figure 5 illustrates transient variations of the flow field variables. The largest change take place at the early times after which the final, steady state value is approached almost asymptotically. It is to be noted that the numerical time marching procedure converges to a proper theoretical steady state answer.

Figure 6. Variation of flow parameters along the length of the nozzle

The numerical simulation results obtained with N=61 is detailed below from Figure 7-10.

• N = 61

Figure 7. Timewise variation of nondimensional mass flow as a function of distance along the nozzle (Non-conservative form)

Figure 8. Timewise variation of nondimensional mass flow as a function of distance along the nozzle (Conservative form)

Figure 9. Timewise variation of density, temperature, pressure and Mach number at the nozzle throat

Figure 10. Variation of flow parameters along the length of the nozzle

Table 1. Comparison of numerical flow field variables with analytical solution

All the flow field variables values are fairly close to the exact analytical solution, however doubling the grid points results in marginally improving the numerical solution. Although as stated it is marginal, hence it is reasonable to state that N=31 provides a grid independent numerical solution. Furthermore, the numerical solution achieved while solving the non-conservative form of governing equations tends to provide flow field variable values close to the analytical solution.

Table 2. Estimation of simulation time

It is inferred from table 2 that the conservative form of governing equations takes longer computational time compared to its counterpart. Moreover, while doubling the number of grid points i.e. N=31 to 61, computational time for solving the conservative form also doubled. Hence, based on the results, the non-conservative form is much faster compared to the conservative form.