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

OBJECTIVE:

To Simulate the isentropic flow through a Quasi 1D subsonic nozzle by deriving both the conservation and non-conservation forms of the governing equations and solve them by implementing the MacCormack's technique using MATLAB.

Problem Statement:

Determine the steady-state temperature distribution for the flow-field variables and to investigate the difference between the two forms of governing equations by comparing their solutions along with the plots stated below.

1. Steady-State distribution of primitive variable inside the nozzle.
2. Time-wise variation of the primitive variable.
3. Variation of Mass rate distribution inside the nozzle at different time steps during the time-marching process.
4. Comparision of Normalized mass flow rate distributions of both forms.

THEORY:

Any numerical solution of isentropic flow though a quasi 1D subsonic-supersonic nozzle is unreasonable, while analytic solution exists for these kind of problems. However we are trying here to establish the broadened capability of CFD in various streams rather in representing the approximated form of anayltic solution. Here, we implement the numerical solution with Conservative and Non-Conservative forms.

Governing Equations are marched in terms of PDE for both forms in order to use Time Marching solution for the problem. These forms represent in numerical aspect only, if we analytic all the methods are same.

Non Conservative Form:

Conservative Form:

Non Dimensional Quantities:

Also as part of the problem setup we represent quantites in terms of non-dimensional quantities.

Nozzle diagram with sections:

Since the problem setup is 1D, we consider only distribution in x-axis. Here Nozzle can be split into two sections

• Convergent Section
• Divergent Section

MacCarmack method:

To discretize the numerical solution we use the technique called MacCarmack method. The MacCormack method is a variation of the two-step Lax-Wendroff Scheme but is much simpler in application.

The application of MacCormack method proceeds in two steps; a predictor step which is followed by a corrector step.

• Predictor step: In the predictor step, a "provisional" value of at time level (denoted by ) is estimated obtained by replacing the spatial and temporal derivatives in the previous first order equation using Forward Difference Scheme. is used for Timestep.
• Corrector step: In the corrector step, the predicted value is corrected according to the equation. Note that the corrector step uses Backward Difference Scheme for spatial derivative. The time-step used in the corrector step is in contrast to the used in the predictor step.

CODE:

close all
clear all
clc

L=3;                    %length of nozzle
n=[31 62 93];           %number of grid points
nt=1:1400;              %number of time steps
cfl = 0.5;              %CFL number
for i = 1:length(n)
x=linspace(0,L,n(i));
[mfr_nc,mfr_throat_nc]=non_conservative(L,n(i),nt,cfl);
[mfr_c,mfr_throat_c]=conservative(L,n(i),nt,cfl);
figure(7)
subplot(2,1,1)
plot(nt,mfr_throat_c,'b',nt,mfr_throat_nc,'r','linewidth',1)
xlabel('Time steps (nt)')
ylabel('Mass flow rate')
legend('Conservative form','Non-Conservative form')
title('Mass flow rate using "Non-Conservative form vs Conservative form" at throat point i=16')

subplot (2,1,2)
plot(x,mfr_c,'b',x,mfr_nc,'r','linewidth',1)
xlabel('Length of nozzle')
ylabel('Mass flow rate')
legend('Conservative form','Non-Conservative form')
title('Mass flow rate using "Non-Conservative form vs Conservative form" in the nozzle')

figure(8)
subplot (2,1,1)
plot(nt,mfr_throat_c,'linewidth',1)
hold on
xlabel('Time steps (nt)')
ylabel('Mass flow rate')
legend('@31 nodes','@62 nodes','@93 nodes')
title('Conservative form: Grid independency test at throat point i=16')
subplot (2,1,2)
plot(nt,mfr_throat_nc,'linewidth',1)
hold on
xlabel('Time steps (nt)')
ylabel('Mass flow rate')
legend('@31 nodes','@62 nodes','@93 nodes')
title('Non-Conservative form: Grid independency test at throat point i=16')
figure(9)
subplot (2,1,1)
plot(x,mfr_c,'linewidth',1)
hold on
xlabel('Length of nozzle')
ylabel('Mass flow rate')
legend('@31 nodes','@62 nodes','@93 nodes')
title('Conservative form: Grid independency test along the nozzle')
subplot (2,1,2)
plot(x,mfr_nc,'linewidth',1)
hold on
xlabel('Length of nozzle')
ylabel('Mass flow rate')
legend('@31 nodes','@62 nodes','@93 nodes')
title('Non-Conservative form: Grid independency test along the nozzle')

end

%Non-Conservative function:
function [mfr_nc,mfr_throat_nc]= non_conservative(L,n,nt,cfl)


x=linspace(0,L,n);
dx=x(2)-x(1);

% initial conditions rho, a, v, T, mfr, M, p
rho = 1-0.3146*x;
T= 1-0.2314*x;
v=(0.1+1.09*x).*T.^0.5;
a= 1+2.2*(x-1.5).^2;
M= v./(T.^0.5);
mfr = rho.*a.*v;
p=rho.*T;

g = 1.4;

for z=1:31
    u=(T(z))^(1/2);
    dt(z) = cfl*dx/(u + v(z));
end
dt=min(dt);

for i = 1:length(nt)
  %storing initial values
    rho_old = rho;
    T_old = T;
    v_old = v;
  for j=2:n-1
  %predictor step
    drho_dx = (rho(j+1)-rho(j))/dx;
    dv_dx = (v(j+1)-v(j))/dx;
    dT_dx = (T(j+1)-T(j))/dx;
    da_dx = (log(a(j+1))-log(a(j)))/dx;
    drho_dt_p(j) = -rho(j)*dv_dx - rho(j)*v(j)*da_dx - v(j)*drho_dx;        %continuty
    dv_dt_p(j) = -v(j)*dv_dx - (1/g)*(dT_dx + (T(j)/rho(j))*drho_dx);       %momentum
    dT_dt_p(j) = -v(j)*dT_dx - (g-1)*T(j)*(dv_dx + v(j)*da_dx);             %energy
  %updating primitives
    v(j) = v_old(j) + dv_dt_p(j)*dt;
    rho(j) = rho_old(j) + drho_dt_p(j)*dt;
    T(j) = T_old(j) + dT_dt_p(j)*dt;
  %corrector step
    drho_dx = (rho(j)-rho(j-1))/dx;
    dv_dx = (v(j)-v(j-1))/dx;
    dT_dx = (T(j)-T(j-1))/dx;
    da_dx = (log(a(j))-log(a(j-1)))/dx; 
    drho_dt_c(j) = -rho(j)*dv_dx - rho(j)*v(j)*da_dx - v(j)*drho_dx;        %continuity    
    dv_dt_c(j) = -v(j)*dv_dx - (1/g)*(dT_dx + (T(j)/rho(j))*drho_dx);       %momentum
    dT_dt_c(j) = -v(j)*dT_dx - (g-1)*T(j)*(dv_dx + v(j)*da_dx);             %energy
  end
  %average time derivatives
    drho_dt_avg = 0.5*(drho_dt_p + drho_dt_c);
    dv_dt_avg = 0.5*(dv_dt_p + dv_dt_c);
    dT_dt_avg = 0.5*(dT_dt_p + dT_dt_c);
  %final corected values
  for k = 2:n-1
    rho(k) = rho_old(k) + drho_dt_avg(k)*dt;
    v(k) = v_old(k) + dv_dt_avg(k)*dt;
    T(k) = T_old(k) + dT_dt_avg(k)*dt;
  end
  %entry boundary conditions
    v(1) = 2*v(2)-v(3);
  %exit boundary conditions
    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);
  %Final update of primitives
    p = rho.*T;
    M= v./T.^0.5;
    error(k) =  abs(v(k)-v(k-1));
    mfr= rho.*a.*v;
  %to find values at throat point
    T_throat(i) = T(16);
    rho_throat(i) = rho(16);
    v_throat(i) = v(16);
    M_throat(i) = M(16);
    p_throat(i) = p(16);
    error_throat(i)=error(16);
    mfr_throat(i)=mfr(16);
    if n==31
for ii=[1:50:201 700]
    if ii==i
       figure(1)
     plot(x,mfr,'linewidth',1)
     hold on
     title ('Non-Conservative form - Variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process')
     xlabel('Length of nozzle')
     ylabel('Mass flow rate')   
    elseif i==1400
     plot(x,mfr,'*r')
     legend('@1','@51','@101','@151','@201','@700','@1400') 
    end
end
end
  
end
if n==31
figure(2)
     plot(x,v,'b',x,p,'r',x,T,'k',x,rho,'m',x,M,'g',x,mfr,'c')
     ylim([-0.1 3.7])
     title (sprintf('Time step = %d',i))
     xlabel('Length of nozzle')
     ylabel('Non dimensional parameters')
     legend('Velocity','Pressure','Temperature','Density','Mach number','Mass flow rate','location','northwest')
     title('Non-Conservative form - Steady-state distribution of primitive variables inside the nozzle')

figure (3)
plot(nt,rho_throat)  
hold on
plot(nt,v_throat,'linewidth',1)   
plot(nt,p_throat,'linewidth',1) 
plot(nt,T_throat,'linewidth',1)    
plot(nt,M_throat,'linewidth',1)
plot(nt,mfr_throat,'linewidth',1)
xlabel('Time steps (nt)')
ylabel('Non dimensional parameters')
legend('Density','Velocity','Pressure','Temperature','Mach number','Mass flow rate')
title('Non-Conservative form - Time-wise variation of the primitive variables at throat point i=16')
end
mfr_nc=mfr;
mfr_throat_nc=mfr_throat;
end

%Conservative function:
function [mfr_c,mfr_throat_c]=conservative(L,n,nt,cfl)

x=linspace(0,L,n);
dx=x(2)-x(1);

% initial conditions rho, a, v, T, M, mfr, p
for ii=1:length(x)
    if (x(ii)<=0.5)
        rho(ii) = 1;
        T(ii) = 1;
    elseif (x(ii)>0.5) & (x(ii)<=1.5)
        rho(ii) = 1-0.366*(x(ii)-0.5);
        T(ii)= 1-0.167*(x(ii)-0.5);
    elseif (x(ii)>1.5) & (x(ii)<=3.5)
        rho(ii) = 0.634-0.3879*(x(ii)-1.5);
        T(ii)= 0.833-0.3507*(x(ii)-1.5);
    end
end
a= 1+2.2*(x-1.5).^2;     %area
g = 1.4;                 %gamma value
% solution vectors u1, u2, u3
u1 = rho.*a;
u2 = 0.59*ones(1,n);
v=u2./u1;
u3 = u1.*((T./(g-1))+(g/2).*v.^2);

M= v./(T.^0.5);          %mach number
mfr = rho.*a.*v;         %mass flow rate
p=rho.*T;                %pressure


for z=1:31               %to find dt value
u=(T(z))^(1/2);
dt(z) = cfl*dx/(u + v(z));
end
dt=min(dt);

for i = 1:length(nt)
    %storing initial values
    u1_old = u1;
    u2_old = u2;
    u3_old = u3;
   
    for j=2:n-1
   %flux vectors
        f1 = u2;
        f2 = ((u2.^2)./u1)+((g-1)/g)*(u3-((g/2)*(u2.^2)./u1));
        f3 = ((g.*u2.*u3)./u1) - ((g*(g-1)/2)*(u2.^3)./(u1.^2));
   %source term
        j1 = (1/g)*rho(j)*T(j)*((a(j+1)-a(j))/dx);        
   %predictor step
        du1_dt_p(j)= -(f1(j+1)-f1(j))/dx;           %continuity   
        du2_dt_p(j)= -(f2(j+1)-f2(j))/dx + j1;      %momentum
        du3_dt_p(j)= -(f3(j+1)-f3(j))/dx;           %energy
   %updating solution vectors          
        u1(j) = u1_old(j) + du1_dt_p(j)*dt;
        u2(j) = u2_old(j) + du2_dt_p(j)*dt;
        u3(j) = u3_old(j) + du3_dt_p(j)*dt;
   %updating primitives
        rho=u1./a;
        T=(g-1).*(u3./u1-(g/2).*(u2./u1).^2);  
   %updating flux terms
        f1=u2;
        f2=u2.^2./u1+((g-1)/g).*(u3-(g/2)*u2.^2./u1);
        f3=g.*u2.*u3./u1-0.5*g*(g-1).*u2.^3./u1.^2;
   %source term for corrector step
        j2 = (1/g)*(rho(j)*T(j)*(a(j)-a(j-1))/dx);    
   %corrector step
        du1_dt_c(j) = -(f1(j)-f1(j-1))/dx;            %continuity   
        du2_dt_c(j) = -(f2(j)-f2(j-1))/dx + j2;       %momentum
        du3_dt_c(j) = -(f3(j)-f3(j-1))/dx;            %energy
    end
    %average time derivatives
        du1_dt_avg = 0.5*(du1_dt_p + du1_dt_c);
        du2_dt_avg = 0.5*(du2_dt_p + du2_dt_c);
        du3_dt_avg = 0.5*(du3_dt_p + du3_dt_c);
    %final corected values
    for k = 2:n-1
        u1(k) = u1_old(k) + du1_dt_avg(k)*dt;
        u2(k) = u2_old(k) + du2_dt_avg(k)*dt;
        u3(k) = u3_old(k) + du3_dt_avg(k)*dt;
    end
    %entry boundary conditions
        u1(1) = a(1);
        u2(1) = 2*u2(2)-u2(3);
        u3(1) = u1(1)*(T(1)/(g-1) + (g*(u2(1)/u1(1))^2)/2);
    %exit boundary conditions
        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 update of primitives
        rho = u1./a;
        v = u2./u1;
        T=(g-1).*(u3./u1-(g/2).*(u2./u1).^2);  
        p = rho.*T;
        M= v./T.^0.5;
        error(k) =  abs(v(k)-v(k-1));
        mfr= u2;
  %to find values at throat point
    T_throat(i) = T(16);
    rho_throat(i) = rho(16);
    v_throat(i) = v(16);
    M_throat(i) = M(16);
    p_throat(i) = p(16);
    error_throat(i)=error(16);
    mfr_throat(i)=mfr(16);  
    
if n==31  
    for ii=[1:50:201 700]
    if ii==i
       figure(4)
     plot(x,mfr,'linewidth',1)
     hold on
     title ('Conservative form - Variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process')
     xlabel('Length of nozzle')
     ylabel('Mass flow rate')   
    elseif i==1400
     plot(x,mfr,'*r')
     legend('@1','@51','@101','@151','@201','@700','@1400') 
    end
end
end
  
end
if n==31
figure(5)
     plot(x,v,'b',x,p,'r',x,T,'k',x,rho,'m',x,M,'g',x,mfr,'c')
     ylim([-0.1 3.7])
     title (sprintf('Time step = %d',i))
     xlabel('Length of nozzle')
     ylabel('Non dimensional parameters')
     legend('Velocity','Pressure','Temperature','Density','Mach number','Mass flow rate','location','northwest')
     title('Conservative form - Steady-state distribution of primitive variables inside the nozzle')

figure (6)
plot(nt,rho_throat)  
hold on
plot(nt,v_throat,'linewidth',1)   
plot(nt,p_throat,'linewidth',1) 
plot(nt,T_throat,'linewidth',1)    
plot(nt,M_throat,'linewidth',1)
plot(nt,mfr_throat,'linewidth',1)
xlabel('Time steps (nt)')
ylabel('Non dimensional parameters')
legend('Density','Velocity','Pressure','Temperature','Mach number','Mass flow rate')
title('Conservative form - Time-wise variation of the primitive variables at throat point i=16')
end
mfr_c=mfr;
mfr_throat_c=mfr_throat;

end

EXPLANATION:

• Code is started with clear all, close all and clc to enable a fresh start without the manipulation of previous values.
• Input parameters are assigned for the variables x, dx, n, gamma, nt and C.
• Non Conservative and Conservative forms are created as separate functions with the neccessary inputs which is detailed below this section.
• tic toc command is used for calculating the time elasped between the functions.
• Plotting is done for the following solutions
  1)Timestep variation at throat area of both forms
  2)Variation of Flow Rate Distribution of both forms
  3)Comparison of Normalized Mass Flow Rate of both forms
• Each plot is computed as subplots of four primitive variables such as Density, Temperature, Pressure and Mach Number at throat as well as Convergent/Divergent section.

Non Conservation Function:

The purpose of the function is to calculate the values of Mass flow rate, Pressure, Mach number, Density, velocity, Temperature, Density at throat, velocity at throat, Temperature at throat, Mass flow rate at throat, Pressure at throat, Mach number at throat by the known equations for all the grid points.

• For this initial boundary conditions are established for Density, Temperature & Velocity and profiles for calculating area & throat.
• Now Outer Time loop is computed for for loop to calculate all the grid points for various formulas inside the loop.
• MacCarmack technique is used to solve the numerical solution of governing equations such as Continuity, Momentum & Energy Equation.
• The technique implies two methods for solving namely Predictor method followed by Corrector method.
• Predictor & Corrector method are assigned within the for loop individually for calculating values.
• The calculated values are averaged in the next step and updating the solution with old values.
• Boundary Conditions are applied within the for loop and subsequently used for calculating the Mass Flow Rate, Mach Number & Pressure and also calculating the values throat section.
• Finally plotting is done for Distance(x) & Mass Flow Rate for differrent Timesteps.
• Program is ended.

Conservation Function:

The purpose of the function is to calculate the values of Mass flow rate, Pressure, Mach number, Density, velocity, Temperature, Density at throat, velocity at throat, Temperature at throat, Mass flow rate at throat, Pressure at throat, Mach number at throat by the known equations for all the grid points.
To determine the initial profile of Density & Temperature, for loop is used inputing the known conditions within an else-if command.

• Initial values of area, throat, velocity & Pressure are calculated using formulas.
• In order to simplify and avoid errors, equation is splitted into Solution Vector(U), Flux Vector(F) and Source Term(J).
• Now Outer Time loop is computed for for loop to calculate all the grid points for various formulas inside the loop.
• MacCarmack technique is used to solve the numerical solution of governing equations such as Continuity, Momentum & Energy Equation.
• The technique implies two methods for solving namely Predictor method followed by Corrector method.
• Predictor & Corrector method are assigned within the for loop individually for calculating values.
• The calculated values are averaged in the next step and updating the solution with old values.
• Boundary Conditions are applied within the for loop and subsequently used for calculating the Mass Flow Rate, Mach Number & Pressure and also calculating the values throat section.
• Finally plotting is done for Distance(x) & Mass Flow Rate for differrent Timesteps.
• Program is ended.

OUTPUTS:

Non-conservative form:

fig(1)

fig(2)

fig(3)

fig(1) shows Steady state distribution of non dimensional primitive variables at the end of time step. These values of flow field variables obliges the exact values shown in John D. anderson book.

fig(2) shows Time wise variation of primitive variables at throat point i =16. Around the time step of 500 the flow field variables obtained stable solution.

fig(3) shows the variation of mass flow rate distribution inside the nozzle at different time steps during the time marching process. The mass flow rate finally settled at 0.58.

Conservative form:

fig(4)

fig(5)

fig(6)

fig(4) shows Steady state distribution of non dimensional primitive variables at the end of time step. These values of flow field variables obliges the exact values shown in John D. anderson book.

fig(5) shows Time wise variation of primitive variables at throat point i =16. Around the time step of 600 the flow field variables obtained stable solution.

fig(6) shows the variation of mass flow rate distribution inside the nozzle at different time steps during the time marching process. The mass flow rate finally settled at 0.56.

Non-conservative form vs conservative form:

fig(7)

fig(7) shows comparision of mass flow rate using both conservative and non-conservative forms.

from subplot 1 : Both the forms reached the stable and exact solution by the end of time steps. But the conservative form obatined the stability at lesser time steps.

from subplot 2 : The difference between both the forms is very less in the nozzle. But the mass flow rate in the conservative form is constant where as in non-conservative form it is varying with 0.582±0.003.

From fig(7) we can conclude that the final solution is not much changing when using conservative and non-conservative forms. Although there is a sligth change at the throat area alone in the conservative form, which is not considerably significant. Though the stability is obtained at lesser time steps the computational time is more in case of conservative form than the other.

Grid independency test for both Non-conservative form & Conservative form:

fig(8)

fig(9)

fig(8) and fig(9) shows grid independence test in both the forms.

fig(8) - mass flow rate varying with time step at throat point - subplot(1) - conservative form & subplot(2) - non-conservative form.

At 3 different grid sizes, though the values are different at the starting of time marching process the solution obatined is stable and almost same by the end of time marching.

fig(9) - mass flow rate in the nozzle - subplot(1) - conservative form & subplot(2) - non-conservative form.

The solution obtained using 3 different grid sizes resulting 3 different solutions within the range of 0.589±0.006. i.e. the difference between them is very small which can be neglected/ignored.

From fig(8) and fig(9) we can conclude that the plot is not much changing when the number of grid points are doubled up and nearly tripled up. Although there is a sligth change at the throat area alone in the conservative form, which is not considerably significant.

GRID DEPENDENCE TEST:

OBSERVATION:

Grid dependency is an important factor to determine whether the problem can be solved faster. So here we increased the Grid points from 31 to 61 and populated them in a sheet.

It is clear from the above table that the values are almost same along all the Primitive Variables for both the forms even when the Grid points(n) is changed from 31 to 61 showcasing that it doesn't get affected by the increase of Grid points.

So we can safely conclude that the above problem is Grid Dependant.

CONCLUSION:

Cumulating all the above results we can evidently conclude that Conservative form converges slightly faster than Non Conservative form offering stable solution along the simulation. Yet both the forms are dependable in calculating the solutions which is almost precise.