Simulation Of A 1D Super-sonic Nozzle Using Macormack Method

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

Objective:

Determine the steady-state temperature distribution for the flow field variables and investigate the difference between the two forms of governing equation 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. Comparison of Normalized mass flow rate distributions of both forms.

Theory:

Any numerical solution of isentropic flow through a quasi 1D subsonic-supersonic nozzle is unreasonable, while an analytical solution exists for these kinds of problems. However, we are trying here is to establish the broadened capacity 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 marched in terms of PDE for both forms to use the Time Marching solution for the problem. These forms represent inn numerical aspect only, if we analytic all the methods are the same.

Non-Conservative Form:

Conservative Form:

Non-Dimensional Quantities:

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

Nozzle diagram with sections:

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

• Convergent Section
• Divergent Section

MacCormack Method:

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

The application of the 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 u at the time level n+1 (denoted by `u_i^(n+1)`) is estimated obtained by replacing the spatial and temporal derivatives in the previous first-order equation using Forward Difference Scheme.`Deltat` is used for Timestep.
• Corrector step: In the corrector step, the predicted value `u_i^(n+1)` is corrected according to the equation. Note that the corrector step uses the Backward Difference Scheme for spatial derivatives. The time-step used in the corrector step is `Deltat`/2 in contrast to the `Deltat` used in the predictor step.

MAIN PROGRAM:

For the code to be modular, the program is split into three sections. They are

• The main program of Quasi 1D subsonic-supersonic nozzle.

                Two functions namely are

• Non Conservative form
• Conservative form.

CODE

MAIN CODE

clear all
close all
clc

% Inputs
n = 31;
x=linspace(0,3,n);
dx=x(2)-x(1);
gamma=1.4;

% Data for the CFL Time Loop and Convergence Controller
C=0.5;   % Courant Number
error=1e9;
error2=1e9;
tol=1e-5;

%Area Profile
A=1+2.2*(x-1.5).^2;
n_throat=find(A-1<0.001); % THe value where A/A* is one(or close)

%--------------------NON-CONSERVATIVE FORM SOLVER---------------
tic 
[rho,v,T,p,M,massflow,rho_throat,T_throat,p_throat,M_throat,k]=non_conservative(n,x,dx,gamma,A,n_throat,C,error,tol);
time_nc =toc;
fprintf('Minimum Number Of Iteration for Convergence(NonConservative Form): %fn',k);
fprintf('Computational Time(NonConservative Form): %f sn', time_nc);

% Transient Properties At The Throat

% Transient Variation of the throat properties
Transient_Throat=figure(2);
subplot(4,1,1);
plot(rho_throat,'m');
title('Timewise variation of the properties at the nozzle throat(Non-Conservative Form)');
legend('Density ratio');
ylabel('rho/rho_{0}');
axis([0 k 0.5 1]);
grid minor;

subplot(4,1,2);
plot(T_throat,'r');
legend('Temperature');
ylabel('T/T_{0}');
axis([0 k 0.6 1.2]);
grid minor;

subplot(4,1,3);
plot(p_throat,'g');
legend('Pressure');
ylabel('p/p_{0}');
axis([0 k 0.4 1]);
grid minor;

subplot(4,1,4);
plot(M_throat,'b');
legend('Mach Number');
ylabel('M');
axis([0 k 0.8 1.3]);
grid minor;
xlabel('Number of iteration');


%STAEDY PROFILES


%Steady profile of the properties at the nozzle

figure(3)
subplot(4,1,1);
plot(x,rho,'m');
title('Properties profile in the 1D Supersonic Nozzle Flow(Non-Conservative Form)'); 
ylabel('rho/rho_{0}');
legend('Density ratio','location','southwest');
grid minor;

subplot(4,1,2);
plot(x,T,'r');
ylabel('T/T_{0}');
legend('Temperature ratio','location','southwest');
grid minor;

subplot(4,1,3);
plot(x,p,'g');
ylabel('p/p_{0}');
legend('pressure ratio','location','southwest');
grid minor;

subplot(4,1,4);
plot(x,M,'b');
ylabel('M');
legend('Mach Number','Location','northwest');
grid minor;


%--------------------------CONSERVATIVE FORM SOLVER------------------------

tic

[rho_con,v_con,T_con,p_con,M_con,massflow_con,rho_throat_con,T_throat_con,p_throat_con,M_throat_con,k2]= conservative(n,x,dx,gamma,A,n_throat,C,error2,tol);
time_c = toc;
fprintf('Minimum Number of Iteration for Convergence (Conservative Form): %fn',k2);
fprintf('Computation Time (Conservative Form): %f sn',time_c);


%TRANSIENT PROPERTIES AT THE THROAT

%Transeint varitation of the throat properties

Transient_Throat_con=figure(5);
subplot(4,1,1);
plot(rho_throat_con,'m');
title('Timewise variation of the properties at the nozzle throat (Conservative Form)');
legend('Density ratio');
ylabel('rho/rho_{0}');
axis([0 k2 0.5 1]);
grid minor;

subplot(4,1,2);
plot(T_throat_con,'r');
legend('Temperature');
ylabel('T/T_{0}');
axis([0 k2 0.6 1.2]);
grid minor;

subplot(4,1,3);
plot(p_throat_con,'g');
legend('pressure');
ylabel('p/p_{0}');
axis([0 k2 0.4 1]);
grid minor;

subplot(4,1,4);
plot(M_throat_con,'b');
legend('Mach number');
ylabel('M');
axis([0 k2 0.8 1.3]);
grid minor;
xlabel('Number of iterations');


%  STEADY PROFILES

%Steady profile of the properties at the nozzle 

figure(6)
subplot(4,1,1);
plot(x,rho_con,'m');
title('Properties Profile in the 1D supersonic nozzle flow(Conservative Form)');
ylabel('rho/rho_{0}');
legend('Density ratio','location','southwest');
grid minor;


subplot(4,1,2);
plot(x,T_con,'r');
ylabel('T/T_{0}');
legend('Temperature ratio', 'location','southwest');
grid minor;

subplot(4,1,3);
plot(x,p_con,'g');
ylabel('p/p_{0}');
legend('pressure ratio','location','southwest');
grid minor;

subplot(4,1,4);
plot(x,M_con,'b');
ylabel('M');
legend('Mach Number','Location','northwest');
grid minor;

% Steady massflow rate profile(Comparisson Non-Conservation vs conservative)

figure(7)
plot(x,massflow)
hold on
plot(x,massflow_con)
title('Steady Mass-Flow Profile')
legend('Non-Conservative','Conservative')
grid minor

Non-Conservation Solver:

function [rho,v,T,p,M,massflow,rho_throat,T_throat,p_throat,M_throat,k]=non_conservative(n,x,dx,gamma,A,n_throat,C,error,tol)


% Non-Conservative Approach

% Calculate initial Profiles
rho =1-0.3146*x;          %Density
T = 1-0.2314*x;           %Temperature
v= (0.1+1.09*x).*T.^0.5;  %Velocity

%Initializing for the time loop
k=0;

%outer Time loop

while error>tol
    k=k+1;
    
    rho_old = rho;
    v_old= v;
    T_old = T;
    a_old=sqrt(T_old);
    massflow_old=rho_old.*v_old.*A;
    
    %CFL time loop control
   TimeStep=C*dx./(a_old+v_old);
   dt=min(TimeStep);
   
   
    % Predictor Step
    for j = 2:n-1
  %      drho_dt_p = -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_f = (v(j+1)-v(j))/dx;
          drhodx_f = (rho(j+1)-rho(j))/dx;
          dTdx_f = (T(j+1)-T(j))/dx;
          dlnAdx_f = (log(A(j+1))-log(A(j)))/dx;

          %Continuity Equation
          drhodt_p(j) = -rho(j)*dvdx_f - rho(j)*v(j)*dlnAdx_f - v(j)*drhodx_f;

          %Momentum Equation
          dvdt_p(j) = -v(j)*dvdx_f -(1/gamma)*(dTdx_f + T(j)/rho(j)*drhodx_f);

          %Energy Equation
          dTdt_p(j) = -v(j)*dTdx_f - (gamma-1)*T(j)*(dvdx_f + v(j)*dlnAdx_f);

          %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
        dvdx_b = (v(j)-v(j-1))/dx;
        drhodx_b = (rho(j)-rho(j-1))/dx;
        dTdx_b = (T(j)-T(j-1))/dx;
        dlnAdx_b = (log(A(j))-log(A(j-1)))/dx;
        
        %Continuity Equation
        
        drhodt_c(j) = -rho(j)*dvdx_b - rho(j)*v(j)*dlnAdx_b - v(j)*drhodx_b;
        
        % Momentum Equation
        dvdt_c(j) =-v(j)*dvdx_b -(1/gamma)*(dTdx_b + T(j)/rho(j)*drhodx_b);
        
        %Energy Equation 
        dTdt_c(j) = -v(j)*dTdx_b - (gamma-1)*T(j)*(dvdx_b + v(j)*dlnAdx_b);
    end
      
% Compute the average time derivative
        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 i = 2:n-1
            rho(i) = rho_old(i) + drhodt_avg(i)*dt;
            v(i) = v_old(i) + dvdt_avg(i)*dt;
            T(i) = T_old(i) + dTdt_avg(i)*dt;
        end
        
        
 % Apply the boundary conditions
 
 % Inlet 
            rho(1)=1; 
            v(1)= 2*v(2)-v(3);
            T(1)=1;
  % 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);
             
 %Calculation of other variables of interest
 
   p=rho.*T;                     %Non Dimensional Pressure  
   M=v./sqrt(T);                 %Mach number
   massflow= rho.*v.*A;          %flow rate
  
% Transient Throat Values
   
   rho_throat(k) = rho(n_throat);
   T_throat(k)   = T(n_throat);
   p_throat(k)   = p(n_throat);
   M_throat(k)   = M(n_throat);
   
   %Check For time convergence
   error =max(abs(massflow -massflow_old));
   
 % We want to know how te massflow rate changes with time before the
   % steady solution
   
  TransientMassFlow= figure(1);
  hold on;
   
   if k == 1
       plot(x,massflow ,'--k');
   else if k == 50
           plot(x ,massflow,'r');
       else if k ==100
               plot(x ,massflow,'g');
           else if k == 150
                   plot(x ,massflow,'c');
               else if k == 200
                       plot(x ,massflow,'b');
                   else if k== 700
                          plot(x, massflow,'-.b','linewidth',1.25);
                          plot(x(end),massflow(end),'k*');
                       end
                   end
               end
           end
       end
   end
end

legend('0Deltat','50Deltat','100Deltat','150Deltat','200Deltat','700Deltat');
title('Transient Non-Diemsional Mss Flow Rates(Non-Conservative Form)');
xlabel('non dimensional x/L','Fontsize',11);
ylabel('rhoVA/rho_{0}V_{0}A_{0}');
axis([0 3 0.2 1.8]);
grid minor;
end
                  

Conservation Solver:

function [rho_con,v_con,T_con,p_con,M_con,massflow_con,rho_throat_con,T_throat_con,p_throat_con,M_throat_con,k2]= conservative(n,x,dx,gamma,A,n_throat,C,error2,tol)

% Calculate Initial Profiles
for i=1:length(x)
    if x(i)>=0  && x(i)<=0.5
        rho_con(i)=1;
        T_con(i) =1;
    else if x(i)>=0.5 && x(i)<=1.5
            
        rho_con(i)= 1-(0.366*(x(i)-0.5));
        T_con(i) = 1 - (0.167*(x(i)-0.5));
    
        else if x(i)>=1.5 && x(i)<=30
                
        rho_con(i) = 0.634 - 0.3879*(x(i)-1.5);
         T_con(i) = 0.833 -0.3507*(x(i)-1.5);
            end
        end
    end
end
v_con =0.59./(rho_con.*A);
p_con=rho_con.*T_con;

%Calculate the initial solution vectors
U1 = rho_con.*A;
U2 = rho_con.*v_con.*A;
U3 = rho_con.*A.*((T_con./(gamma-1))+ (gamma/2)*v_con.^2);

%Initialization for the time loop
k2=0;
% Outer time loop with a convergence criteria
while error2>tol
   %convergence variables
   k2=k2+1;
   massflow_old_con=rho_con.*v_con.*A;

% We store the old Solution Vectors
    U1_old = U1;
    U2_old = U2;
    U3_old = U3;
%Defining the Flux vectors:    
      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);

 % CFL Adaptive Time Loop Control
 a_con=sqrt(T_con);
 TimeStep=C*dx./(a_con+v_con);
 dt=min(TimeStep);
 
 %predictor method
  for j = 2:n-1
      % Defining the sourse term     
     J2(j) = (1/gamma).*rho_con(j).*T_con(j).*((A(j+1)-A(j))/dx);
      
      % Componants of U
      
      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;
      
      % Solution Update 
      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
  
  % Predicted values of the primitive variables and flux vectors
  rho_con = U1./A;
  v_con=U2./U1;
  T_con = (gamma-1)*(U3./U1-(gamma/2)*v_con.^2);
  p_con = rho_con.*T_con;
  
  % Computing new set of values using predicted values
  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);
  
  % Corrector Method
  for j = 2:n-1 
      J2(j) = (1/gamma)*rho_con(j).*T_con(j)*((A(j)-A(j-1))/dx);
      
      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
  
   % Average Time Derivatives
   dU1dt_avg = 0.5*(dU1dt_p + dU1dt_c);
   dU2dt_avg = 0.5*(dU2dt_p + dU2dt_c);
   dU3dt_avg = 0.5*(dU3dt_p + dU3dt_c);
  
  
  % Final Solution Update
  for i =2:n-1
      
      U1(i) = U1_old(i) + dU1dt_avg(i)*dt;
      U2(i) = U2_old(i) + dU2dt_avg(i)*dt;
      U3(i) = U3_old(i) + dU3dt_avg(i)*dt;
  end
  
  % Apply the Boundary Conditions
  
  % Inlet 
  rho_con(1)=1;
  T_con(1)=1;
  
   U1(1) = rho_con(1)*A(1);
   U2(1) = 2*U2(2) - U2(3);
   U3(1) = U1(1)*(T_con(1)/(gamma-1)+(gamma/2)*v_con(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);
   
   % Calculation of the new primitive variables
   rho_con = U1./A;
   v_con = U2./U1;
   T_con = (gamma-1)*(U3./U1-(gamma/2)*v_con.^2);
   p_con=rho_con.*T_con;
   
   %Calculation of other variables of interest
   
   M_con=v_con./sqrt(T_con);         % Mach Munber
   massflow_con=rho_con.*v_con.*A;   %Non - dimensional mass Flow rate
   p_con=rho_con.*T_con;             %Non - dimensional pressure 
  
   
   % Capturing flow field values at throat for each time step
   
   rho_throat_con(k2) = rho_con(n_throat);
   T_throat_con(k2)  = T_con(n_throat);
   p_throat_con(k2)  = p_con(n_throat);
   M_throat_con(k2)   = M_con(n_throat);
   
   %Check for time convergence
   error2 =max(abs(massflow_con - massflow_old_con));
   
   %WE wnat to know how the massflow rate changes with time before the
   %steady solution
   
   TransientMassFlow_con=figure(4);
   hold on;
   
   if k2 == 100
       plot(x,massflow_con,'g');
        else if k2 == 200
           plot(x ,massflow_con,'r');
           else if k2 == 700
               plot(x ,massflow_con,'b','linewidth',2);
               end
            end
   end
end
  legend('100Deltat','200Deltat','700Deltat');
  title('Transient Non Dimensional mass Flow Rtae (Non-Conservative Form)');
  xlabel('non dimensional x/L','Fontsize',11);
  ylabel('rhoVA/rho_{0}V_{0}A_{0}');
  axis([0 3 0.5 0.8])
  grid minor;
end
      

Explanation & Results:

A) NON-CONSERVATION  FORM

A.1 - Timewise variation

     A.1.1 - Throat Variable Change

               This section shows the variation of the main properties of the throat during the iterations. The purpose of this section is to see the stability of the scheme and to notice any strangely high or low results in steady values. The non-dimensional values of properties at the throat are well known and that is why this section was chosen to check these calculations.

Here, the stability of the nonconservation form solver has been demonstrated. The slight disturbance during the first few iterations is soon dumped and the solution becomes stable. Although it will be compared lately, here we want to highlight the coherence of the results, as the Mach =1 is achieved at the throat, as it was supposed to be.

A.1.2).  Mass Flow Profile Evolution Through Time

  Another useful parameter to check the solution convergence is the normalized mass flow rate. Under steady conditions, the mass flow rate of the nozzle should be constant at any of the nodes. This is the error parameter that was used to check the convergence in the time loop.

Here, note that the first dashed curve is the variation of the profile at the initial conditions. The strange-looking is simply the product of initial values given to the density, velocity, and the parabolic variation of the nozzle area ratio. Note that after 50 times steps, the mass flow has considerably changed and it keeps changing radially up to 200c steps. Here we can see that the solution is approximating the final solution. Around 700 iterations were enough to achieve the almost constant mass flow rate profile. Also, note that this convergence was to values close to the analytical result of 0.579

A.1.3) - Time and Iteration Needed

A.2 -Steady Solutions

  A.2.1 - Important variables Steady Profile

once the solver was tested to be stable, the final steady profiles were computed and are shown here.

What we can see in these graphs, is that they follow the tendency expected at a nozzle working under these conditions. The density, temperature, and pressure are dropping through the exit of the nozzle, at the expend of an acceleration of the fluid(increase in the Mach Number). Although not compared here, these results are almost identical to the analytical solution of the problem, listed in the Reference Book.

B) CONSERVATION FORM

 B.1 - Timewise Variation

B.1.1 - Throat Variable Changes

As in the nonconservation solver, this section was used to check for convergence in the numerical scheme. Also, note that the change until the steady solution does not go through the same path. That is because of the different numerical schemes(iteration solver).

We can see that the properties get stable real soon in the simulation. The number of iteration is however higher. This happens because the convergence condition was stated in the mass flow rate instead of these parameters. Again the solution is stable and it seems to achieve promising results.

B.1.2) - Mass Flow Profile Evolution Through Time

Once again, the mass flow rate profile provides us an extra insight into the mechanics of the timewise variation o the flow and its approach to the steady.

Here we can appreciate that the solution converges more softly than that of a nonconservation solver. Even from the beginning, the initial condition leads to a more similar profile to the solution.

B.1.3) - Time And Iteration To Converge

 Here, it can be seen that the computational time needed for the conservation solver is a bit higher than for the non-conservation scheme. However, it does not seem like an important difference that would make one decision for one of them, things like accuracy and stability should be checked first.

Also note that the number of iterations is based on the mass flow rate criteria established, and not the primitive variables. Probably, changing to these last ones would make the conversation scheme need fewer iterations(for what can be seen in the previous section graph)

B.2 - STEADY SOLUTIONS

B.2.1 - Important Variable Steady profile

Once again, the solution achieved is not only stable but also pretty accurate. This profile is almost identical to the one achieved in the nonconservation analysis.

B.2.2 - Srady Mass Flow Profile Comparison

The stability of both schemes is undeniable. Although both seem to be reliable, it is now time to analyze accuracy and time, to decide if there is a more fitting scheme for this application.

Here, the steady mass flow profile is compared. This figure illustrates a general advantage of the conservation form of equations. This form always does a better job of preserving math throughout the flow field. This is mainly because the mass flow itself is one of the variables in the equations. This does not mean that the conservation form is a better fit for the numerical scheme, as it will be seen in the following section.

C. GRID DEPENDENCE TEST

The previous discussion does not necessarily mean the superiority of the conservation form, and further analysis must be done. In this path, the grid dependence test is usually critical to decide whether a problem can be solved faster with one numerical scheme or another or if just can be solved with them.

 This test consists of changing the number of the grid on the mesh to check if the results are dependent on the number of points. Any numerical scheme should have as an objective manage to find grid independence. In this case, we change the number of grids from 31 to 61 and perform the same calculation. The result is shown in the table below, stating that the problem is already independent of the number of computation points taken in the control volume, at least to the precision needed here. Also, this result showed that the non-conservation solver is more accurate at computing the primitive variables.