Challenge

Simulation of quasi-one-dimensional Supersonic Nozzle flow using MacCormack method in MATLAB

SShivanshu Kumar28 Jul 2020

SIMULATION OF QUASI-ONE-DIMENSIONAL SUPERSONIC NOZZLE FLOW USING MACCORMACK METHOD

PROBLEM SPECIFICATION:-

In this project, we will simulate the isentropic flow through a quasi 1D subsonic-supersonic nozzle. We will derive both the conservation and non-conservation forms of the governing equations and solve them using the MacCormack's technique. We need to 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.

Need to show the following plots inside the project report:

Steady-state distribution of primitive variables inside the nozzle.
Transient variation of the primitive variables.
Variation of Mass flow rate distribution inside the nozzle at different time-steps during the time-marching process.
Comparison of Normalized mass flow rate distribution of both forms.

Expected Results:

Perform a grid independence test on the solution and explain all the findings.
The code must be modular. We need to write separate functions to solve both forms.
Discuss all the plots we are provided with observed inferences.

-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

THEORY:

NOTE: Refer CFD book by John D. Anderson for a detailed derivation of MacCormack’s Technique.

MacCormack Technique:

It is an explicit finite-difference technique which is second-order accurate in both space and time.

Let’s consider an unsteady, 2-D, inviscid-flow.

Euler’s equations in ‘non-conservative’ form:

Continuity: `(∂ρ)/(∂t)= -(ρ (∂u)/(∂x)+u (∂ρ)/(∂x)+ρ (∂v)/(∂y)+v (∂ρ)/(∂y))` ………. (1)

x-momentum: `(∂u)/(∂t)= -(u (∂u)/(∂x)+v (∂u)/(∂y)+1/ρ (∂ρ)/(∂x))`………. (2)

y-momentum: `(∂v)/(∂t)= -(u (∂v)/(∂x)+v (∂v)/(∂y)+1/ρ (∂ρ)/(∂y))` ………. (3)

Energy: `(∂e)/(∂t)= -(u (∂e)/(∂x)+v (∂e)/(∂y)+P/ρ (∂u)/(∂x)+P/ρ (∂u)/(∂y))` ………. (4)

In the above equations, we assumed no body forces and no volumetric heat addition, i.e. `f=0` and `q ̇=0`. The above equations are hyperbolic w.r.t. time. We assume that the flow field at the time `(t)` is known and we will identify at `(t+∆t).`

Here we are choosing density, `ρ`; as our dependent flow variable for which we are going to derive the MacCormack method.

Let `ρ_(i,j)^t` denote the density at grid point `(i,j)` at the time `(t)`. Density at the same grid point at the time `(t+∆t)`, is given by:

`ρ_(i,j)^(t+∆t)=ρ_(i,j)^t+((∂ρ)/(∂t))_(av) ∆t` ………. (5)

Where `((∂ρ)/(∂t))_(av)` denotes the mean value of `((∂ρ)/(∂t))` between the times `(t) and (t+∆t).`

Values of `((∂ρ)/(∂t))_(av)` is calculated to preserve second-order accuracy without the need to calculate values of the second-order time derivative; which is the term which involves a lot of algebra; as we can notice in the Lax-Wendroff method. With MacCormack’s method, this algebra is circumvented.

Similar relations are written for the other flow field variables (momentum equations and energy equation):

`u_(i,j)^(t+∆t)=u_(i,j)^t+((∂u)/(∂t))_(av) ∆t` ………. (6)

`v_(i,j)^(t+∆t)=v_(i,j)^t+((∂v)/(∂t))_(av) ∆t` ………. (7)

`e_(i,j)^(t+∆t)=e_(i,j)^t+((∂e)/(∂t))_(av) ∆t` ………. (8)

From equation (5):

Average time-derivative, `((∂ρ)/(∂t))_(av)` is obtained from a Predictor-Corrector philosophy as follows:

Predictor Step:

In the continuity equation, equation (1), discretize the spatial derivatives on the right-hand side using forward difference scheme:

………. (9)

In equation (9), all flow variables at a time `(t)` are known values, i.e., the right-side is known.

Now, obtain a predicted value of density,`(ρ ̅ )^(t+∆t)` from the first two terms of Taylor Series, as follows:

`(ρ ̅ )_(i,j)^(t+∆t)=ρ_(i,j)^t+((∂ρ)/(∂t))_(i,j)^t ∆t` ………. (10)

Similarly,

`(u ̅ )_(i,j)^(t+∆t)=u_(i,j)^t+((∂u)/(∂t))_(i,j)^t ∆t` ………. (11)

`(v ̅ )_(i,j)^(t+∆t)=v_(i,j)^t+((∂v)/(∂t))_(i,j)^t ∆t` ………. (12)

`(e ̅ )_(i,j)^(t+∆t)=e_(i,j)^t+((∂e)/(∂t))_(i,j)^t ∆t` ………. (13)

In equations (10) to (13), numbers for the time derivatives on the right-hand side are obtained from equations, (1) to (4), respectively.

Corrector Step:

In the corrector step, we use the calculated predicted value of the time derivative at the time `(t+∆t)`, `((∂ρ)/(∂t))_(i,j)^(t+∆t)`, by substituting the values of `ρ,u,v,` into the RHS of the continuity equation (1), discretizing the spatial derivatives using the backward difference scheme:

………. (14)

The average value of time derivative of density in equation (5) is obtained from the arithmetic mean of `((∂ρ)/(∂t))_(i,j)^t` obtained in equation (9), and, obtained from equation (14):

………. (15)

This allows us to obtain the final, ‘corrected’ valued of density at the time `(t+∆t)` from equation (5):

`(ρ)_(i,j)^(t+∆t)=ρ_(i,j)^t+((∂ρ)/(∂t))_(av)^t ∆t`

A similar procedure can be applied to calculate momentum and energy.

NOTES:

Predictor-Corrector sequence described previously; yields the value of density at grid points `(i,j)` and at the time `(t+∆t).` This sequence is replaced at all grid points to obtain the density throughout the flow field at the time `(t+∆t)`.
To calculate and at the time `(t+∆t)`, the same technique is used, starting with equations (6) to (8) and utilizing the momentum and energy equations in the form of (2) to (4).
In MacCormack’s technique, the use of forward differencing on the predictor and rearward/backward differencing on the corrector is not sacrosanct; the same order of accuracy is obtained from using rearward/backward differencing on the predictor and forward differencing on the corrector. Indeed a time-marching solution can be carried out by alternating between these two sequences at every time-step if we choose.
MacCormack technique is a second-order accurate method because it has a two-step predictor-corrector sequence. It has the same accuracy as of the Lax-Wendroff method.

Numerical Solution of quasi-one-dimensional supersonic nozzle flow using MacCormack method:-

Assumptions: The flow at the inlet to the nozzle comes from a reservoir where the pressure, temperature, and density are denoted by and, respectively. The cross-sectional area of the reservoir is large (theoretically,`A→∞`), and hence the velocity is very small `(V→0)`.

The flow expands isentropically to supersonic speeds at the nozzle exit, where the outlet pressure, temperature, velocity, and Mach number are denoted by `p_e,T_e,V_e, and M_e` respectively. The flow is locally subsonic in the convergent section of the nozzle, sonic at the throat (minimum area), 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. We’re using an asterisk to denote sonic flow values, we have at the throat `V=V^*=a^*`. Similarly, the sonic flow values of pressure, temperature, and area are denoted by `p^*,T^*, and A^*`, respectively.

In reality, the real nozzle flow is a 2-D flow because, with the area changing as a function of ‘x’, in actuality, there will be flow field variations in both the ‘x’ and ‘y’ directions. This is the real physics of the flow.

The assumption of quasi-one-dimensional flow dictates that the flow properties are functions of ‘x’ only. In this problem, we can’t use Euler equations for inviscid flow because the assumption twists the real physics of the flow.

To ensure that the physical principles (Mass conservation, Newton’s second law, and Energy conservation) are satisfied. We must use the integral form of the governing equations and apply these to a control volume consistent with the quasi-one-dimensional assumption.

NOTE: For a detailed derivation of Governing Equations, refer to Chapter-7, CFD book by John D. Anderson.

Governing Equations for unsteady Quasi-one-dimensional flow:

Continuity:

Conservative form: `∂/(∂t )(ρA)+∂/(∂x )(ρAV)=0` ………. (16)

Non-conservative form: `∂/(∂t )(ρA)+V ∂/(∂x) (ρA)=-ρA (∂V)/(∂x)` ………. (17)

Momentum:

Conservative form: `∂/(∂t) (ρAV)+∂/(∂x) (ρAV^2 )=-A (∂P)/(∂x)` ………. (18)

Non-conservative form: `ρ (∂V)/(∂t)+ρV (∂V)/(∂x)=-R[ρ (∂T)/(∂x)+T (∂ρ)/(∂x)]` ………. (19)

Energy:

Conservative form: `∂/(∂t )(ρeA)+∂/(∂x) (ρeVA)=-P∂/(∂x)(AV)` ………. (20)

Non-conservative form: `ρ (∂e)/(∂t)+ρV (∂e)/(∂x)=-P (∂V)/(∂x)-PV ∂/(∂x)(ln⁡A)` ………. (21)

The reason for obtaining the energy equation in the form of equation (21) is that, for a calorically perfect gas, it leads directly to a form of the energy equation in terms of temperature ‘`T`’. For our solution of the quasi-one-dimensional nozzle flow of a calorically perfect gas, this is a fundamental variable ‘`T`’ and therefore it is convenient to deal with it as the primary dependent variable in the energy equation.

For a calorically perfect gas: `e=C_v T`

Hence, the energy equation for conservative and non-conservative form; equation (20) and equation (21) becomes:

Conservative form: `C_v ∂/(∂t) (ρTA)+∂/(∂x) (ρTVA)=-P ∂/(∂x)(AV)` ………. (22)

Non-conservative form: `ρC_v (∂T)/(∂t)+ρVC_v (∂T)/(∂x)=-P (∂V)/(∂x)-PV ∂/(∂x)(ln⁡A)` ………. (23)

When we look at the continuity, momentum, and energy equations. We can see that; there are three equations and four unknowns/variables (`ρ,V,P,T`).

The pressure (`P`) can be eliminated from these equations by using the equation of state:

`P=ρRT` ………. (24)

`(∂P)/(∂x)=R[ρ (∂T)/(∂x)+T (∂ρ)/(∂x)]` ………. (25)

Substituting equation (24) and the derivative term, equation (24); in momentum equations (18) & (19) and energy equations (22) & (23):

Momentum:

Conservative form: `∂/(∂t) (ρAV)+∂/(∂x) (ρAV^2 )=-A[R(ρ (∂T)/(∂x)+T (∂ρ)/(∂x))]` ………. (26)

Non-conservative form: `ρ (∂V)/(∂t)+ρV (∂V)/(∂x)=-R[ρ (∂T)/(∂x)+T (∂ρ)/(∂x)]` ………. (27)

Energy:

Conservative form: `C_v ∂/(∂t )(ρTA)+∂/(∂x )(ρTVA)=-ρRT ∂/(∂x)(AV)` ………. (28)

Non-conservative form: `ρC_v (∂T)/(∂t)+ρVC_v (∂T)/(∂x)=-ρRT[(∂V)/(∂x)-V ∂/(∂x )(ln⁡A )]` ………. (29)

Non-Dimensionalize:

Note that the above-written equations are in terms of dimensional variables. This is fine, and many CFD solutions are carried out directly in terms of such dimensional variables. Indeed, this has an added engineering advantage because it gives us the feeling of the real magnitudes of the real physical quantities as the solution progresses.

However, for nozzle flow, the flow field variables are frequently expressed in terms of non-dimensional variables (where the flow variables are referenced to their reservoir values). Non-dimensional variables, `P/P_0 ,T/T_0 , and ρ/ρ_0` vary between 0 and 1; which is an aesthetic advantage when presenting the result.

NOTE: For nozzle flow problems, fluid dynamicists use non-dimensional terms; whereas CFD practitioners use both dimensional and non-dimensional variables as up to their preference. There’s no difference.

Non-dimensional variables: Dimensionless variables are represented by a prime symbol (`'`).

Dimensionless Temperature: `T^'=T/T_0`

Dimensionless Density: `ρ^'=ρ/ρ_0`

Where `T_0` and `ρ_0` represents reservoir temperature and density respectively

Dimensionless Length: `x^'=x/L`

Where `L` represents the length of the nozzle

Speed of sound in the reservoir: `a_0=√(γRT_o )`

Dimensionless Velocity: `V^'=V/a_0`

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

Dimensionless Area: `A^'=A/A^*`=local area/sonic throat area

Returning to our derived continuity, momentum, and energy equations and substituting non-dimensional variables:

NOTE: These will be the final governing equations that we will solve for this particular problem.

From equations (16) and (17):

Continuity:

Conservative form: `∂/(∂t^' ) (ρ'A')+∂/(∂x^' ) (ρ^' A^' V')=0` ………. (30)

Non-conservative form: `(∂ρ^')/(∂t^' )=-P'(∂V^(' ))/(∂x^(' ) )-ρ^' V'∂/(∂x^' ) (ln⁡A^' )-V'(∂ρ^')/(∂x^' )` ………. (31)

Momentum:

Conservative form: `∂/(∂t^' ) (ρ'A'V')+∂/(∂x^' )[ρ^' A^' (V^')^2+(1/γ)P'A']=1/γ P^' (∂A^')/(∂x^' )` ………. (32)

Non-conservative form: `(∂V^')/(∂t^' )=-V^' (∂V^')/(∂x^' )-1/γ [(∂T^')/(∂x^' )+T^'/ρ^' (∂ρ^')/(∂x^' )]` ………. (33)

Energy:

Conservative form: `∂/(∂t^' )[ρ^' (e^'/(γ-1)+γ/2 (V^')^2 )A']+∂/∂x' [ρ^' (e^'/(γ-1)+γ/2 (V^')^2 ) V^' A^'+ρ^' A'V']=0` ………. (34)

Non-conservative form: `(∂T^')/(∂t^' )=-V^' (∂T^')/(∂x^' )-(γ-1) T^' [(∂V^')/(∂x^' )+V^' ∂/(∂x^' ) (ln⁡ A^' )]` ………. (35)

MESHING:

The first grid point, labeled point 1, is assumed to be in the reservoir. Points are evenly distributed along the x-axis, with `∆x` grid spacing. Point (`i`) is simply an arbitrary grid point, with points (`i-1`) and (`i+1`) as the adjacent points.

For this particular problem, we are taking the number of grid points to be 31 and 61 (as to run the grid dependence test).

NOTE: Theoretically, either form of the governing equations (conservative or non-conservative) is a suitable representation of the physical principles of mass conservation, Newton’s second law, and energy conservation. However, in CFD, there are some good numerical reasons to use one form or the other for the solution of certain flow problems. To examine the differences between results obtained from the non-conservative form of the solutions and those obtained from the conservative form. We will first cast the numerical solution of the non-conservative form of governing equations using the MacCormack method and then we’ll move to the numerical solution of a conservative form of governing equations using the MacCormack method.

NON-CONSERVATIVE FORM OF GOVERNING EQUATIONS:

NOTE: To reduce the complexity of the notation, we will drop the use of the prime (`'`) to denote a dimensionless variable. In the below equations all variables are the non-dimensional variables.

Numerical Solution:

Applying MacCormack method for the numerical solution:

Predictor-Step:

For predictor-step we’re discretizing the spatial terms using forward difference scheme:

Predictor-step for continuity equation is given by:

`((∂ρ)/(∂t))_i^t=-ρ_i^t (V_(i+1)^t-V_i^t)/(∆x)-ρ_i^t V_i^t (ln⁡A_(i+1)-ln⁡A_i )/(∆x)-v_i^t (ρ_(i+1)^t-ρ_i^t)/(∆x)` ………. (36)

Similarly, predictor-step for momentum and energy equations is given by:

`((∂V)/(∂t))_i^t=-V_i^t (V_(i+1)^t-V_i^t)/(∆x)-1/γ [(T_(i+1)^t-T_i)/(∆x)+(T_i^t)/(ρ_i^t ) (ρ_(i+1)^t-ρ_i^t)/(∆x)]` ………. (37)

`((∂T)/(∂t))_i^t=-V_i^t (T_(i+1)^t-T_i^t)/(∆x)-(γ-1) T_i^t [(V_(i+1)^t-V_i^t)/(∆x)+V_i^t (ln⁡A_(i+1)-ln⁡A_i )/(∆x)]`

………. (38)

After predicted time-derivative we update our solution (interim solution update):

Predicted values:

`(ρ ̅ )_i^(t+∆t)=ρ_i^t+((∂ρ)/(∂t))_i^t ∆t` ………. (39)

`(V ̅ )_i^(t+∆t)=V_i^t+((∂V)/(∂t))_i^t ∆t` ………. (40)

`(T ̅ )_i^(t+∆t)=T_i^t+((∂T)/(∂t))_i^t ∆t` ………. (41)

In equations (39) to (41), the variables `ρ_i^t,V_i^t, and T_i^t` are known values (initial conditions) at the time ‘`t`’.

With these predicted values, we calculate our corrector-step:

Corrector Step:

For corrector-step, we are using a backward difference scheme to discretize the spatial terms.

Corrector-step for continuity equation is given by:

………. (42)

Similarly, corrector-step for momentum and energy equations is given by:

………. (43)

………. (44)

Average time-derivatives:

………. (45)

………. (46)

And, ………. (47)

This allows us to obtain the final, ‘corrected’ valued of density, momentum, and energy at the time `(t+∆t)`:

`(ρ)_i^(t+∆t)=ρ_i^t+((∂ρ)/(∂t))_(av)^t ∆t`

`(V)_i^(t+∆t)=V_i^t+((∂V)/(∂t))_(av)^t ∆t`

`(T)_i^(t+∆t)=T_i^t+((∂T)/(∂t))_(av)^t ∆t`

Initial Conditions: (at the time `t=0`)

Depending on the values of x-node we’re setting up the values for `ρ,V, and T`:

`ρ=1-0.3146 x`

`V=1-0.2314 x`

`T=(0.1+1.09 x) T^(1/2)`

Nozzle Shape:

For area (`A`), we’re going to set up a profile:

`A=1+2.2 (x-1.5)^2` `0≤x≤3`

Mach number: `M=V/(√T)`

Calculation of Time-Step:

For both non-conservative and conservative form of governing equations we’re going to take, time-step:

`∆t=C (∆x)/(a+V)` ………. (48)

Where `C`=Courant number

`V`=local flow velocity

`a`=local speed of sound

For an explicit numerical solution to be stable: `C<=1` (the simple stability analysis of a linear hyperbolic equation). For this particular problem, we’re taking, `C=0.5`.

Equation (48), is the Courant-Friedrichs-Lewy (CFL) criterion for a one-dimensional flow, along with `C<=1`; it simply states that `∆t` must be less than, or at best equal to, the time taken by a sound wave to move from one grid point to the next.

`∆x` is the same throughout the flow, but both '`V`' and '`a`' are variables. Hence, at a given grid point at a given time step:

`∆t_i^t=C (∆x)/(a_i^t+V_i^t )` ………. (49)

And, at an adjacent grid point:

`∆t_(i+1)^t=C (∆x)/(a_(i+1)^t+V_(i+1)^t )` ………. (50)

Boundary Conditions:

Subsonic Inflow Boundary: (at node 1)

We must allow one variable to float; we choose the velocity ‘`V`’, because on a physical basis we know that the mass flow through the nozzle must be allowed to adjust to the proper steady-state, and allowing ‘`V`’ to float makes the complete sense as part of this adjustment.

`V_1=2 V_2-V_3`

Density and temperature are fixed (independent of time):

`ρ_1=1 and T_1=1`

Supersonic Outflow Boundary: (at node N)

We must allow all flow-field variables to float:

`V_n=2 V_(n-1)-V_(n-2)`

`ρ_n=2 ρ_(n-1)-ρ_(n-2)`

`T_n=2 T_(n-1)-T_(n-2)`

CONSERVATIVE FORM OF GOVERNING EQUATIONS:

Governing equations for quasi-one-dimensional flow in conservative form is given by:

Continuity:

Conservative form: `∂/(∂t^' ) (ρ'A')+∂/(∂x^' ) (ρ^' A^' V')=0` ………. (51)

Momentum:

Conservative form: `∂/(∂t^' ) (ρ'A'V')+∂/(∂x^' )[ρ^' A^' (V^')^2+(1/γ)P'A']=1/γ P^' (∂A^')/(∂x^' )` ………. (52)

Energy:

Conservative form: `∂/(∂t^' )[ρ^' (e^'/(γ-1)+γ/2 (V^')^2 )A']+∂/(∂x') [ρ^' (e^'/(γ-1)+γ/2 (V^')^2 ) V^' A^'+ρ^' A'V']=0` ………. (53)

The governing equations for quasi-one-dimensional flow can be expressed in a more generic form by defining the elements as the solution vector `U`, the flux vector `F`, and the source term `J`, as follows:

`U_1=ρ'A'`

`U_2=ρ'A'V'`

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

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

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

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

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

With these elements, governing equations (51) to (53) can be written as:

Continuity:

`(∂U_1)/(∂t^' )=-(∂F_1)/(∂x^' )` ………. (54)

Momentum:

`(∂U_2)/(∂t^' )=-(∂F_2)/(∂x^' )+J_2` ………. (55)

Energy:

`(∂U_3)/(∂t^' )=-(∂F_3)/(∂x^' )` ………. (56)

In the conservative form of governing equations, the dependent variables (the variables for which we directly obtain numbers) are not the primitive variables. To obtain the primitive variables (`ρ,V,P,T,etc`), we must decode the elements `U_1`, `U_2` and `U_3` as follows:

`ρ^'=U_1/A^'` ………. (57)

`V^'=U_2/U_1` ………. (58)

`T^'=e^'=(γ-1)(U_3/U_1 -γ/2 (V^')^2 )` ………. (59)

`P^'=ρ^' T^'` ………. (60)

Therefore, after we obtain `U_1,U_2, and U_3` at each time step from the numerical solution of equations (54) to (56), we can immediately calculate the corresponding primitive variables at each time step `ρ^',V^',T^', and P^'` from equations (57) to (60).

“Pure” form of the flux terms (`F`):

When we write the computer program with `F_1,F_2, and F_3` expressed in terms of `ρ^',V^',P^', and e^'`; some instabilities develop during the course of the time-marching solution. This instability is due to a lack of “purity” in the conservative form; which eventually causes numerical problems. When the primitive variables are used to construct `F_1,F_2, and F_3`, numerical difficulties occasionally arise. This is somehow connected to the fact that the dependent variables which appear explicitly in equations (54) to (56) are `U_1,U_2, and U_3`– not the primitive variables. For this particular reason, it is best to define `F_1,F_2, and F_3` directly in terms of the dependent variables `U_1,U_2, and U_3` and avoid the use of primitive variables in equations (54) to (56).

We’ll write the governing equations (54) to (56) as “purely” in terms of the elements of the solution vector, i.e., in terms of `U_1,U_2, and U_3` only. Hence, the “pure” form of the Flux vector is given by:

`F_1=U_2`

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

`F_3=γ((U_2 U_3)/U_1 )-γ(y-1)/2 ((U_2^3)/(U_1^2 ))`

`J_2=(γ-1)/γ (U_3-γ/2 (U_2^2)/U_1 ) (∂(ln⁡A^' ))/(∂x^' )`

This is the “pure” form of governing equations in conservative form. We’ll use this pure form in our computer program to get a stable solution and to achieve– convergence to a steady-state.

Initial Conditions: (at time `t=0`)

For `0≤x^'≤0.5 :` `ρ^'=1.0 and T^'=1.0`

For `0.5≤x^'≤1.5 :` `ρ^'=1.0-0.366(x^'-0.5) and T^'=1.0-0.167(x^'-0.5)`

For `1.5≤x^'≤3.0 :` `ρ^'=0.634-0.3879(x^'-1.5) and T^'=0.833-0.3507(x^'-1.5)`

For the initial condition only, let’s assume a constant mass flow through the nozzle and calculate `V'` as:

`V^'=U_2/(ρ^' A^' )=0.59/(ρ^' A^' )`

The value of 0.59 is chosen for `U_2` because it is the closest to the exact analytical value of the steady-state mass flow (i.e. 0.579).

Therefore, the initial condition for `V'` as a function of `x'` is obtained by substituting the variation. Finally, the initial condition for `U_1,U_2, and U_3`, are obtained by substituting the above variations for `ρ^',T', and V'` into the definitions given for solution vector (`U`), namely,

`U_1=ρ^' A'`

`U_2=ρ^' A'V'`

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

Where `e^'=T' and V'` is calculated such that `U_2 =ρ^' A^' V^'=0.59`

Boundary Conditions:

Subsonic Inflow Boundary: (at node 1)

Two properties are held fixed and one is allowed to float– `ρ' and T'` are fixed at inflow boundary, both equal to 1.0, and allowing `V'` to float:

`(U_1)_(i=1) =(ρ^' A^' )_(i=1)=A_(i=1)^'=`fixed value

`(U_2)_(i=1) =2(U_2)_(i=2) -(U_2)_(i=3)`

`(U_3)_(i=1) =ρ^' [e^'/(γ-1)+γ/2 V^'2 ]A'`

Or, `(U_3)_(i=1) =U_1 [T^'/(γ-1)+γ/2 V^'2 ]` `{∵ ρ^' A^'=U_1 and e^'=T'}`

Supersonic Outflow Boundary: (at node N)

All properties are allowed to float.

`(U_1 )_N=2(U_1 )_(N-1)-(U_1 )_(N-1)`

`(U_2 )_N=2(U_2 )_(N-1)-(U_2 )_(N-1)`

`(U_3 )_N=2(U_3 )_(N-1)-(U_3 )_(N-1)`

MATLAB PROGRAM TO SOLVE THE QUASI-ONE-DIMENSIONAL SUPERSONIC NOZZLE FLOW USING MACCORMACK METHOD

NOTE: To solve the problem in MATLAB: first we’ll create the user-defined functions for both non-conservative and conservative forms of governing equations, and then we’ll call them accordingly in our main program.

User-defined function for the Non-conservative form of governing equations:

% User-defined function for quasi-1-D supersonic nozzle flow 
% Non-conservative form of governing equations:-
 
function [rho1,T1,V1,P1,M1,m1,rho_throat1,V_throat1,T_throat1,P_throat1,M_throat1,residual1]= Quasi_1D_nozzle_flow_non_conservative_func(x,dx,A,C,nt,gamma)
 
% Initial Conditions (at time, t=0):-
% Density
rho_at_t_0=1-0.314*x;
 
% Temperature
T_at_t_0=1-0.2314*x;
 
% Velocity
V_at_t_0=(0.1+1.09*x).*T_at_t_0.^(0.5);
 
rho1=rho_at_t_0;
T1=T_at_t_0;
V1=V_at_t_0;
 
% Mass flow rate
m1=zeros(nt,length(x));
 
% Defining area at throat
throat=find(A==1);
 
% Initial condition of variables at the throat section:
rho_throat1=zeros(nt,1); %Density @ throat
T_throat1=rho_throat1; %Temperature @ throat
V_throat1=T_throat1; %Velocity @ throat
P_throat1=V_throat1; %Pressure @ throat
M_throat1=P_throat1; %Mach number @ throat
 
residual1=zeros(nt,3);
 
% Using MacCormack Method:
% Time-Loop
for k=1:nt
    
    % Creating a copy of old values
    rho_old=rho1;
    T_old=T1;
    V_old=V1;
    
    % Stability criterion (CFL condition) 
    % Time-step (dt)
    dt=min(C*dx./(T1.^0.5+V1));
    
    % Predictor-Step
    % Space-Loop
    for i=2:length(x)-1
        
        % Discretizing using forward difference scheme for spatial derivatives
        dVdx=(V1(i+1)-V1(i))/dx;
        dlogAdx=(log(A(i+1))-log(A(i)))/dx;
        drhodx=(rho1(i+1)-rho1(i))/dx;
        dTdx=(T1(i+1)-T1(i))/dx;
        
        % Continuity equation
        drho_dt_p(i)=-rho1(i)*dVdx-rho1(i)*V1(i)*dlogAdx-V1(i)*drhodx;
        
        % X-Momentum equation
        dV_dt_p(i)=-V1(i)*dVdx-(1/gamma)*(dTdx+(T1(i)/rho1(i))*drhodx);
        
        % Energy equation
        dT_dt_p(i)=-V1(i)*dTdx-(gamma-1)*T1(i)*(dVdx+V1(i)*dlogAdx);
        
        % Interim Solution update
        rho1(i)=rho1(i)+drho_dt_p(i)*dt;
        V1(i)=V1(i)+dV_dt_p(i)*dt;
        T1(i)=T1(i)+dT_dt_p(i)*dt;
    end % end of space loop (Predictor-step)
    
    % Corrector-Step
    % Space-Loop
    for i=2:length(x)-1
        
        % Discretizing using backward difference scheme for spatial derivatives
        dVdx=(V1(i)-V1(i-1))/dx;
        dlogAdx=(log(A(i))-log(A(i-1)))/dx;
        drhodx=(rho1(i)-rho1(i-1))/dx;
        dTdx=(T1(i)-T1(i-1))/dx;
        
        % Continuity equation
        drho_dt_c(i)=-rho1(i)*dVdx-rho1(i)*V1(i)*dlogAdx-V1(i)*drhodx;
        
        % X-Momentum equation
        dV_dt_c(i)=-V1(i)*dVdx-(1/gamma)*(dTdx+(T1(i)/rho1(i))*drhodx);
        
        % Energy equation
        dT_dt_c(i)=-V1(i)*dTdx-(gamma-1)*T1(i)*(dVdx+V1(i)*dlogAdx);
            
    end % end of space loop (Corrector-step)
    
    % Computing the average time derivative
    drho_dt=0.5*(drho_dt_p+drho_dt_c);
    dV_dt=0.5*(dV_dt_p+dV_dt_c);
    dT_dt=0.5*(dT_dt_p+dT_dt_c);
    
    residual1(k,1)=max(abs(drho_dt));
    residual1(k,2)=max(abs(dV_dt));
    residual1(k,3)=max(abs(dT_dt));
    
    % Final solution update
    for j=2:length(x)-1
        rho1(j)=rho_old(j)+(drho_dt(j)*dt);
        V1(j)=V_old(j)+(dV_dt(j)*dt);
        T1(j)=T_old(j)+(dT_dt(j)*dt);
    end
    
% Applying Boundary Conditions:-
% Subsonic Inflow Boundary:
V1(1)=2*V1(2)-V1(3); % Velocity at the inlet
 
% Supersonic Outflow Boundary:
V1(end)=2*V1(end-1)-V1(end-2); % Velocity at outlet
rho1(end)=2*rho1(end-1)-rho1(end-2); % Density at outlet
T1(end)=2*T1(end-1)-T1(end-2); % Temperature at outlet
 
% Mass flow rate
m1(k,:)=rho1.*A.*V1; 
 
% Pressure (P)
P1=rho1.*T1;
 
% Mach number
M1=V1./(T1.^0.5);
 
% Values at throat:-
rho_throat1(k)=rho1(throat);
M_throat1(k)=M1(throat);
P_throat1(k)=P1(throat);
T_throat1(k)=T1(throat);
V_throat1(k)=V1(throat);
end %end of time-loop
end %end of function

User-defined function for the Conservative form of governing equations:

% User-defined function for quasi-1-D supersonic nozzle flow 
% Conservative form of governing equations:-
 
function [rho2,V2,T2,P2,M2,m2,rho_throat2,V_throat2,T_throat2,M_throat2,P_throat2,residual2]= Quasi_1D_nozzle_flow_conservative_func(x,dx,A,C,gamma,nt)
 
% Initial Conditions (at time, t=0)
rho_at_t_0=zeros(1,length(x));
T_at_t_0=rho_at_t_0;
V_at_t_0=T_at_t_0;
 
for i=1:length(x)
    if x(i)>=0 && x(i)<0.5
        rho_at_t_0(i)=1;
        T_at_t_0(i)=1;
    
    elseif x(i)>=0.5 && x(i)<1.5
            rho_at_t_0(i)=1-0.366*(x(i)-0.5);
            T_at_t_0(i)=1-0.167*(x(i)-0.5);
    
    elseif x(i)>=1.5 && x(i)<=3.5
                rho_at_t_0(i)=0.634-0.3879*(x(i)-1.5);
                T_at_t_0(i)=0.833-0.3507*(x(i)-1.5);
   
    end
    V_at_t_0(i)=0.59/(rho_at_t_0(i)*A(i));
end
        
rho2=rho_at_t_0;
T2=T_at_t_0;
V2=V_at_t_0;
 
% Mass flow rate
m2=zeros(nt,length(x));
 
 
% Initial Conditions for the solution vectors (U)
U1=rho2.*A;
U2=rho2.*A.*V2;
U3=rho2.*A.*((T2./(gamma-1))+(gamma/2).*(V2.^2));
 
% Defining area at throat
throat=find(A==1);
 
% Initial condition of variables at the throat section:
rho_throat2=zeros(nt,1); %Density @ throat
T_throat2=rho_throat2; %Temperature @ throat
V_throat2=T_throat2; %Velocity @ throat
P_throat2=V_throat2; %Pressure @ throat
M_throat2=P_throat2; %Mach number @ throat
 
residual2=zeros(nt,3);
 
% Using MacCormack method
% Time-Loop 
 for k=1:nt
     
     % Creating a copy of Old Values
     U1_old=U1;
     U2_old=U2;
     U3_old=U3;
     
     % Stability criterion (CFL condition) 
     % Time-step (dt)
     dt=min(C*dx./(T2.^0.5+V2));
     
     % Initial Conditions for flux vectors (F)
     F1=U2;
     F2=(U2.^2)./U1+((gamma-1)/gamma)*(U3-0.5*gamma*(U2.^2)./U1);
     F3=gamma*U2.*U3./U1-0.5*gamma*(gamma-1)*(U2.^3)./(U1.^2);
          
     % Predictor-Step
     % Space-Loop
     for i=2:length(x)-1
         % Discretizing using forward difference scheme for spatial derivatives
         % Source Vector
         J2=(1/gamma)*rho2(i)*T2(i)*((A(i+1)-A(i))/dx);
         
         % Continuity Equation
         dU1_dt_p(i)=-(F1(i+1)-F1(i))/dx;
         
         % Momentum Equation
         dU2_dt_p(i)=-(F2(i+1)-F2(i))/dx+J2;
         
         % Energy Equation
         dU3_dt_p(i)=-(F3(i+1)-F3(i))/dx;
         
         % Interim Solution Update
         
         U1(i)=U1(i)+dU1_dt_p(i)*dt;
         U2(i)=U2(i)+dU2_dt_p(i)*dt;
         U3(i)=U3(i)+dU3_dt_p(i)*dt;
         
     end
     
     % Obtaining the Primitive variables form our solution vector (U1,U2,U3):-
     rho2=U1./A;
     V2=U2./U1;
     T2=(gamma-1)*((U3./U1)-((gamma/2)*(V2.^2)));
     
     % Updating Flux Vectors
     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
     % Space-Loop
     for i=2:length(x)-1
         % Discretizing using backward difference scheme for spatial derivatives
         % Source Vector
         J2=(1/gamma)*rho2(i)*T2(i)*((A(i)-A(i-1))/dx);
         
         % Continuity Equation
         dU1_dt_c(i)=-(F1(i)-F1(i-1))/dx;
         
         % Momentum Equation
         dU2_dt_c(i)=-(F2(i)-F2(i-1))/dx+J2;
         
         % Energy Equation
         dU3_dt_c(i)=-(F3(i)-F3(i-1))/dx;
       
     end
     
     % Computing the average time derivative
     dU1_dt=0.5*(dU1_dt_p+dU1_dt_c);
     dU2_dt=0.5*(dU2_dt_p+dU2_dt_c);
     dU3_dt=0.5*(dU3_dt_p+dU3_dt_c);
     
     residual2(k,1)=max(abs(dU1_dt));
     residual2(k,2)=max(abs(dU2_dt));
     residual2(k,3)=max(abs(dU3_dt));
     
     % Final Solution Update
     for j=2:length(x)-1
         U1(j)=U1_old(j)+dU1_dt(j)*dt;
         U2(j)=U2_old(j)+dU2_dt(j)*dt;
         U3(j)=U3_old(j)+dU3_dt(j)*dt;
     end
     
% Subsonic Inflow Boundary
U1(1)=rho2(1)*A(1);
U2(1)=2*U2(2)-U2(3);
U3(1)=U1(1)*((T2(1)/(gamma-1))+0.5*gamma*V2(1)^2);
 
% Supersonic Outflow Boundary
U1(end)=2*U1(end-1)-U1(end-2);
U2(end)=2*U2(end-1)-U2(end-2);
U3(end)=2*U3(end-1)-U3(end-2);
 
% Final values of different variables:
rho2=U1./A; %Density
V2=U2./U1; %Velocity
T2=(gamma-1)*((U3./U1)-0.5*gamma*V2.^2); %Temperature
M2=V2./(T2.^0.5); %Mach number
P2=rho2.*T2; %Pressure
m2(k,:)=rho2.*A.*V2; %Mass flow rate
 
rho_throat2(k)=rho2(throat);
V_throat2(k)=V2(throat);
T_throat2(k)=T2(throat);
P_throat2(k)=P2(throat);
M_throat2(k)=M2(throat);
 
end %end of the time loop
end %end of function

Main Program:

% Simulation of Quasi-one-dimensional Supersonic Nozzle Flow using MacCormack method
 
clear all
close all
clc
% Defining Geometry and meshing
% number of grid points
nx=31;
% length of domain
L=3;
% meshing
x=linspace(0,L,nx);
% grid spacing
dx=x(2)-x(1);
% Area of nozzle
A= 1+2.2*(x-1.5).^2;
 
gamma=1.4;
nt=1400; %time-steps
C=0.5; %Courant number
 
% Calling function for non-conservative form
tic
[rho1,T1,V1,P1,M1,m1,rho_throat1,V_throat1,T_throat1,P_throat1,M_throat1,residual1]= Quasi_1D_nozzle_flow_non_conservative_func(x,dx,A,C,nt,gamma);
Time_taken_by_non_conservative_form=toc 

% Calling function for conservative form
tic
[rho2,V2,T2,P2,M2,m2,rho_throat2,V_throat2,T_throat2,M_throat2,P_throat2,residual2]= Quasi_1D_nozzle_flow_conservative_func(x,dx,A,C,gamma,nt);
Time_taken_by_conervative_form=toc

% Plots for different variables:
 
% Plots for different flow parameters: Non-Conservative form
figure (1)
subplot (4,1,1)
plot (x,M1,'LineWidth',2)
grid on
grid minor;
ax=gca;
ax.FontSize=14;
ylim ([0,4.5])
yticks (0:1.5:4.5)
ylabel ('Mach number','FontSize',12,'FontWeight','bold')
title ('Variation of different Flow Parameters: Non-Conservative Form of Governing Equations','FontSize',25,'FontWeight','Bold','color','g')
 
subplot (4,1,2)
plot (x,P1,'LineWidth',2)
grid on
grid minor
ax=gca;
ax.FontSize=14;
yticks (0:0.2:1);
ylabel ('Pressure ratio','FontSize',12,'FontWeight','bold')
 
subplot (4,1,3)
plot (x,rho1,'LineWidth',2)
grid on
grid minor
ax=gca;
ax.FontSize=14;
yticks (0:0.2:1)
ylabel ('Density ratio','FontSize',12,'FontWeight','bold')
 
subplot (4,1,4)
plot (x,T1,'LineWidth',2)
grid on
grid minor
ax=gca;
ax.FontSize=14;
ylim ([0,1])
yticks (0:0.2:1)
ylabel ('Temperature ratio','FontSize',12,'FontWeight','bold')
xlabel ('x/L','FontSize',12,'FontWeight','bold')
 
% Plots for different flow parameters: Conservative form
figure (2)
subplot (4,1,1)
plot (x,M2,'LineWidth',2)
grid on
grid minor;
ax=gca;
ax.FontSize=14;
ylim ([0,4.5])
yticks (0:1.5:4.5)
ylabel ('Mach number','FontSize',12,'FontWeight','bold')
title ('Variation of different Flow Parameters: Conservative Form of Governing Equations','FontSize',25,'FontWeight','Bold','color','g')
 
subplot (4,1,2)
plot (x,P2,'LineWidth',2)
grid on
grid minor
ax=gca;
ax.FontSize=14;
yticks (0:0.2:1);
ylabel ('Pressure ratio','FontSize',12,'FontWeight','bold')
 
subplot (4,1,3)
plot (x,rho2,'LineWidth',2)
grid on
grid minor
ax=gca;
ax.FontSize=14;
yticks (0:0.2:1)
ylabel ('Density ratio','FontSize',12,'FontWeight','bold')
 
subplot (4,1,4)
plot (x,T2,'LineWidth',2)
grid on
grid minor
ax=gca;
ax.FontSize=14;
ylim ([0,1])
yticks (0:0.2:1)
ylabel ('Temperature ratio','FontSize',12,'FontWeight','bold')
xlabel ('x/L','FontSize',12,'FontWeight','bold')
 
% Plots for different flow parameters at 'Throat Section': Non-conservative vs. Conservative form
figure (3)
subplot (4,1,1)
plot (1:nt,M_throat1,'LineWidth',2)
hold on
plot (1:nt,M_throat2,'LineWidth',2)
grid on
grid minor
ax=gca;
ax.FontSize=12;
yticks (0:0.2:2)
ylabel ('Mach Number','FontSize',15,'FontWeight','bold')
l_g_d=legend('Non-Conservative','Conservative');
l_g_d.FontSize=15;
title ('Variation of different Flow Parameters at the ''Throat Section'' of Nozzle','FontSize',25,'FontWeight','bold','color','g')
 
subplot (4,1,2)
plot (1:nt,P_throat1,'LineWidth',2)
hold on
plot (1:nt,P_throat2,'LineWidth',2)
grid on
grid minor
ax=gca;
ax.FontSize=12;
yticks (0:0.2:1)
ylabel ('P/P_0','FontSize',15,'FontWeight','bold')
l_g_d=legend('Non-Conservative','Conservative');
l_g_d.FontSize=15;
 
subplot (4,1,3)
plot (1:nt,rho_throat1,'LineWidth',2)
hold on
plot (1:nt,rho_throat2,'LineWidth',2)
grid on
grid minor
ax=gca;
ax.FontSize=12;
yticks (0:0.2:1)
ylabel ('rho/rho_0','FontSize',15,'FontWeight','bold')
l_g_d=legend('Non-Conservative','Conservative');
l_g_d.FontSize=15;
 
subplot (4,1,4)
plot (1:nt,T_throat1,'LineWidth',2)
hold on
plot (1:nt,T_throat2,'LineWidth',2)
grid on
grid minor
ax=gca;
ax.FontSize=12;
yticks (0:0.2:1)
ylabel ('T/T_0','FontSize',15,'FontWeight','bold')
xlabel ('Time Steps','FontSize',15,'FontWeight','bold')
l_g_d=legend('Non-Conservative','Conservative');
l_g_d.FontSize=15;
 
% Plot for Residuals: Non-conservative vs. Conservative 
figure (4)
subplot (2,1,1)
plot (1:nt,residual1,'LineWidth',2)
grid on
grid minor
ax=gca;
ax.FontSize=12;
xlabel ('Time Steps','FontSize',18,'FontWeight','bold')
ylabel ('Residuals','FontSize',20,'FontWeight','bold')
yticks (0:0.2:1.2)
l_g_d=legend('|drho/dt|_{av}','|dV/dt|_{av}','|dT/dt|_{av}');
l_g_d.FontSize=15;
title ('Residuals: Non-Conservative Form','FontSize',25,'FontWeight','bold','color','g')
 
subplot (2,1,2)
plot (1:nt,residual2,'LineWidth',2)
grid on
grid minor
ax=gca;
ax.FontSize=15;
xlabel ('Time Steps','FontSize',18,'FontWeight','bold')
ylabel ('Residuals','FontSize',20,'FontWeight','bold')
yticks (0:0.5:4)
l_g_d=legend('|dU_1/dt|_{av}','|dU_2/dt|_{av}','|dU_3/dt|_{av}');
l_g_d.FontSize=14;
title ('Residuals: Conservative Form','FontSize',25,'FontWeight','bold','color','g')
 
% Plot for Mass Flow Rate variation: Non-conservative form
figure (5)
plot (x,m1(50,:),'-*','LineWidth',2)
hold on
grid on
grid minor
plot (x,m1(100,:),'-x','LineWidth',2)
hold on
plot (x,m1(150,:),'-d','LineWidth',2)
hold on
plot (x,m1(200,:),'-p','LineWidth',2)
hold on
plot (x,m1(500,:),'-o','LineWidth',2)
hold on
plot (x,m1(700,:),'-s','LineWidth',2)
hold on
plot (x,m1(1400,:),'--','LineWidth',2)
xlabel ('x/L','FontSize',15,'FontWeight','bold')
ylabel ('rhoAV/rho_0A^*a_0','FontSize',15,'FontWeight','bold')
l_g_d=legend ('50 Deltat','100 Deltat','150 Deltat','200 Deltat','500 Deltat','700 Deltat','1400 Deltat');
l_g_d.FontSize=14;
title ('Variation of Non-dimensional Mass Flow Rate: Non-conservative Form','FontSize',25,'FontWeight','bold','color','g')
 
% Plot for Mass Flow Rate variation: Conservative form
figure (6)
plot (x,m2(50,:),'-*','LineWidth',2)
hold on
grid on
grid minor
plot (x,m2(100,:),'-x','LineWidth',2)
hold on
plot (x,m2(150,:),'-d','LineWidth',2)
hold on
plot (x,m2(200,:),'-p','LineWidth',2)
hold on
plot (x,m2(500,:),'-o','LineWidth',2)
hold on
plot (x,m2(700,:),'--','LineWidth',2)
xlabel ('x/L','FontSize',15,'FontWeight','bold')
ylabel ('rhoAV/rho_0A^*a_0','FontSize',15,'FontWeight','bold')
l_g_d=legend ('50 Deltat','100 Deltat','150 Deltat','200 Deltat','500 Deltat','700 Deltat');
l_g_d.FontSize=14;
title ('Variation of Non-dimensional Mass Flow Rate: Conservative Form','FontSize',25,'FontWeight','bold','color','g')
 
% Plot for variation of Non-dimensional Mass Flow Rate-Steady State: Non-conservative vs. Conservative form vs. Exact Solution
figure (7)
plot (x,0.579*ones(length(x),1),'--','Linewidth',2)
hold on
grid on
grid minor
plot (x,m1(end,:),'LineWidth',2)
hold on
plot (x,m2(end,:),'LineWidth',2)
ax=gca;
ax.FontSize=12;
xlabel ('x/L','FontSize',15,'FontWeight','bold')
ylabel ('rhoAV/rho_0A^*a_0','FontSize',15,'FontWeight','bold')
legend ('Exact','Non-conservative','Conservative','FontSize',15)
title ('Variation of Non-dimensional Mass Flow Rate: Steady State','FontSize',25,'FontWeight','bold','color','g')

RESULTS & DISCUSSIONS

From the above plots and data, we can notice that, as the Mach number increases (flows reaches from subsonic to the supersonic regime), the primitive variables: Pressure, Density, and Temperature; shows declivity after the sonic region (Mach=1) till the exit of the nozzle. Although not compared here, these results are almost identical to the analytical solution (results validated from the CFD book by John D. Anderson). The steady-state results obtained from the time-marching solution for the governing equations in non-conservative form and conservative form are essentially the same but with a few slight, but notable differences. The converged solution is presented in the above table, which is the results after 1400 time steps. After doing a quick comparison of the results shown in the above table for non-conservative and conservative forms; we can conclude that for this particular problem both forms give the same results. The flow problem studied in the above table is the isentropic, subsonic-supersonic flow through a nozzle, and for such flow, the choice of the form of equations is not important.

From the above plot for variation of different primitive variables w.r.t time steps, we can see that the conservative form, for which we used the modified initial conditions, converges to the steady-state faster than the non-conservative form.

The behavior of residuals for the conservative form is not as good as that for the non-conservative form of the equations. For the non-conservative form, at early times the residuals are approximately on the order of `10^-1`, but they decay to about `10^-8` after 1400 time steps. In contrast, for the conservative form of the equations, the residuals at early times are on the higher-order side but decay only to `10^-5` after 1400 time steps.

The most dramatic difference between the non-conservative and conservative form of governing equations can be seen in the mass flow distribution. Note that for non-conservative form, after 500-time steps (labeled as 500`∆t`) the mass flow is moving toward the steady-state value. After 1400 time steps the flow distribution starts showing consistency and after 5000-time steps (not shown in the figure), the mass flow has converged to a horizontal line; i.e. `ρAV=`constant. Whereas, for the conservative form of governing equations, the transient mass flow rate deviates away from the initial conditions; the result after 100-time steps (labeled as 100`∆t` ) shows a somewhat “humped” distribution. After 200 time steps (labeled as 200`∆t`), the mass flow rate distribution is started becoming more constant, at 500-time steps (labeled as 500`∆t`) the flow distribution start showing consistency, and after 700-time steps (labeled as 700`∆t`), it is almost (but not quite) equal to a constant value. Moreover, it is quite close to the exact analytical value of 0.579.

On comparing the mass flow distribution for non-conservative form with the conservative form, we can see that for conservative form the variations in mass flow are much less severe. We can notice milder behavior in the case of conservative form due to our assumed initial condition of constant mass flow.

In the above plot, we compare the steady-state variations of mass flow obtained with the non-conservative and conservative forms of the governing equations (both after 1400 time steps). Here we see that the steady-state mass flow distribution predicted by the conservative form of the equations is much more satisfactory than that obtained from the non-conservative form, on two accounts:

The conservative form gives a distribution i.e. much closer to being a constant. In contrast, the non-conservative results have a sizeable variation, with some spurious oscillations at both the inflow and outflow boundary. Of course, on a practical basis, these variations are not apparent, and the mass flow essentially appears to be constant.
The steady-state mass flow results obtained with the conservative form are, on the whole, much closer to the exact analytical solution `(ρ^' A^' V^'=0.579).`

The comparison shown in the above figure illustrates a general advantage of the conservative form of the equations. The conservative form does a better job of preserving mass throughout the flow field. It is because the mass flow is one of the dependent variables in the equations–the mass flow is a primary result of these equations. On the other hand, the dependent variables in the non-conservative form of the equations are the primitive variables, and the mass flow obtained as a secondary result. Here, the conservative form of governing equations does a better job of conserving the mass throughout the flow field.

The above discussion does not necessarily establish a definite superiority of the conservative form results over the non-conservative form results. Quite the contrary, let us look at the primitive variables; in particular, density, temperature, pressure, and Mach number at nozzle throat, as tabulated in the above table. The first row gives the exact, analytical results. The second and third-row give the numerical results for the non-conservative and conservative forms, respectively. Note that the non-conservative form results are distinctly closer to the exact values. The second last and last row in the above table gives the non-conservative form and conservative form results, respectively; but with grid points as twice as many as the first two rows (61 in comparison to 31 grid points). A comparison of the last two rows with the first two rows indicates grid independence for the conservative form and non-conservative form results. Note that, by doubling the number of grid points, the steady-state numerical results are slightly closer to the exact, analytical values (but still not close as the non-conservative form results). For all the practical purposes, we have grid independence with 31 grid points (suffice to produce the steady-state results, for all practical purposes).

The above table shows the time taken by each form of the governing equation to solve this problem for two different numbers of grid points. By observing the data, we can say that the non-conservative form of governing equation takes less time to solve this nozzle flow problem (for both the number of grid points). It can be also noticed that as we increase the number of grid points, the value of flow variables approach towards the exact, analytical solution. But this is also followed by the increase in the runtime of the simulation. Since for faster results, fewer grid points are used, but the deviation of the flow variables w.r.t. the exact analytical solution in such a case is more. Hence, by taking a lesser grid point, a compromise should be made between the accuracy of the solution and the runtime of the simulation.

CONCLUSIONS

In summary, for the given flow problem, we can’t establish a logical superiority of the conservative form of the governing equations over the non-conservative form of the governing equations. From all our previous discussions, we can only make the following conclusions:

The conservative form yields a better mass flow distribution. Compared to non-conservative form, the conservative form simply does a better job of conserving mass.
The non-conservative 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-conservative form does a better job than a conservative form.
There is no logical superiority of using either form in terms of the accuracy of the results. The non-conservative form seems to produce slightly more accurate results for the primitive variables, and the conservative form seems to produce slightly more accurate results for the flux variables. The results, in either case, are satisfactory.
While comparing the amount of calculational effort to achieve a solution for non-conservative form and conservative form of governing equations, we note that the solution of the conservative form requires marginally more work. Most of this is because of the need to decode the primitive variables from the flux variables; such decoding is unnecessary when you are solving the non-conservative form.
Hence, it is clear that both the conservative and non-conservative form of governing equations return almost identical steady-state profiles. We can use either of them to solve this particular problem as we are not capturing any shocks in the present problem.

Simulation of quasi-one-dimensional Supersonic Nozzle flow using MacCormack method in MATLAB

Ship work like this
in your own course.

Talk to a career counsellor

Ship work like thisin your own course.

Talk to a career counsellor

Ship work like this
in your own course.