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

Simulation of a 1D supersonic nozzle flow using Maccormack method

The flow at the inlet to the nozzle arrives from the reservoir where the pressure and temperature are denoted by Po and Co. The cross-sectional area of the reservoir is very large and hence the fluid velocity is very low. The flow expands from subsonic to supersonic flow, whereas the velocity of the fluid at the throat is equal to the sonic velocity.

Governing equation

The assumptions made for solving nozzle flow problems is that the flow properties varies only in the x-direction and remains uniform across any given cross-section. This type of problem is called quasi 1D flow problem. The integral form of the governing equations are

Continuity equation

Momentum equation

Energy equations

The integral form of governing equations should be converted into corresponding PDE form such that the PDE can be numerically solved to obtain the primitive variables such as velocity, density, temperature, pressure and so on. A control volume of negligible dimensions(tending to zero) is selected and the above integral form of governing equation is applied to CV and solved to obtain both conservative and non-conservative form of PDE.

Conservative form of governing equation

Continuity equation

Momentum equation

Energy equation

Non-Conservative form of Governing equation

Continuity equation

Momentum equation

Energy equation

, , , , ,,

The prime in the above equations represents the non-dimensional form of governing equation. Non-dimensional form of Navier stokes equation can ease the analysis of the problem and it reduces the number of parameters.

The difference between the conservative and non-conservative form of governing equation is that both forms give the same analytical solution but the different numerical output. Conservative form provides accurate results and it can provide better results in the presence of shocks or discontinuities.

Steps for solving the Nozzle flow problem

Meshing

The given 1D domain is divided into n number of elements. For Nozzle flow problem the length of nozzle is chosen as 3 and the length is divided into 30 elements.

Discretization

The PDE is discretized and solved for the points in the domain. Maccormack method is used for discretizing the equation. Maccormack method gives better results for the hyperbolic equation and it is also called as a shock-capturing method. These methods can provide stable solutions in the presence of shock. Hence they are used in nozzle flow problems. In Maccormack method combination of both forward and backward differencing schemes can be used to solve for the primitive variables. Non-conservative and conservative form of governing equations are solved by using this method

The time step used for the simulation is

The minimum value of is used such that the time marching solution is following the actual unsteady flow variation that exists in nature.

Steps involved in Maccormack method

• Predictor step(forward differencing)
• Solution update
• Corrector step( Backward differencing)
• Determining the average of the time derivative
• Final solution update

Non-conservative form of Governing equation

In the non-conservative form of governing equation the continuity, momentum and energy equations can be discretized directly to obtain the primitive variables such as density, velocity and temperature.

The initial and boundary conditions are

Conservative form of governing equation

In the conservative form of governing equation, the primitive variables are grouped and they are represented in terms of solution vector U, flux vector F and source term J. The resulting equation is solved to obtain updated values of the density, velocity, temperature and so on

Now the governing equation can be written as

To obtain the primitive variables the solution vector U1, U2, U3 must be decoded.

It is better to write F1, F2, F3 in terms of U1, U2, U3 to avoid the occurrence of instabilities during the time marching solution.

The initial and boundary condition for the conservative form of governing equations are

Initial conditions

Matlab code

Main Program

%%Main Program
clear all
close all
clc
%inputs
n=31
x=linspace(0,3,n);
dx=x(2)-x(1);
a=1+(2.2*(x-1.5).^2)
gamma=1.4
nt=1400
%% calculation
% non conservative form
tic
[v_nc,t_nc,rho_nc,mass_flow_rate_nc,mach_nc,p_nc,t_nct,rho_nct,mach_nct,p_nct,mass_flow_rate_nct]=non_conservative_form(x,dx,n,a,gamma,nt)
simu_time_non_conservative=toc
%% conservative form
tic
[v_c,t_c,rho_c,mass_flow_rate_c,mach_c,p_c,t_ct,rho_ct,mach_ct,p_ct,mass_flow_rate_ct]=conservativeform(x,dx,n,a,gamma,nt)
simu_time_conservative=toc
%% Plotting Steady state variation
%Temperature
subplot(4,1,1)
grid on
title('Steady state distribution of primitive variables')
plot(x,t_nc,'k')
hold on
plot(x,t_c,'*')
xlabel('Length of Nozzle')
ylabel('Temperature')
hold on
%density
subplot(4,1,2)
grid on
plot(x,rho_nc,'k')
hold on
plot(x,rho_c,'*')
xlabel('Length of Nozzle')
ylabel('density')
hold on
%mach number
subplot(4,1,3)
grid on
plot(x,mach_nc,'k')
hold on
plot(x,mach_c,'*')
xlabel('Length of Nozzle')
ylabel('mach number')
hold on
%Pressure
subplot(4,1,4)
grid on
plot(x,p_nc,'k')
hold on
plot(x,p_c,'*')
xlabel('Length of Nozzle')
ylabel('Pressure')
legend('non conservative form','conservative form')
hold on       
%% Plotting transient variation
time=1:nt
%Temperature
subplot(4,1,1)
grid on
title('Timewise distribution of primitive variables at Throat')
plot(time,t_nct,'k')
hold on
plot(time,t_ct,'r')
xlabel('Time step')
ylabel('Temperature')
hold on
%%density
subplot(4,1,2)
grid on
title('Timewise distribution of primitive variables at Throat')
plot(time,rho_nct,'k')
hold on
plot(time,rho_ct,'r')
xlabel('Time step')
ylabel('Density ')
hold on
%mach number
subplot(4,1,3)
grid on
title('Timewise distribution of primitive variables at Throat')
plot(time,mach_nct,'k')
hold on
plot(time,mach_ct,'r')
xlabel('Time step')
ylabel('mach number')
hold on
%Pressure
subplot(4,1,4)
grid on
title('Timewise distribution of primitive variables at Throat')
plot(time,p_nct,'k')
hold on
plot(time,p_ct,'r')
xlabel('Time step')
ylabel('Pressure')
hold on
legend('non conservative form','conservative form')
grid on
%% Plotting Normalized mass flow rate
title('Mass flow rate for both conservative and non conservative form')
plot(x,mass_flow_rate_nc,'k')
hold on
plot(x,mass_flow_rate_c,'r')
hold on
grid on
legend('non conservative form','conservative form')
xlabel('Length of nozzle')
ylabel('Mass flow rate')
%% Mass flow rate at throat different time steps
title('Mass flow rate at throat for different time steps') 
plot(time,mass_flow_rate_nct,'k')
hold on
plot(time,mass_flow_rate_ct,'r')
hold on
grid on
legend('non conservative form','conservative form')
%% exact vs analytical
rho_e=[0.995,0.994,0.992,0.990,0.987,0.983,0.978,0.970,0.958,0.942,0.920,0.888,0.844,0.787,0.716,0.634,0.547,0.461,0.382,0.315,0.258,0.213,0.176,0.147,0.124,0.105,0.090,0.078,0.068,0.059,0.052]
plot(x,rho_e,'*')
hold on
plot(x,rho_c,'k')
hold on
plot(x,rho_nc,'r')
xlabel('length of nozzle')
ylabel('Density')
grid on
legend('Exact','conservative','nonconservative')
title('Exact vs Numerical values for density')

Conservative form

function [v_c,t_c,rho_c,mass_flow_rate_c,mach_c,p_c,t_ct,rho_ct,mach_ct,p_ct,mass_flow_rate_ct]=conservativeform(x,dx,n,a,gamma,nt,dt)
for w=1:n
    if x(w)>=0 && x(w)<=0.5
        rho_c(w)=1
        t_c(w)=1
      elseif x(w)>=0.5 && x(w)<=1.5
        rho_c(w)=1-(0.366*(x(w)-0.5))
        t_c(w)=1-(0.167*(x(w)-0.5))
      elseif x(w)>=1.5 && x(w)<=3
        rho_c(w)=0.634-(0.3879*(x(w)-1.5))
        t_c(w)=0.833-(0.3507*(x(w)-1.5))
    end
end
%time step
c=0.5
throat=find(a==1)
% Properties
p_c=rho_c.*t_c
v_c=0.59./(rho_c.*a)
dt= min(c*(dx./((t_c.^0.5)+v_c)))
mass_flow_rate_c=rho_c.*a.*v_c
u1=rho_c.*a
u2= u1.*v_c
u3=((t_c./(gamma-1))+(gamma*0.5.*(v_c.^2))).*u1
mach_c=v_c./(t_c.^0.5)
% initialization
p_init=p_c
v_init=v_c
mass_flow_rate_init=mass_flow_rate_c
u1_init=u1
u2_init=u2
u3_init=u3
mach_init=mach_c
for d=1:nt
% storing values
    u1_old=u1;
    u2_old=u2;
    u3_old=u3;
    F1=u2;
    F2=((u2.^2)./u1)+(((gamma-1)/gamma)*(u3-(((u2.^2)./u1)*0.5*gamma)))
    F3=((gamma.*u2.*u3)./u1)-(0.5*gamma*(gamma-1).*((u2.^3)./(u1.^2)))
for j=2:n-1
    dadx(j)=(log(a(j+1))-log(a(j)))/(dx)
    J(j)=(((gamma-1)/gamma)*(u3(j)-(((u2(j)^2)/u1(j))*0.5*gamma)))*dadx(j)
    % Preditor step
    du1_dt_p(j)=-((F1(j+1)-F1(j))/dx)
    du2_dt_p(j)=-((F2(j+1)-F2(j))/dx)+J(j)
    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
%properties update
v_c=u2./u1
t_c=(gamma-1)*((u3./u1)-(gamma*0.5.*v_c.^2))
rho_c=u1./a
p_c=rho_c.*t_c
F1=u2;
F2=((u2.^2)./u1)+(((gamma-1)/gamma)*(u3-(((u2.^2)./u1)*0.5*gamma)))
F3=((gamma.*u2.*u3)./u1)-(0.5*gamma*(gamma-1).*((u2.^3)./(u1.^2)))
for j=2:n-1
    dadx(j)=(log(a(j))-log(a(j-1)))/(dx)
    J(j)=(((gamma-1)/gamma)*(u3(j)-(((u2(j)^2)/u1(j))*0.5*gamma)))*dadx(j)
    % corrector step
    du1_dt_c(j)=-((F1(j)-F1(j-1))/dx)
    du2_dt_c(j)=-((F2(j)-F2(j-1))/dx)+J(j)
    du3_dt_c(j)=-((F3(j)-F3(j-1))/dx)
    du1_dt(j)=0.5*(du1_dt_p(j)+du1_dt_c(j))
    du2_dt(j)=0.5*(du2_dt_p(j)+du2_dt_c(j))
    du3_dt(j)=0.5*(du3_dt_p(j)+du3_dt_c(j))
    u1(j)=u1_old(j)+(du1_dt(j)*dt)
    u2(j)=u2_old(j)+(du2_dt(j)*dt)
    u3(j)=u3_old(j)+(du3_dt(j)*dt)
end
u1(1)=rho_c(1)*a(1)
u2(1)=2*u2(2)-u2(3)
u3(1)=((t_c(1)/(gamma-1))+(gamma*0.5*(v_c(1)^2)))*u1(1)
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)
%properties update
v_c=u2./u1
t_c=(gamma-1)*((u3./u1)-(gamma*0.5.*v_c.^2))
rho_c=u1./a
p_c=rho_c.*t_c
mass_flow_rate_c=rho_c.*a.*v_c
mach_c=v_c./(t_c.^0.5)
% Throat values
t_ct=t_c(throat)
rho_ct=rho_c(throat)
mach_ct=mach_c(throat)
p_ct=p_c(throat)
mass_flow_rate_ct=mass_flow_rate_c(throat)
%% plotting mass flow rate inside the nozzle for different time step
     if d==50
     plot(x,mass_flow_rate_c,'k')
     hold on
     elseif d==100
     plot(x,mass_flow_rate_c,'*')
     hold on
     elseif d==200
     plot(x,mass_flow_rate_c,'b')
     hold on
     elseif d==600
     plot(x,mass_flow_rate_c,'g')
     hold on
     elseif d==800
     plot(x,mass_flow_rate_c,'o')
     hold on
     elseif d==1400
     plot(x,mass_flow_rate_c,'r')
     hold on
     grid on
     title('Variation of mass flow rate along the nozzle for different time step-conservative form')
     legend('50','100','200','600','800','1400')
     xlabel('length of nozzle')
     ylabel('mass flow rate')
 end

end
end

Non Conservative form

function [v_nc,t_nc,rho_nc,mass_flow_rate_nc,mach_nc,p_nc,t_nct,rho_nct,mach_nct,p_nct,mass_flow_rate_nct]=non_conservative_form(x,dx,n,a,gamma,nt,dt)
throat=find(a==1)
rho_nc=1-0.3146*x
t_nc=1-0.2314*x
v_nc=(0.1+1.09*x).*t_nc.^0.5
%Time step
c=0.5
dt= min(c*(dx./((t_nc.^0.5)+v_nc)))

%% calculating primitive variables by Macormack method
for d=1:nt
rho_old_nc=rho_nc
t_old_nc=t_nc
v_old_nc=v_nc
    for j=2:n-1
        %% Preditor Step (Forward differencing)
        dvdx=(v_nc(j+1)-v_nc(j))/dx
        dtdx=(t_nc(j+1)-t_nc(j))/dx
        drhodx=(rho_nc(j+1)-rho_nc(j))/dx
        dadx=(log(a(j+1))-log(a(j)))/dx
        % Continuity equation
        drhodt_p(j)=-(rho_nc(j)*dvdx)-(rho_nc(j)*v_nc(j)*dadx)-(v_nc(j)*drhodx)
        % Momentum equation
        dvdt_p(j)=-(v_nc(j)*dvdx)-((1/gamma)*(dtdx+(t_nc(j)/rho_nc(j))*drhodx))
        % Energy equation
        dtdt_p(j)=-(v_nc(j)*dtdx)-((gamma-1)*t_nc(j)*(dvdx+(v_nc(j)*dadx)))
    
    v_nc(j)=v_nc(j)+dvdt_p(j)*dt
    t_nc(j)=t_nc(j)+dtdt_p(j)*dt
    rho_nc(j)=rho_nc(j)+drhodt_p(j)*dt
    end
     %% corrector step
    for j=2:n-1
         dvdx=(v_nc(j)-v_nc(j-1))/dx
        dtdx=(t_nc(j)-t_nc(j-1))/dx
        drhodx=(rho_nc(j)-rho_nc(j-1))/dx
        dadx=(log(a(j))-log(a(j-1)))/dx
        % Continuity equation
        drhodt_c(j)=-(rho_nc(j)*dvdx)-(rho_nc(j)*v_nc(j)*dadx)-(v_nc(j)*drhodx)
        % Momentum equation
        dvdt_c(j)=-(v_nc(j)*dvdx)-((1/gamma)*(dtdx+(t_nc(j)/rho_nc(j))*drhodx))
        % Energy equation
        dtdt_c(j)=-(v_nc(j)*dtdx)-((gamma-1)*t_nc(j)*(dvdx+(v_nc(j)*dadx)))
    end
    drhodt=0.5*(drhodt_p+drhodt_c)
    dvdt=0.5*(dvdt_p+dvdt_c)
    dtdt=0.5*(dtdt_p+dtdt_c)
 for j=2:n-1 
       v_nc(j)=v_old_nc(j)+dvdt(j)*dt
    t_nc(j)=t_old_nc(j)+dtdt(j)*dt
    rho_nc(j)=rho_old_nc(j)+drhodt(j)*dt
 end
 % starting point
 v_nc(1)=2*v_nc(2)-v_nc(3)
 t_nc(1)=1
 rho_nc(1)=1
 % end point
 v_nc(n)=2*v_nc(n-1)-v_nc(n-2)
 t_nc(n)=2*t_nc(n-1)-t_nc(n-2)
 rho_nc(n)=2*rho_nc(n-1)-rho_nc(n-2)
mass_flow_rate_nc=rho_nc.*a.*v_nc
mach_nc=v_nc./(t_nc.^0.5)
p_nc=rho_nc.*t_nc

%% Throat values
t_nct=t_nc(throat)
rho_nct=rho_nc(throat)
mach_nct=mach_nc(throat)
p_nct=p_nc(throat)
mass_flow_rate_nct=mass_flow_rate_nc(throat)

%% plotting mass flow rate inside the nozzle for different time step
     if d==50
     plot(x,mass_flow_rate_nc,'k')
     hold on
     elseif d==100
     plot(x,mass_flow_rate_nc,'*')
     hold on
     elseif d==200
     plot(x,mass_flow_rate_nc,'b')
     hold on
     elseif d==600
     plot(x,mass_flow_rate_nc,'g')
     hold on
     elseif d==800
     plot(x,mass_flow_rate_nc,'o')
     hold on
     elseif d==1400
     plot(x,mass_flow_rate_nc,'r')
     hold on
     grid on
     title('Variation of mass flow rate along the nozzle for different time step-Non conservative form')
     legend('50','100','200','600','800','1400')
     xlabel('length of nozzle')
     ylabel('mass flow rate')
 end
 end
 end

Results and discussion

1) Steady-state distribution of primitive variables inside the nozzle

From the above graph, it can be concluded Mach number increases with an increase in the value of x. At the throat, the value of the Mach number is equal to 1. And there is a gradual decrease in the temperature, pressure and density. Both conservative and non-conservative form gives similar results.

At timestep=1400 the numerical solution matches with the exact solution with the least amount of error.

2)Time vice variation of primitive variables

The above graph shows the timewise distribution of primitive variables at the throat

3) Variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process

Above graph represents the variation of mass flow rate along the length of the nozzle at different time steps. The value of mass flow rate obtained from both conservative and non-conservative form shows consistency at the time step of 600 and above.

The total time taken for solving the conservative form is more when compared to non-conservative form.

4)Comparison of Normalized mass flow rate distributions of both forms

It can be observed from the graph that the mass flow rate is preserved and shows stability for conservative form rather than the non-conservative form. This is because the mass flow rate is the primitive variable for the conservative form. But in non-conservative form, the primitive variables are density, velocity and so on. The mass flow rate is calculated separately for non-conservative form.

5) Grid independence test

From the table, it can be concluded that as the number of grid points increases the solution accuracy increases and the non-conservative form show better results than the conservative form.

Conclusion:

The non-conservative form produces accurate results for the primitive variable and the conservative form gives better results for the flux variable.