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

                                                                                                                Date: 14-02-2021

 

SIMULATION OF A 1D SUPER-SONIC NOZZLE FLOW SIMULATION  USING MACORMACK METHOD

 

 

INTRODUCTION: -

 

 We consider steady, isentropic flow inside the convergent - divergent nozzle (CD-nozzle). Generally, this CD-nozzle is used to rise the velocity form subsonic to supersonic. In reverse fashion also we can do like supersonic to subsonic, then this is called CD-diffuser. We have analytical solution for this CD-nozzle in x-direction by assuming other two directions are neglected. We assume that flow properties vary only with ‘x’. That means we can say that uniform flow properties across any given cross- section (Area only depends on ‘x’). That’s why this type of flow is called Quasi 1D-flow.

                                          McCormack’s method is an explicit, finite difference method to find numerical time marching solution. This is advanced method than Lax-Wendorff technique. Compare to Lax-Wendorff method, this McCormack’s method has more accurate numerical solution and simple to programme. We use predictor step and correct step to find solution for each time step.

 

 Here, we conduct Grid independence test to remove the sensitivity on number of grid points.

 

MAIN OBJECTIVE: -

• Derivation for both conservative and non-conservative forms for isentropic & quasi 1D- CD- nozzle.
• Find the numerical solution for both forms by help of McCormack’s method.
• Time wise variation of flow filed variables.
• Steady state distribution of the flow field variables for both forms inside the nozzle.
• Variation of mass flow rate inside the nozzle at different time steps for both forms.
• Grid independence test for both forms.

Non-conservative forms of the governing differential equations derivation: -

• In this form we always consider that our control volume is moving with the fluid flow.

Continuity equation: -

We start with integral form of continuity equation

 

 

                         `del/del_t(intintint rho*dv) +intint rho*v*ds =0`         ---------------1

               where
                                                           ρ - density
                                                           ѵ - control volume
                                                            s – surface area

                    

• Here term1(volume integral) indicates the rate of change of mass inside the control volume and term2(surface integral) indicates the net mass flow (outflow-inflow) .
• Now we apply our continuity equation for infinitesimal control volume(dѵ) and integrate it . In our problem we deal with one dimension so ‘dx’ will be my control volume. You can see this in figure1.
• From fig1 we can write term1 as

                                   `del/del_t (intintint rho*dv) = del/del_t (rho*A*dx)`  -----------2

• In the limit of ‘dx’ control volume become A*dx . so we can say our surface area is ‘A’

Then surface integral becomes

 

 

Fig (1): Control volume for deriving the partial differential equations for unsteady, quasi 1D-flow

 

   ∬ ρ 𝐯. 𝐝𝐬 = −ρ*v*A + ( ρ + dρ)*(v+dv) *(A+dA)  = −ρ*v*A + ρ*v*A + ρ*v*dA + ρ*dv*A + ρ*dv*dA + 𝑑ρ*v*A +𝑑ρ*v*dA + 𝑑ρ*dv*A + 𝑑ρ*dv*dA

➢ Neglect very small elements then our equation becomes

∬ ρ 𝐯. 𝐝𝐬 = ρ*v*dA + ρ*dv*A + ρ𝑑*v*A = d(ρAV) ---3

➢ Substitute equation 2 and 3 in equation1 then we end up with

ꝺ/ꝺt(ρAdx) + d(ρAv) =0
ꝺ/ꝺt(ρA) + (d(ρAv)/dx) =0

➢ We can write d(ρAv)/dx =ꝺ/ꝺx(ρAv) [By taking limits of ‘dx’ goes to ‘0’]

So,
ꝺ/ꝺt(𝛒𝐀) + ꝺ/ꝺx(𝛒𝐀𝐯) = 0 ---- 4

➢ If we expand the above partial differential equation we end up with

ꝺ/ꝺt(𝛒𝐀) +(𝛒𝐀) ∗ ꝺv/ꝺx + (𝛒𝐯) ∗ ꝺA/ꝺx + (𝐀𝐯) ∗ ꝺ𝛒/ꝺx = 0 ---I

➢ This is 1 dimensional continuity governing partial differential equation.
➢ Here we consider unsteady quasi 1D – flow.

Momentum equation: -

• This is derived from fundamental concept of newtons second law.
• Integral form of the x-component of the momentum equation for inviscid flow with no body forces

 ꝺ/ꝺt(∭(ρu) dѵ) + ∬(ρu𝐯). 𝐝𝐬 = - ∬(p. 𝐝𝐬)x ----5
    (Term1)                 (Term2)             (Term3)

➢ From figure1 we can write that  

Term1 = ꝺ/ꝺt (∭ρ udѵ) = ꝺ/ꝺt (ρvAdx) --6

Term2 = ∬(ρu𝐯). 𝐝𝐬 = - ρ`v^2` A +(ρ + dρ) `(v+dv)^2`(A+dA) --7

➢ Form figure2 we can write the resultant of surface force (pressure)

Term3 = ∬(p. 𝐝𝐬)x = -pA +(p+dp)(A+dA)-2p(dA/2) ---8

➢ Substitute equation 6,7 and 8 in equation 5. Then we end up with

ꝺ/ꝺt (ρvAdx) - ρ`v^2` A +(ρ + dρ)`(v+dv)^2`(A+dA) = -pA +(p+dp)(A+dA)-2p(dA/2)

➢ Remove product of differential terms in above equation

ꝺ/ꝺt (ρvAdx) + d(ρ`v^2` A) = -A dp

➢ By taking limits of dx goes to zero then
ꝺ/ꝺt (ρvAdx) + ꝺ/ꝺx(ρ`v^2`A) = -A ꝺp/ꝺx----8

 

 

Fig (2): -The forces in the x direction acting on the control volume

 

➢ Multiply euqation4 with ‘v’ and subtract with equation 8 then we get

ꝺ/ꝺt (ρvAdx) -v*ꝺ/ꝺt(ρA) - v*ꝺ/ꝺx(ρAv) + ꝺ/ꝺx(ρ`v^2`A) = -A ꝺp/ꝺx

➢ From the expansion of derivative terms we end up with

(𝛒𝐀) ꝺv/ꝺt +( 𝛒𝐀𝐯) ꝺv/ꝺx = -A ꝺp/ꝺx ----9

➢ Here we can eliminate pressure terms by using equation of state

P = ρRT

➢ If you derivate above equation we end up with

ꝺp/ꝺx = R(ρ ∗ꝺT/ꝺx + T*ꝺρ/ꝺx)

➢ Substitute ꝺp/ꝺx in equation 9 then
(𝛒𝐀) ꝺv/ꝺt +( 𝛒𝐀𝐯) ꝺv/ꝺx = -R(𝛒 ∗ꝺT/ꝺx + T*ꝺ𝛒/ꝺx) ----II

Energy equation: -

➢ Consider integral form of the energy equation

ꝺ/ꝺt(∭ρ(e +`v^2`/2) dѵ) + ∬ ρ(e +`v^2`/2)𝐯. 𝐝𝐬 = - ∬(pv). 𝐝𝐬

➢ Like continuity and momentum equation, From figure 1 and figure 2 we can write above
equation as follows

ꝺ/ꝺt[ρ (e+`v^2`/2)Adx] - ρ (e+`v^2`/2)VA + (ρ + dρ) (e+de+`(v+dv)^2`/2)(v+dv)(A+dA) = -[-pvA+(p+dp)(v+dv)(A+dA)-2p(pvdA/2)]

➢ Remove product of differential terms in above equation

ꝺ/ꝺt[ρ (e+`v^2`/2)Adx] + d(ρevA) + d(ρ`v^3`A)/2 = -d(pAv)

➢ We can write above equation as follows
ꝺ/ꝺt[ρ (e+`v^2`/2)Adx] + d[ρ (e+`v^2`/2)vA] = -d(pAv)

➢ By taking limits of dx goes to zero then

ꝺ/ꝺt[ρ (e+`v^2`/2)A] + ꝺ/ꝺx[ρ (e+`v^2`/2)vA] = - ꝺ/ꝺx (pAv) ----10

➢ Multiply equation 9 with ‘v’ and subtract from equation10 then we end up with

ꝺ/ꝺt(ρeA) + ꝺ/ꝺx(ρevA) = -p ꝺ(Av)/ꝺx ------11

➢ Multiply equation 4 with ‘e’ and subtract with equation11 then

(ρA) ꝺe/ꝺt + (ρAv) ꝺe/ꝺx = -p ꝺ(Av)/ꝺx

➢ Expand RHS of above equation

(ρ) ꝺe/ꝺt + (ρv) ꝺe/ꝺx = -p ꝺv/ꝺx - p(v/A) ꝺ(A)/ꝺx
(ρ) ꝺe/ꝺt + (ρv) ꝺe/ꝺx = -p ꝺv/ꝺx – (pv) ꝺ(In A)/ꝺx

➢ For calorically perfect gas internal energy ‘e’ = Cv T .Substitute this in above equation

(ρCv) ꝺT/ꝺt + (ρvCv) ꝺT/ꝺx = -p ꝺv/ꝺx – (pv) ꝺ(In A)/ꝺx

➢ Like momentum equation substitute value of ‘ꝺp/ꝺx’ in above equation then we end up with

(𝛒𝐂𝐯) ꝺT/ꝺt + (𝛒vCv) ꝺT/ꝺx = -(𝛒RT) [ ꝺv/ꝺx +v ꝺ(In A)/ꝺx ] --------III

 

Non-Dimensional of Governing differential equations of non-conservative form: -

➢ Because of this non-dimensional of the flow quantities, we don’t have to worry about units and values of primary variables . These values are always between ‘0’ to ‘1’ for    CD -nozzle .
➢ T’ =T/T0 , ρ’ = ρ/ρ0 , x’ = x/L , v’ = v/a0 , t’= t/(L/ a0) and A’ =A/A*  these are the  dimensionless quantities.

where , Speed of sound `a_0 =(gammaRT)^0.5`

Specific heat capacity at constant colume `C_v =R/(gamma-1)`

➢ Apply these dimensionless quantities to equation I,II and III then

Continuity equation: -

ꝺ(ρ’A’)/ꝺt’*(`rho_0`A*/(L/`a_0`)) + ρ’A’ *(ꝺv’/ꝺx) *( `rho_0`A*`a_0`/L) + ρ’v’ *(ꝺA’/ꝺx)*( `rho_0`A*`a_0`/L) +v’A’ *(ꝺρ’/ꝺx)*( `rho_0`A*`a_0`/L) = 0 ---12

➢ A’ is function of x’ only so there is no variation in A’ with respect to time(t’)

So above equation we can written as

ꝺρ’/ꝺt’ = -ρ’*ꝺv’/ꝺx’ - (ρ′v′) ∗ ꝺ (In A’)/ꝺx’ - v′ ∗ ꝺρ′/ꝺx’

Momentum equation: -

-ρ’v’*ꝺv’/ꝺt’(`rho_0a_0^2`/L) + ρ’v’*ꝺv’/ꝺx’(`rho_0a_0^2`/L) = -R(ρ’ꝺT’/ꝺx’+ T’ꝺ ρ’/ꝺx’ )*( `rho_0T_0`/L) -13

We can write R`T_0`/`a_0^2` = 𝛾R`T_0`/𝛾`a_0^2` = `a_0^2`/ 𝛾`a_0^2` = 1/𝛾

So finally, we can write equation13 with above simplification as
ꝺv’/ꝺt’ = -v’* ꝺv’/ꝺx’ – (1/γ)(ꝺT’/ꝺx’ + (T’/ρ′)*ꝺρ’/ꝺx’)

Energy equation: -

ρ’Cv *(ꝺT’/ꝺt’) *(`rho_0a_0T_0`/L) + ρ’Cv*v’ *(ꝺT’/ꝺx’) *( `rho_0a_0T_0`/L) = ρ’RT’[ ꝺv’/ꝺx’+V’*ꝺ(InA’)/ꝺx’]* ( `rho_0a_0T_0`/L) --14

we can write R/Cv = R/R/ (γ − 1) = γ − 1

So finally, we can write equation14 with above simplification as

ꝺT’/ꝺt’ = -v’* ꝺT’/ꝺx’ – (γ − 1)T’ [ ꝺv’/ꝺx’ +v’ ꝺ(In A’)/ꝺx’]

Continuity equation: ꝺρ’/ꝺt’ = -ρ’*ꝺv’/ꝺx’ - (𝛒′𝐯′) ∗ ꝺ (In A’)/ꝺx’ - 𝐯′ ∗ ꝺ𝛒′/ꝺx’ --IV
Momentum equation: ꝺv’/ꝺt’ = -𝐯’* ꝺv’/ꝺx’ – (1/𝛄)(ꝺT’/ꝺx’ + (T’/𝛒′)*ꝺ𝛒’/ꝺx’) --V
Energy equation: ꝺT’/ꝺt’ = -v’* ꝺT’/ꝺx’ – (𝛄 − 𝟏)T’ [ ꝺv’/ꝺx’ +v’ ꝺ(In A’)/ꝺx’] --V

 

CFD-analysis of non-conservation forms by using McCormack’s method: -

 

• For finding numerical solution we have used McCormack’s method here which implicit, finite difference method and used to find time marching solution.
• For implementing this method, we start with dividing of along x-axis of the nozzle into a number of discrete grid points as shown in figure 3.
• ‘N’ is the number of gird points and ‘dx’ is the space between gid points and ‘dx’ is constant always.
• Here we start to find flow variables at grid point ‘i’ by help of the ‘i+1’ and ‘i-1’.

Form equation I,II and III we can write discretized forms of the governing equations as below

 

Fig (3): - Grid point distribution along nozzle

 

Step1(predictor step): -

• Apply forward difference method to find flow variables

(ꝺρ/ꝺt)^t_{Predictor at ‘i’}  =  = -ρ^t_i ((v^t_i+1 - v^t_i )/dx) -ρ^t_i ((InA^t_i+1 - InA^t_i )/dx) -ρ^t_i ((ρ^t_i+1 - ρ^t_i )/dx)

(ꝺv/ꝺt)^t _{Predictor at ‘i’}   = -v^t_i ((v^t_i+1 - v^t_i )/dx)- (1/) ((T^t_i+1 - T^t_i )/dx) + (( T^t_i / ρ^t_i ) *(ρ^t_i+1 - ρ^t_i )/dx)))

(ꝺT/ꝺt)^t _{Predictor at ‘i’}   = -v^t_i ((T^t_i+1 - T^t_i )/dx)- () T^t_i (((v^t_i+1 - v^t_i )/dx) + v^t_i ((InA^t_i+1 - InA^t_i )/dx))

• Find flow variables by help of flow gradients

              ρ^t+dt_i   = ρ^t_i  + ((ꝺρ/ꝺt)^t_{Predictor at ‘i’}) *(dt)

             v^t+dt_i   = v^t_i  + ((ꝺv/ꝺt)^t_{Predictor at ‘i’}) *(dt)

             T^t+dt_i   = T^t_i  + ((ꝺT/ꝺt)^t_{Predictor at ‘i’}) *(dt)

 

Step2(Corrector step): -

• Apply Backward difference method to find flow variables

(ꝺρ/ꝺt)^t+dt_{corrector at ‘i’}  =  = -ρ^t+dt_i ((v^t+dt_i - v^t+dt_i-1 )/dx) -ρ^t+dt_i ((InA^t+dt_i - InA^t+dt_i-1 )/dx) -ρ^t+dt_i ((ρ^t+dt_i - ρ^t+dt_i-1 )/dx)

(ꝺv/ꝺt)^t+dt_{corrector at ‘i’}   = -v^t+dt_i ((v^t+dt_i - v^t+dt_i-1 )/dx)- (1/) ((T^t+dt_i - T^t+dt_i-1 )/dx) + (( T^t+dt_i / ρ^t+dt_i ) *(ρ^t+dt_i - ρ^t+dt_i-1 )/dx)))

(ꝺT/ꝺt)^t+dt _{corrector at ‘i’}   = -v^t+dt_i ((T^t+dt_i - T^t+dt_i-1 )/dx)- () T^t+dt_i (((v^t+dt_i - v^t+dt_i-1 )/dx) + v^t+dt_i ((InA^t+dt_i - InA^t+dt_i-1 )/dx))

Step3 :-

• Averaging the gradients from corrector and predictor steps

(ꝺρ/ꝺt)_avg = 0.5*((ꝺρ/ꝺt)^t_{Predictor at ‘i’}   + (ꝺρ/ꝺt)^t+dt_{corrector at ‘i’}   )

(ꝺv/ꝺt)_avg =  0.5*( (ꝺv/ꝺt)^t _{Predictor at ‘i’}   + (ꝺv/ꝺt)^t+dt_{corrector at ‘i’}   )

(ꝺT/ꝺt)_avg = 0.5*( (ꝺT/ꝺt)^t _{Predictor at ‘i’}   + (ꝺT/ꝺt)^t+dt _{corrector at ‘i’}   )

              ρ^t+dt_i   = ρ^t_i  + ((ꝺρ/ꝺt)^t_avg) *(dt)

             v^t+dt_i   = v^t_i  + ((ꝺv/ꝺt)^t_avg) *(dt)

             T^t+dt_i   = T^t_i  + ((ꝺT/ꝺt)^t_avg) *(dt)

• These are the final corrected values for time step ‘t+dt’ .
• For each time step we repeat the same process for each grid points except boundary points

 

Courant number: -

• Stability of numerical solution is checked by courant number.
• In this particular problem we consider courant number as 0.5.

Here courant number(C) = dt*(a+v)/dx

Where, ‘a’ is speed of sound and ‘v’ is velocity

• For non-dimensional numbers also we end up with same formula and instead of dimensional quantities we use non dimensional quantities for speed of sound and velocity
• Form courant number formula we can find minimum time step value

                              dt = C*dx/(a+v)

Boundary conditions:

 

• In figure 3 first and last gird points are considered as boundary points.
• These values are not varying with time

Density at 1^st node ‘ρ’ = 1

Velocity at 1^st node ‘v’ = v₁ = 2v₂ -v₃

Temperature at 1^st node ‘T’ = 1

• Above mentioned values are for non-dimensional variables.

Density at N^th node ‘ρ’ = 2 ρ _(N-1) - ρ _(N-2)

Velocity at N^th node ‘v’ = v₁ = 2v_(N-1) -v_(N-2)

Temperature at N^th node ‘T’ = 2T_(N-1) -T_(N-2)

• Above mentioned values are for non-dimensional variables.

Initial conditions: -

• These initial conditions are same and fix at t=0 sec
• They will vary for next time steps
• Flow variables at t=0 sec are vary with respect to ‘x’ as we can see below

 Density (ρ) = 1- 0.3146x

Temperature (T)= 1- 0.2314x

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

Area (A) = 1+ 2.2 (x-1.5)²

• These are valid for 0 <= x <= L where ‘ L’ is the length of nozzle.

 

OUTPUT: -

 

 

                             Fig(4)                                                                           Fig(5)

                                           Fig(6)

 

Conclusions from above graphs: -

• From figure 4 we can say that we got converge values of primitive values at 1400 time step. We can see the steady state values in figure 5.
• From figure 6 we can say that mass flow rate distribution by increasing the time step value and we can say that mass flow rate is adjusting by increasing time steps and It became stable .

Conservative forms of the governing differential equations derivation: -

 

From equation 12 we can write dimensionless continuity equation as

ꝺ/ꝺt’(ρ′A′) + ꝺ/ꝺx(ρ′A′v′) = 0 ---15

➢ From equation 13 we can write dimensionless momentum equation as in terms of pressure
term

ꝺ/ꝺt’(ρ′A′v′) + ꝺ/ꝺx’(ρ′A′`v'^2` + (1/ γ)p’A’) = (1/ γ)p’ ꝺ A’/ꝺx’ ----16

➢ From equation 10 we can write dimensionless energy equation as
ꝺ/ꝺt’[ρ′ (e’+`v'^2`2)A’] + ꝺ/ꝺx’[ρ′ (e’+`v'^2`/2)v’A’+ p’A’v’)] = 0 ----17

➢ Let’s consider

U1 = ρ′A′
U2 = ρ′A′v′
U3 = ρ′ (e’/( γ − 1)+ (γ/2) `v'^2`)A’
F1 = ρ′A′v′ ----18
F2 = ρ′A′`v'^2` + (1/ γ)p’A’ -----19
F3 = ρ′ (e’/ ( γ − 1)+ (γ/2) `v'^2`) v’A’ + ρ′A′v′ -----20
J2 = (1/ γ)p’ ꝺ A’/ꝺx’ --21

➢ ‘U’ is called the solution vector = [𝑈1 𝑈2 𝑈3]′
➢ ‘F’ is called the flux vector = [𝐹1 𝐹2 𝐹3]′

➢ Source term ‘J2’ as mention above
Then our final governing equations are in conservation form

Continuity equation:    ꝺU₁/ꝺt’ = -ꝺF₁/ꝺx’    --VII

Momentum equation:   ꝺU₂/ꝺt’ = -ꝺF₂/ꝺx’ + J₂ --VIII

Energy equation:     ꝺU₃/ꝺt’ = -ꝺF₃/ꝺx’   ---IX

  

CFD-analysis of conservation forms by using McCormack’s method: -

• Like we discuss above for non-conservative form, here we have predictive step and corrective step.
• Before we move on to the predictive step we have to find F₁,F₂ and F₃
• We can find F₁,F₂ and F₃ as the function of U₁,U₂ and U₃ like below

    F₁  = U₂     ----22

    F₂  = (U₂²/U₁) +(`(gamma-1)/gamma`)/(U₃ – (`gamma/2`) U₂²/U₁)  -----23

    F₃  =    (U₂ U₃/ U₁) –  (`gamma(gamma-1)/2`)/2 (U₂³/U₁²)    -----24

 

• We use F₁,F₂ and F₃ as inputs and same as non-conservative form we found gradients of primitive variables.
• Use that gradients to find primitive variables for next corrective step.
• In corrective step find gradients of primitive variables and use them to find final primitive variables for first time step.
• Repeat same method to find values for each time step till we attain steady state.

Boundary conditions:

• In figure 3 first and last gird points are considered as boundary points.
• These values are not varying with time

 At 1^st node U2(1) = 2U₂(2) -U₂(3)

• Above mentioned values are for non-dimensional variables.

 At ‘n^th’ node

‘U1’ = 2U1(n-1)-U1(n-2)

‘U2’ = 2U2(n-1)-U2(n-2)

‘U3’ = 2U3(n-1)-U3(n-2)

 

• Above mentioned values are for non-dimensional variables.

Initial conditions: -

• These initial conditions are same and fix at t=0 sec
• They will vary for next time steps
• Flow variables at t=0 sec are vary with respect to ‘x’ as we can see below

For 0 to 0.5 meters

 Density (ρ) = 1 kg/m³

Temperature (T)= 1 k

     For 0.5 to 1.5 meters

Density (ρ) = 1- 0.366(x-0.5)

Temperature (T)= 1- 0.167(x-0.5)

For 1.5 to 3 meters

Density (ρ) = 0.634 - 0.3146(x-1.5)

Temperature (T)= 0.833 - 0.3507(x-1.5)

• Area we can given as

Area (A) = 1+ 2.2 (x-1.5)²

• Considering constant initial mass flow rate which is value of 0.59 kg/sec.

 

OUTPUT: -

 

  

                                              Fig(7)                                                                                         Fig(8)

 

                                                       Fig(9)

 

Conclusion form above figures: -

 

• From figure 7 we can say that we got converge values of primitive values at 1400 time step. We can see the steady state values in figure 8.
• From figure 9 we can say that mass flow rate distribution by increasing the time step value and we can say that mass flow rate is adjusting by increasing time steps and It became stable.

 

Comparison of normalized mass flow rate distribution of both forms: -

 

                                                 fig(10)

 

Conclusion form above figures: -

 

• Form figure 9 we can say that large fluctuations are happening while we use non conservative forms because we use density, velocity and temperature as the primitive variables. While using non conservative forms thre is a chance that solution may blow up.
• But while using conservation form we can say that there is no large fluctuations.

 

Grid independent test: -

 

 

                                     Fig(10)                                                                  Fig(11)

 

Conclusion form above figures: -

• From figure 10 and 11 we can say that at number of node points 31 and 61 we achieve same conserve solution for both forms.
• So instead of using 61 node points we can use 31 node points to find solutions.
• It will save computational time and will give accurate solution.

Main differences observed between both forms: -

• Theoretically and mathematically both forms are same but numerically(CFD), there is
• Regarding stability, conserve form solutions are more stable compare to non-conservative forms. We can see this difference in figure 9. Fluctuations in non-conservative form primitive variables are more compare to conservative form primitive variables.
• Regarding making a programme, conservative forms have complex program compare to non-conservative forms.

 

MATLAB-CODE:-

MAIN CODE:-

%% MAIN PROGRAMME
clc;
clear all;
close all;
%% Geometry Inputs
n=31;
L=3;
x=linspace(0,3,n);
dx=x(2)-x(1);
%% Initial profiles at t='0' sec (Dimentionless quantities) for non conservative form
rho_non = 1-(0.3146.*x);
T_non = 1-(0.2314.*x);
v_non=(0.1+(1.09.*x)).*(T_non.^(0.5));
A = 1+(2.2.*(x-1.5).^2);
mass_flow_rate_non_intial= rho_non.*A.*v_non;
%% Courant number for stabilty 
nt=1400;
c=0.5;
a=(T_non).^(0.5); % Spead of sound in terms of non dimensional form
dt=min((c.*dx)./(a+v_non));
gama=1.4;
%% Non-Conservative form
[rho,v,T,rho_non_throat,v_non_throat,T_non_throat,mass_flow_rate] = non_conservative_macormack(nt,n,gama,x,rho_non,T_non,v_non,A,dt);

%% Convergence criteria of primary variables at throat section w.r.t  number of time steps
figure(1)
subplot(3,1,1)
plot(rho_non_throat,'-m','linewidth',2);
ylabel('Density');
title('Convergence criteria of primary variables at throat section w.r.t  number of time steps');
subplot(3,1,2)
plot(v_non_throat,'-m','linewidth',2);
ylabel('Velocity');
subplot(3,1,3)
plot(T_non_throat,'-m','linewidth',2);
xlabel('Number of time steps');
ylabel('Temperature');
%% Mass flow rate at diffrent time steps for non conservative form
figure(2)
plot(x,mass_flow_rate_non_intial,'linewidth',2 );
hold on;
plot(x,mass_flow_rate(1,:),'linewidth',2 );
hold on;
plot(x,mass_flow_rate(2,:),'linewidth',2 ); %mass flow rate at steady state
hold on;
plot(x,mass_flow_rate(3,:),'linewidth',2 ); %mass flow rate at steady state
hold on;
plot(x,mass_flow_rate(4,:),'linewidth',2 ); %mass flow rate at steady state
plot(x,mass_flow_rate(5,:),'linewidth',2 );
hold on;
plot(x,mass_flow_rate(6,:),'linewidth',2 ); %mass flow rate at steady state
hold on;
xlabel('length of the nozzle in meters');
ylabel('mass flow rate distribution');
title('Variation of mass flow rate distribution inside the nozzle at diffrent timesteps during time marching ');
legend(' 0_t_h time step','50_t_h time step','100_t_h time step','200_t_h time step','300_t_h time step','600_t_h time step',' 1400_t_h time step ');
%% Steady State distribution of primary variables inside the nozzle for non-convergence form
s_s_non_con_pv =[rho' v' T'];
figure(3)
subplot(3,1,1)
plot(x,s_s_non_con_pv(:,1),'-r','linewidth',2);
ylabel('Density');
title('Primary variables inside the nozzle for non-conservative form at steady state');
subplot(3,1,2)
plot(x,s_s_non_con_pv(:,2),'-r','linewidth',2);
ylabel('Velocity');
subplot(3,1,3)
plot(x,s_s_non_con_pv(:,3),'-r','linewidth',2);
xlabel('Length of nozzle in meters');
ylabel('Temperature');

%% Initial profiles at t='0' sec (Dimentionless quantities) for conservative form
[rho_con,T_con] = intial_vlaues_conservative(x,n);
u1 = rho_con.*A;
v_con = 0.59./u1;
u2 = rho_con.*A.*v_con;
mass_flow_rate_con_intial= u2;
u3= rho_con.*A.*(((T_con)/(gama-1))+((gama/2).*v_con.^(2)));
%% Conservative form
[u1_con,u2_con,u3_con,u1_t,u2_t,u3_t,mass_flow_rate_con] = conservative_macormack(u1,u2,u3,A,dx,dt,gama,n(1),nt,T_con,rho_con);

%% Convergence criteria of primary variables at throat section w.r.t  number of time steps
figure(4)
subplot(3,1,1)
plot(u1_t,'-m','linewidth',2);
ylabel('U1');
title('Convergence criteria of primary variables at throat section w.r.t  number of time steps');
subplot(3,1,2)
plot(u2_t,'-m','linewidth',2);
ylabel('U2');
subplot(3,1,3)
plot(u3_t,'-m','linewidth',2);
xlabel('Number of time steps');
ylabel('U3');
%% Steady State distribution of primary variables inside the nozzle for convergence form
s_s_con_pv =[u1_con' u2_con' u3_con'];
figure(5)
subplot(3,1,1)
plot(x,s_s_con_pv(:,1),'-r','linewidth',2);
ylabel('U1');
title('Primary variables inside the nozzle for conservative form at steady state');
subplot(3,1,2)
plot(x,s_s_con_pv(:,2),'-r','linewidth',2);
ylabel('U2');
ylim([0.58 0.59]);
subplot(3,1,3)
plot(x,s_s_con_pv(:,3),'-r','linewidth',2);
xlabel('length of nozzle in meters');
ylabel('U3');
%% Mass flow rate at diffrent time steps for conservative form
figure(6)
plot(x,mass_flow_rate_con_intial,'linewidth',2 );
hold on;
plot(x,mass_flow_rate_con(1,:),'linewidth',2 );
hold on;
plot(x,mass_flow_rate_con(2,:),'linewidth',2 ); %mass flow rate at steady state
hold on;
plot(x,mass_flow_rate_con(3,:),'linewidth',2 ); %mass flow rate at steady state
hold on;
plot(x,mass_flow_rate_con(4,:),'linewidth',2 ); %mass flow rate at steady state
plot(x,mass_flow_rate_con(5,:),'linewidth',2 );
hold on;
plot(x,mass_flow_rate_con(6,:),'linewidth',2 ); %mass flow rate at steady state
hold on;
xlabel('length of the nozzle in meters');
ylabel('mass flow rate distribution');
title('Variation of mass flow rate distribution inside the nozzle at diffrent timesteps during time marching ');
legend(' 0_t_h time step','50_t_h time step','100_t_h time step','200_t_h time step','300_t_h time step','600_t_h time step',' 1400_t_h time step ');

%% comparision between conservative and non conservative form non dimensional mass flow rate
mass_flow_rate_analytical = 0.579.*ones(1,31);
figure(7)
plot(x,mass_flow_rate(6,:),'linewidth',2 );
hold on;
plot(x,u2_con,'linewidth',2 ); %mass flow rate at steady state
hold on;
plot(x,mass_flow_rate_analytical,'linewidth',2 ); %mass flow rate at steady state
legend('Non conservative mass flow rate','conservative mass flow rate','analytical mass flow rate');
xlabel('length of nozzle in meters');
ylabel('Mass flow rate');
title('Variation of mass flow rate inside the nozzle at steady state');
%%  Grid independence test
%Geometry Inputs
n1=61;
L=3;
x1=linspace(0,3,n1);
dx_n1=x(2)-x(1);
%Initial profiles at t='0' sec (Dimentionless quantities) for non conservative form
rho_non_n1 = 1-(0.3146.*x1);
T_non_n1 = 1-(0.2314.*x1);
v_non_n1=(0.1+(1.09.*x1)).*(T_non_n1.^(0.5));
A_n1 = 1+(2.2.*(x1-1.5).^2);
% Courant number for stabilty 
nt=1400;
c=0.5;
a_n1=(T_non_n1).^(0.5); % Spead of sound in terms of non dimensional form
dt_n1=min((c.*dx_n1)./(a_n1+v_non_n1));
gama=1.4;
% Non-Conservative form
[rho_n1,v_n1,T_n1] = non_conservative_macormack(nt,n1,gama,x1,rho_non_n1,T_non_n1,v_non_n1,A_n1,dt_n1);
figure(8)
subplot(3,1,1)
plot(x,rho,'linewidth',2);
title('Grid independent test for non-conservative forms of PDEs');
hold on;
plot(x1,rho_n1,'oy','linewidth',1);
ylabel('Density');
legend('Solution for n=31','Solution for n=61');
subplot(3,1,2)
plot(x,v,'linewidth',2);
hold on;
plot(x1,v_n1,'oy','linewidth',1);
ylabel('Velocity');
subplot(3,1,3)
plot(x,T,'linewidth',2);
hold on;
plot(x1,T_n1,'oy','linewidth',1);
xlabel('Length of the nozzle in meters');
ylabel('Temperature');
% Initial profiles at t='0' sec (Dimentionless quantities) for conservative form
[rho_con_n1,T_con_n1] = intial_vlaues_conservative(x1,n1);
u1_n1 = rho_con_n1.*A_n1;
v_con_n1 = 0.59./u1_n1;
u2_n1 = rho_con_n1.*A_n1.*v_con_n1;
u3_n1= rho_con_n1.*A_n1.*(((T_con_n1)/(gama-1))+((gama/2).*v_con_n1.^(2)));
% Conservative form
[u1_con_n1,u2_con_n1,u3_con_n1] = conservative_macormack(u1_n1,u2_n1,u3_n1,A_n1,dx_n1,dt_n1,gama,n1,nt,T_con_n1,rho_con_n1);
figure(9)
subplot(3,1,1)
plot(x,u1_con,'linewidth',2);
hold on;
title('Grid independent test for conservative forms of PDEs');
plot(x1,u1_con_n1,'oy','linewidth',1);
ylabel('U1');
legend('Solution for n=31','Solution for n=61');
subplot(3,1,2)
plot(x,u2_con,'linewidth',2);
hold on;
plot(x1,u2_con_n1,'oy','linewidth',1);
ylim([0.5 0.6]);
ylabel('U2');
subplot(3,1,3)
plot(x,u3_con,'linewidth',2);
hold on;
plot(x1,u3_con_n1,'oy','linewidth',1);
xlabel('Length of the nozzle in meters');
ylabel('U3');

Function code for non-conservative form:

 

%%function code for non-conservative form
function [rho,v,T,d_non_throat,v_non_throat,T_non_throat,mass_flow_rate] = non_conservative_macormack(nt,n,gama,x,rho,T,v,A,dt)

       
%time loop
for k=1:nt
       rho_old = rho;
       T_old=T;
       v_old=v;
       
  %predictive method
    for j=2:n-1
        dv_dx=(v(j+1)-v(j))/(x(j+1)-x(j));
        dloga_dx = ((log(A(j+1)))-log(A(j)))/(x(j+1)-x(j));
        drho_dx = (rho(j+1)-rho(j))/(x(j+1)-x(j));
        dT_dx = (T(j+1)-T(j))/(x(j+1)-x(j));

     %continuty equation
        drho_dt_p(j) = (-rho(j)*(dv_dx))-(rho(j)*v(j)*dloga_dx)-(v(j)*drho_dx);
     %momentum equation
        dv_dt_p(j) = (-v(j)*(dv_dx))-((1/gama)*((dT_dx)+(((T(j)/rho(j)))*drho_dx)));
     %energy equation
        dT_dt_p(j) = (-v(j)*(dT_dx))-(((gama-1)*T(j))*((dv_dx)+(v(j)*dloga_dx)));
        
        
        rho(j) = rho(j)+((drho_dt_p(j))*dt);
         v(j) = v(j)+((dv_dt_p(j))*dt);
         T(j) = T(j)+((dT_dt_p(j))*dt);
    end
%Corrective method
    for i=2:n-1
        dv_dx=(v(i)-v(i-1))/(x(i)-x(i-1));
        dloga_dx = ((log(A(i)))-log(A(i-1)))/(x(i)-x(i-1));
        drho_dx = (rho(i)-rho(i-1))/(x(i)-x(i-1));
        dT_dx = (T(i)-T(i-1))/(x(i)-x(i-1));
     %continuty equation
        drho_dt_c(i) = (-rho(i)*(dv_dx))-(rho(i)*v(i)*dloga_dx)-(v(i)*drho_dx);
        
     %momentum equation
        dv_dt_c(i) = (-v(i)*(dv_dx))-((1/gama)*((dT_dx)+(((T(i)/rho(i)))*drho_dx)));
     %energy equation
        dT_dt_c(i) = (-v(i)*(dT_dx))-(((gama-1)*T(i))*((dv_dx)+(v(i)*dloga_dx)));
        

    end
    
        drho_dt = 0.5*(drho_dt_p + drho_dt_c);
        dT_dt = 0.5*(dT_dt_p + dT_dt_c);
        dv_dt = 0.5*(dv_dt_p + dv_dt_c);
        
    for i= 2:n-1
        rho(i) = rho_old(i)+((drho_dt(i))*dt);
        v(i) = v_old(i)+((dv_dt(i))*dt);
        T(i) = T_old(i)+((dT_dt(i))*dt);
        
    end
     % Boundary conditions
        
        %inlet
        v(1) = 2*v(2)-v(3);
        %outlet
        rho(n)= 2*rho(n-1)-rho(n-2);
        v(n)= 2*v(n-1)-v(n-2);
        T(n) = 2*T(n-1)-T(n-2);
        %
      
          for  i= 2:n-1
       % variations wrt time 
            d_non_throat(k)= rho(16);
            v_non_throat(k)= v(16);
            T_non_throat(k)= T(16);  
        
          end
    
   
    if(k==50)
         mass_flow_rate(1,:) = rho.*A.*v;
    elseif(k==100)
         mass_flow_rate(2,:) = rho.*A.*v;
    elseif(k==200)
         mass_flow_rate(3,:) = rho.*A.*v;
    elseif(k==300)
         mass_flow_rate(4,:) = rho.*A.*v;
    elseif(k==500)
         mass_flow_rate(5,:) = rho.*A.*v;
    elseif(k==nt)
         mass_flow_rate(6,:) = rho.*A.*v;
    end
   
     %plot(x,rho);
     %title(sprintf('time step = %d',k));
     %pause(0.03); 
     
end
    
    
    

end

 Function code for initial values for conservative form:

function [rho_con,T_con] = intial_vlaues_conservative(x,n)
 
if (n==31)
    rho_in = zeros(1,31);
    T_in = zeros(1,31);

    for i=1:5
    rho_in(i) = 1;
    T_in(i) = 1;
    end
    for i=6:16
    rho_in(i) = 1-0.366*(x(i)-0.5);
    T_in(i) = 1-0.167*(x(i)-0.5);
    end
    for i=17:31
    rho_in(i) = 0.634-0.3879*((x(i))-1.5);
    T_in(i) = 0.833-0.3507*(x(i)-1.5);
    end
    rho_con = rho_in ;
    T_con = T_in ;
    elseif (n==61)
    rho_in = zeros(1,61);
    T_in = zeros(1,61);

    for i=1:11
     rho_in(i) = 1;
    T_in(i) = 1;
    end
    for i=12:31
    rho_in(i) = 1-0.366*(x(i)-0.5);
    T_in(i) = 1-0.167*(x(i)-0.5);
    end
    for i=32:61
    rho_in(i) = 0.634-0.3879*((x(i))-1.5);
    T_in(i) = 0.833-0.3507*(x(i)-1.5);
    end
    rho_con = rho_in ;
    T_con = T_in ;
end
end

Functional code for conservative form:

function [flux1,flux2,flux3,u1_t,u2_t,u3_t,mass_flow_rate_con] = conservative_macormack(u1,u2,u3,A,dx,dt,gama,n,nt,T_con,rho_con)

                
 for k=1:nt
     
       u1_old = u1;
       u2_old = u2;
       u3_old = u3;
    
       f1= u2;
       f2 = ((u2.^(2))./u1) + (((gama-1)/gama).*(u3-((gama/2).*((u2.^(2))./u1))));
       f3 = (gama.*((u2.*u3)./u1)) - ((gama*(gama-1)/2).*(((u2.^(3))./(u1.^(2)))));
      
       
       
%predicted step
    for j= 2:n-1
       
     %continuty equation
        du1_dt_p(j) =-((f1(j+1)-f1(j))/dx);
     %momentum equation
         j2(j) = (1/gama).*rho_con(j).*T_con(j).*((A(j+1)-A(j))/dx);
         du2_dt_p(j) = (-(f2(j+1)-f2(j))/dx)+ j2(j);
      %energy equation
         du3_dt_p(j) = -((f3(j+1)-f3(j))/dx);
     
         u1(j)= u1(j)+((du1_dt_p(j))*dt);
         u2(j)= u2(j)+((du2_dt_p(j))*dt);
         u3(j)= u3(j)+((du3_dt_p(j))*dt);
    end
    
       rho_con = u1./A;
       T_con = (gama-1).*((u3./u1)-((gama/2).*((u2./u1).^(2))));
       f1= u2;
       f2 = ((u2.^(2))./u1) + (((gama-1)/gama).*(u3-((gama/2).*((u2.^(2))./u1))));
       f3 = (gama.*((u2.*u3)./u1)) - ((gama*(gama-1)/2).*(((u2.^(3))./(u1.^(2)))));
    
          
%corrected step
      for j= 2:n-1   
           
     %continuty equation
        du1_dt_c(j) =-((f1(j)-f1(j-1))/dx);
     %momentum equation
         j2(j) = (1/gama).*rho_con(j).*T_con(j).*((A(j)-A(j-1))/dx);
         du2_dt_c(j) = (-(f2(j)-f2(j-1))/dx)+ j2(j);
      %energy equation
         du3_dt_c(j) = -((f3(j)-f3(j-1))/dx);

     end
             %averaging
      
          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);
          
     for i=2:n-1
               u1(i) = u1_old(i) +((du1_dt(i))*(dt));
               u2(i) = u2_old(i) +((du2_dt(i))*(dt));
               u3(i) = u3_old(i) +((du3_dt(i))*(dt));    
              
     end
             
      
      %boundary conditions
     
      %inlet
       u2(1) = 2*u2(2)-u2(3);
      %outlet
        
      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);
      
       rho_con = u1./A;
       v_con = u2./u1;
       T_con = (gama-1).*((u3./u1)-((gama/2).*((u2./u1).^(2))));
      
      for i=2:n-1
               u1_var(k)=u1(16);
               u2_var(k) =u2(16);
               u3_var(k)=u3(16);

      end
      
     if(k==50)
         mass_flow_rate_con(1,:) = u2;
    elseif(k==100)
         mass_flow_rate_con(2,:) = u2;
    elseif(k==200)
         mass_flow_rate_con(3,:) = u2;
    elseif(k==300)
         mass_flow_rate_con(4,:) = u2;
    elseif(k==500)
         mass_flow_rate_con(5,:) = u2;
    elseif(k==nt)
         mass_flow_rate_con(6,:) = u2;
    end
           
           
 end
           
           flux1=u1;
           flux2=u2;
           flux3=u3;
           u1_t = u1_var;
           u2_t = u2_var ;
           u3_t = u3_var;
           
      
      
end