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

Simulation of 1D Supersonic Nozzle Flow Using Macormack Method

Main Code:

%Program to solve 1D quasi linear supersonic nozzle for conservative
%and non-conservative form
clear all
close all
clc

%Number of grids
n=31;

%Length of nozzle
x=linspace(0,3,n);

nt=1400;
m=0.579*ones(n);

%Function to solve the governing equation non conservative form
tic;
[rho_nc,T_nc,P_nc,M_nc,m_nc,time_nc]=nonconservative(n,x,nt)
time_nc=toc;
sprintf(\'Timetaken=%dfs\',time_nc)

%Function to solve the governing equation conservative form
tic;
[rho_c,T_c,P_c,M_c,m_c,time_c]=conservative(n,x,nt)
time_c=toc;
sprintf(\'Timetaken=%dfs\',time_c)

%Plot
figure(4)
plot(x,m_nc,\'r\')
hold on
plot(x,m_c,\'b\')
hold on
plot(x,m,\'--\')
xlabel(\'length of nozzle\')
ylabel(\'mass flow\')
title(\'Comparision of mass flow variation\')
legend(\'Non Conservative Form\',\'Conservative Form\',\'Analytical Solution\')

Function Code for Non Conservative Form:

%Function to solve the governing equation non conservative form
function [rho,T,P,M,m,count]=nonconservative(n,x,nt)

%Discretization of grid
dx=x(2)-x(1);

%Area
A=1+(2.2*(x-1.5).^2);

%error criteria

error=1;
tol=1e-4;
m_old=0;
gamma=1.4;
count=1;
throat=(n+1)/2;
nt_grid=linspace(1,nt,nt);

%cfl number
c=0.5;
%initial profiles
rho=(1-(0.3146*x));
T=(1-(0.2314*x));
v=(0.1+(1.09*x)).*(T.^0.5);
m=0.579*ones(n);
m_initial=rho.*A.*v;

for h=1:n
   dt_x=(c*dx)/(T(h)^0.5+v(h));
end

dt=min(dt_x)

for k=1:nt
    rho_old=rho;
    v_old=v;
    T_old=T;
    
    %Predictor Step
    for i=2:n-1
        drhodx=(rho(i+1)-rho(i))/dx;
        dvdx=(v(i+1)-v(i))/dx;
        dTdx=(T(i+1)-T(i))/dx;
        dlogAdx=(log(A(i+1))-log(A(i)))/dx;
        drhodt_p(i)=-rho(i)*dvdx-rho(i)*v(i)*dlogAdx-v(i)*drhodx;
        dvdt_p(i)=-v(i)*dvdx-(1/gamma)*(dTdx+(T(i)/rho(i))*drhodx);
        dTdt_p(i)=-v(i)*dTdx-(gamma-1)*T(i)*(dvdx+(v(i)*dlogAdx));
        rho(i)=rho_old(i)+drhodt_p(i)*dt;
        v(i)=v_old(i)+dvdt_p(i)*dt;
        T(i)=T_old(i)+dTdt_p(i)*dt;
    end
    
    %corrector step
    
    for j=2:n-1
        drhodx=(rho(j)-rho(j-1))/dx;
        dvdx=(v(j)-v(j-1))/dx;
        dTdx=(T(j)-T(j-1))/dx;
        dlogAdx=(log(A(j))-log(A(j-1)))/dx;
        drhodt_c(j)=-rho(j)*dvdx-rho(j)*v(j)*dlogAdx-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)*dlogAdx));
    end
    
    drhodt=0.5*(drhodt_p+drhodt_c);
    dvdt=0.5*(dvdt_p+dvdt_c);
    dTdt=0.5*(dTdt_p+dTdt_c);
    
    for z=2:n-1
        rho(z)=rho_old(z)+drhodt(z)*dt;
        v(z)=v_old(z)+dvdt(z)*dt;
        T(z)=T_old(z)+dTdt(z)*dt;
    end
    
    v(1)=2*v(2)-v(3);
    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;
    M=v./(T.^0.5);
    m=rho.*v.*A;
    mt(k,:)=rho.*v.*A;
    rho_throat(k)=rho(throat);
    P_throat(k)=P(throat);
    T_throat(k)=T(throat);
    M_throat(k)=M(throat);
    error=max(abs(m-m_old));
    if(error>tol)
        count=count+1;
        m_old=m;
    end
    
end

figure(1)

%mass flow rate at different time steps
plot(x,m_initial,x,mt(50,:),x,mt(100,:),x,mt(200,:),x,mt(500,:),x,mt(700,:),x,m,\'--\');
title(\'mass flow rate at different time steps non conservative form\')
xlabel(\'length of nozzle\')
ylabel(\'mass flow\')
legend(\'time step=0\',\'time step=50\',\'time step=100\',\'time step=200\',\'time step=500\',\'time step=700\',\'Analytical\')

%flow variables plot

figure(2)

subplot(4,1,1)
plot(x,M,\'r\');
xlabel(\'length of nozzle\')
ylabel(\'Mach number (M)\')
title(\'flow variables-non conservative form\')

subplot(4,1,2)
plot(x,P,\'y\')
xlabel(\'length of nozzle\')
ylabel(\'Pressure\')
title(\'flow variables-non conservative form\')

subplot(4,1,3)
plot(x,rho,\'b\')
xlabel(\'length of nozzle\')
ylabel(\'Density\')
title(\'flow variables-non conservative form\')

subplot(4,1,4)
plot(x,T,\'g\')
xlabel(\'length of nozzle\')
ylabel(\'Temperature\')
title(\'flow variables-non conservative form\')

%flow variables at throat section

figure(3)

subplot(4,1,1)
plot(nt_grid,M_throat,\'r\');
xlabel(\'length of nozzle\')
ylabel(\'Mach number (M)\')
title(\'flow variables at throat-non conservative form\')

subplot(4,1,2)
plot(nt_grid,P_throat,\'y\')
xlabel(\'length of nozzle\')
ylabel(\'Pressure\')
title(\'flow variables at throat-non conservative form\')

subplot(4,1,3)
plot(nt_grid,rho_throat,\'b\')
xlabel(\'length of nozzle\')
ylabel(\'Density\')
title(\'flow variables at throat-non conservative form\')

subplot(4,1,4)
plot(nt_grid,T_throat,\'g\')
xlabel(\'length of nozzle\')
ylabel(\'Temperature\')
title(\'flow variables at throat-non conservative form\')

end

1.1) Output For Non Conservative Form:

Calculations done at,

n=31,

nt=1400,

c=0.5

1) Mass flow rate at different time step for non conservative form

Fig 1: Mass flow rate at different time steps for non conservarive form

2) Non dimensional flow variables for Non conservative form

Fig 2: Flow Variables for Non conservative form

3) Non dimensional flow variables at throat section for Non conservative form

Fig 3: Flow Variables at throat section for non conservative form

Function Code for Conservative Form:

%Function to solve the governing equation conservative form
function [rho,T,P,M,m,count]=conservative(n,x,nt)

%Discretization of grid
dx=x(2)-x(1);

%Area
A=1+(2.2*(x-1.5).^2);

%error criteria
error=1;

%tolerance
tol=1e-4;
m_old=0;
m=0.597*ones(n);
gamma=1.4;
count=1;
throat=(n+1)/2;
nt_grid=linspace(1,nt,nt);

%cfl number
c=0.5;
u2i=0.59;

%initial profiles of flow variables

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))
        rho(i)=0.634-(0.3879*(x(i)-1.5));
        T(i)=0.833-(0.35307*(x(i)-1.5));
    end
    
    v(i)=u2i/(rho(i)*A(i));
    u1(i)=rho(i)*A(i);
    u2(i)=rho(i)*A(i)*v(i);
    u3(i)=rho(i)*A(i)*((T(i)/(gamma-1))+((gamma/2)*(v(i)^2)));
    
end

%copy of initial values

u1_inital=u1;
u2_initial=u2;
u3_inital=u3;

%time step calculation

for h=1:n
    dt_x=(c*dx)/(T(h)^0.5+v(h));
end
    dt=min(dt_x);
    
    %time loop
    for k=1:nt
        u1_old=u1;
        u2_old=u2;
        u3_old=u3;
        
        %predictor step
        for i=2:n-1
            F1=u2;
            F2=((u2.^2)./u1)+((gamma-1)/gamma).*(u3-(0.5*gamma.*u2.^2)./u1);
            F3=((gamma.*u2.*u3)./u1)-(0.5*gamma*(gamma-1).*u2.^3)./(u1.^2);
            dF1dx=(F1(i+1)-F1(i))/dx;
            dF2dx=(F2(i+1)-F2(i))/dx;
            dF3dx=(F3(i+1)-F3(i))/dx;
            J2(i)=(1/gamma).*(rho(i).*T(i).*((A(i+1)-A(i))/dx));
            
            %continuity equitation
            du1dt_p(i)=-dF1dx;
            
            %momentum equation
            du2dt_p(i)=-dF2dx+J2(i);
            
            %energy equation
            du3dt_p(i)=-dF3dx;
            
            %updating values of predictor method
            u1(i)=u1(i)+du1dt_p(i)*dt;
            u2(i)=u2(i)+du2dt_p(i)*dt;
            u3(i)=u3(i)+du3dt_p(i)*dt;
            
        end
        
        %updating flow variables
        rho=u1./A;
        v=u2./u1;
        T=(gamma-1)*((u3./u1)-((0.5*gamma)*(v.^2)));
        
        %corrector step
        for j=2:n-1
            F1=u2;
            F2=((u2.^2)./u1)+((gamma-1)/gamma).*(u3-(0.5*gamma.*u2.^2)./u1);
            F3=((gamma.*u2.*u3)./u1)-(0.5*gamma*(gamma-1).*u2.^3)./(u1.^2);
            dF1dx=(F1(j)-F1(j-1))/dx;
            dF2dx=(F2(j)-F2(j-1))/dx;
            dF3dx=(F3(j)-F3(j-1))/dx;
            J2(j)=(1/gamma).*(rho(j).*T(j).*((A(j)-A(j-1))/dx));
            
            %continuity equitation
            du1dt_c(j)=-dF1dx;
            
            %momentum equation
            du2dt_c(j)=-dF2dx+J2(j);
            
            %energy equation
            du3dt_c(j)=-dF3dx;
            
        end
        
        %average values
        du1dt=0.5*(du1dt_p+du1dt_c);
        du2dt=0.5*(du2dt_p+du2dt_c);
        du3dt=0.5*(du3dt_p+du3dt_c);
        
        %updating values of predictor method
        for z=2:n-1
            u1(z)=u1_old(z)+du1dt(z)*dt;
            u2(z)=u2_old(z)+du2dt(z)*dt;
            u3(z)=u3_old(z)+du3dt(z)*dt;
        end
        
        %inlet boundary condition
        u1(1)=rho(1)*A(1);
        u2(1)=2*u2(2)-u2(3);
        u3(1)=u1(1)*(T(1)/(gamma-1)+((gamma/2)*(v(1)^2)));
        
        %outlet boundary condition
        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);
        
        %calculating the non dimensional flow variables
        
        for i=2:n-1
            rho(i)=u1(i)./A(i);
            v(i)=u2(i)./u1(i);
            T(i)=(gamma-1)*((u3(i)./u1(i))-((0.5*gamma)*(v(i).^2)));
        end
        
        %computing pressure,Mach number,mass flow
        P=rho.*T;
        M=v./(T.^0.5);
        m=rho.*v.*A;
        mt(k,:)=rho.*v.*A;
        
        %flow variables at throat
        rho_throat(k)=rho(throat);
        P_throat(k)=P(throat);
        T_throat(k)=T(throat);
        M_throat(k)=M(throat);
        error=max(abs(m-m_old));
        
        if(error>tol)
            count=count+1;
            m_old=m;
        end
        
    end
    
%plot of mass flow rate at different time steps
    
             
figure(5)

plot(x,u2_initial,x,mt(50,:),x,mt(100,:),x,mt(200,:),x,mt(500,:),x,mt(700,:),x,m,\'--\');
title(\'mass flow rate at different time steps conservative form\')
xlabel(\'length of nozzle\')
ylabel(\'mass flow\')
legend(\'time step=0\',\'time step=50\',\'time step=100\',\'time step=200\',\'time step=500\',\'time step=700\',\'Analytical\')

%flow variables plot

figure(6)

subplot(4,1,1)
plot(x,M,\'r\');
xlabel(\'length of nozzle\')
ylabel(\'Mach number (M)\')
title(\'flow variables conservative form\')

subplot(4,1,2)
plot(x,P,\'y\')
xlabel(\'length of nozzle\')
ylabel(\'Pressure\')
title(\'flow variables conservative form\')

subplot(4,1,3)
plot(x,rho,\'b\')
xlabel(\'length of nozzle\')
ylabel(\'Density\')
title(\'flow variables conservative form\')

subplot(4,1,4)
plot(x,T,\'g\')
xlabel(\'length of nozzle\')
ylabel(\'Temperature\')
title(\'flow variables conservative form\')

%fluid properties at throat section

figure(7)

subplot(4,1,1)
plot(nt_grid,M_throat,\'r\');
xlabel(\'length of nozzle\')
ylabel(\'Mach number (M)\')
title(\'flow variables at throat-conservative form\')

subplot(4,1,2)
plot(nt_grid,P_throat,\'y\')
xlabel(\'length of nozzle\')
ylabel(\'Pressure\')
title(\'flow variables at throat-conservative form\')

subplot(4,1,3)
plot(nt_grid,rho_throat,\'b\')
xlabel(\'length of nozzle\')
ylabel(\'Density\')
title(\'flow variables at throat-conservative form\')

subplot(4,1,4)
plot(nt_grid,T_throat,\'g\')
xlabel(\'length of nozzle\')
ylabel(\'Temperature\')
title(\'flow variables at throat-conservative form\')
end

1.2) Output for Conservative Form

Calculations done at,

n=31,

nt=1400,

1) Mass flow rate at different steps

Fig 1: Mass flow rate at different steps for conservative form

2) Non dimensional flow variables for Conservative form

Fig 2: Flow variables for conservative form

3) Non dimensional flow variables at throat section for Conservative form

Fig 3: Flow variables at throat section for Conservative form

4) Comparision of mass flow variations for both forms and analytical solution

Fig 4: Comparision of mass flow variation

Output results after changing number of grid points

n=61,

nt=1400,

c=0.5

1) Output for Non Conservative form

a) Mass flow rate at different time steps for non conservative form

Fig 1: Mass flow rate at different time steps for non conservative form

b) Non dimensional flow variables for non conservative form

Fig 2: Non dimensional flow variables for non conservative form

c) Non dimensional flow variables at throat section for non conservative form

Fig 3: Non dimensional flow variables at throat section for non conservative form

2) Output results for Conservative form

a) Mass flow rate at different time steps for conservative form

Fig 1:Mass flow rate at different time steps for conservative form

b) Non dimensional flow variables for conservative form

Fig 2: Non dimensional flow variables for conservative form

c) Non dimensional flow variables at throat section for conservative form

Fig 3: Non dimensional flow variables at throat section for conservative form

*Comparision of mass flow vairiations for both forms and analytical solution

Fig : Comparision of mass flow variation

Non conservative form:

* From the above plots, it is clear that mass flow rate at time steps 500, 700 and

analytical solution are nearly equal.

* There are few fluctuations at beginning but the solution remains stable

*After 500 time steps the solution becomes stable and is nearer to the exact analytical

solution.

* At the initial time steps (0,50dt,100dt) there are oscillations in value due to the large

values of time derivatives.

* As time steps increases there are no fluctuation and solution gets stable.

Conservative form:

* In this case, there are lesser fluctuations in the solution due to the better initial

conditions applied compared to non conservative form.

* After 500 time steps solution is stable and residue is lower.

* After 700 time steps it is constant and nearly equal to the exact analytical solution

(mass flow rate=0.579) and steady state is achived.

Comparision of conservative and non conservative form:

* The conservative form yields a better mass flow rate with less fluctuation compared to

the non conservative form. The non conservative form generates oscillation at initial time

steps due to the residuals. In conservative form mass flow is one of the dependant

variables whereas in non conservative form the dependants are primitive variables like

pressure, Mach number etc.

Grid dependancy:

* Grid dependance or number of grid points is a important factor in obtaining a stable

numerical solution.

* Grid independent solution is obtained when the solution is not changing as a function

of grid points.

* At 61 grid points we have stable solution and the result improves marginally with

analytical solution, hence the solution is grid independant.

Courant number:

*Courant number is another important factor which determines the stability of given

numerical solution.

* In this problem courant number is taken as c=0.5 since for a hyperbolic 1D linear

equation C<=1.

* Though the equation is non linear in nature the assumption of courant number c=0.5

is good as we obtain a stable solution.

Faster Solution:

* Comparing both the forms non conservative equations takes lesser execuation time

than conservative form.

* As the grid points are doubled time taken for non conservative equation increases as

more calculations involved.

Execution time comparision for both forms: