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

Aim:

To simulate isentropic flow through a quasi 1D subsonic-supersonic nozzle by deriving both the conservative and non-conservative forms of governing equations using MacCormack's method.

1. To plot all the necessary graphs for both the conservative and non conservative forms.

2. Time wise variation of primitive variables.

3. Variation of mass flow rate at different time steps.

4. Comparission of mass flow rate in both the forms.

5. Grid dependency test on the solution.

Theory :

MacCormack method:

This method is a explicit finite difference technique with second order accurate in both space and time. It is a student friendly and easy to program for a hyperbolic PDE. This method is a multi stage method comprising of predictor and corrector steps mainly. The result is averaged so that the solution converges to the accurate value. This method is highly prefferd because of combination of both FDS and BDS and finding its average which makes it more accurate.

Method:

We considered a steady ienstropic flow through a convergent divergent nozzle. The flow at the inlet to the nozzle comes from a reservoir. Since the area of cross section of reservoir is very large we can assume our velocity to be approximately zero. The flow is generaly a 2D flow but we are assuming that there is no flow in along Y direction. Here the flow is subsonic till the throat and later it is enhanced to a supersonic properties. The inlet properties are stagnent properties.

The governing equations:

Conservative form:

1. Continuity equation

2. Momentum equation

3. Energy equation

All the terms here are non-dimentionalized with the stangnant properties.

Now we have to setup the above equations using finite difference expressions as per MacCormacks method.

Initially the equations are divided into solution vector U, flux vector F and source term J.

Hence the primitive variables are defined as

1.First step is predictor time step which will be solved using forward differencing scheme for Flux vector

2. Update solution

3. Corrector step

Here we have to calculate new flux terms and angain repeat the predictor step with backward differencing scheme.

4. Average step

We calculate the average of predictor and corrector step.

5. Solution for correct values at t+δt timme step.

6. we calculate the primitive values at the end which provide real quantities of pressure machnumber temperature and density.

Non-conservative form:

1. Continuity equation

2. Momentum equation

3. Energy equation

All the terms here are non-dimentionalized with the stangnant properties.

Now we have to setup the above equations using finite difference expressions as per MacCormacks method.

1.First step is predictor time step which will be solved useing forward differencing scheme

2. Now solution update

3. Corrector step is same as the predictor step but here backward differencing is carried out.

4. Average of predictor and corrector step

5. Final corrected values for the t+δt time step

Initial conditions:

1. Conservative form:

inlet conditions

outlet condiitions

2. Non-conservative form

inlet conditions

v1 = 2*v2 - v3

outlet conditions

**********************************************MAIN PROGRAM*********************************************

close all
clc


%mesh
n = 61;
x = linspace(0,3,n);
dx = x(2)-x(1);


gamma = 1.4; 
cv = 0.718;

nt = 1400;
A = 1+2.2*(x - 1.5).^2;
throatc = find(A==1);

[rho_nonconserve,T_nonconserve,v_nonconserve,p_nonconserve,M_nonconserve,m_nonconserve,rho_thrnc,T_thrnc,p_thrnc,M_thrnc] = non_conservative_form(n,x,dx,gamma,nt,A);
[rho_conserve,T_conserve,v_conserve,p_conserve,m_conserve,M_conserve,rho_thrc,T_thrc,p_thrc,M_thrc]=conservative_form(n,x,dx,gamma,nt,A,throatc);

rho_v_A_nc = rho_nonconserve.*v_nonconserve.*A;
rho_v_A_c = rho_conserve.*v_conserve.*A;

figure(5)
title('Non dimensional quantities of conservative form n = 61')
subplot(4,1,1)
plot(1:nt,T_thrc,'b');
xlabel('Along the nozzle');
ylabel('Temperature');
legend('Temperature');
subplot(4,1,2)
plot(1:nt,p_thrc,'g');
xlabel('Along the nozzle');
ylabel('Pressure');
legend('Pressure');
subplot(4,1,3)
plot(1:nt,rho_thrc,'r');
xlabel('Along the nozzle');
ylabel('Density');
legend('Density');
subplot(4,1,4)
plot(1:nt,M_thrc,'k');
xlabel('Along the nozzle');
ylabel('Mach No');
legend('Mach no');


figure(6)
title('Non dimensional quantities of non conservative form n = 61')
subplot(4,1,1)
plot(1:nt,T_thrnc,'b');
xlabel('Along the nozzle');
ylabel('Temperature');
legend('Temperature');
subplot(4,1,2)
plot(1:nt,p_thrnc,'g');
xlabel('Along the nozzle');
ylabel('Pressure');
legend('Pressure');
subplot(4,1,3)
plot(1:nt,rho_thrnc,'r');
xlabel('Along the nozzle');
ylabel('Density');
legend('Density');
subplot(4,1,4)
plot(1:nt,M_thrnc,'k');
xlabel('Along the nozzle');
ylabel('Mach No');
legend('Mach no');


figure(7)
plot(x,rho_v_A_nc,'r');
hold on
plot(x,rho_v_A_c,'b');
title('Mass flow rate comparission for n=61')
xlabel('Along the nozzle');
ylabel('Mass flow rate');
legend('mass flow rate conservative form','mass flow rate non conservative');

Non-conservative form function

function [rho_nonconserve,T_nonconserve,v_nonconserve,p_nonconserve,M_nonconserve,m_nonconserve,rho_thrnc,T_thrnc,p_thrnc,M_thrnc] = non_conservative_form(n,x,dx,gamma,nt,A)


%non-conservative form
%initial profiles
rho_nonconserve = 1 - 0.3146*x;
T_nonconserve = 1 - 0.2314*x;
v_nonconserve = (0.1 + 1.09*x).*T_nonconserve.^0.5;
throatnc=find(A==1);
C = 0.5;


for k = 1:nt
    rho_old = rho_nonconserve;
    T_old = T_nonconserve;
    v_old = v_nonconserve;
    dt = min(C*dx./(sqrt(T_nonconserve)+v_nonconserve));
   %predictor method
   for i = 2:n-1
      dvdx = (v_nonconserve(i+1)-v_nonconserve(i))/dx;
      drhodx = (rho_nonconserve(i+1)-rho_nonconserve(i))/dx;
      dTdx = (T_nonconserve(i+1)-T_nonconserve(i))/dx;
      dlogAdx = (log(A(i+1))-log(A(i)))/dx;
      
      drhodt_p(i) = -(rho_nonconserve(i)*dvdx) - (rho_nonconserve(i)*v_nonconserve(i)*dlogAdx) - (v_nonconserve(i)*drhodx);
    
      dvdt_p(i) = -v_nonconserve(i)*dvdx - (1/gamma)*(dTdx + (T_nonconserve/rho_nonconserve)*(drhodx));
       
      dTdt_p(i) = -v_nonconserve(i)*dTdx - (gamma-1)*(T_nonconserve(i))*(dvdx + v_nonconserve(i)*dlogAdx);
      
      rho_nonconserve(i) = rho_nonconserve(i) +drhodt_p(i)*dt;
      v_nonconserve(i) = v_nonconserve(i) +dvdt_p(i)*dt;
      T_nonconserve(i) = T_nonconserve(i) +dTdt_p(i)*dt;
   end
    %corrector method
    for i = 2:n-1
      dvdx = (v_nonconserve(i)-v_nonconserve(i-1))/dx;
      drhodx = (rho_nonconserve(i)-rho_nonconserve(i-1))/dx;
      dTdx = (T_nonconserve(i)-T_nonconserve(i-1))/dx;
      dlogAdx_c = (log(A(i))-log(A(i-1)))/dx;
      
      drhodt_c(i) = -rho_nonconserve(i)*dvdx - rho_nonconserve(i)*v_nonconserve(i)*dlogAdx_c - v_nonconserve(i)*drhodx;
    
      dvdt_c(i) = -v_nonconserve(i)*dvdx - (1/gamma)*(dTdx + (T_nonconserve/rho_nonconserve)*(drhodx));
       
      dTdt_c(i) = -v_nonconserve(i)*dTdx - (gamma-1)*(T_nonconserve(i))*(dvdx + v_nonconserve(i)*dlogAdx_c);
      
    end
   
    drhodt = 0.5*(drhodt_p + drhodt_c);
    dvdt = 0.5*(dvdt_p + dvdt_c);
    dTdt = 0.5*(dTdt_p + dTdt_c);
    
    for i = 2:n-1
        rho_nonconserve(i) = rho_old(i) + drhodt(i)*dt;
        v_nonconserve(i) = v_old(i) + dvdt(i)*dt;
        T_nonconserve(i) = T_old(i) + dTdt(i)*dt;
    end
    
    
    %inlet
    v_nonconserve(1) = 2*v_nonconserve(2) - v_nonconserve(3);
    
    %outlet
    rho_nonconserve(n) = 2*rho_nonconserve(n-1) - rho_nonconserve(n-2);
    v_nonconserve(n) = 2*v_nonconserve(n-1) - v_nonconserve(n-2);
    T_nonconserve(n) = 2*T_nonconserve(n-1) - T_nonconserve(n-2);
    
    m_nonconserve = rho_nonconserve.*v_nonconserve.*A;
    M_nonconserve = v_nonconserve./sqrt(T_nonconserve);
    p_nonconserve = rho_nonconserve.*T_nonconserve;
    
    rho_thrnc(k) = rho_nonconserve(throatnc);
    T_thrnc(k) = T_nonconserve(throatnc); 
    p_thrnc(k) = p_nonconserve(throatnc);
    M_thrnc(k) = M_nonconserve(throatnc);
    
    
        figure(1)
        title('Mass flow rate at various time steps in non_conservative form n = 61')
    if k == 50
        plot(x,m_nonconserve)
        hold on
    elseif k==100
        plot(x,m_nonconserve)
    elseif k==150
        plot(x,m_nonconserve)
    elseif k==200
        plot(x,m_nonconserve)
    elseif k==400
        plot(x,m_nonconserve)
    elseif k==700
        plot(x,m_nonconserve)
    elseif k==1400
        plot(x,m_nonconserve)
        legend('k=50','k=100','k=150','k=200','k=400','k=700','k=1400');
        xlabel('No of iterations');
        ylabel('Mass flow rate');
    end

end

figure(2)
title('Non conservative form')
subplot(4,1,1)
plot(x,v_nonconserve);
xlabel('Along length of nozzle');
ylabel('Velocity');
subplot(4,1,2)
plot(x,M_nonconserve);
xlabel('Along length of nozzle');
ylabel('Mach number');
subplot(4,1,3)
plot(x,T_nonconserve)
xlabel('Along length of nozzle');
ylabel('Temperature');
subplot(4,1,4)
plot(x,rho_nonconserve)
xlabel('Along length of nozzle');
ylabel('Density');
end

Conservative form function

function [rho_conserve,T_conserve,v_conserve,p_conserve,m_conserve,M_conserve,rho_thrc,T_thrc,p_thrc,M_thrc]=conservative_form(n,x,dx,gamma,nt,A,throatc)
%conservative form
%initial variables
for i = 1:n
    if (x(i)>=0 && x(i)<=0.5)
        rho_conserve(i) = 1;
        T_conserve(i) = 1;
    elseif (x(i)>=0.5 && x(i)<=1.5)
        rho_conserve(i) = 1-0.366*(x(i)-0.5);
        T_conserve(i) = 1-0.167*(x(i)-0.5);
    elseif (x(i)>=1.5 && x(i)<=3.5)
        rho_conserve(i) = 0.634-0.3879*(x(i)-1.5);
        T_conserve(i) = 0.833-0.3507*(x(i)-1.5);
    end
end

u1 = rho_conserve.*A;
v_conserve = 0.59./u1;
u2 = u1.*v_conserve;
u3 = u1.*((T_conserve./(gamma-1)) + (gamma*0.5)*(v_conserve.^2));

C = 0.5;

for k = 1:nt
    u1_old = u1;
    u2_old = u2;
    u3_old = u3;
    
    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)));
    dt = min(C.*dx./(sqrt(T_conserve)+v_conserve));
    %predictor step
    for i = 2:n-1
        
        
        df1dx(i) = (f1(i+1) - f1(i))/dx;
        df2dx(i) = (f2(i+1) - f2(i))/dx;
        df3dx(i) = (f3(i+1) - f3(i))/dx;
        
        dadx(i) = (A(i+1)-A(i))/dx;
        j2(i) = (1/gamma)*rho_conserve(i)*T_conserve(i)*dadx(i);
        du1dt_p(i) = -df1dx(i);
        du2dt_p(i) = -df2dx(i) + j2(i);
        du3dt_p(i) = -df3dx(i);
        %solution update
        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
    
    rho_conserve = u1./A;
    v_conserve = u2./u1;
    T_conserve = (gamma - 1)*((u3./u1) - ((gamma*(v_conserve.^2))/2));
    
    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 i = 2:n-1
        
        
        df1dx(i) = (f1(i) - f1(i-1))/dx;
        df2dx(i) = (f2(i) - f2(i-1))/dx;
        df3dx(i) = (f3(i) - f3(i-1))/dx;
        
        dadx(i) = (A(i)-A(i-1))/dx;
        j2(i) = (1/gamma)*rho_conserve(i)*T_conserve(i)*dadx(i);
        
        du1dt_c(i) = -df1dx(i);
        du2dt_c(i) = -df2dx(i) + j2(i);
        du3dt_c(i) = -df3dx(i);
    end
    
    du1dt=0.5*(du1dt_p+du1dt_c);
    du2dt=0.5*(du2dt_p+du2dt_c);
    du3dt=0.5*(du3dt_p+du3dt_c);
    
    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
        
    %boundary conditions
    %inlet
    u1(1) = rho_conserve(1)*A(1);
    u2(1) = 2*u2(2) - u2(3);
    u3(1) = u1(1)*((T_conserve(1)/(gamma-1)) + (gamma*0.5)*(v_conserve(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_conserve = u1./A;
    v_conserve = u2./u1;
    T_conserve = (gamma - 1)*((u3./u1) - ((gamma*(v_conserve.^2))/2));
    m_conserve = rho_conserve.*v_conserve.*A;
    M_conserve = v_conserve./sqrt(T_conserve);
    p_conserve = rho_conserve.*T_conserve;

    rho_thrc(k) = rho_conserve(throatc);
    T_thrc(k) = T_conserve(throatc); 
    p_thrc(k) = p_conserve(throatc);
    M_thrc(k) = M_conserve(throatc);
            figure(3)
            title('Mass flow rate at various time steps in conservative form for n = 61')
    if k == 100
        plot(x,m_conserve)
        hold on
    elseif k==200
        plot(x,m_conserve)
    elseif k==400
        plot(x,m_conserve)
    elseif k==700
        plot(x,m_conserve)
    elseif k==1400
        plot(x,m_conserve)    
    legend('k=100','k=200','k=400','k=700','k=1400');
    xlabel('No of iterations');
    ylabel('Mass flow rate');
    end
    
end

figure(4)
title('Conservative form');
subplot(4,1,1)
plot(x,v_conserve);
xlabel('Along length of nozzle');
ylabel('Velocity');
subplot(4,1,2)
plot(x,M_conserve);
xlabel('Along length of nozzle');
ylabel('Mach number');
subplot(4,1,3)
plot(x,T_conserve)
xlabel('Along length of nozzle');
ylabel('Temperature');
subplot(4,1,4)
plot(x,rho_conserve)
xlabel('Along length of nozzle');
ylabel('Density');

end

Resullts obtained:

1.Primitive variables

• Non - conservative form result for n = 31

• Conservative form results for n = 31

Here we can clearly observe that both conservation and non-conservation method have same plot with minute difference but the trend and the values are approximately same. We can conclude that by the end of 1400 iterations both the forms provide similar results which is very close to exact solution.

2. Convergence of solution to true value at the throat

• Non-conservative form n = 31

• Non-conservative form n = 61

• Conservative form n = 31

• Conservative form n = 61

Grid dependency test results:

		Temperature	Pressure	Density	Mach no
n = 31	Non-con	0.8488	0.5625	0.6627	0.9484
n = 61	Non-con	0.8476	0.5603	0.6612	0.9491
n = 31	Conserve	0.8411	0.5488	0.6524	0.9779
n = 61	Conserve	0.8352	0.5331	0.6382	0.9993

From here we can say that conservative form is way closer to the true solution and with higher grid values we get a good result. But conservation method is taking more number of iterations to arrive at a stable result.

3. Mass flow variation with respect to iterations

• Conservation form

Here we can see that the at 100th iteration the value is highly fluctuating and at 200th iteration it started to settle. At 400 it just almost started to become stable and we can see that by 700th iteration solution already is constant and stabiliezed.

• Non-conservation form

In the same case as above the solution was highly disturebed and by 700 it started to settle and stabilized.

From both the methods we can observe that conservative form is way better for mass flow rate stabilization and it converges a bit slower than non conservative form.We have already observed that conservative method is unstablte at lower grid points. Now let us also take a look at mass flow rate for higher grid rates.

• Conservative form n = 61

• Non-conservative form

Here we can clearly observe that the non conservative form is highly fluctuating for n = 61 and it can be seen that the solution is not settled in both the forms when no of nodes are high. This depicts that the increase in elements number makes graph more clear about the variation in the trend.

4. Mass flow rate comparission

The mass flow rate in the conservative form is smooth and consistent while in the case of non-conservative form it is unstable. Since the mass flow rate is itself considered as a vector which is a variable in conservative form which is why it is stabilised while in non-conservative form it is derived from other variables.

Here we can see that non conservative form is more close to the true value even without the direct dependency of mass flow rate. It is also more clear that as the grid size is smaller the convergence is more fast and accurate to the real value.

Conclusion:

1. Non conservative form converges faster when the elements number is increased when compared to non-conservative form.

2. Conservative form is stable but not really as accurate as non-conservative form.

3. Mass flow rate is more stable and less fluctuating in the case of conservative form.

4. The primitive variables like rho, T, v, M are approximately the same at the end of time marching for both the methods.

5. When the elements number is increased the conservative form takes more iterations to converge.