Week 7 - Simulation of a 1D Supersonic nozzle flow simulation using Macormack Method

AIM

To simulate the isentropic flow through a quasi 1D subsonic-supersonic nozzle by deriving both conservative and non-conservative forms of governing equation and solve them by using MacCormack's technique using MATLAB.

 

OBJECTIVE

• Need to determine the steady-state temperature distribution for the flow-field variables and to investigate the difference between the two forms of governing equations by comparing their solutions.
• To plot steady-state distribution of primitive variables inside the nozzle.
• To plot time-wise variation of the primitive variables.
• To plot variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process.
• To plot comparison of Normalized mass flow rate distributions of both forms.

 

INTRODUCTION

The flow of fluids through nozzles is an essential phenomenon in many industries, including aerospace, chemical engineering, and propulsion systems. Understanding the flow behavior is crucial for designing and optimizing these systems. In this regard, numerical simulations have emerged as powerful tools to predict the flow characteristics. In this project, we aim to simulate the isentropic flow through a quasi 1D subsonic-supersonic nozzle. The simulation will involve deriving both the conservation and non-conservation forms of the governing equations and solving them using the MacCormack's technique.

Additionally, we will determine the steady-state temperature distribution for the flow-field variables and investigate the difference between the two forms of governing equations by comparing their solutions. The report will present various plots, including the steady-state distribution of primitive variables inside the nozzle, time-wise variation of the primitive variables, variation of mass flow rate distribution inside the nozzle at different time steps during the time-marching process, and a comparison of normalized mass flow rate distributions of both forms. We will also perform a grid independence test on the solution and explain all our findings. Finally, we will discuss all the plots provided with observed inferences.

 

 

THEORY

 

Quasi one-dimensional flow

            Quasi one-dimensional flow refers to the fluid flow through a duct or channel, where the flow velocity and properties are assumed to vary only in one dimension, typically the axial direction. The flow is considered to be uniform in the transverse direction, and all flow properties, such as density, pressure, temperature, and velocity, are assumed to be functions of only one independent variable, usually the axial distance.

 

The assumption of quasi-one-dimensional flow simplifies the mathematical model used to analyze the fluid flow and enables the use of one-dimensional equations of motion. These equations describe the conservation of mass, momentum, and energy for the fluid flowing through the channel.

Quasi one-dimensional flow is widely used in the analysis of internal combustion engines, gas turbines, and aerospace propulsion systems. It provides a practical approach to analyze the fluid flow, and the results obtained can be used for the design and optimization of these systems.

 

Convergent Divergent Nozzle

            A nozzle is a type of device designed to control the direction or characteristics of fluid flow at the exit of a pipe (or such system). A nozzle is often a pipe or tube of varying cross-sectional area, thus can be used to direct or modify flow of a fluid.     

In the convergent section, the cross-sectional area of the nozzle decreases gradually, which increases the velocity of the fluid and decreases the pressure. The fluid then enters the throat section, where the cross-sectional area of the nozzle is the smallest. At this point, the fluid reaches its maximum velocity, known as the critical velocity.

In the divergent section, the cross-sectional area of the nozzle increases, which causes the velocity of the fluid to decrease while the pressure increases. This acceleration of the fluid in the nozzle converts the thermal energy of the fluid into kinetic energy, which is used to propel objects, such as rockets, at high speeds.

Convergent-divergent nozzles are widely used in aerospace propulsion systems, such as jet engines and rocket engines. They are also used in industrial applications, such as in the manufacture of plastics, where they are used to extrude molten plastic at high speeds. The design of convergent-divergent nozzles is critical to their performance, and their efficiency is determined by the ratio of the area of the throat to the area of the exit section.

The flow at the inlet of the nozzle comes from a reservoir where pressure and temperature are denoted as by  `ρ_0` and `T_0`  respectively. The cross sectional area of the reservoir is large (theoretically, A→∞) and hence the velocity is very small (V→0). Thus p₀ and T₀ are the stagnation values, or total pressure and temperature respectively. The flow expands isentropically to supersonic speeds at the nozzle exits, where the exit pressure, temperature, velocity, and Mach number are denoted by Pe, Te, Ve and Me respectively. The flow is locally subsonic in the convergent section, sonic at the throat section and supersonic at the divergent section. The sonic flow (M=1) at the throat means that the local velocity at this location is equal to the local speed of sound. Using an asterisk to denote sonic flow values, we have at the throat                 V = V∗ = a∗ Similarly, the sonic flow values of pressure and temperature are denoted by P∗ and T∗, respectively.

            Although the area of nozzle changes as a function of distance along the nozzle x and therefore in reality the flow field is two dimensional (xy), we make an assumption that the flow properties vary only with x, this tendency to assume uniform flow properties across any given cross-section. Such flow is defined as Quasi-one-dimensional flow.

            In order to solve quasi-one-dimensional problem, we need PDE or Integral form of equations for CFD analysis. These PDE or Integral form of the governing equations are in both conservative and non-conservative form.

            Before, writing the governing equations in conservative and non-conservative form, we need to non-dimensionalize the primitive variables for removing the constraints of dimensions to analyse the numerical solution. It will certainly ease the method of analysis for a given problem.

 

 

Mach number

     It is a dimensionless quantity representing the ratio of flow velocity to the local speed of the sound.

`M=u/c`

Where,

     M - Mach Number

     u - Flow Velocity

     c - Velocity Sound

 

            Supersonic and subsonic flow refer to the speed faster or slower than the speed of the sound. Anything going faster than the speed of the sound 343.2m/s, is supersonic flow. Anything going slower than the speed of the sound is called as the subsonic.

 

Convergent Section: M < 1 Subsonic flow

Throat Section: M = 1 Sonic flow

Divergent Section: M > 1 Supersonic flow

 

Governing Equation

     The Governing Equation for the Quasi One-Dimensional, Steady, Inviscid, Isentropic Flow is as follows,

Continuity Equation: `ρ_(1.) V_(1.) A_(1=) ρ_(2.) V_(2.) A_1`

Momentum Equation: `p_(1.) A_1+ρ_(1.) V_1^2.A_1+∫_(A_2)^(A_1)p.dA=p_(2.) A_(2+) ρ_(2.) V_2^2.A_2`

Energy Equation:  `h_1+(V_1^2)/2=h_2+(V_2^2)/2`

     We assume that flow is quasi one-dimensional flow i.e., at a given cross-sectional area A, all flow variables vary primarily along one direction, say x, the flow properties are uniform across this section. But in reality the flow field is two-dimensional (the flow varies in the two-dimensional x-y space).

 

Conservative and Non-Conservative Form

            In fluid dynamics, the governing equations that describe the behaviour of a fluid flow can be expressed in two forms: the conservative form and the non-conservative form.

The conservative form of the governing equations expresses the conservation of mass, momentum, and energy in terms of a set of conservation equations. These equations describe the rate of change of the fluid properties, such as density, velocity, and temperature, in terms of their gradients in space and time. The conservative form is derived by applying the principle of conservation of mass, momentum, and energy to a control volume in the fluid flow.

The non-conservative form of the governing equations expresses the same physical laws as the conservative form but in a different mathematical form. It describes the change in fluid properties in terms of the fluxes of these properties across the boundaries of a control volume. The non-conservative form is derived by applying the principle of conservation of mass, momentum, and energy to a differential element of the fluid flow.

Both the conservative and non-conservative forms of the governing equations are valid and can be used to solve fluid dynamics problems. However, the choice of the form depends on the specific application and the numerical method used to solve the equations. For example, the conservative form is often used in numerical methods that require a conservation law, while the non-conservative form is often used in methods that rely on the characteristics of the flow.

 

Non Dimensional Quantities

 

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

Temperature `T'=T/T_0`

Density  `ρ'=ρ/ρ_0`

Non - dimensional distance `X'=X/L`

Non dimensional velocity `V'=V/a_0` whereas a₀ is the local sonic velocity  `a_0=√(γRT_0 )`

Non-dimensional time  `t^'=t/((L/a_0 ) )` where  `L/a_0` is characteristic time scale.

Non dimensional area `A^'=A/A^**`

 

Conservative form :

    In the conservative form of equations, the variables that are being differentiated does not exist as coefficients of any derivative in the equation.

 

Continuity equation : 

`(∂(ρ'A'))/(∂t')+(∂(ρ'A'V'))/(∂x')=0`

 

Energy equation:

`(∂ρ^' (T^'/(γ-1)+γ/2 (V^' )^2 ) A^')/(∂t^' )+∂(ρ^' (T^'/(γ-1)+γ/2 (V^' )^2 ) A^' V^'+p^' A^' V^' )/(∂x^' )=0`

 

The above is the generic form of non-dimensional conservation form of the continuity, momentum and energy equations for quasi-one-dimensional flow. Let us define the above equations in a simple generic form, which will be used in calculations instead of the above rigid and complex form. Elements of the solutions vector U, the flux vector F, and the source term J as follows,

`U_1=ρ^' A^'`

`U_2=ρ^' A^' V^'`

`U_3=ρ^' (T^'/(γ-1)+γ/2 (V^' )^2 ) A^'`

`F_1=ρ^' A^' V^'`

`F_2=ρ^' A^('(V^' )^2 )+1/γ p^' A^'`

`F_3=ρ^' (T^'/(γ-1)+γ/2 (V^' )^2 ) A^' V^'+p^' A^' V^'`

`J_2=1/γ p'(∂A^')/(∂x')`

Before setting up numerical solution for conservative form, we must keep in mind that in the conservation form of the governing equations, the dependent variables (the variables for which we directly obtain numbers) are not primitive variables. In order to calculate the primitive variables, we need to decode the dependent variables and update our primitive variables in further steps.

 

 

Non Conservative form :

      The non-conservative form has its differentiated variables as coefficients of some derivatives in the equation.

Continuity equation :

`∂ρ'/(∂t')=ρ'∂V'/(∂x')-ρ'V'(∂(lnA))/(∂x')-V'∂ρ'/(∂x')`

 

Momentum equation :

`∂V'/(∂t')=V'(∂V')/(∂x')-1/γ((∂T')/(∂x')+T'ρ'(∂ρ')/(∂x'))`

 

Energy equation :

`∂T'/(∂t')=V'(∂T')/(∂x')-(γ-1)T'[(∂V')/(∂x')+V'(∂(lnA'))/(∂x')]`

 

 

MacCormack’s technique:

MacCormack’s technique is an explicit finite difference technique that is second-order accurate. It is a variant of the Lax-Wendroff method but a lot simpler in the application.

MacCormack’s is a two-step predictor-corrector sequence with forward differences for the spatial derivatives in the predictor step and the backward differences for the spatial derivatives in the corrector step. The exact method is to first obtain the forward differences for the spatial derivatives in the predictor step for all the grid points in a particular time step and then to obtain the predicted values for all the primitive variables from the first two terms of a Taylor series. Then, the backward difference of the previously obtained predicted values is obtained in the corrector step.

Once the predictor-corrector steps are obtained, the final values for all the primitive variables for the equations being solved are obtained by taking the average of the values obtained from both the predictor and the corrector steps.

We are calculating both the predictor and the corrector values and then taking the average of the two so that the final solution obtained is second-order accurate. Otherwise, with either just the predictor step or with just the corrector step, the solution would just be only first-order accurate.

 

Steps involved in MacCormack's method

• Predictor step (Replace the spatial derivatives with on the right hand side with forward differencing)
• Solution update
• Corrector step (Replace the spatial derivatives with on the right hand side with backward differencing)
• Determining the average of the time derivative
• Final solution update 

 

Predictor Step

In predictor step we replace the spatial derivatives with on the right hand side with forward differencing. In order to get clear idea , if we consider continuity equation 

`((∂ρ)/(∂t))_(i,j)^t=-(ρ_(i,j)^t.(u_(i+1,j)^t-u_(i,j)^t)/dx+u_(i,j)^t.(ρ_(i+1,j)^t-ρ_(i,j)^t)/dx)`

Here we can see that on right hand side we replace spatial derivatives with forward differencing. 

 

Corrector step

In corrector step , we first obtain a predicted value of the time derivative at time t+Δt by substituting the predictive values of rho, u and v into the right side of the  equation . Then , we replace the spatial derivatives with backward differencing.

`((∂ρ)/(∂t))_(i,j)^(t+Δt)=-(ρ_(i,j)^(t+Δt).(u_(i,j)^(t+Δt)-u_(i-1,j)^(t+Δt))/dx+u_(i,j)^(t+Δt).(ρ_(i,j)^(t+Δt)-ρ_(i-1,j)^(t+Δt))/dx)`

The average value is obtained from the average of predicted and corrector step. 

`((∂ρ)/(∂t))_(avg)=0.5[((∂ρ)/(∂t))_p+((∂ρ)/(∂t))_c ]`

Here p stands for predictor step and c for corrector step in MacCormack's method.

 

 

BOUNDARY CONDITIONS

 

• In order to examine the boundary conditions, we consider two infinitely weak Mach waves at the boundary points, both propagating  with speed of sound a . One propagating towards the left direction and other moving towards the right direction, both  named as Left Running characteristic and right running characteristic.
• 

 

Subsonic Inflow Boundary Condition -

From the above figure we can see that at the Subsonic Inlet boundary, left running characteristic would be going towards left and the flow coming through inlet which is moving towards right, and will have the subsonic speed. The relative velocity of the left characteristic will be `-(a_1- V_1)` and relative velocity for the right running characteristic will be `(a_1+V_1)`. Hence, we can see here that the left running characteristic running characteristic is propagating out of the nozzle region i.e. outside the domain of influence and right running characteristic propagates inside the domain.

Therefore, as a conclusion, at the subsonic inflow boundary, we most stipulate the value of Two Fixed Flow-Field Variables that are independent of time, whereas the value of One other variable must be allowed to Float.

So, for the subsonic inflow boundary, we choose that the Density and the Energy Terms to be stipulated. And the Velocity term to be allowed to float with respect to time. The Velocity at the inlet boundary can be estimated by Linear Extrapolation, from the velocities obtained at inner nodes at every time step.

 

Supersonic Outflow Boundary Condition -

Now, by applying the same concept at the outflow boundary. The fluid flow is flowing towards right direction same as above, however now the speed of the stream is supersonic. So, the relative velocity of left characteristic `(V_N -a_N)`  and right characteristic (`(V_N +a_N)` , both dominated dominated by the speed of the fluid flow. Hence, both of the characteristics are propagating towards right direction, and i.e. out of the Domain.

Therefore, at the Supersonic Outflow Boundary, since both the characteristic and the stream is going of the domain, So there are no Flow-Field Variables which require their values to be stipulates at the outflow boundary, Hence all variables must be allowed to Float at this boundary.

Here, at the outflow boundary all variables must be allowed to float, therefore we will calculate the variables using the Linear Extrapolation at every time Step. 

 

Time step calculation 

     The governing system of equations is hyperbolic with respect to time. So there is stability constraint on this system which is known as Courant-Friedrich-Lewy (CFL) condition. The Courant number is a simple stability analysis of a linear hyperbolic equations that gives result CFL≤1" for an explicit numerical solution to be stable.

     CFL condition can be given as,

     CFL condition = (any quantity)*(time-step with which information would travel)/(element step size)

     So, in our case, the quantity which is going to travel is velocity, or in other sense velocity is the quantity which will be travelling with certain time from one node to other. So here velocity of the fluid at a point in the flow is coupled with local speed of sound. The quantity can be expressed to be as a+v

`(a+v^' ).(Δt)/(Δx)`

     We will assume the Courant number as C = 0.5, accordingly will calculate the Δt. In addition to it, we will calculate  `(Δt)_i^t`

     at all the grid points, i = 1 to i = n, and then we'll choose the minimum value of Δt in our calculations.

     After considering all these governing equations and their respective pre-requisites to solve the problem, then we will move towards to the MATLAB code for execution.

 

Boundary Conditions

Initial Conditions

 

MATLAB CODING

 

Main Code

clear all
close all
clc
 
%Inputs
n = 31;
x = linspace(0,3,n);
dx = x(2)-x(1);
%Courant number
C = 0.5;
nt = 1400;
gamma = 1.4;
 
 
tic;
[rho_NC,t_NC,v_NC,a_NC,pressure_NC,rho_throat_NC,temp_throat_NC, pressure_throat_NC,mass_flow_NC,mass_flow_throat_NC,mass_flow_time_NC] = non_conservative_form(n,x,dx,gamma,nt,C);
Time_taken = toc;
fprintf('Time elapsed for non conservative form = %4.2f secondn', Time_taken);
 
tic;
[rho_C,t_C,p_C,v_C,a_C,rho_throat_C,temp_throat_C, pressure_throat_C,mass_flow_C,mass_flow_throat_C,mass_flow_time_C] = conservative_form(n,x,dx,gamma,nt,C);
Time_taken_conservative = toc;
fprintf('Time elapsed for conservative form = %4.2f secondn', Time_taken_conservative);
 
 
% Mass flow rate comparison between conservative and Non-conservative forms
figure(7)
plot(x,mass_flow_NC,'linewidth',1.5)
hold on
plot(x,mass_flow_C,'linewidth',1.5)
title('Mass flow rate comparison','FontSize',12)
xlabel('Distance along Nozzle (Non-Dimensional)');
ylabel('Mass flow rate-(Non-Dimensional)');
grid on
legend('Non-conservative','conservative')

Code Explanation

• We are going to solve this 1 D quasi nozzle flow problem in both forms Non-conservations and conservations using Maccormick's technique.
• First of all as usual, clear all , close all and clc has used to clear the workspace variables, figure window and command window. Then , given the necessary inputs for this problem such as n, x, nt, gamma, etc.. 
• Next line , we have called the function to solve this problem with Non-conservation and Conservation governing equations. We have written two different function code for the same.
• Also to calculate the amount of time taken to execute the solution, tic and toc has used.  

 

Non Conservative Function Code

function [rho_NC,t_NC,v_NC,a_NC,pressure_NC,rho_throat_NC,temp_throat_NC, pressure_throat_NC,mass_flow_NC,mass_flow_throat_NC,mass_flow_time_NC] = non_conservative_form(n,x,dx,gamma,nt,C)
  
%Calculate initial profiles
rho_NC = 1-0.31468*x;
t_NC = 1-0.2314*x; %t for temperature
v_NC = (0.1+1.09*x).*t_NC.^0.5;
 
%Area
a_NC = 1+2.2*(x-1.5).^2;
throat_loc_NC = find(a_NC==1);
 
 %Outer time loop
for k = 1:nt
 
    dt_NC = min(C.*(dx./(sqrt(t_NC)+(v_NC))));
 
 
    rho_old_NC = rho_NC;
    v_old_NC = v_NC;
    t_old_NC = t_NC;
 
    %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_NC(j+1))-log(a_NC(j)))/dx;
 
        %Continuity equation
        drhodt_p(j) = -rho_NC(j)*dvdx-rho_NC(j)*v_NC(j)*dlogadx-v_NC(j)*drhodx;
 
        %Momentum equation
        dvdt_p(j) = -v_NC(j)*dvdx - (1/gamma)*(dtdx +(t_NC(j)/rho_NC(j))*drhodx);
 
        %Energy equation
        dtdt_p(j) = -v_NC(j)*dtdx-(gamma-1)*t_NC(j)*(dvdx+v_NC(j)*dlogadx);
 
        %Solution update
        rho_NC(j) = rho_NC(j)+drhodt_p(j)*dt_NC;
        v_NC(j) = v_NC(j)+dvdt_p(j)*dt_NC;
        t_NC(j) = t_NC(j)+dtdt_p(j)*dt_NC;
 
    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_NC(j))-log(a_NC(j-1)))/dx;
 
        %Continuity equation
        drhodt_c(j) = -rho_NC(j)*dvdx-rho_NC(j)*v_NC(j)*dlogadx-v_NC(j)*drhodx;
 
        %Momentum equation
        dvdt_c(j) = -v_NC(j)*dvdx - (1/gamma)*(dtdx +(t_NC(j)/rho_NC(j))*drhodx);
 
        %Energy equation
        dtdt_c(j) = -v_NC(j)*dtdx-(gamma-1)*t_NC(j)*(dvdx+v_NC(j)*dlogadx);
 
     end
 
    %Compute average time derivative to obtain the final solution
    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_NC(i) = rho_old_NC(i) + drhodt(i)*dt_NC;
        v_NC(i) = v_old_NC(i) + dvdt(i)*dt_NC;
        t_NC(i) = t_old_NC(i) + dtdt(i)*dt_NC;
    end
 
    %Apply boundary condition manually
    %Inlet boundary condition
    v_NC(1) = 2*v_NC(2) - v_NC(3);
 
    %Outlet boundary condition
    rho_NC(n) = 2*rho_NC(n-1)-rho_NC(n-2);
    v_NC(n) = 2*v_NC(n-1)-v_NC(n-2);
    t_NC(n) = 2*t_NC(n-1)-t_NC(n-2);
 
    pressure_NC = rho_NC.*t_NC;
    mass_flow_NC=rho_NC.*a_NC.*v_NC;
 
    %Primitive variables at location throat
 
    pressure_throat_NC(k)=pressure_NC(throat_loc_NC);
    rho_throat_NC(k)=rho_NC(throat_loc_NC);
    temp_throat_NC(k)=t_NC(throat_loc_NC);
    Velocity_throat_NC(k)=v_NC(throat_loc_NC);
    mass_flow_throat_NC(k) = mass_flow_NC(throat_loc_NC);
 
    mass_flow_time_NC(k,:)  = mass_flow_NC;
    pressure_time_NC(k,:) = pressure_NC;
  
end
 
 %Plotting of the primitive variables along the domain
figure (1)
subplot (4,1,1)
plot(x,v_NC./sqrt(t_NC),'LineWidth',1.5,'Color','r')
xlabel('x')
ylabel('Mach')
legend('Mach number')
title('Primitive variables along the length of nozzle -Non-Conservative Form','FontSize',12)
 
subplot (4,1,2)
plot(x,pressure_NC,'LineWidth',1.5,'Color','k')
xlabel('x')
ylabel('Pressure')
legend('Pressure')
 
subplot (4,1,3)
plot(x,rho_NC,'LineWidth',1.5,'Color','g')
xlabel('x')
ylabel('Density')
legend('Density')
 
subplot (4,1,4)
plot(x,t_NC,'LineWidth',1.5,'Color','b')
xlabel('x')
ylabel('Temperature')
legend('Temperature')
 
%Variation of primitive variables at throat location with time steps
figure(2)
subplot(4,1,1)
plot(1:nt,rho_throat_NC,'LineWidth',1.5,'Color','r')
xlabel('Time step')
ylabel('Density throat')
legend('Density at throat')
title('Time-step Variation of primitive variables at throat area- Non conservative form','FontSize',12)
 
subplot(4,1,2)
plot(1:nt,temp_throat_NC,'LineWidth',1.5,'Color','k')
xlabel('Time step')
ylabel('Temp throat')
legend('Temperature at throat')
 
subplot(4,1,3)
plot(1:nt,pressure_throat_NC,'LineWidth',1.5,'Color','g')
xlabel('Time step')
ylabel('Pressure throat')
legend('Pressure at throat')
 
subplot(4,1,4)
plot(1:nt,Velocity_throat_NC./sqrt(temp_throat_NC),'LineWidth',1.5,'Color','b')
xlabel('Time step')
ylabel('Mach')
legend('Mach number at throat')
 
%Plotting Mass flow rate at different time steps along the domain
figure(3)
plot(x,mass_flow_time_NC(50,:),'LineWidth',1.5,'Color','g')
xlabel('Legth of the domain')
ylabel('Mass flow rate')
title('Mass flow rate at different time steps- Non-conservative Form','FontSize',12)
hold on
plot(x,mass_flow_time_NC(100,:),'LineWidth',1.5,'Color','r')
plot(x,mass_flow_time_NC(150,:),'LineWidth',1.5,'Color','k')
plot(x,mass_flow_time_NC(200,:),'LineWidth',1.5,'Color','m')
plot(x,mass_flow_time_NC(700,:),'LineWidth',1.5,'Color','b')
grid on
legend('50','100','150','200','700')
 
 
end

Code explanation

• Calculated the initial profile and area of the domain. And then we have find the location of the throat area where Area = 1. 
• Outer time loop has started and it will execute till number of time steps whatever we have given. Within this time loop , using Maccormick's method we have calculated the governing equations for predictor and then for corrector method from node number 2 to (n-1) because we know the boundary conditions. 
• Within this time loop calculation of the  value of time step size has done.
• After the corrector method, we have taken the average of time derivative of Predictor and Corrector step in order to obtain the final solution. 
• Next, we have updated the final solution for the last time and assigned the inlet and outlet boundary conditions.
• After that, calculated the flow filed variables throughout the nozzle and also at the throat location. 
• And then outer time loop comes to an end and from the next line plotting has been started.

 

Conservative Function Code

function [rho_C,t_C,p_C,v_C,a_C,rho_throat_C,temp_throat_C, pressure_throat_C,mass_flow_C,mass_flow_throat_C,mass_flow_time_C] = conservative_form(n,x,dx,gamma,nt,C)
 
 
%Calculate initial profiles
for i = 1:n
 
    if (x(i) >= 0 && x(i) <= 0.5)
        rho_C(i) = 1;
        t_C(i) = 1;
    elseif (x(i) > 0.5 && x(i) <= 1.5)
        rho_C(i) = 1 - (0.366*(x(i) - 0.5));
        t_C(i) = 1 - (0.167*(x(i) - 0.5));
    elseif (x(i) > 1.5 && x(i) <= 3.5)
        rho_C(i) = 0.634 - (0.3879*(x(i) - 1.5));
        t_C(i) = 0.833 - (0.3507*(x(i) - 1.5));
    end
 
end
 
 %Area & Throat profiles
a_C = 1 + (2.2*(x - 1.5).^2);
throat_loc_C = find(a_C==1);
 
%Initial profile fot Velocity and Pressure
v_C = 0.59./(rho_C.*a_C);
p_C = rho_C.*t_C;
 
%Calculate the Initial Solution Vectors(U)
U1 = rho_C.*a_C;
U2 = rho_C.*a_C.*v_C;
U3 = rho_C.*a_C.*(t_C./(gamma - 1)) + ((gamma/2).*v_C.^2);
 
%-------Macromack Method----%
% ----Conservation Equation------------%
 
%Outer time loop
for k = 1:nt
 
    dt_C = min(C.*dx./(sqrt(t_C)+v_C)); %Calculating the minimum time steps size value
 
    U1_old = U1;
    U2_old = U2;
    U3_old = U3;
 
    %Compute of flux vectors F1,F2,F3
    F1 = U2;
    F2 = ((U2.^2)./U1)+((gamma-1)/gamma)*(U3-(gamma/2)*((U2.^2)./U1));
    F3 = ((gamma*U2.*U3)./U1)- (gamma*(gamma-1)*((U2.^3)./(2*U1.^2)));
 
    %Predictor step
    for j = 2:n-1
        %Source terms in momentum equation
        j2(j) = (1/gamma)*rho_C(j)*t_C(j)*((a_C(j+1)-a_C(j))/dx);
 
        % Continuity, Momentum & Energy Equations
        dU1_dt_p(j) = -((F1(j+1)-F1(j))/dx);
        dU2_dt_p(j) = -((F2(j+1)-F2(j))/dx)+j2(j);
        dU3_dt_p(j) = -((F3(j+1)-F3(j))/dx);
 
        %Solution Updating
        U1(j) = U1(j)+dU1_dt_p(j)*dt_C;
        U2(j) = U2(j)+dU2_dt_p(j)*dt_C;
        U3(j) = U3(j)+dU3_dt_p(j)*dt_C;
    end
 
    %Calculation of primitive variables rho,V,T
    rho_C = U1./a_C;
    v_C = U2./U1;
    t_C = (gamma-1)*((U3./U1)-(0.7*v_C.^2));
    %p_C = rho_C.*t_C;
 
    %Updating flux vectors
    F1 = U2;
    F2 = ((U2.^2)./U1)+((gamma-1)/gamma)*(U3-((gamma*0.5*(U2.^2))./(U1)));
    F3 = ((gamma*U2.*U3)./U1)-((gamma*0.5*(gamma-1)*(U2.^3))./(U1.^2));
 
    %Corrector step
    for j = 2:n-1
        %Source terms in momentum equation
        j2_step(j) = (1/gamma)*rho_C(j)*t_C(j)*((a_C(j)-a_C(j-1))/dx);
 
        % Continuity, Momentum & Energy Equations
        dU1_dt_c(j) = -((F1(j)-F1(j-1))/dx);
        dU2_dt_c(j) = -((F2(j)-F2(j-1))/dx)+j2_step(j);
        dU3_dt_c(j) = -((F3(j)-F3(j-1))/dx);
    end
 
    %Compute average of predictor and corrector step
    dU1_dt_avg = 0.5*(dU1_dt_p + dU1_dt_c);
    dU2_dt_avg = 0.5*(dU2_dt_p + dU2_dt_c);
    dU3_dt_avg = 0.5*(dU3_dt_p + dU3_dt_c);
 
    %Final solution update
    for m = 2:n-1
        U1(m) = U1_old(m) + dU1_dt_avg(m)*dt_C;
        U2(m) = U2_old(m) + dU2_dt_avg(m)*dt_C;
        U3(m) = U3_old(m) + dU3_dt_avg(m)*dt_C;
    end
 
    %Applying Boundary conditions
    %Inlet boundary conditions(Node = 1)
    U1(1) = rho_C(1)*a_C(1);
    U2(1) = 2*U2(2)-U2(3);
    v_C(1) = U2(1)./U1(1);
    U3(1) = U1(1)*((t_C(1)/(gamma-1))+((gamma/2)*(v_C(1)).^2));
 
    %Outlet boundary condition(Node = n)
    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 update of primitive values
    rho_C = U1./a_C;
    v_C = U2./U1;
    t_C = (gamma-1)*((U3./U1-(gamma/2)*(v_C).^2));
    %Mass flow rate calculation
    mass_flow_C = rho_C.*a_C.*v_C;
    %p_C = rho_C.*t_C;
 
    %Primitive variables at location throat
    pressure_throat_C(k) = p_C(throat_loc_C);
    rho_throat_C(k) = rho_C(throat_loc_C);
    temp_throat_C(k) = t_C(throat_loc_C);
    Velocity_throat_C(k)=v_C(throat_loc_C);
    mass_flow_throat_C(k) = mass_flow_C(throat_loc_C);
 
    mass_flow_time_C(k,:) = mass_flow_C; %Storing mass flow rate at every time steps
 
end  %Time loop end
 
%Plotting of the primitive variables along the nozzle
figure (4)
subplot (4,1,1)
plot(x,v_C./sqrt(t_C),'LineWidth',1.5,'Color','r')
xlabel('x')
ylabel('Mach')
legend('Mach number')
title('Primitive variables along length of nozzle -Conservative Form','FontSize',12)
 
subplot (4,1,2)
plot(x,p_C,'LineWidth',1.5,'Color','k')
xlabel('x')
ylabel('Pressure')
legend('Pressure')
 
subplot (4,1,3)
plot(x,rho_C,'LineWidth',1.5,'Color','g')
xlabel('x')
ylabel('Density')
legend('Density')
 
subplot (4,1,4)
plot(x,t_C,'LineWidth',1.5,'Color','b')
xlabel('x')
ylabel('Temperature')
legend('Temperature')
 
%Variation of primitive variables at throat location with time steps
figure(5)
subplot(4,1,1)
plot(1:nt,rho_throat_C,'LineWidth',1.5,'Color','r')
xlabel('Time step')
ylabel('Density throat')
legend('Density at throat')
title('Time-step Variation of primitive variables at throat area-Conservative_form','FontSize',12)
 
subplot(4,1,2)
plot(1:nt,temp_throat_C,'LineWidth',1.5,'Color','k')
xlabel('Time step')
ylabel('Temp throat')
legend('Temperature at throat')
 
subplot(4,1,3)
plot(1:nt,pressure_throat_C,'LineWidth',1.5,'Color','g')
xlabel('Time step')
ylabel('Pressure throat')
legend('Pressure at throat')
 
subplot(4,1,4)
plot(1:nt,Velocity_throat_C./sqrt(temp_throat_C),'LineWidth',1.5,'Color','b')
xlabel('Time step')
ylabel('Mach')
legend('Mach number at throat')
 
%Plotting Mass flow rate at different time steps along the domain
figure(6)
plot(x,mass_flow_time_C(50,:),'Color','g')
xlabel('Length of the domain')
ylabel('Mass flow rate')
title('Mass flow rate at different time steps-Conservative Form','FontSize',12)
hold on
plot(x,mass_flow_time_C(100,:),'LineWidth',1.5,'Color','r')
plot(x,mass_flow_time_C(150,:),'LineWidth',1.5,'Color','k')
plot(x,mass_flow_time_C(200,:),'LineWidth',1.5,'Color','m')
plot(x,mass_flow_time_C(700,:),'LineWidth',1.5,'Color','b')
grid on
legend('50','100','150','200','700')
 
 end

 

Analytical Solution Coding

%analytical solution for 1D supersonic nozzle flow
clear all
close all
clc
 
%number of grid points
n=31;
x=linspace(0,3,n);
 
%constant
gamma=1.4;
 
%initialize mach number
m=0;
 
%area variation
a=1+(2.2*(x-1.5).^2);
b=a*-1;
 
%Finding mach number
%inlet upto throat
for i=1:16
    residual=1;%difference between LHS and RHS assume one
    while residual>1e-9%condition till 1e-9 difference(residual)
        m=m+1e-4;
        residual=((1/m^2)*((2/(gamma+1))*(1+((gamma-1)/2)*m^2))^((gamma+1)/(gamma-1)))-a(i)^2;%RHS - LHS
    end
    mach_no(i)=m;
end
%throat to outlet
for i=17:31
    residual=1;%difference between LHS and RHS assume one
    while residual>1e-9
        m=m+1e-4;
        residual=((1/m^2)*((2/(gamma+1))*(1+((gamma-1)/2)*m^2))^((gamma+1)/(gamma-1)))-a(i)^2;%RHS - LHS
        residual=residual*-1;%since increasing area,Rhs higher
    end
    mach_no(i)=m;
end
 
%Finding temperature,pressure,density,velocity analytical formula
temp=(1+((gamma-1)/2).*mach_no.^2).^(-1);
rho=(1+((gamma-1)/2).*mach_no.^2).^(-1/(gamma-1));
pressure=rho.*temp;
v=mach_no.*(temp.^0.5);
mass=rho.*a.*v;
 
%plotting
%nozzle outline
figure(1)
%subplot(5,1,1);
plot(x,a);
hold on;
plot(x,b);
plot([x(16) x(16)],[-5 5],'--','color','k');
 
 
figure(2)
%mach number plot
subplot(5,1,1);
plot(x,mach_no);
hold on;
% plot([0 x(16)],[mach_no(16) mach_no(16)],'--','color','k');
% plot([x(16) x(16)],[0 mach_no(16)],'--','color','k');
% plot([0 x(31)],[mach_no(31) mach_no(31)],'--','color','k');
% yticks([mach_no(1) mach_no(16) mach_no(31)]);
% yticklabels({mach_no(1),mach_no(16),mach_no(31)});
ylabel('Mach');
 
%pressure plot
subplot(5,1,2);
plot(x,pressure);
hold on;
% plot([0 x(16)],[pressure(16) pressure(16)],'--','color','k');
% plot([x(16) x(16)],[0 pressure(16)],'--','color','k');
% plot([0 x(31)],[pressure(31) pressure(31)],'--','color','k');
% yticks([pressure(31) pressure(16) pressure(1)]);
% yticklabels({pressure(31),pressure(16),pressure(1)});
ylabel('Pr');
 
 
%temperature plot
subplot(5,1,3);
plot(x,temp);
hold on;
% plot([0 x(16)],[temp(16) temp(16)],'--','color','k');
% plot([x(16) x(16)],[0 temp(16)],'--','color','k');
% plot([0 x(31)],[temp(31) temp(31)],'--','color','k');
% yticks([temp(31) temp(16) temp(1)]);
% yticklabels({temp(31),temp(16),temp(1)});
ylabel('Temp');
 
 
%density plot
subplot(5,1,4);
plot(x,rho);
hold on;
% plot([0 x(16)],[rho(16) rho(16)],'--','color','k');
% plot([x(16) x(16)],[0 rho(16)],'--','color','k');
% plot([0 x(31)],[rho(31) rho(31)],'--','color','k');
% yticks([rho(31) rho(16) rho(1)]);
% yticklabels({rho(31),rho(16),rho(1)});%'rho_e/rho_o'
ylabel('Density');
 
%velocity plot
subplot(5,1,5);
plot(x,v);
hold on;
% plot([0 x(16)],[v(16) v(16)],'--','color','k');
% plot([x(16) x(16)],[0 v(16)],'--','color','k');
% plot([0 x(31)],[v(31) v(31)],'--','color','k');
% yticks([v(1) v(16) v(31)]);
% yticklabels({v(1),v(16),v(31)});
xlabel('x');
ylabel('Vel');

 

 

RESULTS

 

1. Flow field variables along the length of Nozzle 

Non-conservative form

Conservative form 

In the above figures we can see how the variables change over the non dimensionalized distance. We can observe that at the throat section i.e., at x’ = 1.5 the Mach number seems to increase indicating the flow is changing from subsonic to supersonic over the nozzle, The pressure is compensated for an increase in velocity. The variation in flow distribution throughout the length of the nozzle is the same for both forms.

 

2. Time wise variations of flow field variables at nozzle throat (Area=1)

 

Non-conservative form

 

Conservative form 

From the above figures we can say that both conservative and non-conservative forms have deviations in the beginning of the time marching, but the conservative form had relatively less deviations compared to non-conservative flow. Over the time both have converged into a stable solution.

 

3. Variations of mass flow rate distributions through the nozzle  at different times

 

Non-conservative form

Conservative form 

 

From the above figures we can observe that the mass flowrate tends to stabilise with increase in timestep it is completely stable for both the forms, also we can observe that the Non conservative mass flow rate has more deviation than the conservative timestep this might be due to the reason that the control volume is moving in conservative.

 

4.Mass flow variations obtained from the non-conservations and conservations form of governing equations 

 

From the above graph it is observed that the mass flow rate for both the process are very much near to analytical solution. The Non Conservative mass flow rate is deviating more than the conservative mass flow rate hence in a study about mass flow rate the conservative form is much more reliable than the non-conservative form.

 

Grid Independence test

       In general, when we perform any calculations using numerical methods, we calculate our governing equations across fixed node points. These node points essentially form a grid/mesh. As the number of nodes increases, the spatial difference between the points decreases and the solution obtained is more accurate. But this in turn, increases the overall computation time. Hence, it is very important to arrive at a value for the number of nodes or basically, the mesh density for our physical domain under consideration so that we can estimate the solution to be within the acceptable range of accuracy and save our computation time by not considering mesh with higher densities. 

     The code was run again with 61 node points and the above table gives us a picture of how the solutions obtained compare with the actual analytical solution. The values are that of the primitive variables in the throat section of our domain. The analytical were obtained with reference from the textbook "Computational Fluid Dynamics by John D Anderson."

      It becomes quite obvious that a higher number of nodes lead to a solution that is more accurate. But, the solutions obtained from the mesh with 31 nodes were also reasonably accurate. Hence, one can decide if one wants a higher accuracy that would take more computational time or stop at a lower mesh density with a slight compromise on the solution accuracy. 

 

CONCLUSION

• The non-conservation form seems to produce slightly more accurate results for the primitive variables, and the conservation form seems to produce slightly more accurate results for the flux variables. The results in either case are certainly satisfactory.
• The conservation form yields a better mass flow distribution. The conservation form simply does a better job of conserving mass.
• The non-conservation form leads to smaller residuals. The amount by which the residuals decay is often used as an index of "quality" of the numerical algorithm. In this sense, the non-conservation form does a better job.
• Comparing the amount of calculational effort to achieve a solution, we note that the solution of the conservation form requires marginally more work. Most of this is due to the need to decode the primitive variables from the flux variables; such decoding is not necessary when we are solving the non-conservation form.
• Conservative form is much more accurate and reliable than non-conservative form. In order for us to achieve more accurate results with the non-conservative the timesteps must be increased.