Week 7 - Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method
Aim: Simulation of a 1D Super-sonic nozzle flow simulation using Macormack Method.
Objective:
1. Steady-state distribution of primitive variables inside the nozzle.
2. Time-wise variation of the primitive variables.
3. Variation of Mass flow rate distribution inside the nozzle at different time steps during the time-marching process.
4. Comparison of Normalized mass flow rate distributions of both forms.
Theory:
MacCormack technique:
It is an explicit finite difference technique that is second-order accurate in both space and time.
Let's consider an unsteady,2-d inviscid flow
Eulers equations in non-conservative form:
Continuity:`(δρ)/(δt)=-(ρ(δu)/(δx)+u(δρ)/(δx)+ρ(δv)/(δv)+v(δρ)/(δy))`------1
x-momentum:`(δu)/(δt)=-(u(δu)/(δx)+u(δu)/(δy)+1/ρ.(δρ)/(δx))` ----------2
y-momentum:`(δv)/(δt)=-(u(δv)/(δx)+v(δv)/(δy)+1/ρ.(δρ)/(δy))` ------------3
Energy: `(δu)/(δt)=-(u(δe)/(δx)+u(δe)/(δy)+P/ρ(δu)/(δx)+P/ρ.(δu)/(δy))` --------4
In the above equation, we assumed nobody forces and no volumetric heat addition i.e f=0 and q=0 the above equation is
hyperbolic w.r.t time we assume that the flow field at the time (t) is known and we will identify at (t+Deltat).
Here we are choosing density rho as our dependent flow variable for which we are going to derive the MacCormack method,
Let (ρ^t)i, j denotes the density at grid points (i,j) at the time(t) density at the same grid point at the time (t+Deltat) is given by,
`(ρ^(t+ Δt))i, j = (ρ^t)i, j + ((δρ)/(δt))av Δt`---------(5)
Where denotes (δρ/δt)av the mean value of (δρ/δt) between the times(t) and (t+Deltat)
values of (δρ/δt)av are calculated to preserve second-order accuracy without the need to calculate values of the second-order
time derivative which is the term that involves a lot of algebra.
Similar relations are written for the other flow field variables
`(u^t+ Δ t)i, j = (u^t)i, j+(δu)/(δt)av Δt`--------------(6)
`(v^t+ Δ t)i, j = (v^t)i, j+((δv)/(δt))av Δt`--------------(7)
`(e^t+ Δ t)i, j = (e^t)i, j+((δe)/(δt))av Δt`--------------(8)
From equation (5),
The average time derivative (δρ/δt )av is obtained from a predictor and corrector method.
Predictor step:
In the continuity equation (1) discretize the spatial derivatives on the right-hand side using a forward difference scheme.
The predicted value of density (rho)t+DEltat from the first two terns of the Taylor series as follows,
`(ρ^t+ Δ t)i, j = (ρ^t)i, j+((δρ)/(δt))av Δt`--------------(9)
similarly,
`(u^t+ Δ t)i, j = (u^t)i, j+(δu)/(δt)av Δt`--------------(10)
`(v^t+ Δ t)i, j = (V^t)i, j+((δv)/(δt))av Δt`--------------(11)
`(e^t+ Δ t)i, j = (e^t)i, j+((δe)/(δt))av Δt`--------------(12)
In the above equations, 9 to 12 numbers for the time derivative on the right-hand side are obtained from the equation.
Corrector step:
In the corrector step, we use the calculated predicted value of the time derivative at the time (t + Δt)((δρ/δt)^(t− Δ t))(i, j)δρ
by substituting the values of rho,u,v into the RHS of the continuity equation(1) discretizing spatial derivatives by using backward
difference scheme
-------(13)
The average value of the time derivative of density in equation (5) is obtained from the arithmetic mean of (δρ/δt)^t))i,j
obtained in equation (9) and (delta rho)/(deltat)t+Deltat at (i,j) obtained from equation 13
this allows us to obtain the final corrected value of density at the time (t+deltat) from equation 5
(ρ^(t+ Δt))i, j = (ρ^t)i, j + (δρ/δt )av Δt
A similar procedure can be applied to calculate momentum and energy
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 the cross-sectional area of the reservoir is large (A.infinity)and hence velocity is very small (V>0).
We consider the steady, isentropic flow through a convergent-divergent nozzle as sketched in The flow at the inlet to the nozzle
comes from a reservoir where the pressure and temperature are denoted by p0, and To, respectively. The Thcross-sectional area of
the reservoir is large ( theoretically, A >infinity ), and hence the velocity is very small ( V> O). Thus, Po and T0 are the stagnation
values,_ or total
pressure and total temperature, respectively. The flow expands 1sentrop1cally. to supersonic speeds at the nozzle exit, where
the exit pressure, temperature, velocity, and Mach number are denoted by Pe, Te, V., and Me, respectively. The flow IS locally
subsonic in the convergent section if the nozzle, some at the throat (minimum area), and supersonic at the d1vergent section.
flow
(M = 1) at the throat means that the local velocity at this location 1s 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. The area of the sonic throat is denoted by A*. We assume that at a given
section, where the cross-sectional area is A, the flow properties are uniform across that section. Hence, although the area of the
nozzle changes as a function of distance along with the nozzle, x, and therefore in reality the flow field is two-dimensional (the flow
varies in. the two-dimensional .ry space), we make the assumption that the flow properties-vary only with x; this is tantamount
to assuming uniform flow properties across any given cross-section. Such flow is defined as quasi-one-dimensional flow. The
governing continuity, momentum, and energy equations for this quasi-one-dimensional, steady, isentropic flow can be
expressed, respectively, as
Continuity:
Conservative form:
`δ/(δt)(ρA) +δ/(δx)(ρAV )`=0---------------(15)
Non conservative form:
`δ/(δt)(ρA) +δ/(δx)(ρAV )`+V((δ)/(δx))(ρA)=-ρA(δV)/(δx)`--------------(16)
Momentum:
`δ/(δt)(ρAV)+δ/(δx)(ρAV^2)=-A(δP)/(δx)-----------------(17)
Non conservative form:
`ρ(δV)/(δt)+ρV(δV)/(δx)=-R[ρ(δT)/(δx)-T(δρ)/(δx)]-----------(18)
Energy:
conservative:
`δ/(δt)(ρAV)+δ/(δx)(ρeAV)=-Pδ/(δx)(AV)`------------(19)
Non-conservative:
`ρ (δe)/(δt)+ρV(δe)/(δx)=-P(δV)/(δx)-PV(δV)/(δx)(lnA)---------(20)
The reason for obtaining the energy equation in the form of equation 21 is that for a calorically perfect gas it leads directly to
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 is a fundamental variable T and therefore it is convenient to deal with it as the primary dependent
variable in the energy section.
For a calorically perfect gas e= CvT
hence the energy equation of conservative and non-conservative form from equation 19 and 20 becomes
conservative f or m: `CV δ/(δt)(ρTA) +δ/(δx) (ρTV A) = − P δ/(δt)(AV )`-------------------(21)
Nonconservative: `ρCv(δT)/(δt) + ρV (δT)/(δx)= − P (δV)/(δx)− PV (δ)/(δx) (ln A)`-----------------(22)
Momentum
conservative:
`δ/(δt)(ρAV ) +δ/(δx) (ρAV^2) = − A[R(ρ (δT)/(δx) + T (δP)/(δx)` − − − − − − − − − − − (23)
Nonconservative:
`ρCv(δT)/(δt) + ρV(δT)/(δx) = − P(δV)/(δx) − PV (δ)/(δx)(ln A)`----------------------------------(24)
Momentum:
conservtaive:
`δ/(δt)(ρAV ) +δ/(δx) (ρA V^2) = − A[R(ρ(δT)/(δx) + T(δP)/(δx ))]`-----------------------------(25)
nonconservative f or m
`ρ (δV)/(δt) +ρV (δV)/(δx) ( = − [R(ρ(δT)/(δx) - T(δP)/(δx )]`-----------------------------(26)
Non-dimensional variables:
Dimensionless temperature T' =T/To
Dimensionless density rho'=rho/rho0
where To and rho0 represent reservoir temperature and density respectively.
dimensionless length x'=x/L
where L represents the length of the nozzle.
speed pf sound in the reservoir ao=squareroot of gammaRTo.
dimensionless velocity V'=V/a0
dimensionaless time t'=T/L/ao.
dimensionless area A'=A/A local area /sonic throat area.
Returning to our derived continuity momentum and energy equations and substituting non-dimensional variables.
Meshing:
The first grid point labeled point 1 is assumed to be in the reservoir points are evenly distributed along the x-axis with Deltax
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)
Average time derivatives:
Initial conditions(at the time =0)
Depending on the values of x-node we are setting up the values for rho, V, and T.
rho=1-0.3146x.
V=1-0.2314x
T=(0.1=1.09x)T^0.5.
Nozzle shape:
For area(A), we are going to set up a profile.
A=1+2.2(x-1.5)^2 0<x<3.
Mach number:
`M=V/(root T)`
Calculation of time step:
For both non-conservative and conservative forms of governing equation, we are going to take time step.
`Δt=c((Δx)/a + v)`
where C=courant number
V=local flow velocity
a=locl speed of sound.
for an explicit numerical solution to be stable C<1, here we are taking C=0.5.
Equation 40 is the CFL criterion for a one-dimensional flow along with C<1 it simply states that delta 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
Boundary conditions:
subsonic inflow boundary (at node1 )
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 complete sense as part of this
adjustment
V1=2V2-V3
Density and temperature are fixed (independent of time)
rho1=1 and T1=1
Supersonic outflow boundary(at node N)
we must allow all flow field variables to float.
Vn=2Vn-1-Vn-2
rhon=2rhon-1-rhon-2
Tn=2Tn-1-Tn-2
Conservative form of governing equations:
Governing equations for quasi-one-dimensional flow in conservative form is given by
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.
In the conservative form of governing equations the dependent variables to obtain the primitive variables(rho, V, P, T, etc...) we
must decode the elements U1, U2, U3 as follows
ρ'=(U1)/A
V'=(U2)/(U1)
`T'=e'-(γ − 1)((U3)/(U2))-γ/2(V')^2`
`P'=rho'T`
Therefore after we obtain U1, U2, and U3 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 rho', V', T' and P' from equations 57 to 60.
The pure form of flux vector is given by;
This is the pure form of governing equation in conservative form.
Initial conditions:
for 0<x'0.5: rho'=1.0 and T' =1.0
for 0.5<1.5: rho'=1.0-0.366(x'-0.5) and T'=1.0-0.167(x'-0.5)
for 1.5<x'<3.0: rho'=0.634-0.3879(x'-1.5) and T'=0.833-0.3507(x'-1.5) .
for the initial condition only,lets assume a constant mass flow through the nozzle and calculate V as,
`V'=(U2)/(ρ' A')=(0.59)/(ρ' A')`
The value of 0.59 is chosen for U2 because it is the closest to the exact analytical value of the steady state mass flow
(i.e)=0.57
Boundary conditions:
Subsonic inflow boundary (at node 1)
Two properties are held fixed and one is allowed to float rho and T are fixed at inflow boundary both equal to 1 and allowing V
to float.
