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

OBJECTIVE OF THE CHALLENGE:

Objective of the challenge is to simulate the isentropic flow through a 1D subsonic-supersonic nozzle by deriving the governing (continuity,momentum and energy) equations in both conserved and non conserved form using Maccormack's technique.

Governing equations in non conservative form :

continuity equation (non-dimensional):

Momentum equation(non-dimensional):

Energy equation(non-dimensional):

all the primitive variables like ρ',v',A',T,'x', in the above equations are non-dimensionalised ,which are calculated by dividing them with the values at inlet which is considered as reservior

Conservative form of equations:

continuity equation:

=0

momentum equation:

Energy equation:

Maccormack's technique:

It is a widely used discretization technique to obtain numerical solution for hyperbolic Partial Differential Equtions

In Maccormack's method solution is found in two steps namely predictor step and corrector step explicitely. In predictor step an initial solution is found by using the initial boundary conditions by using forward differencing scheme and in the corrector step solution found in the predictor step is used to find an updated solution and then mean of both the solutions(otained in predictor and corrector steps) gives the final solution.

For non-conservative form no changes in the equations have been made to perform simulations but for conservative form governing equations are not solved directly instead , flux terms are written in the form of solution vectors inturn these solution vectors are function of primitive variables .

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));

j2=(1/gamma)*rho_c(k)*t_c(k)*dloga_dx;

where u1,u2,u3 are functions of primitive variables

u1=rho_c.*a;
u2=rho_c.*a.*v_c;
v_c=u2./u1;
u3=rho_c.*((t_c./(gamma-1))+(gamma/2).*v_c.^2).*a;

with the above terms governing equations in conserved form can be written as

After this setup the whole process is same as followed in non-conserved form. Finally after corrector step's solution update primitive variables are extracted from the Solution terms u1,u2,u3.

Code for conservative_function:

function [rho_c,v_c,t_c,pr_c,rho_thr,v_thr,t_thr,m_c,time_c,pr_thr]=conservative_form_quassi_1D(n,x,dx,gamma,c,dt)
clear all
close all
clc

% grid points
n=31;
x=linspace(0,3,n);
% mesh
dx=x(2)-x(1);
% area
L = 3;
a = 1 +2.2*((x-(L/2)).^2);


% initial values of primitive variables
for i = 1:n
    if x(i)>= 0 && x(i)<=L/6
        rho_c(i) = 1;
        t_c(i) = 1;
    end
    if x(i)>=L/6 && x(i)<=L/2
        rho_c(i) = 1 -0.366*(x(i)-0.5);
        t_c(i) = 1 -0.167*(x(i)-0.5);
    end
    if x(i)>=L/2
        rho_c(i) = 0.634-0.3879*(x(i)-1.5);
        t_c(i) = 0.833-0.3507*(x(i)-1.5);
    end
end
  
gamma=1.4;
c=0.5;
% defining solution vector variables
v_c = 0.59./(rho_c.*a);
u1 = rho_c.*a;
u2 = u1.*v_c;
u3 = u1.*((t_c/(gamma-1)) + (gamma*0.5*(v_c.^2)));

%throat values
thr=find(a==1);
%dt = 0.0267;

% predictor step
for z=1:1400
  
dt = min(c*(dx./(sqrt(t_c)+v_c)));
time_c=dt*1400;
  
  u1_old=u1;
  u2_old=u2;
  u3_old=u3;
  
  
  % defining flux terms
    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));
 
    for k=2:n-1
        
        df1dx(k)=(f1(k+1)-f1(k))/dx;
        df2dx(k)=(f2(k+1)-f2(k))/dx;
        df3dx(k)=(f3(k+1)-f3(k))/dx;
    
        dloga_dx=(a(k+1)-a(k))/dx;
        
        j2=(1/gamma)*rho_c(k)*t_c(k)*dloga_dx; 
   
        du1dt_p(k)=-df1dx(k);
        du2dt_p(k)=-df2dx(k)+j2;
        du3dt_p(k)=-df3dx(k);
   
        % solution update
        u1(k)=u1(k)+du1dt_p(k)*dt;
        u2(k)=u2(k)+du2dt_p(k)*dt;
        u3(k)=u3(k)+du3dt_p(k)*dt;
      end
      
    
 
 % extracting primitive variables from solution vector
 rho_c=u1./a;
 v_c=u2./u1;
 t_c=(gamma-1)*((u3./u1)-(gamma/2)*(v_c.^2));
 pr=rho_c.*t_c;

 
  % defining flux terms
    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 step
    for k=2:n-1
        
        df1dx(k)=(f1(k)-f1(k-1))/dx;
        df2dx(k)=(f2(k)-f2(k-1))/dx;
        df3dx(k)=(f3(k)-f3(k-1))/dx;
   
        dloga_dx=(a(k)-a(k-1))/dx;
        j2=(1/gamma)*rho_c(k)*t_c(k)*dloga_dx; 
        
        du1dt_c(k)=-df1dx(k);
        du2dt_c(k)=-df2dx(k)+j2;
        du3dt_c(k)=-df3dx(k);
      end
      
   
          %finding average values
          du1dt_avg=0.5*(du1dt_p+du1dt_c);
          du2dt_avg=0.5*(du2dt_p+du2dt_c);
          du3dt_avg=0.5*(du3dt_p+du3dt_c);
for j=2:n-1      
         % solution update  
         u1(j)=u1_old(j)+du1dt_avg(j)*dt;
         u2(j)=u2_old(j)+du2dt_avg(j)*dt;
         u3(j)=u3_old(j)+du3dt_avg(j)*dt;
 end
 
   
 % boundary conditions
% inlet

% extracting primitive variables 


rho_c = u1./a;
v_c = u2./u1;
t_c = (gamma-1)*((u3./u1) - ((gamma/2)*(v_c.^2)));
pr_c=rho_c.*t_c;


    
% Boundary Nodes
% Inlet
u1(1) = rho_c(1)*a(1);
u2(1) = 2*u2(2) - u2(3);
u3(1) = u1(1)*((t_c(1)/(gamma-1))+(0.5*gamma*(v_c(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);

% mass flow rate
m_c=(rho_c.*a.*v_c);
 
 % throat values for different time step
 rho_thr(z)=rho_c(thr);
 v_thr(z)=v_c(thr);
 t_thr(z)=t_c(thr);
 pr_thr(z)=pr(thr);
 
 
 %plot for mass flow rate
    figure(4)
   if z==1
    hold on
     plot(x,m_c,'r','linewidth',3)
     hold on 
   elseif z==50
     plot(x,m_c,'y','linewidth',3)
     hold  on
           elseif z==100
     plot(x,m_c,'g','linewidth',3)
     hold on
                   elseif z==200
     plot(x,m_c,'c','linewidth',3)
     hold on
                          elseif z==700
     plot(x,m_c,'--b')
    
     hold on
     xlabel('time steps')

     ylabel('mass flow rate ')
     title('mass flow rate--time wise variation along the length of the nozzle')
     legend('0dt','50dt','100dt','200dt','700dt')
     
   end  
  

  end

endfunction

  

Code for non- conservative function:

function [rho,v,t,m_n,pr,a,rho_th,v_th,t_th,pr_th,time_n]=Non_conservative_form(x,n,gamma,dx,c);
clear all
close all
clc

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

% defining input variables
rho= 1-0.3146.*x;
t=1-0.2314.*x;
v=(0.1 + 1.09.*x).*t.^0.5;

% area
a= 1+2.2.*(x - 1.5).^2;

gamma=1.4
c=0.5
nt=1400
% throat value
th=find(a==1)
for i=1:nt
dt=min(c.*(dx./(sqrt(t)+v)));
time_n=dt*nt;
 
rho_old=rho;
v_old=v;
t_old=t;

for j=2:n-1
  % Non-conservative form
   %predictor step
   drhodx=(rho(j+1)-rho(j))/dx;
   dvdx=(v(j+1)-v(j))/dx;
   dtdx=(t(j+1)-t(j))/dx;
   dlog_a_dx=(log(a(j+1))-log(a(j)))/dx;
   
   drhodt_p(j)=-rho(j)*(dvdx)-rho(j)*v(j)*(dlog_a_dx)-v(j)*drhodx;
   dvdt_p(j)=-v(j)*(dvdx)-(1/gamma)*(dtdx+(t(j)/rho(j))*drhodx);
   dtdt_p(j)=-v(j)*(dtdx)-(gamma-1)*t(j)*(dvdx+v(j)*dlog_a_dx);
   
   %solution update 
   rho(j)=rho(j)+drhodt_p(j)*dt;
   v(j)=v(j)+dvdt_p(j)*dt;
   t(j)=t(j)+dtdt_p(j)*dt;
end
    
 for j=2:n-1 
   %corrector step
   drhodx=(rho(j)-rho(j-1))/dx;
   dvdx=(v(j)-v(j-1))/dx;
   dtdx=(t(j)-t(j-1))/dx;
   dlog_a_dx=(log(a(j))-log(a(j-1)))/dx;
   
   drhodt_c(j)=-rho(j)*(dvdx)-rho(j)*v(j)*(dlog_a_dx)-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)*dlog_a_dx);
 end
 
 
   %computing avg values
   drhodt_avg=0.5*(drhodt_p+drhodt_c);
   dvdt_avg=0.5*(dvdt_p+dvdt_c);
   dtdt_avg=0.5*(dtdt_p+dtdt_c);
   
    % final solution update
 for k=2:n-1
       rho(k)=rho_old(k)+drhodt_avg(k)*dt;
       v(k)=v_old(k)+dvdt_avg(k)*dt;
       t(k)=t_old(k)+dtdt_avg(k)*dt;
  end   
       
  %appying boundary conditions
  % inlet
  v(1)=2*v(2)-v(3);

  %outlet
  rho(n)=2*rho(n-1)-rho(n-2);
  v(n)=2*v(n-1)-v(n-2);
  t(n)=2*t(n-1)-t(n-2);
  
  % mass flow rate
  m_n=rho.*a.*v;
  % pressure ratio
  pr=rho.*t;
  
  % thaoat values 
  rho_th(i)=rho(th);
  v_th(i)=v(th);
  t_th(i)=t(th);
  pr_th(i)=pr(th);
  
  %plot for mass flow rate
   figure(2)
  if i==1
   hold on
    plot(x,m_n,'r','linewidth',3)
    hold on 
  elseif i==50
    plot(x,m_n,'b','linewidth',3)
    hold  on
          elseif i==100
    plot(x,m_n,'g','linewidth',3)
    hold on
                  elseif i==200
    plot(x,m_n,'c','linewidth',3)
    hold on
                         elseif i==700
    plot(x,m_n,'--b','linewidth',3)
    hold on
                              
    xlabel('time steps')
    ylabel('mass flow rate ')
    title('mass flow rate--time wise variation along the length of nozzle')
    legend('0dt','50dt','100dt','200dt','700dt')
    title('Non -conservative mass flow rate at different time steps ')
    
                     
  end
    
end
 
endfunction
 plot(x,m_n)
 

 

Main code:

clear all
close all
clc
  
 n=31;
 nt=1400
 c=0.5;
 x=linspace(0,3,n);
 dx=x(2)-x(1);
 gamma=1.4
 % non conseravtive form
 [rho,v,t,m_n,pr,a,rho_th,v_th,t_th,pr_th,time_n]=Non_conservative_form(x,n,gamma,dx,c);
 
 
 %conservative form
 [rho_c,v_c,t_c,pr_c,rho_thr,v_thr,t_thr,m_c,time_c,pr_thr]=conservative_form_quassi_1D(n,x,dx,gamma,c)

 % primitive variables at steady state non-conservative form
 figure(1)
 
 %density
 subplot(4,1,1)
 plot(x,rho)
 xlabel('length of the nozzle')
 ylabel('density')
 title('primitive variables at steady state')
 
 %velocity
 subplot(4,1,2)
 plot(x,v)
 xlabel('length of the nozzle')
 ylabel('velocity')
 
 %temperature
 subplot(4,1,3)
 plot(x,t)
 xlabel('length of the nozzle')
 ylabel('temperature')
 
 %pressure ratio
 subplot(4,1,4)
 plot(x,pr)
 xlabel('length of the nozzle')
 ylabel('pressure ratio')
 
  % primitive variables at steady state conservative form
  figure(5)
 
 %density
 subplot(4,1,1)
 plot(x,rho_c)
 xlabel('length of the nozzle')
 ylabel('density')
 title('primitive variables at steady state')
 
 %velocity
 subplot(4,1,2)
 plot(x,v_c)
 xlabel('length of the nozzle')
 ylabel('velocity')
 
 %temperature
 subplot(4,1,3)
 plot(x,t_c)
 xlabel('length of the nozzle')
 ylabel('temperature')
 
 %pressure ratio
 subplot(4,1,4)
 plot(x,pr_c)
 xlabel('length of the nozzle')
 ylabel('pressure ratio')

 
 % variation of primitive variables with respect to time
 
 figure(6)
 
 %density
 subplot(4,1,1)
 plot(rho_th)
 xlabel('Time steps')
 ylabel('density')
 title('variation of primitive variables- time marching(non-conservative)')
 
 %velocity
 subplot(4,1,2)
 plot(v_th)
 xlabel('Time steps')
 ylabel('velocity')
 
 %temperature
 subplot(4,1,3)
 plot(t_th)
 xlabel('Time steps')
 ylabel('temperature')
 
 %pressure ratio
 subplot(4,1,4)
 plot(pr_th)
 xlabel('rime steps')
 ylabel('pressure ratio')

 
 % variation of primitive variables with respect to time conservative
 figure(3)
 
 %density
 subplot(4,1,1)
 plot(rho_thr)
 xlabel('Time steps')
 ylabel('density')
 title('variation of primitive variables- time marching(conservative)')
 
 %velocity
 subplot(4,1,2)
 plot(v_thr)
 xlabel('Time steps')
 ylabel('velocity')
 
 %temperature
 subplot(4,1,3)
 plot(t_thr)
 xlabel('Time steps')
 ylabel('temperature')
 
 %pressure ratio
 subplot(4,1,4)
 plot(pr_thr)
 xlabel('Time steps')
 ylabel('pressure ratio')
 
 % normalized mass flow rate
 figure(7)
 plot(x,m_n)
 hold on 
 plot(x,m_c)
 xlabel('length 0f nozzle')
 ylabel('mass flow')
title('normalized mass flow rate')
legend('non-conservative','conservative')


 

Results and Discussion:

Distribution of primitive variables through the length of the nozzle at steady state.( Non-conservative form):

Distribution of primitive variables through the length of the nozzle at steady state. (conservative form):

from the above plots ,all the primitive variables except velocity decreases gradually in the form of hyperbolic curves. Inlet conditions are assumed to be a reservoir where temperature and density are assumed to be unity and constant at inlet and velocity is allowed to float and both the plots are alike.

Variation of primitive variables with respect to time (time marching ) non-conservative:

Variation of primitive variables with respect to time (time marching ) conservative:

Above plot says that ,at initial time steps variables were oscillating and are unsteady but got settled after 600 time steps, non-conservative form has less oscillations when compared to conservative form .

These initial oscillations are due to a phenomena well explained with a 'rubber band' example, when a rubber band is stretched to its maximum elastic limit and released it goes most compression in initial stages and very less at later stages the same effect is seen here. Initially due to some sudden contractions and expansions oscillations in the variables is seen.

Mass flow along the length of the nozzle at different time steps non-conservative form :

Number of grids n=31

n=61

The above plots are for mass flow along the length of the nozzle ,initial curve which is coloured red is for initial conditions assumed hence we named it as '0dt' which is in the form of a 'sine' wave. Initial conditions are assumed in such a way that density and velocity are hyperbolic in nature multiplied with area which is parabolic.

In time marching ,mass flow slowly started to converge and moving towards stable conditions and at 200 time steps it is close to constant mass flow rate but after reaching 700 time steps it is constant through out the length of the nozzle exactly equal to the analytical value 0.579 calculated at throat using primitive variables .

But in real situations it is never going to happen as we have fixed density and temperature to be constant resulting in velocity very close to zero but not actually.
second plot refers to mass flow rate for 61 grid points which is discussed in grid independency section.

Mass flow along the length of the nozzle at different time steps conservative form :

N=31

n=61

Plots for mass flow along the length for conserved form shows that , the initial conditions are chosen in such a way which has very less instabilities as the time marching continues instabilities develop in the flow and after reaching 700 time steps it converges to a constant value which is again 0.579 ,calculated at throat.

Normalized mass flow:

N=31

N=61

The above plots shows the advantage of conservative form of equations over non-conservative form in preserving mass in the flow filed which does not conclude that conservative form of equations gives better results. Better preservation of mass in conserved form is because mass flow itself is a dependent variable in the equations where as in non conservative form primary variables are dependent variables and mass flow is obtained from them secondarily. By this we can justify that why conservative form equations are labelled so.

Grid independency test:

FORM	density ratio	temperature ratio	pressure ratio	velocity ratio
Exact solutions calculated at throat	0.634	0.833	0.528	1.000
Non-conservative(n=61)	0.637	0.8359	0.5326	0.9137
conservative(n=61)	0.638	0.8351	0.533	0.9132
Non-conservative(n=31)	0.638	0.836	0.534	0.914
conservative(n=31)	0.650	0.840	0.546	0.899

Grid independency test:

Grid independence is very seriously considered in CFD. In general we simulate solutions for n grid points over the flow field until you reach a solution, if you double the grid points to 2n then you might get an different solution with different values it is because you are decreasing the size of grid or mesh by increasing the grid points and you have to continue it till you reach a solution which is insensitive to grid size .

In our case by increasing the grid points to 61 from initial size of 31 grids we have obtained a solutions which are tabulated above. we haven't achieved better results even after doubling the grid points .Hence our simulation is Grid independent.

conclusion:

1.From the above results and discussion, we can say that there is no superiority of conservative form over non-conservative form of equations, but conservative form gives better mass conservation and in the non-conservative form primitive variables are preserved well.

2.Complexity of solving conservative form is slightly more when compared to non-conservative form because in conservative form primitive variables are not directly solved but have to be decoded from the flux terms.