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

AIM

The aim of this project is to simulate one dimensional super-sonic flow using Macormack method

OBJECTIVES

The objectives of the project are to

• write a program to solve the conservative and the non-conservative form of equations
• perform a comparison study between the normalised mass flow rate of both the forms of equations
• find the number of iterations for the convergence of the solutions

ONE DIMENSIONAL SUPERSONIC FLOW OVER A CONVERGENT-DIVERGENT NOZZLE

A nozzle is a device which is used to control the characteristics of fluid flow as it exits a closed chamber. Now in the above figure we can see that we are seeing a convergent-divergent nozzle in which we have assumed that the flow parameters are only changing in the X direction only. Such a condition where the parameters change only in one direction is known as "Quasi One Dimensional Flow. The flow at the inlet comes from a reservoir with large area of cross-section and we can also see that the flow is in subsonic conditions at the convergent part, as we move towards divergent section the flow turns out to be in supersonic state.

NON CONSERVATIVE FORM OF EQUATIONS (MACORMACK METHOD)

 The Non- conservative form of the governing equations indicates that the control volume is moving. So lets define the equations and we are to take a note that all the variables defined are in the non-dimensional form

Continuity equation:- `frac[delrho'][delt']=-rho'frac[delV'][delx']-rho'V'frac[del(ln A')][delx']-V'frac[delrho'][delx']`

Momentum Equation:- `frac[delV'][delt']=-Vfrac[delV'][delx']-1/gamma(frac[delT'][delx']+frac[T'][rho']frac[delrho'][delx'])`

Energy Equation:- `frac[delT'][delt']=-Vfrac[delT'][delx']-(gamma-1)T'[frac[delV'][delx']+V'frac[del(ln A')][delx']]`

To solve the above equations we have to apply Macormack method to form the difference equations. Macormack method is a numerical technique which uses explicit scheme of Finite Difference method to solve the equations. So, the first step would be to apply the predictor step in Macormack technique.  The MacCormack method is an explicit finite difference method that is 2nd order accurate in both time and space. The predictor step we use forward differencing scheme for the first order time derivative terms. We have to calculate the the flow variables from the previous time step. In the corrector step we use bacward difference scheme for the time derivative terms and then compute the corrected values of the flow field. Finally we find the average values of both the predictor and corrector step for steady state values. 

Predictor step

`(frac[delrho][delt])^t=-rho_i frac[V_(i+1)^t-V_i^t][trianglex]-rho_i^tV_i^tfrac[ln A_(i+1)-lnA_i][trianglex]-V_i^tfrac[rho_(i+1)^t-rho_i^t][trianglex]`

`(frac[delV][delt])_i^t=-V_i^tfrac[V_(i+1)^t-V^t][trianglex]-1/gamma(frac[T_(i+1)^t-T_i^t][trianglex]+frac[T_i^t][rho_i^t]frac[rho_(i+1)^t-rho_i^t][trianglex])`

`(frac[delT][delt])_i^t=-V_i^tfrac[T_(i+1)^t-T_i^t][trianglex]-(gamma-1)T_i^t(frac[V_(i+1)^t-V_I^t][trianglex]+V_i^tfrac[ln A_(i+1)-ln A_i][trianglex])`

 Now, we obtain predicted values of `rho`,V,T as follows

`barrho_i^(t+trianglet)=rho_i^t+(frac[delrho][delt])_i^t trianglet`

`barV_i^(t+trianglet)=V_i^t+(frac[delV][delt])_i^t trianglet`

`barT_i^(t+trianglet)=T_i^t+(frac[delT][delt])_i^t trianglet`

Corrector Step

`(frac[delrho][delt])_i^(t+trianglet)=-rho_i frac[V_(i+1)^(t+trianglet)-V_i^(t+trianglet)][trianglex]-rho_i^(t+trianglet)V_i^(t+trianglet)frac[ln A_(i+1)-lnA_i][trianglex]-V_i^(t+trianglet)frac[rho_(i+1)^(t+trianglet)-rho_i^(t+trianglet)][trianglex]`

`(frac[delV][delt])_i^(t+trianglet)=-V_i^(t+trianglet)frac[V_(i+1)^(t+trianglet)-V^(t+trianglet)][trianglex]-1/gamma(frac[T_(i+1)^(t+trianglet)-T_i^(t+trianglet)][trianglex]+frac[T_i^(t+trianglet)][rho_i^(t+trianglet)]frac[rho_(i+1)^(t+trianglet)-rho_i^(t+trianglet)][trianglex])`

`(frac[delT][delt])_i^(t+trianglet)=-V_i^(t+trianglet)frac[T_(i+1)^(t+trianglet)-T_i^(t+trianglet)][trianglex]-(gamma-1)T_i^(t+trianglet)(frac[V_(i+1)^(t+trianglet)-V_i^(t+trianglet)][trianglex]+V_i^(t+trianglet)frac[ln A_(i+1)-ln A_i][trianglex])`

After the predictor and corrector step we need to compute the average of the predictor and corrector derivatives

`(frac[delrho][delt])_(av)=0.5[(frac[delrho][delt])_i^t+(frac[delrho][delt])_i^(t+trianglet]]`

`(frac[delV][delt])_(av)=0.5[(frac[delV][delt])_i^t+(frac[delV][delt])_i^(t+trianglet]]`

`(frac[delT][delt])_(av)=0.5[(frac[delT][delt])_i^t+(frac[delT][delt])_i^(t+trianglet]]`

MATLAB CODE 

Main Program

clear all;
close all;
clc;

n=31;
x=linspace(0,3,n);
m=0.579*ones(n);

tic;
% solving equations in the non-conservation form
[rho_nc,T_nc,p_nc,M_nc,m_nc,n_nc] = non_conservative_form;
time_nc=toc;

tic;
% solving equations in the consevtion form
[rho_c,T_c,p_c,M_c,m_c,n_c] = conservative_form;
time_c=toc;

% Comparing the normalized mass flow rate between conservation & non-conservation forms
figure(7)
plot(x,m_nc,x,m_c,x,m,'g');
xlabel('Length');
ylabel('Mass flow rate');
title(sprintf('Comparison of mass flow variations of conservative and non-cionservative form'));
legend('Non-conservation form','Conservation form','Analytical solution');

FUNCTION CODE FOR NON CONSERVATIVE FORM

function [rho,t,p,M,m,count]= non_conservative_form

    % initialization
    n=31;
    nt=1400;
    x = linspace(0,3,n);
    dx = x(2)-x(1); 
    A = 1 + 2.2*(x-1.5).^2;
    throat = (n+1)/2;   
    nt_x = linspace(1,nt,nt);
    gamma = 1.4;
    error = 1; 
    count = 1; 
    m_old = 0;
    C=0.5; 
    tol=1e-4;

    % Initial conditions
    rho = 1 - 0.3146*x;
    t = 1 - 0.2314*x;
    v = (0.1 + 1.09.*x).*t.^(1/2);
    m= 0.579*ones(n);
    m_initial=rho.*A.*v;

    % loop to calculate dt
    for k = 1:n
        dt_x(k) = C*dx/(t(k)^(1/2)+v(k));
    end
    dt = min(dt_x);

    %outer time loop
    for j = 1:nt
        
        % updating old values
        rho_old = rho;
        v_old = v;
        t_old = t;
        % predictor step
        for i = 2:n-1
            drho_dx = (rho(i+1) - rho(i))/dx;
            dv_dx = (v(i+1) - v(i))/dx;
            dt_dx = (t(i+1) - t(i))/dx;
            dlogA_dx = (log(A(i+1)) - log(A(i)))/dx;

            drho_dtp(i) = -rho(i)*dv_dx - rho(i)*v(i)*dlogA_dx - v(i)*drho_dx; % continuity equation
            dv_dtp(i) = -v(i)*dv_dx - (1/gamma)*(dt_dx+(t(i)/rho(i))*drho_dx); % momentum equation
            dt_dtp(i) = -v(i)*dt_dx - (gamma-1)*t(i)*(dv_dx + v(i)*dlogA_dx);% energy equation

            % predictor step update
            rho(i) = rho_old(i) + drho_dtp(i)*dt;
            v(i) = v_old(i) + dv_dtp(i)*dt;
            t(i) = t_old(i) + dt_dtp(i)*dt;
        end

        % corrector step
        for i = 2:n-1
            drho_dx = (rho(i) - rho(i-1))/dx;
            dv_dx = (v(i) - v(i-1))/dx;
            dt_dx = (t(i) - t(i-1))/dx;
            dlogA_dx = (log(A(i))-log(A(i-1)))/dx;

            drho_dtc(i) = - rho(i)*dv_dx - rho(i)*v(i)*dlogA_dx - v(i)*drho_dx; % continuity equation
            dv_dtc(i) = - v(i)*dv_dx - (1/gamma)*(dt_dx + (t(i)/rho(i))*drho_dx); % momentum equation
            dt_dtc(i) = - v(i)*dt_dx - (gamma-1)*t(i)*(dv_dx + v(i)*dlogA_dx); % energy equation
        end

        % average value of time derivatives
        drho_dt = 0.5*(drho_dtp + drho_dtc);
        dv_dt = 0.5*(dv_dtp + dv_dtc);
        dt_dt = 0.5*(dt_dtp + dt_dtc);
        
        % corrector step update
        for i = 2:n-1            
            rho(i) = rho_old(i) + drho_dt(i)*dt;
            v(i) = v_old(i) + dv_dt(i)*dt;
            t(i) = t_old(i) + dt_dt(i)*dt;
        end

        % Inlet boundary condition
        v(1) = 2*v(2) - v(3);

        % Outlet boundary 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);

        % Calculating pressure, mass flow rate and Mach number
        p = rho.*t;
        M = v./t.^0.5;
        m = rho.*v.*A;
        m_t(j,:) = rho.*v.*A; 
        
        % Capturing flow field variables at the throat
        rho_throat(j) = rho(throat);
        T_throat(j) = t(throat);
        p_throat(j) = p(throat);
        M_throat(j) = M(throat);

        % convergence
        error = max(abs(m-m_old));
        if error>tol
            count = count + 1;            
            m_old = m;
        end

    end
    
    % plotting the solutions
    figure(1);
    title(sprintf('Flow-field variables across nozzle in Non-conservation form n No. of time steps = %g',nt));
    subplot(4,1,1);
    plot(x,M,'b','linewidth',2);
    xlabel('Length');
    ylabel('Mach no');
    title(sprintf('Mach number'));
    subplot(4,1,2);
    plot(x,rho,'g','linewidth',2);
    xlabel('Length');
    ylabel('Density');
    title(sprintf('Non-dimensional Density'));
    subplot(4,1,3);
    plot(x,t,'m','linewidth',2);
    xlabel('Length');
    ylabel('Temperature');
    title(sprintf('Non-dimensional Temperature'));
    subplot(4,1,4);
    plot(x,p,'r','linewidth',2);
    xlabel('Length');
    ylabel('Pressure');
    title(sprintf('Non-dimensional Pressure'));

    % Variarion of flow field variables at nozzle throat
    figure(2)
    plot(nt_x,rho_throat,'b','linewidth',2);
    hold on
    plot(nt_x,T_throat,'r','linewidth',2);
    plot(nt_x,p_throat,'g','linewidth',2);
    plot(nt_x,M_throat,'y','linewidth',2);
    title(sprintf(['Timewise variation of Flow-field variables at nozzle throat Non-conservation form']));
    xlabel('Number of time steps');
    legend('Density','Temperature','Pressure','Mach number');

    % Mass flow variation across nozzle at different time steps
    figure(3)
    plot(x,m_initial,x,m_t(50,:),x,m_t(100,:),x,m_t(150,:),x,m_t(200,:),x,m_t(700,:),x,m,'r');
    xlabel('Non-dimensional distance through nozzle');
    ylabel('Non-dimensional mass flow');
    legend('0Deltat','50Deltat','100Deltat','150Deltat','200Deltat','700Deltat','Analytical solution');
    title(sprintf('Non-dimensional mass flow rate across the nozzle Non-conservation form'));
end

When we go through the function code we can see that the first we define a user function by using function command. The subsequent lines of the code we define some of the parameters such as 

• Define mesh (dx) = x(2)-x(1)
• Area of nozzle (a) = 1+2.2(x-1.5)^2
• Ratio of specific heat(gamma) = 1.4
• CFL Number (C) = 0.5
• Density(rho) = 1-0.3146x
• Temperature= 1-0.2314x
• Velocity(v) = (0.1+1.09x)T^0.5
• mass flow rate (m_initial) = rho*a*v

So when we enter the time loop we use the for command and define the time loop. Next in the predictor step we define the spacial loop for the inner nodes and write the continuity, momentum and energy equations by defining some standard derivative terms like dvdx, drhodx,dtdx,dlogadx. In the same fashion the corrector step is also defined and then we write down the equations to compute the average derivatives. Next we define the boundary conditions at the inlet and outlet of the nozzle and compute the flow parameters at the throat section. We then define the convergence criteria for the mass flow rate and then plot the following graphs

• Mass flow rate versus length of nozzle
• timewise variations of flow parameters
• flow behaviour throughout the length of nozzle.

RESULTS AND DISCUSSIONS

Qualitative aspects of flow field parameters

The above plot shows how the flow parameters inside the convergent divergent nozzle is taking place and we can infer that the flow parameters changes inside the nozzle configuration. The variation of pressure, density, and temperature shows there is a decrease in the plot which shows us that the nozzle pressure, density and temperature is dropping throughout the length of the nozzle. The Mach number increases because the flow is sub sonic at the inlet boundary condition. As we progress along the length of convergent- divergent nozzle the flow has become a supersonic flow towards the outlet boundary conditions. The normalised pressure ratio across the nozzle exerts a force on the gas to accelerate it through the nozzle. The variation of the normalised flow parameters can be expressed as a function of the Mach number which are defined by

Pressure Ratio`p/p_o=(1+frac[gamma-1][2]M^2)^(frac[-gamma][gamma-1`

Density Ratio `rho/rho_o=(1+frac[gamma-1][2]M^2)^(frac[-1][gamma-1`

 Temperature Ratio `T/T_o=(1+frac[gamma-1][2]M^2)^-1`

Timewise variations of flow field parameters

Now we have to have a interest in the timewise variations of the flow field parameters, so to give a basic idea which shows us the variation of density, temperature pressure and mach number at the nozzle throat plotted with the number of time steps. The initial condition for the flow field parameters start from zero and the time step ends at 1400. With increasing time as we progress from left to right we can see that the time marching plots where the time is increasing from 0 and ends at 1400. We can see in the plots that the changes are taking place in early stages of the time and as we move onto to the right hand side of the plots we are able to see that the solution has reached the steady state value of the flow field. The flow field variables are driven by a stronger potential which changes rapidly. At the later stage as the steady state is approached, the changes become much more smaller as the time progresses. we can also see that there is no artificial viscosity has been added for the simulation.

Mass flow rate plot versus length of the nozzle

The insight for another concept of fluid mechanics for timewise variation of the flow is provided by the mass flow variation as shown in the above plot. There are 7 different curves which suggests the behaviour of the mass flow rate which is coming into the nozzle. The 7 curves shows us different time step in the time-marching solution. The first curve which is labelled as "time step=0" it shows the initial conditions at which the flow occurs. Next we have the time step=50 which indicates after 50 time steps the mass flow variation has changed by a certain amount. After 100 time steps the mass flow distribution has changed radically. As and when we have increased the number of time steps from 150, 200, 500 we are able to see that the behaviour of the mass flow begins to settle down and when the time step is reached to 700 we can see that the mass flow variation is a straight horizontal line in the graph. This indicates that it has reached a steady state value throughout the nozzle.

CONSERVATIVE FORM OF EQUATIONS IN (MACORMACK METHOD)

The Non- conservative form of the governing equations indicates that the control volume is moving. So lets define the equations and we are to take a note that all the variables defined are in the non-dimensional form

Continuity Equation `frac[del(rho'A')][delt']+frac[del(rho'A'V')][delx]=0`

Momentum Equation `frac[del(rho'A'V')][delt']+frac[del[rho'A'V'+(1/gamma)p'A']][delx']=1/gamma p'frac[delA'][delx']`

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

the above equations are conservation form of continuity, momentum and energy equations for quasi 1D flow. Now let us define the elements of solution vector U, the flux Vector F and source term J as follows

`U_1=rho'A'`

`U_2=rho'A'V'`

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

`F_1=rho'A'V'`

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

`F_3=rho'(frac[e'][gamma-1]+gamma/2V'^2)V'A'+rho'A'V'`

`J_2=1/gammap"frac[delA'][delx']`

With these updates we can write the equations of the conservation form as

`frac[delU_1][delt']=-frac[delF_1][delx']`

`frac[delU_2][delt']=-frac[delF_2][delx']+J_2`

`frac[delU_3][delt']=-frac[delF_3][delx']`

The above equations represent the conservative form of quasi 1D flow and these are the equations we wish to solve. Now, we may note that the numerical solution will give us values for U1, U2, and U3 which is known as solution vector. In order to obtain the primitive variables we must decode the solution vector as follows.

`rho'=U_1/A'`

`V'=U_2/U_1`

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

`p'=rho'T'`

Now we have to compute the pure form of the flux terms in order to express primitive variables solution from the solution vectors and the flux terms 

`F_1=U_2`

`F_2=U_2^2/U_1+frac[gamma-1][gamma](U_3-frac[gammaU_2^2][2U_1])`

`F_3=gammafrac[U_2U_3][U_1]-frac[gamma(gamma-1)][2]frac[U_2^3][U_1^2]`

`J_2=[gamma-1]/gamma(U_3-frac[gammaU_2^2][2U_1])frac[del(ln A')][delx']`

Predictor Step

`((del U_1)/(del t'))_i^t = -(F_(1(i+1))^t - F_(1(i))^t)/(triangle x')`

`((del U_2)/(del t'))_i^t = (F_(2(i+1))^t - F_(2(i))^t)/(triangle x')+J_2`

`((del U_3)/(del t'))_i^t = -(F_(3(i+1))^t - F_(3(i))^t)/(triangle x')`

Corrector Step

Now the solution vector for the must be updated from the above equations for the next time step

`bar(U)_(1(i))^(t+Delta t) = ((del U_1)/(del t'))_i^t .triangle t`

`bar(U)_(2(i))^(t+Delta t) = ((del U_2)/(del t'))_i^t .triangle t`

`(bar(del U_3)/(del t'))_i^(t+triangle t) = -(bar(F)_(3(i))^(t+triangle t) - bar(F)_(3(i-1))^(t+triangle t))/(triangle x')`

After the predictor and corrector step we calculate the average of the derivatives calculated from the predictor and corrector step.

`((del U_1)/(del t'))_(av)=1/2[((del U_1)/(del t'))_i^t+(bar(del U_1)/(del t'))_i^(t+trianglet)]`

`((del U_2)/(del t'))_(av)=1/2[((del U_2)/(del t'))_i^t+(bar(del U_2)/(del t'))_i^(t+triangle t))`

`((del U_3)/(del t'))_(av)=1/2[((del U_3)/(del t'))_i^t+(bar(del U_3)/(del t'))_i^(t+triangle t)]`

FUNCTION CODE FOR CONSERVATIVE FORM

function [rho,t,p,M,m,count]= conservative_form
    
    % Initialization
    n=31;
    nt=1400;
    C=0.5;
    tol=1e-4;
    x = linspace(0,3,n);
    dx = x(2)-x(1);
    A = 1 + 2.2*(x-1.5).^2;
    throat = (n+1)/2; 
    nt_x = linspace(1,nt,nt);
    gamma = 1.4;
    error = 1;
    count = 1;
    m_old = 0;
    m = 0.579*ones(n);

    % Initial conditions
    u2_i = 0.59;  
    for i= 1:n
        if ((x(i)>=0) && (x(i)<=0.5))
            rho(i) = 1;
            t(i) = 1;
        elseif ((x(i)<1.5) && (x(i)>0.5))
            rho(i) = 1 - 0.366*(x(i)-0.5);
            t(i) = 1 - 0.167*(x(i)-0.5);
        elseif ((x(i)>=1.5) && (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
        v(i) = u2_i/(rho(i)*A(i));
        u1(i) = rho(i)*A(i);
        u2(i) = rho(i)*A(i)*v(i);
        u3(i) = rho(i)*A(i)*(t(i)/(gamma-1) + (gamma/2)*v(i)^2);
    end
      
    % copy of initial values
    u1_initial = u1;
    u2_initial = u2;
    u3_initial = u3;

    % loop to calculate dt
    for k = 1:n
        dt_x(k) = C*dx/(t(k)^(1/2)+v(k));
    end
    dt = min(dt_x);

    % time loop
    for j = 1:nt
        u1_old = u1;
        u2_old = u2;
        u3_old = u3;
        
        % predictor step
        for i = 2:n-1
            % computing flux and source terms
            f1_i = u2(i);
            f1_ip1 = u2(i+1);
            f2_i = u2(i)^2/u1(i) + ((gamma-1)/gamma)*(u3(i) - (gamma/2)*u2(i)^2/u1(i));
            f2_ip1 = u2(i+1)^2/u1(i+1) + ((gamma-1)/gamma)*(u3(i+1) - (gamma/2)*u2(i+1)^2/u1(i+1));
            j2  = (1/gamma)*rho(i)*t(i)*((A(i+1) - A(i))/dx);
            f3_i = gamma*u2(i)*u3(i)/u1(i) - (gamma*(gamma-1)/2)*u2(i)^3/u1(i)^2;
            f3_ip1 = gamma*u2(i+1)*u3(i+1)/u1(i+1) - (gamma*(gamma-1)/2)*u2(i+1)^3/u1(i+1)^2;

            du1_dtp(i) = -(f1_ip1 - f1_i)/dx;    % continuity eq      
            du2_dtp(i) = j2 - (f2_ip1 - f2_i)/dx;  % momentum eq
            du3_dtp(i) = -(f3_ip1 - f3_i)/dx;    % energy eq

            % updating the predicted values
            u1(i) = u1(i) + du1_dtp(i)*dt;
            u2(i) = u2(i) + du2_dtp(i)*dt;
            u3(i) = u3(i) +du3_dtp(i)*dt;
        end

        % Updating primitive variables after predictor step
        rho = u1./A;
        v = u2./u1;
        t = (gamma-1)*((u3./u1) - (gamma/2)*v.^2);
    
        % corrector step
        for i = 2:n-1
            f1_i = u2(i);
            f1_ic1 = u2(i-1);
            f2_i = u2(i)^2/u1(i) + ((gamma-1)/gamma)*(u3(i) - (gamma/2)*u2(i)^2/u1(i));
            f2_ic1 = u2(i-1)^2/u1(i-1) + ((gamma-1)/gamma)*(u3(i-1) - (gamma/2)*u2(i-1)^2/u1(i-1));
            j2  = (1/gamma)*rho(i)*t(i)*((A(i) - A(i-1))/dx);
            f3_i = gamma*u2(i)*u3(i)/u1(i) - (gamma*(gamma-1)/2)*u2(i)^3/u1(i)^2;
            f3_ic1 = gamma*u2(i-1)*u3(i-1)/u1(i-1) - (gamma*(gamma-1)/2)*u2(i-1)^3/u1(i-1)^2;

            du1_dtc(i) = -(f1_i - f1_ic1)/dx;        
            du2_dtc(i) = j2 - (f2_i - f2_ic1)/dx;        
            du3_dtc(i) = -(f3_i - f3_ic1)/dx;
        end

        %average value
        du1_dt = 0.5*(du1_dtp + du1_dtc);
        du2_dt = 0.5*(du2_dtp + du2_dtc);
        du3_dt = 0.5*(du3_dtp + du3_dtc);

        for i = 2:n-1
            % computing corrected values
            u1(i) = u1_old(i) + du1_dt(i)*dt;
            u2(i) = u2_old(i) + du2_dt(i)*dt;
            u3(i) = u3_old(i) + du3_dt(i)*dt;
        end
        
        % Inlet BC
        u1(1) = rho(1)*A(1);
        u2(1) = 2*u2(2) - u2(3);
        v(1) = u2(1)/(rho(1)*A(1));
        u3(1) = u1(1)*(t(1)/(gamma-1) + (gamma/2)*v(1)^2);

        % Outlet 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);
        
        % computing the values of the primitive variables
        for i= 1:n
            rho(i) = u1(i)/A(i);
            v(i) = u2(i)/u1(i);
            t(i) = (gamma-1)*(u3(i)/u1(i) - gamma/2*v(i)^2);
        end

        % Pressure, mass flow rate and Mach number
        p = rho.*t;
        M = v./t.^0.5;
        m = rho.*v.*A;
        m_t(j,:) = rho.*v.*A; 
        
        % Flow field variables at the throat
        rho_throat(j) = rho(throat);
        T_throat(j) = t(throat);
        p_throat(j) = p(throat);
        M_throat(j) = M(throat);

        % Convergence
        error = max(abs(m-m_old));
        if error>tol
            count = count + 1;            
            m_old = m;
        end
    end
    
    % plotting steady state solutions of flow-field variables
    figure(4);
    title(sprintf('Flow-field variables across nozzle Conservation form n No. of time steps = %g',nt));
    subplot(4,1,1);
    plot(x,M,'b','linewidth',2);
    xlabel('Length');
    ylabel('Mach number');
    title(sprintf('Mach number'));
    subplot(4,1,2);
    plot(x,rho,'g','linewidth',2);
    xlabel('Length');
    ylabel('Density');
    title(sprintf('Non-dimensional Density'));
    subplot(4,1,3);
    plot(x,t,'m','linewidth',2);
    xlabel('Length');
    ylabel('Temperature');
    title(sprintf('Non-dimensional Temperature'));
    subplot(4,1,4);
    plot(x,p,'r','linewidth',2);
    xlabel('Length');
    ylabel('Pressure');
    title(sprintf('Non-dimensional Pressure'));

    %Variarion of flow field variables at nozzle throat
    figure(5)
    plot(nt_x,rho_throat,'b','linewidth',2);
    hold on
    plot(nt_x,T_throat,'r','linewidth',2);
    plot(nt_x,p_throat,'g','linewidth',2);
    plot(nt_x,M_throat,'y','linewidth',2);
    title(sprintf('Timewise variation of Flow-field variables at nozzle throat conservative form'));
    xlabel('Number of time steps');
    legend('Density','Temperature','Pressure','Mach number');

    % Mass flow variation across nozzle at different time steps
    figure(6)
    plot(x,u2_initial,'--',x,m_t(50,:),x,m_t(100,:),x,m_t(150,:),x,m_t(200,:),x,m_t(700,:),x,m,':');
    xlabel('Non-dimensional distance through nozzle');
    ylabel('Non-dimensional mass flow');
    legend('0Deltat','50Deltat','100Deltat','150Deltat','200Deltat','700Deltat','Analytical solution');
    title(sprintf('Mass flow rate across the nozzle Conservation form'));
end

When we go through the function code we can see that the first we define a user function by using function command. The subsequent lines of the code we define some of the parameters such as 

• Define mesh (dx) = x(2)-x(1)
• Area of nozzle (a) = 1+2.2(x-1.5)^2
• Ratio of specific heat(gamma) = 1.4
• CFL Number (C) = 0.5
• Density(rho) = 1-0.3146x
• Temperature= 1-0.2314x
• Velocity(v) = (0.1+1.09x)T^0.5
• mass flow rate (m_initial) = rho*a*v

So when we enter the time loop we use the for command and define the time loop. Next in the predictor step we define the spacial loop for the inner nodes and write the continuity, momentum and energy equations by defining some standard derivative terms like solution vector U, flux vector F and the source term J. In the same fashion the corrector step is also defined and then we write down the equations to compute the average derivatives. Next we define the boundary conditions at the inlet and outlet of the nozzle and compute the flow parameters at the throat section. We then define the convergence criteria for the mass flow rate and then plot the following graphs

• Mass flow rate versus length of nozzle
• timewise variations of flow parameters
• flow behaviour throughout the length of nozzle.

RESULTS AND DISCUSSIONS

Qualitative aspects of flow field parameters

The above plot shows how the flow parameters inside the convergent divergent nozzle is taking place and we can infer that the flow parameters changes inside the nozzle configuration. The variation of pressure, density, and temperature shows there is a decreases in the plot which shows us that the nozzle pressure, density and temperature is dropping throughout the length of the nozzle. The Mach number increases because the flow is sub sonic at the inlet boundary condition. As we progress along the length of convergent- divergent nozzle the flow has become a supersonic flow towards the outlet boundary conditions. The normalised pressure ratio across the nozzle exerts a force on the gas to accelerate it through the nozzle. The variation of the normalised flow parameters can be expressed as a function of the Mach number which are defined by

Pressure Ratio``p/p_o=(1+frac[gamma-1][2]M^2)^(frac[-gamma][gamma-1``

Density Ratio `rho/rho_o=(1+frac[gamma-1][2]M^2)^(frac[-1][gamma-1`

 Temperature Ratio `T/T_o=(1+frac[gamma-1][2]M^2)^-1`

Timewise variations of flow field parameters

Now we have to have a interest in the timewise variations of the flow field parameters, so to give a basic idea which shows us the variation of density, temperature pressure and mach number at the nozzle throat plotted with the number of time steps. The initial condition for the flow field parameters start from zero and the time step ends at 1400. With increasing time as we progress from left to right we can see that the time marching plots where the time is increasing from 0 and ends at 1400. We can see in the plots that the changes are taking place in early stages of the time and as we move onto to the right hand side of the plots we are able to see that the solution has reached the steady state value of the flow field. The flow field variables are driven by a stronger potential which changes rapidly. At the later stage as the steady state is approached, the changes become much more smaller as the time progresses. we can also see that there is no artificial viscosity has been added for the simulation.

Mass flow rate plot versus length of the nozzle

The insight for another concept of fluid mechanics for timewise variation of the flow is provided by the mass flow variation as shown in the above plot. There are 7 different curves which suggests the behaviour of the mass flow rate which is coming into the nozzle. The 7 curves shows us different time step in the time-marching solution. The first curve which is labelled as "time step=0" it shows the initial conditions at which the flow occurs. Next we have the time step=50 which indicates after 50 time steps the mass flow variation has changed by a certain amount. After 100 time steps the mass flow distribution has some kind of a humped distribution . As and when we have increased the number of time steps from 150, 200, 500 we are able to see that the behaviour of the mass flow begins to settle down and when the time step is reached to 700 we can see that the mass flow variation is a straight horizontal line in the graph. This indicates that it has reached a steady state value throughout the nozzle.

COMPARITIVE STUDY BETWEEN THE CONSERVATIVE AND NON CONSERVATIVE FORM FOR NON DIMENSIONAL MASS FLOW RATE

Now to perform the comparison study lets go back to those plots which shows us the mass flow rates of both the forms of governing equations. At time step zero both the plots shows us that the plots for mass flow distribution starts of at the same point as it has but as and when we increase the number of time steps the distributions converge to steady state value. In the non conservative form the distribution of mass flow is quite scattered and in that of conservative form the scattering of the plot is lot organised. The conservation does a better plot of preserving the mass throughout the flow field because mass itself is a dependent variable in the equation. Dependent variables on the non conservative form of equations are primitive variables and the mass flow is obtained only as a secondary result. The mass flow rate plots are less fluctuating in Conservative form and are closer to one another as compared to that of Non-conservative form.

CONCLUSION

We can draw some of the conclusions from the simulation to be as follows

1. The mass flow rate plots is less scattered in the conservative form and are closer to one another compared to that of non-conservative form where the plot is more scattered in nature.
2. The conservative form tends to be more accurate in the numerical discretization since the plot is close to the analytical solution and the non conservative part tends to move away from the exact solution and the scattering is more in non-conservative form.
3. Now for the time step of 1400 we can see that in the timewise variations plots in both the form of equations. The non conservative form takes approximately 550 time steps to reach steady state while the conservative form takes approximately 500 time step to reach steady state. Hence we can infer that conservative form takes less number of iterations for convergent solution.
4. Time taken for the non-conservative form of equation(11.367 seconds) is less than the conservative form(15.505 seconds). Non-conservative form is faster but as the time step changes there are drastic changes in the mass flow rate distribution compared to conservative form which shows no drastic changes in the plots.
5. The non conservative form is faster than conservative form
6. Conservative form of equation involves more complex calculations and hence simulation is slower. But the accuracy is better than non-conservative form.