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

Theory:

Any numerical solution of isentropic flow through a quasi 1D subsonic-supersonic nozzle is unreasonable, while an analytic solution exists for this kind of problem. However, we are trying here to establish the broad-banded capability of CFD in various streams rather than in representing the approximated form of the analytic solution. Here we implement the numerical solution with conservative and non-conservative forms.

Governing equations are matched in terms of PDE in both forms in order to use a time marching solution for the problem. These forms represent in numeric aspect only, if we analytic all the methods are same.

Non-conservative form:

Conservative form:

Non-dimensional quintiles:

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

Nozzle diagram with sections:

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

1. Convergent section
2. Divergent section

McCormack method:

To discretize the numerical solution we use the technique called the MacCormack method. This method is a variation of the two-step lax-wend-off scheme but is much simpler in the application.

The application of the McCormack method proceeds in two steps, a predictor step which is followed by a corrector step.

1. Predictor step: in this step, a provisional value of u at time level n+1 is estimated and obtained by replacing the spatial and temporal derivatives in the previous first-order equation using the following difference scheme delta t is used as a time step.
2. Corrector method: in the corrector step the predicted value n+1 is corrected according to the equation. Note that the corrector step uses a backward differences scheme. The time step used in the corrector method is delta t/2 in contrast to the delta t used in the predictor step.

Main programming:

In order for the code to be modular, the

program is split into three sections.

clear all
close all
clc
%inputs
% o of grid poits
n=31;
x=linspace(0,3,n);
dx=x(2)-x(1);
gamma=1.4;
%no of time step
nt=1400;
%courant number
CFL=0.5;

%non coservative form
tic;
[mass_flow_rate_nc, pressure_nc, mach_number_nc, rho_nc, v_nc, t_nc, rho_throat_nc, v_throat_nc, t_throat_nc, mass_flow_rate_throat_nc, pressure_throat_nc, mach_number_throat_nc]=non_conservative_form(x,dx,n,gamma,nt,CFL);
time_in_non_conservative_form=toc;
hold on

%conservative form function
tic;
[mass_flow_rate_c, pressure_c, mach_number_c, rho_c, v_c, t_c, rho_throat_c, v_throat_c, t_throat_c, mass_flow_rate_throat_c, pressure_throat_c, mach_number_throat_c]=conservative_form(x, dx, n, gamma, nt, CFL);
time_in_conservative_form=toc;
fprintf('time elapsed for non conservative form: 0.4%f',time_in_conservative_form);
hold on

%plotting
%time step variation of properties at throat area in non conservative form

figure(3)
%density
subplot(4,1,1)
plot(linspace(1,nt,nt),rho_throat_nc,'color','m')
ylabel('density')
legend({'density at throat'});
axis([0 1400 0 1.3]);
grid minor;
title('time step variation at throat area - non conservativa form')


%temparature
subplot(4,1,2)
plot(linspace(1,nt,nt),t_throat_nc,'color','c')
ylabel('temparature')
legend({'Temparature at throat'});
axis([0 1400 0.6 1]);
grid minor;

%pressure
subplot(4,1,3)
plot(linspace(1,nt,nt),pressure_throat_nc,'color','g')
ylabel('pressure')
legend({'pressure at throat'});
axis([0 1400 0.4 1.1]);
grid minor;
%mach number

subplot(4,1,4)
plot(linspace(1,nt,nt),mach_number_throat_nc,'color','y')
ylabel('time steps')
legend({'mach number at throat'});
axis([0 1400 0.6 1.4]);
grid minor;

%steady state simulation for primitve value for non conservative form

figure(4)

%density
subplot(4,1,1)
plot(x,rho_nc,'color','c')
ylabel({'density'})
legend('density')
axis([0 3 0 1])
grid minor;
title('variation of floe rate distribution- non conservative form')

%temparature
subplot(4,1,2)
plot(x,t_nc,'color','m')
ylabel({'temparature'})
legend('temparature')
axis([0 3 0 1])
grid minor;

%pressure
subplot(4,1,3)
plot(x,pressure_nc,'color','g')
ylabel({'pressure'})
legend('pressure')
axis([0 3 0 1])
grid minor;


%mach number
subplot(4,1,4)
plot(x,mach_number_nc,'color','y')
ylabel({'mach number'})
xlabel('nozzle length')
legend('mach number')
axis([0 3 0 1])
grid minor;

%plotting for conservative form
%time step variation of properties at throat area in conservative form
figure(5)

%density

subplot(4,1,1)
plot(rho_throat_c,'color','c')
ylabel({'density'})
legend('density at throat')
axis([0 1400 0 1])
grid minor;
title('timestep variation at throat area- conservative form')

%temparature
subplot(4,1,2)
plot(t_throat_c,'color','m')
ylabel({'temparature'})
legend('temparature at throat')
axis([0 1400 0.6 1])
grid minor;

%pressure
subplot(4,1,3)
plot(pressure_throat_c,'color','g')
ylabel({'pressure'})
legend('pressure at throat')
axis([0 1400 0.2 0.8])
grid minor;


%mach number
subplot(4,1,4)
plot(mach_number_throat_c,'color','y')
ylabel({'mach number'})
xlabel('time steps')
legend('mach number at throat')
axis([0 1400 0.6 1.4])
grid minor;

%steady state simulation for primitve values for conservative form
figure(6)
%density
subplot(4,1,1)
plot(x,rho_c,'color','c')
ylabel({'density'})
legend('density')
axis([0 3 0 1])
grid minor;
title('variation of flow rate distribution- conservative form')

%temparature
subplot(4,1,2)
plot(x,t_c,'color','m')
ylabel({'temparature'})
legend('temparature')
axis([0 3 0 1])
grid minor;

%pressure
subplot(4,1,3)
plot(x,pressure_c,'color','g')
ylabel({'pressure'})
legend('pressure')
axis([0 3 0 1])
grid minor;


%mach number
subplot(4,1,4)
plot(x,mach_number_nc,'color','y')
ylabel({'mach number'})
xlabel('nozzle length')
legend('mach number')
axis([0 3 0 1])
grid minor;

%comparision of narmalised mass flow rate of both
figure(7)
hold on
plot(x,mass_flow_rate_nc,'r','linewidth',1.5)
hold on;
plot(x,mass_flow_rate_c,'b','linewidth',1.5)
hold on
grid on

legend('non conservative', 'conservative')
xlabel('length of ht domain')
ylabel('mass flow rate')
title('comparision of normalised floe rate of bboth forms')

Explanation:

1. As a basic start, we start clear all, close all clc, command to start the page.
2. Input parameters are assigned for the variables x, dx, n, gamma,nt, CFL.
3. Nonconservative, conservative forms are created as a separate function with the necessary inputs as instructed in the video.
4. Tic and toc command is used for calculating the time elapsed.
5. Plotting is done for the following solutions

Time step variation at throat area of both forms

Variation of flow rate distribution of both forms

Comparison of the normalized mass flow rate of both forms

6. Each plot is computed as subplots of four primitive variables such as density, temperature, pressure, and Mach number at the throat as well as convergent and divergent sections.

function [mass_flow_rate_nc, pressure_nc, mach_number_nc, rho_nc, v_nc, t_nc, rho_throat_nc, v_throat_nc, t_throat_nc, mass_flow_rate_throat_nc, pressure_throat_nc, mach_number_throat_nc]=non_conservative_form(x,dx,n,gamma,nt,CFL);
%initial boundary condition
rho_nc=1-0.3146*x;
t_nc=1-0.2341*x;
v_nc=(0.1+1.09*x).*t_nc.^0.5;
%area and throat profile
a=1+2.2*(x-1.5).^2;
throat=find(a==1);

%outer time loop
for k=1:nt
  %time-step control
  dt=min(CFL.*dx./(sqrt(t_nc)+v_nc));
  %saving a copy of old values
  rho_old=rho_nc;
  t_old=t_nc;
  v_old=v_nc;
  
  %predictor methods
  for j=2:n-1
    %simplifying the equation spatial terms
    dv_dx=(v_nc(j+1)-v_nc(j))/dx;
    drho_dx=(rho_nc(j+1)-rho_nc(j))/dx;
    dlog_a_dx=(log(a(j+1))-log(a(j)))/dx;
    dt_dx=(t_nc(j+1)-t_nc(j))/dx;
    
    %continuty equation
    drho_dt_p(j)=-rho_nc(j)*dv_dx-rho_nc(j)*v_nc(j)*dlog_a_dx-v_nc(j)*drho_dx;
    
    %momentum equation
    dv_dt_p(j)=-v_nc(j)*dv_dx-(1/gamma)*(dt_dx+(t_nc(j)*drho_dx/rho_nc(j)));
    
    %energy equation
    dT_dt_p(j)=-v_nc(j)*dt_dx-(gamma-1)*t_nc(j)*(dv_dx+v_nc(j)*dlog_a_dx);
    
    %solution update
    
    rho_nc(j)=rho_nc(j)+drho_dt_p(j)*dt;
    v_nc(j)=v_nc(j)+dv_dt_p(j)*dt;
    t_nc(j)=t_nc(j)+dT_dt_p(j)*dt;
    
end   
%corrector methods
  for j=2:n-1
    %simplifying the equation spatial terms
    dv_dx=(v_nc(j)-v_nc(j-1))/dx;
    drho_dx=(rho_nc(j)-rho_nc(j-1))/dx;
    dt_dx=(t_nc(j)-t_nc(j-1))/dx;
    dlog_a_dx=(log(a(j))-log(a(j-1)))/dx;
    
    %continuty equation
    drho_dt_c(j)=-rho_nc(j)*dv_dx-rho_nc(j)*v_nc(j)*dlog_a_dx-v_nc(j)*drho_dx;
    
    %momentum equation
    dv_dt_c(j)=-v_nc(j)*dv_dx-(1/gamma)*(dt_dx+(t_nc(j)*drho_dx/rho_nc(j)));
    
    %energy equation
    dT_dt_c(j)=-v_nc(j)*dt_dx-(gamma-1)*t_nc(j)*(dv_dx+v_nc(j)*dlog_a_dx);
     
end
%computing average time derivative
drho_dt=0.5*(drho_dt_p+drho_dt_c);
dv_dt=0.5*(dv_dt_p+dv_dt_c);
dT_dt=0.5*(dT_dt_p+dT_dt_c);

%final solution update
for i=2:n-1
rho_nc(i)=rho_old(i)+drho_dt(i)*dt;
v_nc(i)=v_old(i)+dv_dt(i)*dt;
t_nc(i)=t_old(i)+dT_dt(i)*dt;
end
%applying boundary comditions
%at inlet
rho_nc(1)=1;
v_nc(1)=2*v_nc(2)-v_nc(3);
t_nc(1)=1;

%at out let
rho_nc(n)=2*rho_nc(n-1)-rho_nc(n-2);
v_nc(n)=2*v_nc(n-1)-v_nc(n-2);
t_nc(n)=2*t_nc(n-1)-t_nc(n-2);



%defining mass fllow rate pressure and mach number
mass_flow_rate_nc=rho_nc.*a.*v_nc;
pressure_nc=rho_nc.*t_nc;
mach_number_nc=(v_nc./sqrt(t_nc));

%calculating variables at throat
rho_throat_nc(k)= rho_nc(throat);
t_throat_nc(k)=t_nc(throat);
v_throat_nc(k)=v_nc(throat);
mass_flow_rate_throat_nc(k)=mass_flow_rate_nc(throat);
pressure_throat_nc(k)=pressure_nc(throat);
mach_number_throat_nc(k)=mach_number_nc(throat);

%plotting
figure(2)
if k==50
  plot(x,mass_flow_rate_nc,'m','linewidth',1.25);
  hold on
  elseif k==100
  plot(x,mass_flow_rate_nc,'c','linewidth',1.25);
  hold on
  elseif k==200
  plot(x,mass_flow_rate_nc,'g','linewidth',1.25);
  hold on
  elseif k==400
  plot(x,mass_flow_rate_nc,'y','linewidth',1.25);
  hold on
  elseif k==800
  plot(x,mass_flow_rate_nc,'k','linewidth',1.25);
  hold on 
end
% labeling
title('mass flow rate at different time steps- non conservative form');
xlabel('distance');
ylabel('mass flow rate');
legend({'50^t^h timestep','100^th timestep','200^th timestep','400th timestep','800th timestep'})
axis([0 3 0 2]);
grid on

end
endfunction%

Explanation:

1. The purpose of the function is to calculate the values of mass flow rate, pressure, Mach number, density, velocity, temperature, density at the throat, the velocity at the throat, temp at the throat, the mass flow rate at the throat, the pressure at the throat, Mach number at the throat by the known equation for all the grid points.
2. For this initial boundary conditions are established for density, temperature, and velocity and profiles for calculating area and throat.
3. Now outer temp loop is computed for the loop to calculate all the grid points for various formulas inside the loop.
4. McCormack technique is used to solve the numerical solution of governing equations such as continuity, momentum, and energy equations.
5. The technique implies two methods for solving namely predictor and then corrector method.
6. Predictor and corrected method are assigned within the for loop individually for calculating values.
7. Boundary conditions are applied within the for loop and subsequently used for calculating the values throat section.
8. Finally plotting is done for distance(x) and mass flow rate difference time steps.

function [mass_flow_rate_c, pressure_c, mach_number_c, rho_c, v_c, t_c, rho_throat_c, v_throat_c, t_throat_c, mass_flow_rate_throat_c, pressure_throat_c, mach_number_throat_c]=conservative_form(x, dx, n, gamma, nt, CFL);
  
  %determining initial profile density and tempdir
  for i=1:n
    if (x(i)>=0 && x(i)<=0.5)
      rho_c(i)=1;
      t_c(i)=1;
    elseif (x(i)>=0.5 && x(i)<=1.5)
      rho_c(i)=1-0.366*(x(i)-0.5);
      t_c(i)=1-0.167*(x(i)-0.5);
    elseif (x(i)>= 1.5 && x(i)<=3)
      rho_c(i)=0.634-0.3879*(x(i)-1.5);
      t_c(i)=0.833-0.3507*(x(i)-1.5);
    endif
  endfor
  %%area & throat profiles
  a=1+2.2*(x-1.5).^2;
  throat=find(a==1);
  
  %initial profile for velocity and pressure\
  v_c=0.59./(rho_c.*a);
  pressure_c=rho_c.*t_c;
  
  %calclating the initial solution vectors(u)
  u1=rho_c.*a;
  u2=rho_c.*a.*v_c;
  u3=rho_c.*a.*(t_c./(gamma-1))+((gamma/2).*v_c.^2);
  %outer time loop
  for k =1:nt
    % time step control
    dt= min(CFL.*dx./(sqrt(t_c)+v_c));
    
    %solving a copy of old value
    u1_old=u1;
    u2_old=u2;
    u3_old=u3;
    %calculating the flux vector(F)
    f1=u2;
    f2=((u2.^2)./u1)+((gamma-1)/gamma)*(u3-(gamma/2)*((u2.^2)./(u1)));
    f3=((gamma*u2.*u3)./u1)-(gamma*(gamma-1)*((u2.^3)./(2*(u1.^2))));
     %predictor method
for j=2:n-1
      %calculating the source term(j)
      j2(j)=(1/gamma)*rho_c(j)*t_c(j)*((a(j+1)-a(j))/dx);
      
      %continuty equation
      du1_dt_p(j)=-(f1(j+1)-f1(j))/dx;
      
      %momentum equation
      du2_dt_p(j)=-((f2(j+1)-f2(j))/dx)+j2(j);
      
      %energyequation
      du3_dt_p(j)=-(f3(j+1)-f3(j))/dx;
      
      %updating the new values predictor
      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
%%updating the flux vector
 f1=u2;
 f2=((u2.^2)./u1)+((gamma-1)/gamma)*(u3-((gamma*0.5*(u2.^2))./(u1)));
 f3=((gamma*u2.*u3)./u1)-((gamma*0.5*(gamma-1)*(u2.^3))./(u1.^2));
 %corrector method
for j=2:n-1
   %updating the source term(j)
   j2(j)=(1/gamma)*rho_c(j)*t_c(j)*((a(j)-a(j-1))/dx);
  
  %continuty equation
      du1_dt_c(j)=-((f1(j)-f1(j-1))/dx);
      
      %momentum equation
      du2_dt_c(j)=-((f2(j)-f2(j-1))/dx)+j2(j);
      
      %energyequation
      du3_dt_c(j)=-((f3(j)-f3(j-1))/dx);
end      
   %computing the avg time derivativ
du1_dt=0.5*(du1_dt_c+du1_dt_p);
du2_dt=0.5*(du2_dt_c+du2_dt_p);
du3_dt=0.5*(du3_dt_c+du3_dt_p);

%final solution update
for i=2:n-1
  
 u1(i)=u1_old(i)+(du1_dt(i)*dt);
 u2(i)=u2_old(i)+(du2_dt(i)*dt);
 u3(i)=u3_old(i)+(du3_dt(i)*dt);
end
%applying boundary conditions
%at inlet
u1(1)=rho_c(1)*a(1);
u2(1)=2*u2(2)-u2(3);
v_c(1)=u2(1)./u1(1)
u3(1)=u1(1)*((t_c(1)/(gamma-1))+((gamma/2)*(v_c(1)).^2));

%at  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 primitive value update
rho_c=u1./a;
v_c=u2./u1;
t_c=(gamma-1)*((u3./u1-(gamma/2)*(v_c).^2));

%defining mass fllow rate pressure and mach number
mass_flow_rate_c=rho_c.*a.*v_c;
pressure_c=rho_c.*t_c;
mach_number_c=(v_c./sqrt(t_c));

rho_throat_c(k)=rho_c(throat);
t_throat_c(k)=t_c(throat);
v_throat_c(k)=v_c(throat);
mass_flow_rate_throat_c(k)=mass_flow_rate_c(throat);
pressure_throat_c(k)=pressure_c(throat);
mach_number_throat_c(k)=mach_number_c(throat);

%plotting
figure(1)
if k==50
  plot(x,mass_flow_rate_c,'m','linewidth',1.25);
  hold on
  elseif k==100
  plot(x,mass_flow_rate_c,'c','linewidth',1.25);
  hold on
  elseif k==200
  plot(x,mass_flow_rate_c,'g','linewidth',1.25);
  hold on
  elseif k==400
  plot(x,mass_flow_rate_c,'y','linewidth',1.25);
  hold on
  elseif k==800
  plot(x,mass_flow_rate_c,'k','linewidth',1.25);
  hold on 
end
% labeling
title('mass flow rate at different time steps- conservative form')
xlabel('distance')
ylabel('mass flow rate')
legend({'50^t^h timestep','100^th timestep','200^th timestep','400th timestep','800th timestep'})
axis([0 3 0.4 1])
grid on



end
endfunction
%

Explanation:

1. The purpose of the function is to calculate the values of mass flow rate, pressure, Mach number, density, velocity, temperature, density at the throat, the velocity at the throat, temp at the throat, the mass flow rate at the throat, the pressure at the throat, Mach number at the throat by the known equation for all the grid points.
2. To determine the initial profile and density and temp for loop is used input the known conditions within an else if condition.
3. Initial values of area, throat, velocity, and pressure are calculated using formulas.
4. In order to simplify and avoid errors, equations are splatted into solution vector(U), flux vector(F), and source team(J).
5. Now outer time loop is computed for the loop to calculate all the grid points for various formulas inside the loop.
6. McCormack technique is used to solve the numerical solution of governing equations such as continuity, momentum, and energy equation.
7. The technique implies two methods for solving namely the predictor method followed by the corrector method.
8. Predictor and corrector methods are assigned within the for loop individually for calculating values.
9. The calculated values are avg in the next step and updating the solution with old values.
10. Boundary conditions are applied within the for loop and subsequently used for calculating h=the mass flow rate, Mach number, and pressure and also calculating the values throat section.
11. Finally plotting is done for distance x and mass flow rate different time steps.

outputs:

time step variation at primitive variable at throat area:

Observations:

From the above plots, we can conclude that initial variations observed in non-conservative form for a certain period of time make it consume a higher amount of time to achieve a stable solution when compared to a conservative form

variation in flow rate distribution:

Observation:

Here we can observe that the above plots will say that both the forms are in the same phase giving exact same solution for all the values.

mass flow rate at different time steps:

Observation:

Above plots shows that mass flow rate different time steps along the distance oscillates initially, however the solution becomes stable after the 400^th time step for the both forms attaining a constant values

comparison plot for a normalized mass flow rate of both forms:

Observations:

Even though the final solution of both the forms is almost similar, the non-conservative form oscillates much when compared to the conservative form. The letter also gives an almost stable solution along the length. By this, we can conclude that the conservative form attains stability in a lesser number of iterations with the earliest time.

Conclusion:

Cumulating all the above results we can evidently conclude that conservative form converges slightly faster than non-conservative form offering stable solution along with the simulation. Yet both the forms are dependable in calculating the solution which is almost précised.