Simulation of Quasi 1-D Supersonic nozzle flow using Maccormack method

Objectives of the Project:

(1) Solve the Conservative and Non-conservative form of governing equation of Quasi 1-D Supersonic nozzle flow using the Maccormack method.

(2) Compare the results of the Conservative and Non-conservative form.

(3) Perform Grid dependence test.

Theory:

The Maccormack method-

This method is employed to solve the governing equations of quasi 1D supersonic nozzle flow rather than the differencing schemes as we have studied earlier. Maccormack method is an explicit finite difference technique which is second order accurate in both time and space.

According to the maccormack method, the value of variable at grid point i,j at time level t+dt is given as follows-

where the average time derivative is the mean value between time t and t+dt. The average time derivative is obtained from the predictor-corrector philosophy.

Predictor Step-

In the predictor step, the spatial terms on the right-hand side of the equation are replaced by forward differences. The RHS of the equation is known as all the values of variables are for time-step 't'. Next, the value of the variable on the LHS is calculated by employing the time marching solution and for time-step 't+dt' value of the variable at all grid points is obtained.

Consider the following equation for example-

This is the predicted value of variable which will be used in the corrector step to evaluate the corrected value of variable.

Corrector Step-

In the corrector step, the spatial terms on the right-hand side are replaced by backward differences. The value of variables in RHS are taken which are obtained in the predicted step for time-step 't+dt'. Then the corrected time derivative is calculated.

Now, the predicted and corrected value of time derivative is used to calculate the average time derivative which in turn gives the value of the variable at a given point for a particular time-step. The same steps are repeated until the values of variables are obtained for all grid points for a given number of time-steps.

Governing equations for Quasi 1-D Supersonic nozzle flow-

The conservative and non-conservative form of governing equations are derived using the integral form of continuity, momentum and energy equation.

The derivation is given in the book- Computational Fluid Dynamics by John Anderson Jr.

(1) Non-conservative Form-

All the variables are the non-dimensional values. These equations are solved using Maccormack method and the steady-state values of all variables i.e. Temperature, Velocity and density are obtained. From these variables, the Mach number and pressure is obtained.

Calculation of Time-step-

In order for the solution to be stable, the concept of Courant flow number was introduced in previous projects. For a linear Hyperbolic equation, the value of CFL<1 will give a stable solution.

For this project, CFL=0.5 and value of time-step is calculated as follows-

where, C is the Courant number, V is the local flow velocity at a point and a is the local speed of sound. Value of dt at all grid points is calculated and the minimum of them is chosen for the simulation.

Boundary Conditions-

At the inlet:

Velocity at the inlet is changing with time, while the density and temperature are constant with respect to time.

Velocity at the inlet is calculated by the method of linear extrapolation. Values of velocity at node 2 and 3 are used for the same.

At the Outlet:

All the variables at the outlet are kept changing with respect to time and linear extrapolation is used to calculate the variables at the outlet.

Initial Conditions-

Area of the nozzle will remain constant throughout the time of simulation.

The flow-field variables vary with respect to time. At t=0, their values are given as-

(2) Conservative Form-

The governing equations are expressed in a generic form which consists of the solution vector U, the flux vector F and the source term J.

The equations for Quasi 1-D supersonic flow can be expressed in a similar generic form as follows-

The equations which are solved using Maccormack method are-

Further, the flux vectors and source terms are expressed in terms of solution vectors. Also, to find the field variables, they are written in terms of solution vectors.

Also,

Calculation of time step is the same as that of the procedure followed for non-conservation form.

Boundary conditions-

At inlet:

U2 and U3 are changing with respect to time whereas U1 is kept as constant.

At outlet:

Initial conditions-

Area of the nozzle is constant throughout the runtime of simulation and its value is same as that taken for non-conservative form.

Initial velocity is given as,

At t=0, values of density and temperature throughout the nozzle are given as follows,

MATLAB Script:

Non-conservative form of governing equation

%% 1D Supersonic Nozzle Flow (Non-Conservative form) using Maccormack method

clc
clear all
close all

% Length Domain
L=3;
n=31;
x=linspace(0,L,n);
dx=x(2)-x(1);
gamma=1.4;

% Initial Conditions
rho_initial=1-(0.3146*x);
T_initial=1-(0.2314*x);
u_initial=(0.1+(1.09*x)).*(T_initial.^0.5);
A=1+(2.2*((x-1.5).^2));

% grid point for throat
th=find(A==1);

% Field variables
rho=rho_initial;
T=T_initial;
u=u_initial;

% Time Domain (iterations carried out=1400)
CFL=0.5;
for z=1:n
    dt_ini(z)=CFL*(dx/(((T_initial(z))^0.5)+(u_initial(z))));
end
dt=min(dt_ini);
n_itr=1400;
t_sim=(n_itr-1)*dt;
t=0:dt:t_sim;

%% Maccormack method
for k=1:length(t)
    rho_old=rho;
    T_old=T;
    u_old=u;
    
    for i=2:n-1
        % Defining spatial derivatives according to forward difference scheme
        dudx=(u(i+1)-u(i))/dx;
        dlnAdx=(log(A(i+1))-log(A(i)))/dx;
        drhodx=(rho(i+1)-rho(i))/dx;
        dTdx=(T(i+1)-T(i))/dx;
    
        % Predictor step
        drhodt_p(i)=-(rho(i)*dudx)-((rho(i)*u(i))*dlnAdx)-(u(i)*drhodx);
        dudt_p(i)=-(u(i)*dudx)-((1/gamma)*(dTdx+((T(i)/rho(i))*drhodx)));
        dTdt_p(i)=-(u(i)*dTdx)-((gamma-1)*T(i)*(dudx+(u(i)*dlnAdx)));
        
        % Step update
        rho(i)=rho_old(i)+(drhodt_p(i)*dt);
        T(i)=T_old(i)+(dTdt_p(i)*dt);
        u(i)=u_old(i)+(dudt_p(i)*dt);        
    end
    
    for j=2:n-1
        % Defining spatial derivatives according to backward difference scheme
        dudx=(u(j)-u(j-1))/dx;
        dlnAdx=(log(A(j))-log(A(j-1)))/dx;
        drhodx=(rho(j)-rho(j-1))/dx;
        dTdx=(T(j)-T(j-1))/dx;
    
        % Corrector step
        drhodt_c(j)=-(rho(j)*dudx)-((rho(j)*u(j))*dlnAdx)-(u(j)*drhodx);
        dudt_c(j)=-(u(j)*dudx)-((1/gamma)*(dTdx+((T(j)/rho(j))*drhodx)));
        dTdt_c(j)=-(u(j)*dTdx)-((gamma-1)*T(j)*(dudx+(u(j)*dlnAdx)));
    end
    
    for m=2:n-1
        % Final value for the particular time step
        drhodt_avg(m)=(0.5*dt*(drhodt_p(m)+drhodt_c(m)));
        dTdt_avg(m)=(0.5*dt*(dTdt_p(m)+dTdt_c(m)));
        dudt_avg(m)=0.5*dt*(dudt_p(m)+dudt_c(m));
        
        rho(m)=rho_old(m)+drhodt_avg(m);
        T(m)=T_old(m)+dTdt_avg(m);
        u(m)=u_old(m)+dudt_avg(m);
    end
    
    % At the boundary (by interpolation method)
    % Inlet
    u(1)=(2*u(2))-u(3);
    % Outlet
    u(n)=(2*u(n-1))-u(n-2);
    T(n)=(2*T(n-1))-T(n-2);
    rho(n)=(2*rho(n-1))-rho(n-2);
    % Pressure ratio
    p=rho.*T;
    % Mach number
    M=u./(T.^0.5);
    % Mass flow rate
    mass=rho.*u.*A;
    
    % Field variables at throat section
    rho_throat(k)=rho(th);
    T_throat(k)=T(th);
    p_throat(k)=p(th);
    M_throat(k)=M(th);
    
    % Throat section residual values
    rho_res(k)=drhodt_avg(th);
    T_res(k)=dTdt_avg(th);
    u_res(k)=dudt_avg(th);

    % Plot of mass flow rate at different time instances
    plot_mass(x,k,mass)
  
end

% Plots at throat section
plot_throat_2(rho_throat,p_throat,T_throat,M_throat)
plot_throat_3(rho_res,T_res,u_res)

% Plots at Steady-state
plot_ss(x,M,p,rho,T)

Functions to plot the graphs -

Mass flow rate at different time instances

function[]=plot_mass(x,k,mass)

figure(1)
hold on
if (k==1)
    plot(x,mass)
    elseif (k==51)
        plot(x,mass,'g')
        elseif (k==101)
            plot(x,mass,'r')
            elseif (k==151)
                plot(x,mass,'c')
                elseif (k==201)
                    plot(x,mass,'y')
                    elseif (k==701)
                        plot(x,mass,'black')
                        legend('0dt','50dt','100dt','150dt','200dt','700dt','Location','best')
                        xlabel('Non-dimensional distance through nozzle')
                        ylabel('Non-dimensional mass flow rate')
                        title('Distribution of Mass flow rate in the nozzle at different time instances')
                        grid on
end
end

Field variables at throat section i.e. i=16

function[]=plot_throat_2(rho_throat,p_throat,T_throat,M_throat)

% Exact analytical solution at throat at steady state
rho_exact=0.634*ones(1,1400);
T_exact=0.833*ones(1,1400);
p_exact=0.528*ones(1,1400);
M_exact=1*ones(1,1400);

% Plots for convergence (at throat section)
figure(2)
plot(rho_throat)
hold on 
plot(rho_exact,'r')
xlabel('Number of iterations')
ylabel('Density')
legend('Numerical solution','Exact solution','Location','best')
title('Density at Throat section')
    
figure(3)
plot(T_throat)
hold on 
plot(T_exact,'r')
xlabel('Number of iterations')
ylabel('Temperature')
legend('Numerical solution','Exact solution','Location','best')
title('Temperature at Throat section')
    
figure(4)
plot(p_throat)
hold on 
plot(p_exact,'r')
xlabel('Number of iterations')
ylabel('Pressure')
legend('Numerical solution','Exact solution','Location','best')
title('Velocity at Throat section')

figure(5)
plot(M_throat)
hold on 
plot(M_exact,'r')
xlabel('Number of iterations')
ylabel('Mach number')
legend('Numerical solution','Exact solution','Location','best')
title('Mach number at Throat section')

end

Residual at throat section i.e at i=16

function[]=plot_throat_3(rho_res,T_res,u_res)

% Plots for Residual (at throat section)
figure(6)
plot(rho_res)
hold on
plot(T_res,'g')
hold on
plot(u_res,'r')
xlabel('Number of iterations')
ylabel('Residual')
title('Residual at Throat section')
legend('Density','Temperature','Velocity','Location','best')

end

Steady-state value of field variables

function[]=plot_ss(x,M,p,rho,T)

% Plots at Steady-state
figure(7)
subplot(4,1,1)
plot(x,M)
xlabel('Length of nozzle')
ylabel('Mach number')
grid minor
subplot(4,1,2)
plot(x,p)
xlabel('Length of nozzle')
ylabel('Pressure ratio')
grid minor
subplot(4,1,3)
plot(x,rho)
xlabel('Length of nozzle')
ylabel('Density ratio')
grid minor
subplot(4,1,4)
plot(x,T)
xlabel('Length of nozzle')
ylabel('Temperature ratio')
grid minor
suptitle('Qualitative aspects of Quasi 1-D Supersonic nozzle flow')

end

Results for Non-conservative form

This plot shows the steady-state values of field variables throughout the nozzle. As the nozzle is a convergent-divergent type, the Pressure decreases and the Mach number increases along the length of the nozzle. The profile obtained is in accordance with the profile given in the book by Anderson Jr.

The above four plots show the variation in the value of field variables as the simulation advances. It is observed that after some 600-700 iterations, the solution converges. Although the final solution does not exactly match with the exact analytical solution, the result obtained is satisfactory.

The above plot shows the residual term for each variable as the number of iterations increases. It can be seen that the residual is very small like 1e-6 or nearly equal to zero after 600-650 iterations. Similarly, from the convergence plot, it is seen that solution converges after roughly 650 iterations for all field variables.

This plot shows the mass flow rate in the nozzle at different time instances. After 700dt, the mass flow rate is nearly constant throughout the nozzle.

MATLAB Script:

Conservative form of governing equation

%% 1D Supersonic Nozzle Flow (Conservative form) using Maccormack method

clc
clear all
close all

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

%% Initial conditions (primitive variables)
% density and temperature
for i=1:n
if (x(i)<=0.5)
    rho_initial(i)=1;
    T_initial(i)=1;
elseif (x(i)>0.5) && (x(i)<=1.5)
        rho_initial(i)=1-(0.366*(x(i)-0.5));
        T_initial(i)=1-(0.167*(x(i)-0.5));
else
    rho_initial(i)=0.634-(0.3879*(x(i)-1.5));
    T_initial(i)=0.833-(0.3507*(x(i)-1.5));
end
end
% area
A=1+(2.2*((x-1.5).^2));
% grid point for throat
th=find(A==1);
% velocity
u_initial=0.59./(rho_initial.*A);

rho=rho_initial;
T=T_initial;
u=u_initial;

%% Time Domain (iterations carried out=1400)
CFL=0.5;
for z=1:n
    dt_ini(z)=CFL*(dx/(((T_initial(z))^0.5)+(u_initial(z))));
end
dt=min(dt_ini);
n_itr=1400;
t_sim=(n_itr-1)*dt;
t=0:dt:t_sim;

%% Initial conditions (Solution matrix)
U1_initial=rho_initial.*A;
U2_initial=rho_initial.*A.*u_initial;
U3_initial=rho_initial.*A.*((T_initial/(gamma-1))+((gamma/2)*(u_initial.^2)));

U1=U1_initial;
U2=U2_initial;
U3=U3_initial;

%% Maccormack method
for j=1:length(t)
    U1_old=U1;
    U2_old=U2;
    U3_old=U3;
    
    % Flux variables for predictor step
    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)));

    for k=2:n-1
        dAdx=((A(k+1))-(A(k)))/dx;
        J2=(1/gamma)*rho(k)*T(k)*dAdx;
        
        % Predictor step
        dU1dt_p(k)=-((F1(k+1)-F1(k))/dx);
        dU2dt_p(k)=-((F2(k+1)-F2(k))/dx)+J2;
        dU3dt_p(k)=-((F3(k+1)-F3(k))/dx);
        
        % Step update
        U1(k)=(dU1dt_p(k)*dt)+U1_old(k);
        U2(k)=(dU2dt_p(k)*dt)+U2_old(k);
        U3(k)=(dU3dt_p(k)*dt)+U3_old(k);
    end
    
    % At the boundary (by interpolation method)
    % Inlet
    U2(1)=(2*U2(2))-U2(3);
    U3(1)=U1(1)*((T(1)/(gamma-1))+((gamma/2)*(u(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);
    
    % Values of primitive variables for corrector step
    rho=U1./A;
    u=U2./U1;
    T=(gamma-1)*((U3./U1)-((gamma/2)*(u.^2)));
    
    % Flux variables for corrector step
    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)));
       
    for m=2:n-1
        dAdx=((A(m))-(A(m-1)))/dx;
        J2=(1/gamma)*rho(m)*T(m)*dAdx;
        
        % Corrector step
        dU1dt_c(m)=-((F1(m)-F1(m-1))/dx);
        dU2dt_c(m)=-((F2(m)-F2(m-1))/dx)+J2;
        dU3dt_c(m)=-((F3(m)-F3(m-1))/dx);
    end
    
    for p=2:n-1
        % Average values
        dU1dt_avg(p)=(0.5*dt*(dU1dt_p(p)+dU1dt_c(p)));
        dU2dt_avg(p)=(0.5*dt*(dU2dt_p(p)+dU2dt_c(p)));
        dU3dt_avg(p)=(0.5*dt*(dU3dt_p(p)+dU3dt_c(p)));
        
        % Final values for a particular time-step
        U1(p)=dU1dt_avg(p)+U1_old(p);
        U2(p)=dU2dt_avg(p)+U2_old(p);
        U3(p)=dU3dt_avg(p)+U3_old(p);
    end
    
    % At the boundary (by interpolation method)
    % Inlet
    U2(1)=(2*U2(2))-U2(3);
    U3(1)=U1(1)*((T(1)/(gamma-1))+((gamma/2)*(u(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);
    
    % Final value of primitive variables
    rho=U1./A;
    u=U2./U1;
    T=(gamma-1)*((U3./U1)-((gamma/2)*(u.^2)));
    % Pressure ratio
    p=rho.*T;
    % Mach number
    M=u./(T.^0.5);
    % Mass flow rate
    mass=rho.*u.*A;
    
    % Field variables at throat section
    rho_throat(j)=rho(th);
    T_throat(j)=T(th);
    p_throat(j)=p(th);
    M_throat(j)=M(th);
    
    % Throat section residual values
    U1_res(j)=dU1dt_avg(th);
    U2_res(j)=dU2dt_avg(th);
    U3_res(j)=dU3dt_avg(th);
    
    % Plot of mass flow rate at different time instances
    plot_c_mass(x,j,mass)
   
end

% Plots at throat section
plot_c_throat2(rho_throat,p_throat,T_throat,M_throat)
plot_c_throat3(U1_res,U2_res,U3_res)

% Plots at Steady-state
plot_ss(x,M,p,rho,T)

Functions to plot the graphs -

Mass flow rate at different time instances

function[]=plot_c_mass(x,j,mass)

figure(1)
hold on
if (j==1)
    plot(x,mass)
    elseif (j==51)
        plot(x,mass,'g')
        elseif (j==101)
            plot(x,mass,'r')
            elseif (j==151)
                plot(x,mass,'c')
                elseif (j==201)
                    plot(x,mass,'y')
                    elseif (j==701)
                        plot(x,mass,'black')
                        axis([0 3 0 1.8])
                        legend('0dt','50dt','100dt','150dt','200dt','700dt','Location','best')
                        xlabel('Non-dimensional distance through nozzle')
                        ylabel('Non-dimensional mass flow rate')
                        title('Distribution of Mass flow rate in the nozzle at different time instances')
                        grid on
end
end

Field variables at throat section i.e. i=16

function[]=plot_c_throat2(rho_throat,p_throat,T_throat,M_throat)

% Exact analytical solution at throat at steady state
rho_exact=0.634*ones(1,1400);
T_exact=0.833*ones(1,1400);
p_exact=0.528*ones(1,1400);
M_exact=1*ones(1,1400);

% Plots for convergence (at throat section)
figure(2)
plot(rho_throat)
hold on 
plot(rho_exact,'r')
axis([0 1400 0.5 1])
xlabel('Number of iterations')
ylabel('Density')
legend('Numerical solution','Exact solution','Location','best')
title('Density at Throat section')

figure(3)
plot(T_throat)
hold on 
plot(T_exact,'r')
axis([0 1400 0.65 1])
xlabel('Number of iterations')
ylabel('Temperature')
legend('Numerical solution','Exact solution','Location','best')
title('Temperature at Throat section')

figure(4)
plot(p_throat)
hold on 
plot(p_exact,'r')
axis([0 1400 0 1])
xlabel('Number of iterations')
ylabel('Pressure')
legend('Numerical solution','Exact solution','Location','best')
title('Pressure at Throat section')

figure(5)
plot(M_throat)
hold on 
plot(M_exact,'r')
axis([0 1400 0.8 1.8])
xlabel('Number of iterations')
ylabel('Mach number')
legend('Numerical solution','Exact solution','Location','best')
title('Mach number at Throat section')

end

Residual at throat section i.e at i=16

function[]=plot_c_throat3(U1_res,U2_res,U3_res)

% Plots for Residual (at throat section)
figure(6)
plot(U1_res)
hold on
plot(U2_res,'g')
hold on
plot(U3_res,'r')
xlabel('Number of iterations')
ylabel('Residual')
title('Residual at Throat section')
legend('U1','U2','U3','Location','best')

end

Steady-state value of field variables- Same program as used for non-conservative form

Results for Conservative form

The steady state profile for field variables is in accordance with the plots obtained by the analytical method. It does not differ much from the profiles showed above for non-conservative form.

The above four plots show the variation in the value of field variables as the simulation advances. It is observed that the solution converges after 900 iterations for all field variables. Although the final solution does not exactly match with the exact analytical solution, the result obtained is satisfactory. But the steady-state result of non-conservative form is more closer to the actual value than that of conservative form.

The above plot shows the residual term for each variable as the number of iterations increases. The residual is nearly zero after 900 iterations. The rate of decay of residuals is better for the non-conservative form.

This plot shows the mass flow rate in the nozzle at different time instances. After 700dt, the mass flow rate is nearly constant throughout the nozzle. Mass is conserved better using the conservative form of governing equation.

Grid dependence test:

The results are checked when the number of grid points is varied in the simulation. Here, two cases are studied- 31 grid points and 61 grid points. The results obtained are compared with the exact solution. The values shown in the table are for grid point i=16, i.e. at the throat section.

From the table, it is seen that as we increase the number of grid points, the value of field variables approach towards the exact analytical solution. But this is followed by an increase in the runtime of the simulation. So for faster results, less number of grid points are used, but the deviation of the variable values in such case is more. Hence a compromise should be made between the accuracy of the solution and runtime of simulation.

Conclusions:

From the plots obtained by solving the conservative and non-conservative form of governing equations, the following conclusions are made-

(1) For the same number of grid points, non-conservative form yields better results than the conservative form.

(2) The accuracy of the solution increases as we increase the number of grid points.

(3) The mass flow distribution is better for conservative form. It simply conserves the mass better than the non-conservative form.

(4) Solving the non-conservative form is easier as it involves less mathematical complexity.

(5) The number of iterations required for convergence is less for non-conservative form.

(6) Also, the behavior of residuals for the conservative form is not as good than that of non-conservative form. The rate of decay of residuals is faster for non-conservative form.