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

Objective:

The main object of this project is to write separate functions for conservative and non-conservative forms

• To figure out the minimum number of cycles for which the simulation should be run in order for convergence
• Comparsion the normalized mass flow rate between the conservative forms and non-conservative forms
• To find which method is faster.

Assumptions:

• It is assumed that the flow inside the nozzle is isentropic

• The stream flow is considered Quasi-1D as properties vary along the x-axis and not the y-axis

• Compressible flow

• No consideration of friction and heat transfer

• Considerable pressure ratio is maintained between chocked pressure and back ressure to allow expansion of the flow and avoid waves.

Problem statement

A nozzle is a device that is used to give the direction to the gases coming out of the combustion chamber. The nozzle is a tube which has a capacity to convert the thermo-chemical energy generated in the combustion chamber into kinetic energy.

The nozzle converts the low velocity, high pressure, high-temperature gas in the combustion chamber into high-velocity gas of lower pressure and low temperature. A convergent-divergent nozzle is used if the nozzle pressure ratio is high. High-performance engines in supersonic aircraft generally incorporate some form of a convergent-divergent nozzle.

In the present day, there is a huge development in Aerospace Engineering for various prospects. Extensive research is being carried out in fields like civil and defense prospects.

The virtualization is one of the major developments in the field of research, which has revolutionized Aerospace engineering, along with all other branches. The computational techniques are being used widely for getting better results, which are close to experimental techniques.The flow through a convergent-divergent nozzle is one of the benchmark problems used for modeling the compressible flow through computational fluid dynamics.

Governing Equations

The governing equations of a mathematical model describe how the values of the unknown variables (i.e. the dependent variables) change when one or more of the known (i.e. independent) variables change.

Non Conservation form:

In this Non-conservation flow, a small control volume is taken inside the pipe flow but with a small change, the control volume here is moving with the flow instead of being fixed as we saw in conservative flow.

Conservation form:

In this conservation flow, when a fluid is flowing through a pipe a small control volume is taken in some area inside that pipe and this control volume is fixed So that the equation can be solved to find the properties of the flow-through that points this is known as Convesative flow.

The equations were solved using Macormack method.

Matlab code:

clear all
close all
clc


%Inputs 
% Number of nodes 
n=31;

%Mesh
x=linspace(0,3,n);

%Step size
dx=x(2)-x(1);

%gamma value
gamma = 1.4;

%Time step 
nt = 5000;
c1 = 0.5;
c2 = 0.5;
tol = 1e-6;
mass_tol1 = 1e-3;
mass_tol2 = 1e-2;

%nonconservative

[rho1,v1,T1,total_time1,netmf1,err_v1,err_rho1,err_T1] = nonconserv(n,x,dx,gamma,nt,c1,tol,mass_tol1);

%conservative

[rho2,v2,T2,total_time2,netmf2,err_v2,err_rho2,err_T2] = conserv(n,x,dx,gamma,nt,c2,tol,mass_tol2);

%Velocity plots 
figure(1)
plot(x,v1,\'b\',x,v2,\'r\')
xlabel(\'Non-dimensional Distance\')
ylabel(\'Non-dimensional Velocity\')
legend(\'Non-Conservative\',\'Conservative\')
grid on


%Density plots 
figure(2)
plot(x,rho1,\'b\',x,rho2,\'r\')
xlabel(\'Non-dimensional Distance\')
ylabel(\'Non-dimensional Density\')
legend(\'Non-Conservative\',\'Conservative\')
grid on


%Temperature plots 
figure(3)
plot(x,T1,\'b\',x,T2,\'r\')
xlabel(\'Non-dimensional Distance\')
ylabel(\'Non-dimensional Temperature\')
legend(\'Non-Conservative\',\'Conservative\')
grid on

%Net Mass Flow Rate Plot
figure(4)
subplot(2,1,1)
plot(netmf1)
xlabel(\'Time Step\')
ylabel(\'Net Mass Flow Rate\')
title(\'Non-Conservative Form\')
grid on
subplot(2,1,2)
plot(netmf2)
xlabel(\'Time Step\')
ylabel(\'Net Mass Flow Rate\')
title(\'Conservative Form\')
grid on


%Error in Velocity with time
figure(5)
subplot(2,1,1)
plot(err_v1)
xlabel(\'Time Step\')
ylabel(\'Error in Velocity\')
title(\'Non-Conservative Form\')
grid on
subplot(2,1,2)
plot(err_v2)
xlabel(\'Time Step\')
ylabel(\'Error in Velocity\')
title(\'Conservative Form\')
grid on


%Error in Density with time
figure(6)
subplot(2,1,1)
plot(err_rho1)
axis([0 inf -0.01 0.02])
xlabel(\'Time Step\')
ylabel(\'Error in Density\')
title(\'Non-Conservative Form\')
grid on
subplot(2,1,2)
plot(err_rho2)
axis([0 inf -0.005 0.015])
xlabel(\'Time Step\')
ylabel(\'Error in Density\')
title(\'Conservative Form\')
grid on


%Error in Temperature with time
figure(7)
subplot(2,1,1)
plot(err_T1)
axis([0 inf -0.005 0.01])
xlabel(\'Time Step\')
ylabel(\'Error in Temperature\')
title(\'Non-Conservative Form\')
grid on
subplot(2,1,2)
plot(err_T2)
axis([0 inf -0.01 0.02])
xlabel(\'Time Step\')
ylabel(\'Error in Temperature\')
title(\'Conservative Form\')
grid on

Function for Non-Conservation form:

function [rho,v,T,k,netmf,err_v,err_rho,err_T] = nonconserv(n,x,dx,gamma,nt,c,tol,mass_tol)

% Intial Profile
% Density
rho = (1 - (0.3146*x));  
% Tempaerature
T = (1 - 0.2314*x);
% Velocity
v = (0.1 + (1.09*x)).*(T.^0.5);
% Area
a = (1 + 2.2*(x-1.5).^2)


%Outer time Loop 
for k = 1: nt
    
    dt = min(abs(c*dx./(realsqrt(T)+v)));
    
    %copying old variables
    rho_old = rho;
    v_old = v;
    T_old = T;
    



     %Predictor method
     for j = 2:n-1
         
         dvdx = (v(j+1) - v(j))/dx;
         
         drhodx =(rho(j+1) - rho(j))/dx;
         
         dlogadx = (log(a(j+1)) - log(a(j)))/dx;
         
         dTdx = (T(j+1) - T(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
         v(j) = v(j) + dvdt_p(j)*dt;
         rho(j) = rho(j) + drhodt_p(j)*dt;
         T(j) = T(j) + dTdt_p(j)*dt;
         
     end
     
     %corrector method
     for j = 2:n-1
         
         dvdx = (v(j) - v(j-1))/dx;
         
         drhodx = (rho(j) - rho(j-1))/dx;
         
         dlogadx = (log(a(j)) - log(a(j-1)))/dx;
         
         dTdx = (T(j) - T(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
     
     %computing average time derivatives
     
     drhodt = 0.5*(drhodt_p + drhodt_c)
     
     dvdt = 0.5*(dvdt_p + dvdt_c)
     
     dTdt = 0.5*(dTdt_p + dTdt_c)
     
     %Final Update
          
       for j = 2:n-1
           
           v(j) = v_old(j) + dvdt(j)*dt;
           
           rho(j) = rho_old(j) + drhodt(j)*dt;
           
           T(j) = T_old(j) + dTdt(j)*dt
           
           
       end
       
       
       
       %Applying 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);
       
       %Mass flow rate
       
       inflow = rho(1)*v(1)*a(1);
       
       outflow = rho(n)*v(n)*a(n);
       
       netmf(k) = (outflow-inflow)
       
       %Error tolerance
     err_v(k) = max(v - v_old);
     err_rho(k) = max(rho - rho_old);
     err_T(k) = max(T-T_old);
     if(err_v(k) < tol && err_rho(k) < tol && err_T(k) < tol && netmf(k))
       
     end
end

Function for Conservation form:

function [rho,v,T,k,netmf,err_v,err_rho,err_T] = conserv(n,x,dx,gamma,nt,c,tol,mass_tol)

%Intial Profile
%Density
rho = (1 - (0.3146*x));  
%Tempaerature
T = (1 - 0.2312*x);
%Area
a = 1+2.2*(x-1.5).^2;
%Velocity
v = 0.59./(rho.*a);

% Derivatives
U1 = rho.*a;

U2 = rho.*a.*v;

U3 = rho.*(T/(gamma-1)+(gamma/2)*v.^2).*a;

% Outer Time loop
for k = 1:nt
    
     % storing old values
     U1_old = U1;
     U2_old = U2;
     U3_old = U3;
     
     rho_old = U1_old./a
     v_old =  U2_old./(rho.*a);
     T_old =  (gamma-1)*((U3_old./U1_old) - (gamma/2)*v.^2);
     
     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);

     
     
     for j = 2:n-1
         J2(j) = (1/gamma)*rho(j)*T(j)*(a(j+1) - a(j))/dx;
     end
      
     for i = 1:n
         deltat(i)=(c*(dx/((T(i)^0.5+v(i)))));
     end
     
     dt = min(deltat)
     
     %Predictor method
     for j = 2:n-1
         %Continuity Equation
         dU1dt_p(j) = -(F1(j+1)-F1(j))/dx;
         
         %momentum Equation
         dU2dt_p(j) = -(F2(j+1)-F2(j))/dx + J2(j);
         
         %Energy Equation
          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
     
         %Calculating primitive variables 
         rho = U1./a;
         v = U2./(rho.*a);
         T = (gamma-1)*((U3./U1) - (gamma/2)*v.^2);
     
         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);
         
         for j = 2:n-1
             
            J2(j) = (1/gamma)*rho(j)*T(j)*(a(j) - a(j-1))/dx;
             
         end
         
         %Corrector Method
         
     for j = 2:n-1
         %Continuity Equation
         dU1dt_c(j) = -(F1(j)-F1(j-1))/dx;
         
         %momentum Equation
         dU2dt_c(j) = -(F2(j)-F2(j-1))/dx + J2(j);
         
         %Energy Equation
         dU3dt_c(j) = -(F3(j)-F3(j-1))/dx;
         
     end
        
     %computing average time Derivative
     
     dU1dt = 0.5*(dU1dt_p + dU1dt_c);
     
     dU2dt = 0.5*(dU2dt_p + dU2dt_c);
     
     dU3dt = 0.5*(dU3dt_p + dU3dt_c);
     
     
     %Final update
     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
     
     %Applying Boundary Conditions
     %inlet 
     U2(1) = 2*U2(2) - U2(3)
     
     %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);
     
     
     %Recalculating Primitive variables
     
     rho = U1./a
     v =  U2./(rho.*a);
     T =  (gamma-1)*((U3./U1) - (gamma/2)*v.^2);

     
     %mass flow rate
     
        inflow  = rho(1)*v(1)*a(1);
     
        outflow = rho(n)*v(n)*a(n);
     
         netmf(k) = outflow-inflow;
     
     %Error tolerance
     err_v(k) = max(v - v_old);
     err_rho(k) = max(rho - rho_old);
     err_T(k) = max(T-T_old);
     if(err_v(k) < tol && err_rho(k) < tol && err_T(k) < tol && netmf(k))
     
     end
end

     
     

Results:

The simulation was run for a maximum of 5000-time steps, with CFL = 0.5 and n=31 the tolerance was kept to 1e-6, where the change in rho, v, and T had to fall under the tolerance.Net Mass Flow Rate Tolerance has also been set at 1e-3 for nonconservation and 1e-2 for conservation form.

Non-Conservation Form converged in 676-time steps, by converged it has reached a stable solution.

Conservation Form converged in 4515-time steps, By converged it has reached a stable solution.

Veolcity Plot:

Density Plot:

Temperture Plot:

It has been Observed that conservation form solves in less number of iterations than non- conservative form. However, the net mass flow rate hasn\'t stabilized in the current conditions.

Net Mass Flow Rate:

Net Change in velocity:

Net change in Density:

Net change in Temperature:

Grid dependence:

Even if the number of grid points is doubled and the solution varies only marginally then our original calculation with the initial no of grid points is said to be Grid - independent.

Based on the estimated order of accuracy and the allowable tolerance range this element is purely based on the circumstance.

In our case we assumed n=31, So let check whether it is grid-independent or not initially, to start the discussion it would be better if we plot the graphs with respect to the modified values of n i.e n=61.

when modifying the values of \'n\' the number of time steps should also be modified so that the elapsed non-dimensional time is same for both cases by doing so, the grid independence test can be carried out by comparing the curves obtained for both cases at a particular non-dimensional time.

n=61

nt = 6000

Net Mass Flow Rate:

Net Change in velocity:

Net change in Density:

Net change in Temperature:

Which method is faster

It has been Observed that conservation form solves in less number of iterations than the non- conservative form.