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

AIM : 

 Simulation of isentropic flow through a quasi 1D subsonic-supersonic nozzle using  MacCormack's technique and to determine the steady-state temperature distribution for the flow-field variables and investigate the difference between the two forms of governing equations(conservative and non-conservative ) by comparing their solutions.

OBJECTIVE :

• To find the required plots of primitive variables 
• Perform Grid independency test

   

• To find out which form of governing equations results in faster convergence.
•  To compare Normalized mass flow rate distributions of both forms with exact solution and to find the time taken among both of them, to see which one is faster 

THEORY :

The problem along with its initial and boundary conditions has been taken from the book "Computational Fluid Dynamics - The Basics with Applications by John D Anderson" , Chapter 7

We consider a steady, isentropic flow through a convergent-divergent nozzle. The flow at the inlet to the nozzle comes from a reservoir with a large cross-sectional area ( A tends to infinity), and hence the velocity is very small ( V tends to 0). The flow expands isentropically to supersonic speed at the nozzle exit. It is subsonic in the convergent section, sonic at the throat and supersonic in the divergent section. The flow properties are assumed to vary along only the x-direction although the area is two-dimensional. Such flow is called quasi-one-dimensional flow.

Governing equations in non-conservative form:

Continuity Equation:

`A(delrho)/(delt)+V(del(rhoA))/(delx)=-rhoA(delV)/(delx)`

Momentum Equation:

`rho(delV)/(delt)+rhoV(delV)/(delx)=-(delp)/(delx)`

Energy Equation:

`rho(dele)/(delt)+rhoV(dele)/(delx)=-p(delV)/(delx)-pV(del(lnA))/(delx)`

Generally, for nozzle flows, the flow field variables are expressed as non-dimensional variables. This is done by dividing all the dimensional flow variables by their respective reservoir values.

`implies`  the non-dimensional flow variables are as follows:

`x^'=x/L , rho^'=rho/rho_0, T^'=T/T_0, A^'=A/(A^ast),a_0=sqrt(gammaRT_0),V^'=V/a_0,t^'=t/(L/a_0)`

where, 

`rho_0,T_0,a_0` are the density, temperature and speed of sound at the reservoir.

`L,A^ast` are the nozzle length and throat area respectively.

We have to now simplify and express the above non-conservative form using non-dimensional variables.

Governing equations in Conservative form:

Continuity Equation:

`(delrhoA)/(delt)+(del(rhoAV))/(delx)=0`

Momentum Equation:

`(del(rhoAV))/(delt)+(delrhoAV^2)/(delx)=-A(delp)/(delx)`

Energy Equation:

`(del[rho(e+V^2/2)A])/(delt)+(del[rho(e+V^2/2)VA])/(delx)=-(del(pAV))/(delx)`

Again we have to now simplify and express the above conservative form using non-dimensional variables.

Let U be the solutions vector, F be the flux vector and J be the source term, such that

`U_1=rho^'A^'`

`U_2=rho^'A^'V^'`

`U_3=rho^'(e^'/(gamma-1)+gamma/2*(V^')^2)A^'`

`F_1=rho^'A^'V^'`

`F_2=rho^'A^'(V^')^2+1/gammap^'A^'`

`F_3=rho^'(e^'/(gamma-1)+gamma/2*(V^')^2)V^'A^'+p^'A^'V^'`

`J_2=1/gammap^'(delA^')/(delx^')`

Backsubstitute these values into the equations to non-dimensional conservative form(pure form)

Finally to get primitive variables(`rho,V,T,p`) from these dependent variables in numerical solution(as our answer will be directly in terms of `U_1,U_2,U_3` in steps of time,hence called solution vector), we have to decode them as follows:

`rho^'=U_1/A^' , V^'=U_2/U_1 , T^'=e^'=(gamma-1)*(U_3/U_2-gamma/2*(V^')^2) , p^'=rho^'T^'`

Discretization:

MacCormack technique has been employed in both the governing equations .It is a explicit finite-difference technique which is second-order accuratein both space and time. It uses predictor-corrector setup, in which the predictor step computes the first order time derivatives using the forward difference for spatial derivatives by substituting previous state values of the flow variables. It then updates the flow field variables. These updated values are then used to compute corrector step time derivatives using the backward difference approach. Finally, the corrected solutions for the flow variables are obtained by using the average value of time derivatives computed in the predictor and the corrector steps. These steps are repeated until a steady state is reached or until predefined time steps.

It is relatively easier to apply as compared to the Lax-Wendroff appproach since there is no need to evaluate 2nd time derivative. 

CODE :

clear all
close all
clc
%main program
%inputs
L=3;
n=31;
x=linspace(0,3,n);
dx=x(2)-x(1);
gamma=1.4;
throat=(n-1)/2;
a=1+2.2*(x-1.5).^2;
c=0.5;
%time steps
nt=1400;
m_analytical = 0.579*ones(n);

tic
% solving governing equations in the non-conservation form
[m]=non_cons_form(x,dx,n,c,a,nt,gamma,throat);

% solution by conservative form :
tic
[m_c]= cons_form(x,dx,n,c,a,nt,gamma,throat);

% comparison of mass flow rates by conservative and non_cons method and comparison with exact form:
figure(7)
hold on
plot(x,m,'b');
plot(x,m_c,'r');
plot(x,m_analytical,'y');
ylabel('m^.');
xlabel('Domain Length (L)')
legend('non-conservative form (m^.)','conservative form (m^.)','Exact solution');
title('Mass flowrate at 1D supersonic nozzle flow:(61 gridpoints)');


%function used in main_program_macormack_method

function [m]=non_cons_form(x,dx,n,c,a,nt,gamma,throat)

%initialize
%calculate initial conditions
rho=1-0.3146*x;
t=1-0.2314*x; %temperature
v=(0.1+1.09*x).*t.^0.5;

%outer time loop
for k=1:nt;
    
    rho_old=rho;
    v_old=v;
    t_old=t;
    %dt_array=c.*dx./(sqrt(t)+v);
   dt=min(c.*dx./(sqrt(t)+v));
    
    
    %predictor method
    
    for j=2:n-1
        %drho_dt=-rho(j)*(v(j+1)-v(j))/dx -rho(j)*v(j)*(log(a(j+1))-log(a(j)))/dx -v(j)*(rho(j+1)-rho(j))/dx;
        dvdx=(v(j+1)-v(j))/dx;
        drhodx=(rho(j+1)-rho(j))/dx;
        dtdx=(t(j)-t(j-1))/dx;
        dlogadx=(log(a(j+1))-log(a(j)))/dx;
        
        %continuity equation
        drhodt_p(j)=-rho(j)*dvdx-rho(j)*v(j)*dlogadx-v(j)*drhodx;
        
        %momentum equation
        dvdt_p(j)=-v(j)*dvdx-(1/gamma)*(dtdx+(t(j)/rho(j))*drhodx);
        
        %energy equation
        dtdt_p(j)=-v(j)*dtdx-(gamma-1)*t(j)*(dvdx +v(j)*dlogadx); 
    
        %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
    
     %corrector method
 
    for j=2:n-1
        %drho_dt=-rho(j)*(v(j+1)-v(j))/dx -rho(j)*v(j)*(log(a(j+1))-log(a(j)))/dx -v(j)*(rho(j+1)-rho(j))/dx;
        dvdx=(v(j)-v(j-1))/dx;
        drhodx=(rho(j)-rho(j-1))/dx;
        dlogadx=(log(a(j))-log(a(j-1)))/dx;
        
        %continuity equation
        drhodt_c(j)=-rho(j)*dvdx-rho(j)*v(j)*dlogadx-v(j)*drhodx;
        
        %momentum equation
        dvdt_c(j)=-v(j)*dvdx-(1/gamma)*(dtdx+(t(j)/rho(j))*drhodx);
        
        %energy equation
        dtdt_c(j)=-v(j)*dtdx-(gamma-1)*t(j)*(dvdx +v(j)*dlogadx); 
    
    end
    %compute the average time derivative
    
    drhodt=0.5*(drhodt_p+drhodt_c);
    dvdt=0.5*(dvdt_p+dvdt_c);
    dtdt=0.5*(dtdt_p+dtdt_c);
    
    %final solution update
    
    for i=2:n-1
    rho(i)=rho_old(i)+drhodt(i)*dt;
    v(i)=v_old(i)+dvdt(i)*dt;
    t(i)=t_old(i)+dtdt(i)*dt; 
    end

    %apply the boundary conditions
    %inlet
    v(1)=2*v(2)-v(3);
    
    %outlet
    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);
    
    % defining the values for press, mass flow rate and mac no.:
p=rho.*t; % non-dimensional press
m=rho.*a.*v; % non-dimensional mass flow rate
M=v./sqrt(t); % mach no.
mass_flowrate(k,:)=rho.*a.*v;
% calculation of quantities at the throat point of the nozzle
t_throat(k)=t(throat);
p_throat(k)=p(throat);
rho_throat(k)=rho(throat);
M_throat(k)=M(throat);
    
end
time_elapsed_non_conservative=toc

figure(1)

subplot(4,1,1);
plot(x,rho,'r');
ylabel('rho');
legend('density');
title('Qualitative aspects of 1D supersonic nozzle flow: non-conservative form (dt=0.0224) adaptive timestep');
subplot(4,1,2);
plot(x,t,'g');
ylabel('T');
legend('temperature');

subplot(4,1,3);
plot(x,M,'c');
ylabel('M');
xlabel('Domain Length (L)')
legend('Mach no.');


subplot(4,1,4);
plot(x,p,'m');
ylabel('P');
legend('pressure');
xlabel('Domain Length (L)')

 

figure(2)
plot(x,mass_flowrate(1,:),'--',x,mass_flowrate(50,:),x,mass_flowrate(100,:),x,mass_flowrate(150,:),x,mass_flowrate(200,:),x,mass_flowrate(700,:),x,mass_flowrate(1400,:),':')

xlabel('Domain Length (L)')
ylabel('m^.');
title('variation of mass flowrate at different timestep(non-conservative)');
legend('1st iteration','50th iteration','100th iteration','150th iteration','200th iteration','700th iteration','1400th iteration')

figure(3)

subplot(4,1,1);
plot(rho_throat,'r');
ylabel('rho');
legend('density');
title('Time-wise variation of the primitive variables at nozzle throat: non-conservative form (dt=0.0224) adaptive timestep');
subplot(4,1,2);
plot(t_throat,'g');
ylabel('T');
legend('temperature');

subplot(4,1,3);
plot(M_throat,'c');
ylabel('M');
legend('Mach no.');

subplot(4,1,4);
plot(p_throat,'m');
ylabel('P');
legend('pressure');
xlabel('No. of iterations')
end

%function used in main_program_macormack_method

function [m_c]=cons_form(x,dx,n,c,a,nt,gamma,throat)
  
  % initialize step
  % initializing density and temp
  for i=1:length(x)
if x(i)>=0 && x(i)<0.5
  rho_c(i)=1;
  temp_c(i)=1;
else if x(i)>=0.5 && x(i)<1.5
  rho_c(i)=1-0.366*(x(i)-0.5);
  temp_c(i)=1-0.167*(x(i)-0.5);
else if x(i)>=1.5 && x(i)<=3
  rho_c(i)=0.634-0.3879*(x(i)-1.5);
  temp_c(i)=0.833-0.3507*(x(i)-1.5);
end
end
end
end
% initialize vel and press
v_c=0.59./(rho_c.*a);
p_c=rho_c.*temp_c;

% initialise solution vectors
U1=rho_c.*a;
U2=rho_c.*a.*v_c;
U3=rho_c.*((temp_c./(gamma-1))+(gamma/2)*v_c.^2).*a; 

% time iteration loop
for k=1:nt
  
  % presurving the old solution vectors
  U1_old=U1;
  U2_old=U2;
  U3_old=U3;
  
  % adaptive time step control
  dt=min(c.*dx./(sqrt(temp_c)+v_c));
  
  % defining the flux vectors
  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);
  
  % solution by macormack method
  for j=2:n-1
    % defining source term
    J2(j)=(1/gamma).*rho_c(j).*temp_c(j).*((a(j+1)-a(j))/dx);
    %predictor step
    dU1dt_p(j)=-(F1(j+1)-F1(j))/dx;
    dU2dt_p(j)=-(F2(j+1)-F2(j))/dx + J2(j);
    dU3dt_p(j)=-(F3(j+1)-F3(j))/dx;
    % updating the new values of solution vectors
    U1(j)=U1(j)+dU1dt_p(j)*dt;
    U2(j)=U2(j)+dU2dt_p(j)*dt;
    U3(j)=U3(j)+dU3dt_p(j)*dt;
  end
  % updating the primitive variables and flux vectors
  rho_c=U1./a;
  v_c=U2./U1;
  temp_c=(gamma-1)*(U3./U1-(gamma/2)*v_c.^2);
  p_c=rho_c.*temp_c;
  
  F1=U2;
  F2=(U2.^2./U1)+((gamma-1)/gamma)*(U3-(gamma*U2.^2)./(2*U1));
  F3=(gamma*U2.*U3./U1)-(gamma*(gamma-1)/2)*(U2.^3./U1.^2);
  
  for j=2:n-1
    % updating the source term
    J2(j)=(1/gamma)*rho_c(j).*temp_c(j)*((a(j)-a(j-1))/dx);
    
    % corrector step
    dU1dt_c(j)=-(F1(j)-F1(j-1))/dx;
    dU2dt_c(j)=-((F2(j)-F2(j-1))/dx)+J2(j);
    dU3dt_c(j)=-(F3(j)-F3(j-1))/dx;
  end
  % updating the solution vectors after corrective step
  % averaging the predictive and corrector step
  dU1dt=0.5*(dU1dt_p+dU1dt_c);
  dU2dt=0.5*(dU2dt_p+dU2dt_c);
  dU3dt=0.5*(dU3dt_p+dU3dt_c);
  
  % calculating new solution vectors
  for j=2:n-1
    U1(j)=U1_old(j)+dU1dt(j)*dt;
    U2(j)=U2_old(j)+dU2dt(j)*dt;
    U3(j)=U3_old(j)+dU3dt(j)*dt;
  end
  % apply BC
  % left side/inlet conditions
  U1(1)=rho_c(1).*a(1); 
  U2(1)=2*U2(2)-U2(3);
  U3(1)=U1(1)*(temp_c(1)/(gamma-1)+(gamma/2)*v_c(1).^2);
  
  % right side/outlet 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 solution of the primitive variables
  rho_c=U1./a;
  v_c=U2./U1;
  temp_c=(gamma-1)*(U3./U1-(gamma/2)*v_c.^2);
  p_c=rho_c.*temp_c;
  % define the values of mach no. and mass flow rate
  M_c=v_c./sqrt(temp_c);
  m_c=rho_c.*a.*v_c;
  
  % calculating values of specific quantities at throat
  M_throat_c(k)=M_c(throat);   % mach no
  p_throat_c(k)=p_c(throat);   % press ratio
  temp_throat_c(k)=temp_c(throat);  % non-dimensional temp
  rho_throat_c(k)=rho_c(throat);   % density ratio
  mass_flowrate(k,:)=rho_c.*a.*v_c;
 
  
end
 time_elapsed_conservative=toc



% qualitative aspects of 1D flow in conservative form:
figure(4);

subplot(4,1,1);
plot(x,rho_c,'g');
title('Qualitative aspects of 1D supersonic nozzle flow: cons form(dt=0.0207) adaptive timestep');
ylabel('rho');
legend('density');

subplot(4,1,2);
plot(x,temp_c,'m');
ylabel('T');
legend('temperature');

subplot(4,1,3)
plot(x,M_c,'b');
ylabel('M');
legend('Mach no');

subplot(4,1,4);
plot(x,p_c,'r');
ylabel('p');
legend('pressure');
xlabel('Domain Length (L)')



% comparison of non dimensional mass flow rates at different times
figure(5)
plot(x,mass_flowrate(1,:),'--',x,mass_flowrate(50,:),x,mass_flowrate(100,:),x,mass_flowrate(150,:),x,mass_flowrate(200,:),x,mass_flowrate(700,:),x,mass_flowrate(1400,:),':')

xlabel('Domain Length (L)')
ylabel('m^.');
title('variation of mass flowrate at different timestep(conservative form)');
legend('1st iteration','50th iteration','100th iteration','150th iteration','200th iteration','700th iteration','1400th iteration')

% time wise variation of density,temp,press and mac no. for the cons form:
figure(6)

subplot(4,1,1)
plot(rho_throat_c,'b')
title('Time-wise variation of the primitive variables at nozzle throat:conservative form ');
legend('density');
ylabel('rho');

subplot(4,1,2);
plot(temp_throat_c,'r');
legend('temperature');
ylabel('T');

subplot(4,1,3);
plot(M_throat_c,'m');
legend('mach no');
ylabel('M');
axis([0 1400 0.6 1.4]);

subplot(4,1,4);
plot(p_throat_c,'g');
legend('pressure');
ylabel('P');
xlabel('No. of iterations')

end

OUTPUT:

time_elapsed_non_conservative =

    0.0927


time_elapsed_conservative =

    0.1019

INFERENCE :

We can see that by 700th iteration in non-conservative mass flow rate, it has converged into steady state result (faster one)(even though totally 1400 iterations are present)during time-marching process, but the conservation form yields better mass flow distribution .

The time taken for conservative form( 0.1019 sec) is more than non-conservative form(0.0927 sec)(faster method) ,this could be attributed to the fact that conservative form requires more work(calculations) as we need to decode the primitive variables from the flux variables(dependent variables) , such decoding is not required in  non-conservative form.

Grid Independence Test has been performed with 61 gridpoints plot(last pic of mass flowrate comparison),(both forms are grid independent) we can see that even though the no. of gridpoints has been doubled there is only marginal  improvement in accuracy of answer to exact solution , hence we can say this(61 nodes) is largely a time consumption which can be avoided (as usually more time translates to increased cost/money -unless accuracy is the most important factor )  .Generally speaking truncation error is the root cause of Grid independence

The main noticeable difference among the 2 forms of governing equations is in normalized mass flow rate  distribution plot where the conservative form has an edge ove the other.

Overall there is no superiority of either forms in terms of accuracy of results as non-conservation form produces slightly better accuracy of results for primitive variables , whereas conservation form produces slightly better accuracy of results for flux variables, although are satisfactory enough.

Non-conservative form leads to smaller residuals ,the amount by which residual decay is often used as an index of "quality" of numerical algorithm. In this sense it is slightly better(and lesser time too as a conclusion).