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

AIM: To simulate the 1-D super sonic nozzle flow using Macormack Method

A comparison is made between conservative and non conservative form of governing equation to check which one is faster and how the parameters changes along the length of the nozzle in both the forms, a separate program is written to check at which time step the method converges.

Theory:

We consider the steady state, isentropic flow through convergent-divergent nozzle as shown below.A de Laval nozzle ( convergent-divergent nozzle ) is a tube that is pinched in the middle, making a carefully balanced, asymmetric hourglass shape. It is used to accelerate a hot, preassurised gas passing through it to a higher supersonic speed in the axial (thrust) direction, by converting the heat energy of the flow into kinetic energy. Because of this, the nozzle is widely used in some types of steam turbines and rocket engine nozzles. It also sees use in supersonic jet engine. The flow at inlet to the nozzle comes from a reservoir. The cross-sectional area of the reservoir is large and hence the velocity is very small. Thus pressure and temperature at reservoir are stagnation values or total pressure and total temperatures, respectively. The flow expands isentropically to super-sonic speeds at nozzle exit. The flow is locally subsonic in the convergent section of the nozzle and sonic at the throat and supersonic at the divergent section. The sonic flow is M=1 , where M is mach no. means that local velocity at this location is equal to local speed of sound.

The area of nozzle changes as a function of distance along the nozzle (x), we make the assumption that the flow properties vary only with x, to assume uniform flow properties accross any given cross section. Such flow is defined as quasi-one-Dimensional.

Non- conservative form of the governing equations:

Conservative form of the governiing equations:

note :- The above equations are to be solved using the Maccormacks Techinique. Explanations of the above equations ais mentioned in J.D. Andersons text book CFD The Basics with Applications.

MacCormack method:

The MacCormack method is a variation of the two-step Lax–Wendroff scheme but is much simpler in application. To illustrate the algorithm, consider the following first order hyperbolic equation

The application of MacCormack method to the above equation proceeds in two steps; a predictor step which is followed by a corrector step.

Expected Results :-

1. Steady-state distribution of primitive variables inside the nozzle
2. Time-wise variation of the primitive variables
3. Variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process
4. Comparison of Normalized mass flow rate distributions of both form.
5. Perform a grid independence test on the solution

:-

clear all;
close all;
clc;



% Inputs required
L = 3;
n = 31 ;
x = linspace(0,3,n);
nt =1400;
dx = L/(n-1);
 a = 1+2.2*(x-1.5).^2;  % area profile - Independent of time ;
 c = 0.5 ; % courant number 
 gamma = 1.4 ; % specific heat ratio for air with the assumption that it is ideal
 throat = (n-1)/2;

%  Numerical solution for non conservative form of goevrning equations 
tic 
[rho,v,t,p,rho_th,t_th,p_th,m,M,M_th] = non_conservative(x,dx,n,nt,gamma,a,throat,c);
time_elapsed = toc;
fprintf('Time taken for computation  by Non-Conservative form: %0.4f',time_elapsed);

% Density , pressure, temperature & Mach number for Non-Conservative form :

figure(2);

subplot(4,1,1);
plot(rho_th, 'b');
title ('vatiation of parameters wrt time steps (Non conservative form)');
ylabel('Density');
xlabel('Time steps');
axis([0 1400 0 1.3]);


subplot(4,1,2);
plot(p_th,'g');
xlabel('Time steps');
ylabel('Preassure');
axis([0 1400 0.4 1.1]);


subplot(4,1,3);
plot(t_th,'r');
xlabel('Time steps');
ylabel('Temprature');
axis([0 1400 0.6 1]);


subplot(4,1,4);
plot(M_th,'m');
xlabel('Time steps');
ylabel('Mach Number');
axis([0 1400 0.6 1.4]);




%Spacial quantification of P,Rho,T,M  in Non-Conservative Form:
 figure(3);
subplot(4,1,1);
plot(x,M,'r');
title('Variation of parameters wrt length of the nozzle  (Non-Conservative form)');
ylabel('Mach Number');
xlabel('Length of the nozzle');
axis([0 3 0 4]);


subplot(4,1,2);
plot(x,p,'b');
ylabel('Pressure');
xlabel('Length of the nozzle');
axis([0 3 0 1.01]);


subplot(4,1,3);
plot(x,rho,'m');
ylabel('Density');
xlabel('Length of the nozzle');
axis([0 3 0 1.01]);


subplot(4,1,4);
plot(x,t,'g');
ylabel('Temprature');
axis([0 3 0 1.01]);
xlabel('Length of the nozzle');

% % Numerical solution for conservative form of governing equations 
tic 
[rho_c,v_c,t_c,p_c,rho_th_c,t_th_c,p_th_c,m_c,M_c,M_th_c] =conservative(x,dx,n,nt,gamma,a,throat,c)
time_elapsed_conservative = toc;
fprintf('Time taken for computation  by Conservative form: %0.4f',time_elapsed_conservative);

%density , pressure, temperature & Mach number for conservative form :

figure(5)
subplot(4,1,1);
plot(rho_th_c,'g');
title('Vatiation of parameters wrt time steps ( conservative form)');
ylabel('Density');
axis([0 1400 0 1.3]);
xlabel('Time steps');



subplot(4,1,2);
plot(p_th_c,'b');
ylabel('Pressure');
axis([0 1400 0.4 1.1]);
xlabel('Time steps');



subplot(4,1,3);
plot(t_th_c,'m');
ylabel('Temprature');
axis([0 1400 0.6 1]);
xlabel('Time steps');



subplot(4,1,4);
plot(M_th_c,'r');
ylabel('Mach Number');
axis([0 1400 0.6 1.4]);
xlabel('Time steps');




%Spacial quantification of P,Rho,T,M  in Conservative Form:
 figure(6);


subplot(4,1,1);
plot(x,M_c,'r');
title('Vatiation of parameters wrt Length of the nozzle( conservative form)');
ylabel('Mach Number');
xlabel('Length of the nozzle');
axis([0 3 0 4])


subplot(4,1,2);
plot(x,p_c,'b');
xlabel('Length of the nozzle');
ylabel('Pressure');
axis([0 3 0 1.01])


subplot(4,1,3);
plot(x,rho_c,'m');
ylabel('Density');
axis([0 3 0 1.01])
xlabel('Length of the nozzle');


subplot(4,1,4);
plot(x,t_c,'g');
ylabel('Time');
axis([0 3 0 1.01])
xlabel('Length of the nozzle');


% Mass flow rate comparison for both conservative and non consevative forms
 figure(7);
hold on;
plot(x, m, 'b');
plot(x, m_c, 'r');
line([0 3], [0.579 0.579]);%Analytical solution 

xlabel('Length of the nozzle');
ylabel('Normalized mass flow rate');
title('Dimensionless mass flow comparison ');
legend('non-conservative','conservative');

Conservative form :-

function [rho_c,v_c,t_c,p_c,rho_th_c,t_th_c,p_th_c,m_c,M_c,M_th_c] = conservative(x,dx,n,nt,gamma,a,throat,c)

% Initial conditions 
for i = 1:length(x)
    if x(i) >= 0 && x(i) < 0.5
        rho_c(i) = 1;
        t_c(i) = 1;
    else if x(i) >= 0.5 && x(i) < 1.5
            rho_c(i) = 1.0 -0.366*(x(i) - 0.5);
            t_c(i) = 1.0 - 0.167*(x(i) - 0.5);
    else if x(i) >= 1.5 && x(i) <= 3.0
                rho_c(i) = 0.634 - 0.3879*(x(i) - 1.5);
                t_c(i) = 0.833 - 0.3507*(x(i) - 1.5);
        end
        end
   
  
  end
end
% Initial velocity and pressure:
v_c = 0.59./(rho_c.*a);
p_c = rho_c.*t_c;

% Solution vectors:
U1 = rho_c.*a;
U2 = rho_c.*a.*v_c;
U3 = rho_c.*a.*((t_c./(gamma - 1)) + (gamma/2)*v_c.^2);

% Time Iteration Loop
for k = 1:nt
    
    % Preserving the old solution vectors:
    U1_old = U1;
    U2_old = U2;
    U3_old = U3;
    
    % Adaptive Time-step control:
    dt = min(c.*dx./(sqrt(t_c) + v_c));
    
    % 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);
    
    %% Solution by Mc-Cormack Method:
    for j = 2:n-1
        
        % Defining the source term:
        J2(j) = (1/gamma).*rho_c(j).*t_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;
    t_c = (gamma -1)*(U3./U1 - (gamma/2)*v_c.^2);
    p_c = rho_c.*t_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).*t_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 corrector step:
    % Averaging the predictor and corrector steps:
    dU1dt = 0.5*(dU1dt_p + dU1dt_c);
    dU2dt = 0.5*(dU2dt_p + dU2dt_c);
    dU3dt = 0.5*(dU3dt_p + dU3dt_c);
    
    % Calculating the 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
    % Applying BC:
    %  Inlet conditions
    
    U2(1) = 2*U2(2) - U2(3);
    
    % 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 updation of the primitive variables:
    rho_c = U1./a;
    v_c = U2./U1;
    t_c = (gamma -1)*(U3./U1 - (gamma/2)*v_c.^2);
    p_c = rho_c.*t_c;
    %%
    % Defining the Values of mach Number and mass flow rate:
    M_c = v_c./sqrt(t_c);
    m_c = rho_c.*a.*v_c;
    
    % Calculating values of specific quantites at the throat:
    M_th_c(k) = M_c(throat);  % Mach Number
    p_th_c(k) = p_c(throat);  % Pressure Ratio
    t_th_c(k) = t_c(throat);  % Non-dimensional Temperature
    rho_th_c(k) = rho_c(throat);  % Density Ratio
    
    % Comparison of non-dimensional mass-flow rates at different times:
     figure(4);
    hold on;
    if k == 100
       plot(x,m_c,'b','linewidth',1.75);
     else if k == 200
             plot(x,m_c,'g','linewidth',1.75);
     else if k == 700
             plot(x,m_c,'r','linewidth',1.75);
        end
        end
        end
    end
xlabel('Length of the nozzle');
ylabel('Conservative mass flow rate');
legend('100 Dt','200 Dt','700 Dt');
axis([0 3 0.4 1.6]);


end

Non-conservative form :-

function [rho,v,t,p,rho_th,t_th,p_th,m,M,M_th] = non_conservative(x,dx,n,nt,gamma,a,throat,c)

% Using MacCormack's method to solve non- conservative form of 1-D supersonic nozzle equation. 
rho = 1- 0.3146*x;
t = 1 - 0.2314*x;
v = (0.1+1.09*x).*t.^0.5;


% Outer time loop 

for k = 1:nt 
    
rho_old =rho;
t_old =t;
v_old =v;
 %Adaptive time-step control
    dt = min(c.*dx./(sqrt(t)+v));

    
    %     predictor method 

for j=2:n-1
    dvdx =( v(j+1)-v(j))/dx;
    drhodx =( rho(j+1)-rho(j))/dx;
    dtdx =( t(j+1)-t(j))/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
    dvdx =( v(j)-v(j-1))/dx;
    drhodx =( rho(j)-rho(j-1))/dx;
    dtdx =( t(j)-t(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 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
% BC
% Inlet Boudary conditions 
   v(1) = 2*v(2) -v(3);
 rho(1) = 1;
   t(1) = 1;
% Outlet Boudary conditions 
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 pressure, mass flow rate and mach number:
    p = rho.*t; %Non-Dimensional Pressure
    m = rho.*a.*v; %Non-Dimensional Mass flow rate
    M = v./sqrt(t); %Mach Number
    
    %Calculation of Quantities at the Throat point of the nozzle:
    t_th(k) = t(throat);
    p_th(k) = p(throat);
    rho_th(k) = rho(throat);
    M_th(k) = M(throat);
    
    % Comparison of non-dimensional mass-flow rates at different times:
    figure(1);
    hold on;
    if k == 50
       plot(x,m,'b','linewidth',1.75);
    else if k == 100
            plot(x,m,'g','linewidth',1.75);
    else if k == 150
            plot(x,m,'r','linewidth',1.75);
    else if k == 200
            plot(x,m,'m','linewidth',1.75);
    else if k == 700
            plot(x,m,'-','linewidth',1.75);
            plot(x(end),m(end),'b*');
        end
        end
        end
        end
        end
    end
xlabel('Length of the nozzle');
ylabel('Non- Conservative mass flow rate');
legend('50 Dt','100 Dt','150 Dt','200 Dt','700 Dt');
axis([0 3 0 2.5]);

end     

Resulting graphs:-

When the grid value is 31.

1. Steady-state distribution of primitive variables inside the nozzle

a)Primitive variables variying inside the nozzle for non-conservative form.-

b)Primitive variables variying inside the nozzle for conservative form.-

2. Time-wise variation of the primitive variables

a) Primitive variables variying inside the nozzle for non-conservative form.-

b) Primitive variables variying inside the nozzle for conservative form.-

3. Variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process

a) Non- conservative normalized mass flow rate at specific time steps -

b) Conservative normalized mass flow rate at specific time steps-

4. Comparison of Normalized mass flow rate distributions of both forms

Variation of normalized mass flow rates of both conservative and non-conservative form along the length of the nozzle is given below.

Performing a grid independence test on the solution taking n = 61.

When the grid value is 61.

1. Steady-state distribution of primitive variables inside the nozzle

a)Primitive variables variying inside the nozzle for non-conservative form.-

b)Primitive variables variying inside the nozzle for conservative form.-

2. Time-wise variation of the primitive variables

a) Primitive variables variying inside the nozzle for non-conservative form.-

b) Primitive variables variying inside the nozzle for conservative form.-

3. Variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process

a) Non- conservative normalized mass flow rate at specific time steps -

b) Conservative normalized mass flow rate at specific time steps-

4. Comparison of Normalized mass flow rate distributions of both forms

Variation of normalized mass flow rates of both conservative and non-conservative form along the length of the nozzle is given below.

Results and observation for grid dependence test-

for n = 31		Pressure	Density	Temperature	Mach Number	Time Elapsed
Conservative		0.5487	0.6524	0.841	0.978	11.714 sec
Non-conservative		0.5342	0.6387	0.8365	0.9994	15.518 sec
Exact		0.528	0.634	0.833	1	-
for n= 61
Conservative		0.534	0.6391	0.8355	0.9979	22.712 sec
Non-conservative		0.5327	0.6377	0.8354	0.9997	34.677 sec
Exact		0.528	0.634	0.8331	1	-

• It is observed that conservative form is faster than non-conservative form.
• There is not much diffrence when it comes to increase in the no of mesh points.
• n = 61 takes significantly larger amount of time almost double when compared to n=31.
• Mass flow rate in non-conservative forms show more number of fluctuations when compared to conservative form.
• Accuracy of the solutions remains almost same in both forms.

Conclusion:

Mac-cormacks method as described above, because a two step predictor-corrector sequence is used with forward diffrencing on predictor and backward diffrencing on corrector it is a second order accurate technique. It is much easier to apply.

Conservative form of the equations produces comparitively stable equations and takes less time to compute, but both forms produce equally accurate solutions.