Solving Quasi 1D Supersonic Flow through a Convergent Divergent Nozzle using MATLAB

Abstract : Set up of a time marching, finite difference solution for the Quasi 1D subsonic supersonic nozzle flow. Non-dimensionalizing the Governing flow equations and setting them up with relevant initial and boundary conditions and solving using the Maccormacks technique. Simulating the solution for convergence and comparing the normalized mass flow rate between Non-conservative and Conservative forms and making conclusions. Performing grid dependence test. 

 

I. Introduction : 

We consider the steady state, isentropic flow through convergent-divergent nozzle as shown below. 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. 

 

II. Project Objectives : 

1. Write matlab code to Simulate and solve 1D subsonic supersonic nozzle flow equations using MacCormacks method.
2. Implementing both Conservative and Non conservative forms of governing equations and making seperate functions both forms.
3. To find out minimum number of cycles for which the simulation should be run to achieve convergence.
4. Compare the normalized mass flow rate between Conservative and Non-conservative forms
5. Timewise variations of density, temperature, velocity and mach no at nozzle throat.
6. Distribution of Mass flow through nozzle at different times during time-marching approach to steady state. 
7. Perform grid dependence test.

 

III. Governing Flow equations (Non-Dimensional variables) : 

Defining the Non-Dimensional variables as, 

Temperature : `T' = T/T_0`   

Density : `rho' = rho/rho_0`

Dimensionless length : `x' = x/L`

Speed of sound : `a_o = sqrt (lambdaRT_0)`

Velocity : `v' = v/a_0`

Time : `t' = t/(L/a_0)`

Area : `A' = A/A^(throat area)`

 

Introducing these Dimensionless variables in the governing flow equations, we get : 

III. A) For Non-conservative form :

 Continuity equation :  `(delrho')/(delt') = -rho'((delV')/(delx')) - rho'V'((del(lnA'))/(delx')) - V'((delrho')/(delx'))`

 

Momentum equation : `((delV')/(delt')) = -V'((delV')/(delx'))-1/gamma ((delT')/(delx')+(T')/(rho') (delrho')/(delx') )`

 

Energy equation :  `((delT')/(delt')) = -V'((delT')/(delx')) - (gamma-1)T'[(delV')/(delx')+V'((del(lnA'))/(delx'))]`

 

III. B) For Conservative form : 

Continuity equation : `(del(rho'A'))/(delt') + (del(rho'A'V'))/(delx')`

 

Momentum equation : `(del(rho'A'V'))/(delt') + (del[rho'A'V'^2 + (1/gamma)p'A'])/(delx') = 1/gamma p'(delA')/(delx')`

 

Energy equation : `(del[rho'((e')/(gamma-1) + lambda/2 V'^2)A'])/(delt') + (del[rho'((e')/(gamma-1) + lambda/2 V'^2)V'A' + p'A'V'])/(delx') = 0`

 Now, the generic form of equation can be expressed as : 

Solution vector : U

Flux term : F

Source term : J 

 

Defining the elements of Solution vector U

`U1 = rho'A'`

`U2 = rho'A'V'`

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

 

Pure form of Flux terms :

F1 = U1 

`F_2 = U_2^2/(U_1) + (gamma-1)/gamma (U_3 - gamma/2 (U_2^2/(U_1)))`

`F_3 = gamma (U_2U_3)/U_1 - (gamma(gamma-1))/2 U_2^3/U_1^2`

Source term : 

`J_2 = (gamma-1)/gamma (U_3 - gamma/2 U_2^2/U_1) (del(lnA'))/(delx')`

 

Primitive Variables : 

We need primitive variables to decode the values after calculating the Solutions vector (U1,U2,U3) which is obtained at each time step

`rho' = U_1/A'`

`V' = U_2/U_1`

`T' = e' = gamma-1 (U_3/U_1 - gamma/2 V'^2)`

`p' = rho'T'`

 

IV. Now the Continuity, Momentum and Energy equations for Conservative form :       (respectively) : 

`(delU_1)/(delt') = - (delF_1)/(delx')`

 

`(delU_2)/(delt') = - (delF_2)/(delx') + J2`

 

`(delU_3)/(delt') = - (delF_3)/(delx')`

 

V. Time step :  

A stability constraint exists in the system. The present governing equations are nonlinear PDE's . And an exact stability criterion an only be viewed as general guidance for present nonlinear problem. However it is a Good guidance. 

The courant friedrichs lowry (CFL) condition for time step of 1 D flow , in Non dimensional quantities is given as: 

`Deltat = min(C (Deltax)/(a + V))`

where, C = 0.5 (assumed in the problem)

a is local speed of sound 

V is local flow velocity at a point in the flow.

`Deltat`is the increment in non dimensional time and `Deltax` is increment in nondimensional space. 

 

• For full derivation refer : Computational fluid dynamics - Basics with applications by John D anderson

The problem Specification's has been taken from Chapter 7 of this book, and the code is ran for 1400 time steps as stated in there. The results obtained are directly compared and analyzed from the program created. 

 

VI. Matlab Programs  : 

Main program : 

This program calls the non conservaitive and conservative forms and plots all the results and graphs. 

%% Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
% First function called is of Non-conservative form and second one is of
% Conservative form.
% Ref - (CFD - basics with applications )John D anderson chap 7.

% code by  AADIL SHAIKH
% date :   17/9/2019

clear all
close all
clc

%% inputs
n = 31;
x = linspace(0,3,n);
dx = x(2) - x(1);
throat = (n-1)/2 + 1;
gamma = 1.4;
c = 0.5;

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

%time steps
nt = 1400;

%%
tic
[mass_flow_nc,rho_nc,v_nc,t_nc,M_nc,P_nc,rho_t_nc,t_t_nc,P_t_nc,M_t_nc]  = non_conservative_form(n,x,dx,throat,gamma,c,a,nt);
Nc_time = toc;

tic
[u1,u2,u3,mass_flow_c,rho,v,t,M,P,rho_t,t_t,P_t,M_t] = Conservative_form(n,x,dx,throat,gamma,c,a,nt);
C_time = toc;

%%
for y = 1:n
    analytical_massflow(y) = 0.579;
end
figure(7)
plot(x,mass_flow_nc,'linewidth',1.5)
hold on
plot(x,mass_flow_c,'linewidth',1.5)
plot (x, analytical_massflow,'linewidth',1.5)
xlabel('Non dimensional distance across nozzle')
ylabel('Mass flow variations')
legend('Non conservation form','Conservation form','Exact Solution','location','best')
saveas(gcf,'mass flow compare.png')

 

 VI. A) Non-Conservative form : 

function [mass_flow_nc,rho_nc,v_nc,t_nc,M_nc,P_nc,rho_t_nc,t_t_nc,P_t_nc,M_t_nc] = non_conservative_form(n,x,dx,throat,gamma,c,a,nt)

% calculate initial profiles
rho_nc = 1-0.3146*x;
t_nc = 1-0.2314*x;                  % t = temp
v_nc = (0.1+1.09*x).*t_nc.^0.5;

%% outer time loop
for k = 1:nt
    
    %updating old values
    rho_old = rho_nc;
    v_old = v_nc;
    t_old = t_nc;
    
    % time step
    A = (t_nc.^0.5);
    dt = c.*(dx ./((A)+v_nc));
    dt = min(dt);
    
    % predictor method
    for j = 2:n-1
        dvdx = (v_nc(j+1)-v_nc(j))/dx;
        drhodx = (rho_nc(j+1)-rho_nc(j))/dx;
        dtdx = (t_nc(j+1)-t_nc(j))/dx;
        dlogadx = (log(a(j+1))-log(a(j)))/dx;
        
        % Continuity equation
        drhodt_p(j) = -rho_nc(j)*dvdx - rho_nc(j)*v_nc(j)*dlogadx - v_nc(j)*drhodx;
        
        % Momentum eqn
        dvdt_p(j) = -v_nc(j)*dvdx - (1/gamma)*(dtdx + (t_nc(j)/rho_nc(j))*drhodx);
        
        % Energy eqn
        dtdt_p(j) = -v_nc(j)*dtdx - (gamma - 1)*t_nc(j)*(dvdx + v_nc(j)*dlogadx);
        
        % Soln update
        rho_nc(j) = rho_nc(j) + drhodt_p(j)*dt;
        v_nc(j) = v_nc(j) + dvdt_p(j)*dt;
        t_nc(j) = t_nc(j) + dtdt_p(j)*dt;
    end
    
    % Corrector method
    for j = 2:n-1
        dvdx = (v_nc(j)-v_nc(j-1))/dx;
        drhodx = (rho_nc(j)-rho_nc(j-1))/dx;
        dtdx = (t_nc(j)-t_nc(j-1))/dx;
        dlogadx = (log(a(j))-log(a(j-1)))/dx;
        
        % Continuity equation
        drhodt_c(j) = -rho_nc(j)*dvdx - rho_nc(j)*v_nc(j)*dlogadx - v_nc(j)*drhodx;
        
        % Momentum eqn
        dvdt_c(j) = -v_nc(j)*dvdx - (1/gamma)*(dtdx + (t_nc(j)/rho_nc(j))*drhodx);
        
        % Energy eqn
        dtdt_c(j) = -v_nc(j)*dtdx - (gamma - 1)*t_nc(j)*(dvdx + v_nc(j)*dlogadx);
        
    end
    
    % Compute avg time derivative
    
    drhodt = 0.5*(drhodt_p + drhodt_c);
    dvdt = 0.5*(dvdt_p + dvdt_c);
    dtdt = 0.5*(dtdt_p + dtdt_c);
    
    % Final soln update
    
    for i = 2:n-1
        rho_nc(i) = rho_old(i) + drhodt(i)*dt;
        v_nc(i) = v_old(i) + dvdt(i)*dt;
        t_nc(i) = t_old(i) + dtdt(i)*dt;
    end
    
    %% Apply boundary conditions
    %     Inlet
    v_nc(1) = 2*v_nc(2) - v_nc(3);
    
    % outlet
    v_nc(n) = 2*v_nc(n-1) - v_nc(n-2);
    rho_nc(n) = 2*rho_nc(n-1) - rho_nc(n-2);
    t_nc(n) = 2*t_nc(n-1) - t_nc(n-2);
    
    
    M_nc = v_nc./(t_nc.^0.5);  % Mach no
    P_nc = rho_nc.*t_nc;       % pressure
    mass_flow_nc = rho_nc.*v_nc.*a;
    
    
    rho_t_nc(k) = rho_nc(throat);
    t_t_nc(k) = t_nc(throat);
    P_t_nc(k) = P_nc(throat);
    M_t_nc(k) = M_nc(throat);
    
    %%
    figure(1);
    hold on
    
    if k == 50
        plot(x,mass_flow_nc,'color','g','linewidth',1.5)
    elseif k == 100
        plot(x,mass_flow_nc,'color','m','linewidth',1.5)
    elseif k == 150
        plot(x,mass_flow_nc,'color','b','linewidth',1.5)
    elseif k == 200
        plot(x,mass_flow_nc,'color','r','linewidth',1.5)
    elseif k == 700
        plot(x,mass_flow_nc,'color','c','linewidth',1.5)
    elseif k == 1400
        plot(x,mass_flow_nc,'o','color','bla','linewidth',1.5)
        
        xlabel('Distance through nozzle(x)-(Nondimensional)')
        ylabel('Mass Flow (Nondimensional)')
        grid minor
        legend('50 Deltat','100 Deltat','150 Deltat','200 Deltat','700 Deltat','1400 Deltat')
        title({['Non Conservative form'];['Instantaneous distributions of Mass flow through nozzle at different times during time-marching approach to steady state']})
    end
    
    
    
end

hold on
figure(2)
X = linspace(1,nt,nt);
subplot(2,1,1)
plot(X,rho_t_nc)
xlabel('no. of time steps')
ylabel('Density')
grid on

title({['Non Conservative form'];['Timewise variations of Density & Temperature at nozzle throat']})

subplot(2,1,2)
plot(X,t_t_nc)
xlabel('no. of time steps')
ylabel('Temperature')
grid on

figure(3)
subplot(2,1,1)
plot(X,P_t_nc)
xlabel('no. of time steps')
ylabel('Pressure')
grid on

title({['Non Conservative form'];['Timewise variations of Pressure & Mach no. at nozzle throat']})

subplot(2,1,2)
plot(X,M_t_nc)
xlabel('no. of time steps')
ylabel('Mach number')
grid on


end

 

VI.B ) Conservative Form : 

function [u1,u2,u3,mass_flow_c,rho,v,t,M,P,rho_t,t_t,P_t,M_t] = Conservative_form(n,x,dx,throat,gamma,c,a,nt)

% calculating initial profiles .

for i = 1:n
    if (0 <= x(i)) && (x(i) <= 0.5)
        rho(i) = 1;
        t(i) = 1;
    elseif (0.5 <= x(i)) && (x(i) <= 1.5)
        rho(i) = 1 - 0.366*(x(i) - 0.5);
        t(i) = 1 - 0.167*(x(i) - 0.5);
    elseif (1.5 <= x(i)) && (x(i) <= 3.5)
        rho(i) = 0.634 - 0.3879*(x(i) - 1.5);
        t(i) = 0.833 - 0.3507*(x(i) - 1.5);
    end
end

v = 0.59 ./(rho.*a);

u1 = rho.*a;
u2 = rho.*a.*v;
u3 = rho.*((t/(gamma - 1)) + (gamma/2).*v.^2).*a;
P = rho.*t;

for k = 1:nt
    
    %updating old values
    u1_old = u1;
    u2_old = u2;
    u3_old = u3;
    
    
    % time step
    A = (t.^0.5);
    dt = c.*(dx ./((A)+v));
    dt = min(dt);
    
    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);
    
    %predictor step
    for j = 2:n-1
        
        %j2 = ((gamma-1)/gamma)*(u3(j) - ((gamma/2)*(((u2(j)).^2)./u1(j))))*((log(a(j+1))-log(a(j)))/dx);
        j2(j) = (1/gamma)*rho(j).*t(j).*((a(j+1) - a(j))/dx);
        
        % continuity eqn
        du1dt_p(j) = -((f1(j+1)-f1(j))/dx);
        
        % momentum eqn
        du2dt_p(j) = -((f2(j+1) - f2(j))/dx) + j2(j);
        
        % energy eqn
        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
    
    %updating primitive variables from new u1 u2 u3 values
    rho = (u1)./(a);
    v = u2./u1;
    t = (gamma - 1)*(((u3)./(u1)) - (gamma/2)*((v).^2));
    
    % updating predicted flux term valies from new updated u1 u2 u3 values
    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);
    
    
    % Corrector step
    for j = 2:n-1
        
        %j2 = ((gamma-1)/gamma)*(u3(j) - ((gamma/2)*(((u2(j)).^2)./u1(j))))*((log(a(j))-log(a(j-1)))/dx);
        j2(j) = (1/gamma)*rho(j).*t(j).*((a(j) - a(j-1))/dx);
        
        % continuity eqn
        du1dt_c(j) = -((f1(j)-f1(j-1))/dx);
        
        % momentum eqn
        du2dt_c(j) = -((f2(j) - f2(j-1))/dx) + j2(j);
        
        % energy eqn
        du3dt_c(j) = -((f3(j) - f3(j-1))/dx);
    end
    
    % Compute avg time derivative
    du1dt = 0.5*(du1dt_p + du1dt_c );
    du2dt = 0.5*(du2dt_p + du2dt_c );
    du3dt = 0.5*(du3dt_p + du3dt_c );
    
    % final soln update
    for i = 2:n-1
        u1(i) = u1_old(i) + du1dt(i)*dt;
        u2(i) = u2_old(i) + du2dt(i)*dt;
        u3(i) = u3_old(i) + du3dt(i)*dt;
    end
    
    % Applying boundary conditions
    
    % inflow bc
    u1(1) = rho(1)*a(1);
    u2(1) = 2*u2(2) - u2(3);
    u3(1) = u1(1)*((t(1)./(gamma-1)) + ((gamma/2)*v(1).^2));
    
    %outflow bc
    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);
    
    % corrected values of primitive variables
    rho = (u1)./(a);
    v = u2./u1;
    t = (gamma - 1)*(((u3)./(u1)) - (gamma/2)*((v).^2));
    
    % solutions
    M = v./((t).^0.5);  % Mach no
    P = rho.*t;       % pressure
    mass_flow_c = rho.*v.*a;
    
    rho_t(k) = rho(throat);
    t_t(k) = t(throat);
    P_t(k) = P(throat);
    M_t(k) = M(throat);
    
    figure(4)
    hold on
    
    if k == 50
        plot(x,mass_flow_c,'color','g','linewidth',1.5)
    elseif k == 100
        plot(x,mass_flow_c,'color','m','linewidth',1.5)
    elseif k == 150
        plot(x,mass_flow_c,'color','b','linewidth',1.5)
    elseif k == 200
        plot(x,mass_flow_c,'color','r','linewidth',1.5)
    elseif k == 700
        plot(x,mass_flow_c,'color','c','linewidth',1.5)
    elseif k == 1400
        plot(x,mass_flow_c,'o','color','bla','linewidth',1.5)
        
        xlabel('Distance through nozzle(x)-(Nondimensional)')
        ylabel('Mass Flow (Nondimensional)')
        grid minor
        legend('50 Deltat','100 Deltat','150 Deltat','200 Deltat','700 Deltat','1400 Deltat')
        title({['Conservative form'];['Instantaneous distributions of Mass flow through nozzle at different times during time-marching approach to steady state']})
    end
    
end


hold on
figure(5)
X = linspace(1,nt,nt);
subplot(2,1,1)
plot(X,rho_t)
xlabel('no. of time steps')
ylabel('Density')
grid on

title({['Conservative form'];['Timewise variations of Density & Temperature at nozzle throat']})

subplot(2,1,2)
plot(X,t_t)
xlabel('no. of time steps')
ylabel('Temperature')
grid on

figure(6)
subplot(2,1,1)
plot(X,P_t)
xlabel('no. of time steps')
ylabel('Pressure')
grid on

title({['Conservative form'];['Timewise variations of Pressure & Mach no. at nozzle throat']})

subplot(2,1,2)
plot(X,M_t)
xlabel('no. of time steps')
ylabel('Mach number')
grid on



end

 

 The MacCormacks technique is implemented in both the Non-conservative and Conservative functons.

 

VII. Results : 

VII. A) Timewise variations of density, temperature, velocity and Mach no. at nozzle throat (31 grid)

1. Non-conservative form : 

               

                

 

2. Conservative form :

                

                  

Initial disturbances occur in both the form's variatons but,

We can observe that the Timewise variations in the conservative form is less as compared to the Non conservative form. 

 

VII. B)  Instantaneous distributions of Mass flow through nozzle at different times during time marching approach to steady state (31 grid)

1. Non conservative form

 

2. Conservative form : 

The mass flow distributions through the nozzle vary alot in both forms, however in Conservative solution the variation tends towards the steady state approach abit more than the non-conservative solution. It is a point to note that the initial conditions were  different for both forms.

 

VII. C) Comparison of mass flows between Conservative and Non-conservative forms (31 grid):

                   

Observing the plot above, mass flows of both solutions are very close to the analytical value 0.579 as given in the book. however, The conservation form gives a distribution much closer to being a constant. In contrast, the non conservation results have(on a magnified scale) a sizeable variation, with some spurious oscillations at both inflow and outflow boundaries. 

This shows that the Conservation form does a better job of preserving the mass throughout the flow field. whereas in Nonconservation form the massflow is actually a secondary result due to equations used to obtain above results. 

 

VII. D) Grid dependence test : 

We will increase the grid size from 31 to 61 to see if the solution is dependent or independent of grid size.  

              

All the values are fairly close to the Exact analytical solution however doubling the grid points did improve the numerical solution marginally. We notice that the non-conservative form values are a tad-bit more close to analytical solution than the conservation form.

The same can be said is true for all locations within the nozzle. in otherwords, the steady state numerical solutions are essentially same and we can therefore conclude that our original 31 grid points solution is essentially grid-independent for all practical purposes. 

The grid independent solution does not agree exactly with the analytical results but it is close enough and have great accuracy. 

   

VII. E) Minimum number of cycles for which the simulation should be run to achieve convergence and its Computational time : 

        

Inspite the the fact that conservation form requires more tedious calculations and work it surprisingly converged faster than the Non-conservation form and Minimum cycles for which the simulation should be run displayed. These values are approximated from the graphs 

 

VII. F) Similar Plots for 61 Grid points 

1. Nonconservation form                        

                       

2. Conservation form :

                       

                        

VII. G)  Instantaneous distributions of Mass flow through nozzle at different times during time marching approach to steady state (61 grid)

1. Nonconservative form : 

 

2. Conservative form : 

 

VIII. Conclusions : 

1. There is no Superiority of either form of governing equations over the other for the given flow problem in terms of accuracy of results.  The non-conservation form seems to produce slightly more accurate results for primitive variables and conservation form seems to produce sightly more accurate results for flux variables. 

2. The conservation form yields a better mass flow distribution and does a better job of conserving mass. 

3. The non-conservation form leads to small residuals. 

4. The Conservative method came out to be faster in this case. 

5. The results in either case are certainly satisfactory.

6. The solution is grid independent . 

 

IX : References 

1. Computational fluid dynamics- basics with applications by John D anderson

 

                                                             The End 

keywords - MATLAB, MATLAB-BASICS, NUMERICAL-ANALYSIS, CFD, MACCORMACK, NOZZLE FLOW