MATLAB - Major Project - Simulation of a 1D Super-sonic nozzle flow using MacCormack Technique

Problem Statement: -

Simulate A Quasi 1D Subsonic-Supersonic Nozzle using Conservative and Non-Conservative Forms of Governing Equations using MacCormack Technique.

Desired Objectives: -

1. Derive These governing Equations and Non-Dimensionalise them.

2. Use appropriate Boundary Conditions

3. Compare The Variables in Conservative and Non-Conservative Forms.

4. Compare the mass flow distribution throughout the nozzle at different time Steps

Reference: -

The data for this problem can be taken from Computational Fluid Dynamics, author John D. Anderson, Jr., (1995).

(The data taken from Computational Fluid Dynamics, author John D. Anderson, Jr., (1995) is used for study purposes only and with no intention to infringe the Copyright of the original Author and/or the Publication)

Solution: -

1. AIM: -

A. To Simulate Quasi 1D Subsonic-Supersonic Isentropic Nozzle using MacCormack's Technique in Conservative and Non-Conservative Forms.

B. i. Plot Time-Wise Variation of Primitive Variables.

ii. Plot Mass flow rate Variation at different Time-step within the Nozzle.

iii. Plot Steady-State Variation of Primitive Variables.

C. Comparision of Normalized Mass Flow.

2. Given and Assumed: -

A. The Symbols used for the Variables are as follows: -

i. Density `-> rho`

ii. Velocity(Along x) `-> V`

iii. Temperature `-> T`

iv. Pressure `-> P`

v. Internal Energy `-> e`

vi. Mass Flow rate `-> dot(m)`

vii. Specific Heat Ratio `-> gamma`

viii. Specific Heat at Constant Volume `-> c_v`

ix. Area of Cross-Section(Along x) `-> A`

x. Speed of Sound `-> a`

B. Known Thermodynamic and Fluid Equations used here: -

i. Internal Energy`(e) = c_vT`

ii. Speed of Sound`(a) = sqrt(gammaRT)`

iii. P = rhoRT`

iv. `c_v = R/(gamma - 1)`

C. Convention Used to denote Different Regions and Type of Variables(Let any Property be `phi`)

i. To Denote the Property at Reservior `-> phi_o`

ii. To Denote the Property at Throat `-> phi^("*")`

iii. To Denote the Property at Exit `-> phi_e`

iv. To Denote Non-Dimensional Property `-> phi^'`

D. Subsonic-Supersonic Nozzle with `A/A^("*") = 1 + 2.2(x' -1.5)^2` and `"Length"(L) = 3`

Fig. 1 Diagrammatic Representation of Subsonic-Supersonic Nozzle(This image is inspired/Re-created from Computational Fluid Dynamics, author John D. Anderson, Jr., (1995) for study purpose only and with no intention to infringe the Copyright of the original Author and/or the Publication)

E. Initial Conditions of the Subsonic-Supersonic Nozzle: -

i. For Non-Conservative Form: -

a. `rho^' = 1 - 0.3146*x^' {x^' in [0,3]}`

b. `T^' = 1 - 0.2314*x^' {x^' in [0,3]}`

c. `V^' = (0.1 +1.09*x^')(T^')^(1/2) {x^' in [0,3]}`

ii. For Conservative Form: -

a. `rho^' = { (1,x in[0,0.5)),(1 - 0.366(x^' - 0.5),x in[0.5,1.5]),(0.634 - 0.3879(x^' - 1.5),x in(1.5,3]):}`

b. `T^' = { (1,x in[0,0.5)),(1 - 0.167(x^' - 0.5),x in[0.5,1.5]),(0.833 - 0.3507(x^' - 1.5),x in(1.5,3]):}`

c. `V^' = dot(m)^'/(rho^'A^')`

d. `dot(m)^' = 0.590`

F. The Analytical mass Flow Rate`dot(m)_("Analytical")^' = 0.579`

3. Deriving Governing Equations: -

To Derive the Governing Equations a small section within the Subsonic-Supersonic Nozzle is Considered as Shown Below in Fig. 2

Fig. 2 Diagrammatic Representation of Control Volume in Subsonic-Supersonic Nozzle(This image is inspired/Re-created from Computational Fluid Dynamics, author John D. Anderson, Jr., (1995) for study purpose only and with no intention to infringe the Copyright of the original Author and/or the Publication)

A. Conservative Form of Governing Equation: -

i. Continuity Equation(mass Conservation): -

In the above Control Volume(Orange Shaded Region in Fig. 2) the Total Mass Exchange happens in two ways: -

a. Net mass entering the control Volume can be calculated as: -

Mass Entering from the Left Side`(dot(m)_(Left)) = rhoAV`

Mass Leaving from the Right Side`(dot(m)_(Right)) = rhoAV + (del(rhoAV))/(delx)`

Hence, Net Mass decrease due to mass leaving`(dot(m)_(Net)) =(dot(m)_(Right)) - (dot(m)_(Left))`

`=> dot(m)_(Net) = rhoAV + (del(rhoAV))/(delx) - rhoAV`

`=> dot(m)_(Net) = (del(rhoAV))/(delx)`

b. The Net Change of mass change within the Control Volume: -

Net Mass decrease within the control Volume`(dot(m)_(time)) = -(del(rhoA))/(delt)`

As the Net decrease in mass within the Control Volume is happening due to the Net mass leaving

Hence,

`dot(m)_(Net) = dot(m)_(time)`

`=> (del(rhoAV))/(delx) = - (del(rhoA))/(delt)`

`=> (del(rhoA))/(delt) + (del(rhoAV))/(delx) = 0`

ii. Momentum Equation: -

Assuming the Control Volume Moving Along the Flow in the x direction and for an Inviscid Fluid flow the Forces acting on the Control Volume are shown in Fig. 3

Fig. 3 Diagrammatic Representation of Forces on Control Volume in Subsonic-Supersonic Nozzle(This image is inspired/Re-created from Computational Fluid Dynamics, author John D. Anderson, Jr., (1995) for study purpose only and with no intention to infringe the Copyright of the original Author and/or the Publication)

a. Net Force Acting on the Control Volume as shown in Fig. 3 is as follows

`F_("net") = PA + 2(P(dA)/2) - (P + dP)(A + dA)`

`=> F_("net") = -(AdP + dPdA)`

As `dPdA` is a very small Quantity and so it can be neglected

`=> F_("net") = - AdP" " (eq. 2)`

Now,

As `d(PA) = PdA + AdP`

`=> AdP = d(PA) - PdA" "(eq. 3)`

Now,

Putting `(eq. 3)` Values in `(eq. 2)`

`=> F_("net") = -(d(PA) - PdA)" " (eq. 4)`

b. The Effect of Net Force Applied on the Control Volume in the x direction is given by `ma_x`

Effect of Net Force = `ma_x`

As the Control Volume is moving with the fluid,

`=> ma_x = rhoAdx(DV)/(Dt)`

`=> ma_x = (rhoAdx)((delV)/(delt) + V(delV)/(delx))`

`=> ma_x = (rhoA(delV)/(delt) + rhoAV(delV)/(delx))dx" "(eq. 5)`

Now,

As `(del(rhoAV))/(delt) = V(del(rhoA))/(delt) + rhoA(del(V))/(delt)`

`=> rhoA(del(V))/(delt) = (del(rhoAV))/(delt) - V(del(rhoA))/(delt)" "(eq. 6)`

And

As `(del(rhoAV^2))/(delx) = V(del(rhoAV))/(delx) + rhoAV(del(V))/(delx)`

`=> rhoAV(del(V))/(delx) = (del(rhoAV^2))/(delx) - V(del(rhoAV))/(delx)" "(eq. 7)`

Now,

Putting Values of `(eq. 6)` and `(eq. 7)` in `(eq. 5)`

`ma_x = ((del(rhoAV))/(delt) - V(del(rhoA))/(delt) + (del(rhoAV^2))/(delx) - V(del(rhoAV))/(delx))dx`

`=> ma_x = ((del(rhoAV))/(delt) + (del(rhoAV^2))/(delx) - V(del(rhoA))/(delt) - V(del(rhoAV))/(delx))dx`

`=> ma_x = ((del(rhoAV))/(delt) + (del(rhoAV^2))/(delx) - V((del(rhoA))/(delt) + (del(rhoAV))/(delx)))dx`

Now,

As from `(eq. 1)` the value of `(del(rhoA))/(delt) + (del(rhoAV))/(delx) = 0`

`=> ma_x = ((del(rhoAV))/(delt) + (del(rhoAV^2))/(delx))dx" "(eq. 8)`

As the Net Force acting on the System is causing the net Effect on the System

So,

`ma_x = F_("net")`

`=> ((del(rhoAV))/(delt) + (del(rhoAV^2))/(delx))dx = -(d(PA) - PdA)`

Dividing by `dx`

`=> (del(rhoAV))/(delt) + (del(rhoAV^2))/(delx) = -(del(PA))/(delx) + P(delA)/(delx)`

`=> (del(rhoAV))/(delt) + (del(rhoAV^2 + PA))/(delx) = P(delA)/(delx)`

iii. Energy Equation(Energy Conservation): -

Assuming the Control Volume shown in Fig. 3 is moving with the Fluid in the x direction and For an Inviscid Fluid the change of total Energy per unit time of the given fluid due to external effect and internal effect is given as follows: -

a. The Net Energy Added to the Control Volume per unit time is given by: -

`dot(E)_("external") = PAV + 2(PV(dA)/2) - (P + dP)(A + dA)(V + dV)`

`=> dot(E)_("external") = PAV + 2(PV(dA)/2) - (PAV + AVdP + PVdA + VdPdA + PAdV + AdPdV + PdAdV + dPdAdV)`

As `dPdAdV` is very small quantities and so it can be neglected

`=> dot(E)_("external") = - AVdP - VdPdA - PAdV - PdAdV - AdPdV`

`=> dot(E)_("external") = - AVdP - PAdV - (dA)VdP - (dA)PdV - AdPdV`

`=> dot(E)_("external") = - d(PAV)" "(eq. 10)`

b. The Net Rate of increase of energy of the control fluid moving with the flow: -

`dot(E)_("internal") = ubrace(m(De)/(Dt))_("Internal Energy") + ubrace(m(D(V^2)/2)/(Dt))_("Kinetic Energy")`

`=> dot(E)_("internal") = rhoAdx((D(e + (V^2)/2))/(Dt))`

`=> dot(E)_("internal") = rhoA((del(e + (V^2)/2))/(delt) + V(del(e + (V^2)/2))/(delx))dx`

`=> dot(E)_("internal") = (rhoA(del(e + (V^2)/2))/(delt) + rhoAV(del(e + (V^2)/2))/(delx))dx" "(eq. 11)`

Now,

As `(del(rhoA(e + V^2/2)))/(delt) = rhoA(del(e + V^2/2))/(delt) + (e + V^2/2)(del(rhoA))/(delt)`

`=> rhoA(del(e + V^2/2))/(delt) = (del(rhoA(e + V^2/2)))/(delt) - (e + V^2/2)(del(rhoA))/(delt)" "(eq. 12)`

And,

As `(del(rhoAV(e+V^2/2)))/(delx) = rhoAV(del(e+V^2/2))/(delx) + (e+V^2/2)(del(rhoAV))/(delx)`

`=> rhoAV(del(e+V^2/2))/(delx) = (del(rhoAV(e+V^2/2)))/(delx) - (e+V^2/2)(del(rhoAV))/(delx)" "(eq. 13)`

Now,

Putting `(eq. 12)` and `(eq. 13)` in `(eq. 11)`

`dot(E)_("internal") = ((del(rhoA(e + V^2/2)))/(delt) - (e + V^2/2)(del(rhoA))/(delt) + (del(rhoAV(e+V^2/2)))/(delx) - (e+V^2/2)(del(rhoAV))/(delx))dx`

`=> dot(E)_("internal") = ((del(rhoA(e + V^2/2)))/(delt) + (del(rhoAV(e+V^2/2)))/(delx) - (e + V^2/2)((del(rhoA))/(delt) + (del(rhoAV))/(delx)))dx`

Now,

As from `(eq. 1)` the Value of `(del(rhoA))/(delt) + (del(rhoAV))/(delx) = 0`

`=> dot(E)_("internal") = ((del(rhoA(e + V^2/2)))/(delt) + (del(rhoAV(e+V^2/2)))/(delx))dx`

As the Internal increase of total Energy of the Control Volume is happening due to the External Energy

So,

`dot(E)_("internal") = dot(E)_("external")`

`=> ((del(rhoA(e + V^2/2)))/(delt) + (del(rhoAV(e+V^2/2)))/(delx))dx = - d(PAV)`

Dividing by `dx`

`=> (del(rhoA(e + V^2/2)))/(delt) + (del(rhoAV(e+V^2/2)))/(delx) = - (del(PAV))/(delx)`

`=> (del(rhoA(e + V^2/2)))/(delt) + (del(rhoAV(e+V^2/2) + PAV))/(delx) = 0`

iv. All Governing Equations in Conservative Form: -

The Derivations done above are the governing Equation for the Subsonic-Supersonic Isentropic Quasi 1D Nozzle Flow and these are used to predict(Study) the behavior of the Variables like `rho, V, T, P` within the nozzle. The following are the Governing Equations for Subsonic-Supersonic Quasi 1D Nozzle in Conservative form: -

B. Non-Conservative Form of Governing Equation: -

i. Continuity Equation(mass Conservation): -

In the Section(Conservative Form of Governing Equation) the Continuity Equation is derived, now as the Conservative and Non-Conservative forms of Equations are Mathematically similar and so here `(eq. 1)` is being used to derive the Non-Conservative form of Continuity Equation

As per `(eq. 1)`

`(del(rhoA))/(delt) + (del(rhoAV))/(delx) = 0`

`=> (del(rhoA))/(delt) = - (del(rhoAV))/(delx)`

`=> (del(rhoA))/(delt) = - [rhoA(del(V))/(delx) + V(del(rhoA))/(delx)]`

`=> (del(rhoA))/(delt) = - [rhoA(del(V))/(delx) + rhoV(del(A))/(delx) + AV(del(rho))/(delx)]`

As the Area of the Cross-Section of Nozzle doesn't change with time `(=> A = f^n(x) only)`

Hence,

`=> A(del(rho))/(delt) = - [rhoA(del(V))/(delx) + rhoV(del(A))/(delx) + AV(del(rho))/(delx)]`

Dividing both sides with `A`

`=> (del(rho))/(delt) = - [rho(del(V))/(delx) + (rhoV)/A(del(A))/(delx) + V(del(rho))/(delx)]" "(eq. 15)`

And,

As `d(lnA) = (dA)/A" "(eq. 16)`

Now,

Putting `(eq. 16)` in `(eq. 15)`

`=> (del(rho))/(delt) = - [V(del(rho))/(delx) + rho(del(V))/(delx) + rhoV(del(lnA))/(delx)]`

ii. Momentum Equation: -

In the Section(Conservative Form of Governing Equation) the Momentum Equation is derived, now as the Conservative and Non-Conservative forms of Equations are Mathematically similar and so here `(eq. 9)` is being used to derive the Non-Conservative form of Momentum Equation

As per `(eq. 9)`

`(del(rhoAV))/(delt) + (del(rhoAV^2 + PA))/(delx) = P(delA)/(delx)`

`=> (del(rhoAV))/(delt) + (del(rhoAV^2))/(delx) = - (del(PA))/(delx) + P(delA)/(delx)" "(eq. 18)`

Now,

As `(del(rhoAV))/(delt) = rhoA(del(V))/(delt) + V(del(rhoA))/(delt)" "(eq. 19)`

And,

As `(del(rhoAV^2))/(delx) = rhoAV(del(V))/(delx) + V(del(rhoAV))/(delx)" "(eq. 20)`

Now,

Putting `(eq. 19)` and `(eq. 20)` in `(eq. 18)`

`=> rhoA(del(V))/(delt) + V(del(rhoA))/(delt) + rhoAV(del(V))/(delx) + V(del(rhoAV))/(delx) = - (del(PA))/(delx) + P(delA)/(delx)`

`=> rhoA(del(V))/(delt) + rhoAV(del(V))/(delx) + V[(del(rhoA))/(delt) + (del(rhoAV))/(delx)] = - [(del(PA))/(delx) - P(delA)/(delx)]`

As from `(eq. 1)` the Value of `(del(rhoA))/(delt) + (del(rhoAV))/(delx) = 0`

`=> rhoA(del(V))/(delt) + rhoAV(del(V))/(delx) = - [(del(PA))/(delx) - P(delA)/(delx)]" "(eq. 21)`

Now,

As `(del(PA))/(delx) = A(delP)/(delx) + P(delA)/(delx)`

`=> (del(PA))/(delx) - P(delA)/(delx) = A(delP)/(delx)" "(eq. 22)`

Now,

Putting the Value of `(eq. 22)` in `(eq. 21)`

`=> rhoA(del(V))/(delt) + rhoAV(del(V))/(delx) = - A(delP)/(delx)" "(eq. 23)`

Dividing Both Sides with `A`

`=> rho(del(V))/(delt) + rhoV(del(V))/(delx) = - (delP)/(delx)`

`=> rho(del(V))/(delt) = - [rhoV(del(V))/(delx) + (delP)/(delx)]`

iii. Energy Equation(Energy Conservation): -

In the Section(Conservative Form of Governing Equation) the Energy Conservation Equation is derived, now as the Conservative and Non-Conservative forms of Equations are Mathematically similar and so here `(eq. 14)` is being used to derive the Non-Conservative form of Energy Conservation Equation

As per `(eq. 14)`

`(del(rhoA(e + V^2/2)))/(delt) + (del(rhoAV(e+V^2/2) + PAV))/(delx) = 0`

`=> (del(rhoA(e + V^2/2)))/(delt) + (del(rhoAV(e+V^2/2)))/(delx) = - (del( PAV))/(delx)`

`=> (del(rhoAe))/(delt) + (del(rhoAV^2/2))/(delt)+ (del(rhoAVe))/(delx) + (del(rhoAV(V^2/2)))/(delx) = - (del( PAV))/(delx)`

`=> (del(rhoAe))/(delt) + (del(rhoAVe))/(delx) + (del(rhoA(V^2/2)))/(delt) + (del(rhoAV(V^2/2)))/(delx) = - (del( PAV))/(delx)" "(eq.25)`

Now,

As `(del(rhoAe))/(delt) = rhoA(del(e))/(delt) + e(del(rhoA))/(delt)" "(eq. 26)`

And,

As `(del(rhoAVe))/(delx) = rhoAV(del(e))/(delx) + e(del(rhoAV))/(delx)" "(eq. 27)`

Now,

Putting Values of `(eq. 26)` and `(eq. 27)` in `(eq. 25)`

`=> rhoA(del(e))/(delt) + e(del(rhoA))/(delt) + rhoAV(del(e))/(delx) + e(del(rhoAV))/(delx) + (del(rhoA(V^2/2)))/(delt) + (del(rhoAV(V^2/2)))/(delx) = - (del( PAV))/(delx)`

`=> rhoA(del(e))/(delt) + rhoAV(del(e))/(delx) + e[(del(rhoA))/(delt) + (del(rhoAV))/(delx)] + (del(rhoA(V^2/2)))/(delt) + (del(rhoAV(V^2/2)))/(delx) = - (del( PAV))/(delx)`

Now,

As from `(eq. 1)` the Value of `(del(rhoA))/(delt) + (del(rhoAV))/(delx) = 0`

`=> rhoA(del(e))/(delt) + rhoAV(del(e))/(delx) + (del(rhoA(V^2/2)))/(delt) + (del(rhoAV(V^2/2)))/(delx) = - (del( PAV))/(delx)" "(eq. 28)`

Now,

As `(del(rhoA(V^2/2)))/(delt) = rhoA(del(V^2/2))/(delt) + (V^2/2)(del(rhoA))/(delt)" "(eq. 29)`

And,

As `(del(rhoAV(V^2/2)))/(delx) = rhoAV(del(V^2/2))/(delx) + (V^2/2)(del(rhoAV))/(delx)" "(eq. 30)`

Now,

Putting the Values of `(eq. 29)` and `(eq. 30)` in `(eq. 28)`

`=> rhoA(del(e))/(delt) + rhoAV(del(e))/(delx) + rhoA(del(V^2/2))/(delt) + (V^2/2)(del(rhoA))/(delt) + rhoAV(del(V^2/2))/(delx) + (V^2/2)(del(rhoAV))/(delx) = - (del( PAV))/(delx)`

`=> rhoA(del(e))/(delt) + rhoAV(del(e))/(delx) + rhoA(del(V^2/2))/(delt) + rhoAV(del(V^2/2))/(delx) + (V^2/2)[(del(rhoA))/(delt) + (del(rhoAV))/(delx)] = - (del( PAV))/(delx)`

Now,

As from `(eq. 1)` the Value of `(del(rhoA))/(delt) + (del(rhoAV))/(delx) = 0`

`=> rhoA(del(e))/(delt) + rhoAV(del(e))/(delx) + rhoA(del(V^2/2))/(delt) + rhoAV(del(V^2/2))/(delx) = - (del( PAV))/(delx)`

`=> rhoA(del(e))/(delt) + rhoAV(del(e))/(delx) + rhoAV(del(V))/(delt) + rhoAV^2(del(V))/(delx) = - (del( PAV))/(delx)`

`=> rhoA(del(e))/(delt) + rhoAV(del(e))/(delx) + V[rhoA(del(V))/(delt) + rhoAV(del(V))/(delx)] = - (del( PAV))/(delx)`

Now,

As from `(eq. 23)` the Value of `rhoA(del(V))/(delt) + rhoAV(del(V))/(delx) = - A(delP)/(delx)`

`=> rhoA(del(e))/(delt) + rhoAV(del(e))/(delx) + V[- A(delP)/(delx)] = - (del( PAV))/(delx)`

`=> rhoA(del(e))/(delt) + rhoAV(del(e))/(delx) - AV(delP)/(delx) = - (del( PAV))/(delx)`

`=> rhoA(del(e))/(delt) + rhoAV(del(e))/(delx) = - [(del( PAV))/(delx) - AV(delP)/(delx)]" "(eq. 31)`

Now,

As `(del(PAV))/(delx) = AV(del(P))/(delx) + PA(del(V))/(delx) + PV(del(A))/(delx)`

`=> (del(PAV))/(delx) - AV(del(P))/(delx) = PA(del(V))/(delx) + PV(del(A))/(delx)" "(eq. 32)`

Now,

Putting Values of `(eq. 32)` in `(eq. 31)`

`=> rhoA(del(e))/(delt) + rhoAV(del(e))/(delx) = - [PA(del(V))/(delx) + PV(del(A))/(delx)]`

Now,

Dividing `A` on both Sides

`=> rho(del(e))/(delt) + rhoV(del(e))/(delx) = - [P(del(V))/(delx) + PV/A(del(A))/(delx)]`

Now,

As from `(eq. 16)` value of `(dA)/A = d(lnA)`

`=> rho(del(e))/(delt) + rhoV(del(e))/(delx) = - P(del(V))/(delx) - PV(del(lnA))/(delx)`

`=> rho(del(e))/(delt) = - [rhoV(del(e))/(delx) + P(del(V))/(delx) + PV(del(lnA))/(delx)]`

iv. All Governing Equations in Non-Conservative Form: -

The Derivations done above are the governing Equation for the Subsonic-Supersonic Isentropic Quasi 1D Nozzle Flow and these are used to predict(Study) the behavior of the Variables like `rho, V, T, P` within the nozzle. The following are the Governing Equations for Subsonic-Supersonic Quasi 1D Nozzle in Non-Conservative form: -

C. Non-dimensionalizing Governing Equations: -

Non-Dimensionalizing is a process of converting each Variable into a dimension-less Value. The benefits of this are listed below: -

i. It helps to generate a scalable form of Governing Equation which helps in testing by setting a certain scale for prototyping and developing the real-size machine by setting the scale

ii. It helps to Standardize the Equations

iii. The equations after Non-Dimensionalizing look aesthetic

For Non-Dimensionalizing each term has to be converted into a ratio w.r.t. some standard Values. Here the Non-Dimensional Terms are as follows: -

a. `rho^' = rho/rho_o`

b. `V^' = V/a_o`

c. `T^' = T/T_o`

d. `P^' = P/P_o`

e. `A^' = A/A^("*")`

f. `e^' = e/e_o`

g. `m^' = m/m_o`

h. `t^' = t/(L/a_o)`

i. `x^' = x/L`

i. Non-Dimensionalizing Conservative Form: -

a. Continuity Equation: -

`(del(rhoA))/(delt) + (del(rhoAV))/(delx) = 0`

Applying the Non-Dimensional terms in the above Equation

`=> (del(rho^'A^'))/(delt^')((rho_oA^("*"))/(L/a_o)) + (del(rho^'A^'V^'))/(delx^')((rho_oA^("*")a_o)/(L)) = 0`

Now,

Dividing both sides by `((rho_oA^("*")a_o)/(L))`

`=> (del(rho^'A^'))/(delt^') + (del(rho^'A^'V^'))/(delx^') = 0`

b. Momentum Equation: -

`(del(rhoAV))/(delt) + (del(rhoAV^2 + PA))/(delx) = P(delA)/(delx)`

Applying the Non-Dimensional terms in the above Equation

`=> (del(rho^'A^'V^'))/(delt^')((rho_oA^("*")a_o)/(L/a_o)) + (del(rho^'A^'(V^')^2(rho_oA^("*")a_o^2) + P^'A^'(P_oA^("*"))))/(delx^'(L)) = P^'(delA^')/(delx^')((P_oA^("*"))/(L))`

`=> (del(rho^'A^'V^'))/(delt^')((rho_oA^("*")a_o^2)/(L)) + (del(rho^'A^'(V^')^2(rho_oA^("*")a_o^2) + P^'A^'(P_oA^("*"))))/(delx^'(L)) = P^'(delA^')/(delx^')((P_oA^("*"))/(L))`

Now,

As `P_o = rho_oRT_o` and `a_o = sqrt(gammaRT_o)`

`=> (del(rho^'A^'V^'))/(delt^')((rho_oA^("*")(gammaRT_o))/(L)) + (del(rho^'A^'(V^')^2(rho_oA^("*")(gammaRT_o)) + P^'A^'((rho_oRT_o)A^("*"))))/(delx^'(L)) = P^'(delA^')/(delx^')(((rho_oRT_o)A^("*"))/(L))`

`=> (del(rho^'A^'V^'))/(delt^')((gammarho_oA^("*")RT_o)/(L)) + (del(rho^'A^'(V^')^2gamma + P^'A^'))/(delx^')((rho_oA^("*")RT_o)/(L)) = P^'(delA^')/(delx^')((rho_oRT_oA^("*"))/(L))`

Now,

Dividing both sides by `((rho_oRT_oA^("*")gamma)/(L))`

`=> (del(rho^'A^'V^'))/(delt^') + (del(rho^'A^'(V^')^2 + (1/gamma)P^'A^'))/(delx^') = (1/gamma)P^'(delA^')/(delx^')`

c. Energy Conservation Equation: -

`(del(rhoA(e + V^2/2)))/(delt) + (del(rhoAV(e+V^2/2) + PAV))/(delx) = 0`

Applying the Non-Dimensional terms in the above Equation

`=> (del(rho^'A^'(e^'(e_o) + (V^')^2/2(a_o^2))))/(delt^')((rho_oA^("*"))/(L/a_o)) + (del(rho^'A^'V^'(e^'(e_o)+(V^')^2/2(a_o^2)) + P^'A^'V^'(P_o)))/(delx^')((A^("*")a_o)/(L)) = 0`

Now,

Dividing both sides by `(A^("*")a_o)/(L)`

`=> (del(rho^'A^'(e^'(e_o) + (V^')^2/2(a_o^2))))/(delt^')(rho_o) + (del(rho^'A^'V^'(e^'(e_o)+(V^')^2/2(a_o^2))(rho_o) + P^'A^'V^'(P_o)))/(delx^') = 0`

Now,

As `P_o = rho_oRT_o`, `a_o = sqrt(gammaRT_o)` and `e_o = c_vT_o`

`=> (del(rho^'A^'(e^'(c_vT_o) + (V^')^2/2(gammaRT_o))))/(delt^')(rho_o) + (del(rho^'A^'V^'(e^'(c_vT_o)+(V^')^2/2(gammaRT_o))(rho_o) + P^'A^'V^'(rho_oRT_o)))/(delx^') = 0`

And,

As `c_v = R/(gamma - 1)`

`=> (del(rho^'A^'(e^'((RT_o)/(gamma - 1)) + (V^')^2/2(gammaRT_o))))/(delt^')(rho_o) + (del(rho^'A^'V^'(e^'((RT_o)/(gamma - 1))+(V^')^2/2(gammaRT_o))(rho_o) + P^'A^'V^'(rho_oRT_o)))/(delx^') = 0`

`=> (del(rho^'A^'((e^')/(gamma - 1) + (gamma)(V^')^2/2)))/(delt^')(rho_oRT_o) + (del(rho^'A^'V^'((e^')/(gamma - 1) + (gamma)(V^')^2/2) + P^'A^'V^'))/(delx^')(rho_oRT_o) = 0`

Now,

Dividing both sides by `rho_oRT_o`

`=> (del(rho^'A^'((e^')/(gamma - 1) + (gamma/2)(V^')^2)))/(delt^') + (del(rho^'A^'V^'((e^')/(gamma - 1) + (gamma/2)(V^')^2) + P^'A^'V^'))/(delx^') = 0`

d. All Non-Dimensional Governing Equations in Conservative Form: -

The conversion done above is for the governing equations of Subsonic-Supersonic Quasi 1D Nozzle Flow and the equations thus developed help to predict(study) and scaled size of Subsonic-Supersonic Nozzle. The following are the Non-Dimensional Subsonic-Supersonic Quasi 1D Nozzle flow equation in Conservative form: -

ii. Non-Dimensionalizing Non-Conservative Form: -

a. Continuity Equation: -

`(del(rho))/(delt) = - [V(del(rho))/(delx) + rho(del(V))/(delx) + rhoV(del(lnA))/(delx)]`

Applying the Non-Dimensional terms in the above Equation

`=> (del(rho^'))/(delt^')(rho_o/(L/a_o)) = - [V^'(del(rho^'))/(delx^') + rho^'(del(V^'))/(delx^') + rho^'V^'(del(ln(A^'A^("*"))))/(delx^')](rho_o/(L/a_o))`

`=> (del(rho^'))/(delt^')(rho_o/(L/a_o)) = - [V^'(del(rho^'))/(delx^') + rho^'(del(V^'))/(delx^') + rho^'V^'(del(ln(A^') + ln(A^("*"))))/(delx^')](rho_o/(L/a_o))`

`=> (del(rho^'))/(delt^')(rho_o/(L/a_o)) = - [V^'(del(rho^'))/(delx^') + rho^'(del(V^'))/(delx^') + rho^'V^'(del(lnA^'))/(delx^') + rho^'V^'(del(lnA^("*")))/(delx^')](rho_o/(L/a_o))`

Now,

As `A^("*") = "constant" => (del(lnA^("*")))/(delx^') = 0`

`=> (del(rho^'))/(delt^')(rho_o/(L/a_o)) = - [V^'(del(rho^'))/(delx^') + rho^'(del(V^'))/(delx^') + rho^'V^'(del(lnA^'))/(delx^')](rho_o/(L/a_o))`

Now,

Dividing both sides with `(rho_oa_o)/L`

`=> (del(rho^'))/(delt^') = - [V^'(del(rho^'))/(delx^') + rho^'(del(V^'))/(delx^') + rho^'V^'(del(lnA^'))/(delx^')]`

b. Momentum Equation: -

`rho(del(V))/(delt) = - [rhoV(del(V))/(delx) + (delP)/(delx)]`

Applying the Non-Dimensional terms in the above Equation

`=> rho^'(del(V^'))/(delt^')((rho_oa_o)/(L/a_o)) = - [rho^'V^'(del(V^'))/(delx^')((rho_oa_o^2)/(L)) + (delP^')/(delx^')(P_o/L)]`

Now,

As `P_o = rho_oRT_o` and `a_o = sqrt(gammaRT_o)`

`=> rho^'(del(V^'))/(delt^')((rho_o(gammaRT_o))/(L)) = - [rho^'V^'(del(V^'))/(delx^')((rho_o(gammaRT_o))/(L)) + (delP^')/(delx^')((rho_oRT_o)/L)]`

`=> rho^'(del(V^'))/(delt^')((gammarho_oRT_o)/L) = - [rho^'V^'(del(V^'))/(delx^')(gamma) + (delP^')/(delx^')]((rho_oRT_o)/L)`

Now,

Dividing both sides with `(gammarho_oRT_o)/L`

`=> rho^'(del(V^'))/(delt^') = - [rho^'V^'(del(V^'))/(delx^') + (1/(gamma))(delP^')/(delx^')]`

c. Energy Conservation Equation: -

`rho(del(e))/(delt) = - [rhoV(del(e))/(delx) + P(del(V))/(delx) + PV(del(lnA))/(delx)]`

Applying the Non-Dimensional terms in the above Equation

`=> rho^'(del(e^'))/(delt^')((rho_oe_o)/(L/a_o)) = - [rho^'V^'(del(e^'))/(delx^')((rho_oa_oe_o)/(L)) + P^'(del(V^'))/(delx^')((P_oa_o)/(L)) + P^'V^'(del(lnA^'A^("*")))/(delx^')((P_oa_o)/(L))]`

`=> rho^'(del(e^'))/(delt^')((rho_oe_oa_o)/(L)) = - [rho^'V^'(del(e^'))/(delx^')((rho_oa_oe_o)/(L)) + P^'(del(V^'))/(delx^')((P_oa_o)/(L)) + P^'V^'(del(lnA^' + lnA^("*")))/(delx^')((P_oa_o)/(L))]`

`=> rho^'(del(e^'))/(delt^')((rho_oe_oa_o)/(L)) = - [rho^'V^'(del(e^'))/(delx^')((rho_oa_oe_o)/(L)) + P^'(del(V^'))/(delx^')((P_oa_o)/(L)) + P^'V^'((del(lnA^'))/(delx^') + (del(lnA^("*")))/(delx^'))((P_oa_o)/(L))]`

Now,

As `A^("*") = "constant" => (del(lnA^("*")))/(delx^') = 0`

`=> rho^'(del(e^'))/(delt^')((rho_oe_oa_o)/(L)) = - [rho^'V^'(del(e^'))/(delx^')((rho_oa_oe_o)/(L)) + P^'(del(V^'))/(delx^')((P_oa_o)/(L)) + P^'V^'(del(lnA^'))/(delx^')((P_oa_o)/(L))]`

Now,

As `P_o = rho_oRT_o` and `e_o = c_vT_o`

`=> rho^'(del(e^'))/(delt^')((rho_oa_oc_vT_o)/(L)) = - [rho^'V^'(del(e^'))/(delx^')((rho_oa_oc_vT_o)/(L)) + P^'(del(V^'))/(delx^')((rho_oRT_oa_o)/(L)) + P^'V^'(del(lnA^'))/(delx^')((rho_oRT_oa_o)/(L))]`

Now,

As `c_v = R/(gamma -1)`

`=> rho^'(del(e^'))/(delt^')((rho_oa_oRT_o)/(L(gamma - 1))) = - [rho^'V^'(del(e^'))/(delx^')((rho_oa_oRT_o)/(L(gamma - 1))) + P^'(del(V^'))/(delx^')((rho_oRT_oa_o)/(L)) + P^'V^'(del(lnA^'))/(delx^')((rho_oRT_oa_o)/(L))]`

Now,

Dividing both sides with `((rho_oa_oRT_o)/(L(gamma - 1)))`

`=> rho^'(del(e^'))/(delt^') = - [rho^'V^'(del(e^'))/(delx^') + (gamma - 1)P^'(del(V^'))/(delx^') + (gamma - 1)P^'V^'(del(lnA^'))/(delx^')]`

d. All Non-Dimensional Governing Equations in Non-Conservative Form: -

The conversion done above is for the governing equations of Subsonic-Supersonic Quasi 1D Nozzle Flow and the equations thus developed help to predict(study) and scaled size of Subsonic-Supersonic Nozzle. The following are the Non-Dimensional Subsonic-Supersonic Quasi 1D Nozzle flow equation in Non-Conservative form: -

4. Methods used for Computation: -

A. MacCormack's Technique: -

This is an Explicit predictor-corrector Method that has wide applications within CFD and In the current study, the same method is being used. For Explanation, the Non-Conservative Form of Equation is used

i. Predictor Step: -

In this step, the Time derivative predictor value is calculated by applying Forward differencing to the Space Derivative for each equation of Governing Equation

a. Calculating Predictor Time Derivative: -

Continuity Equation in Non-Conservative form: -

`((del(rho))/(delt))_i^p = - [V_i^t(rho_(i+1)^t - rho_i^t)/(trianglex) + rho_i^t(V_(i+1)^t - V_i^t)/(trianglex) + rho_i^tV_i^t(lnA_(i+1)^t - lnA_i^t)/(trianglex)]`

Similarly,

`((del(V))/(delt))_i^p` and `((del(e))/(delt))_i^p` are calculated

b. Calculating Predictor Variables: -

`bar(rho)^(t+trianglet) = rho_i ^t + ((del(rho))/(delt))_i^p trianglet`

Similarly,

`bar(V)^(t+trianglet)` and `bar(e)^(t+trianglet)`

And using Equations to calculate `bar(P)^(t+trianglet)` and `bar(T)^(t+trianglet)`

ii. Corrector Step: -

In this step, the Time derivative Corrector value is calculated by applying Backward differencing to the Space Derivative for each equation of Governing Equation

a. Calculating Corrector Time Derivative: -

Continuity Equation in Non-Conservative form: -

`((del(rho))/(delt))_i^c = - [bar(V)_i^(t+trianglet)(bar(rho)_(i)^(t+trianglet) - bar(rho)_(i-1)^(t+trianglet))/(trianglex) + bar(rho)_i^(t+trianglet)(bar(V)_(i)^(t+trianglet) - bar(V)_(i -1)^(t+trianglet))/(trianglex) + bar(rho)_i^(t+trianglet)bar(V)_i^(t+trianglet)(lnbar(A)_(i)^(t+trianglet) - lnbar(A)_(i-1)^(t+trianglet))/(trianglex)]`

Similarly,

`((del(V))/(delt))_i^c, ((del(e))/(delt))_i^c,((del(T))/(delt))_i^c "and"((del(P))/(delt))_i^c` is calculated

iii. Calculating Average Time Derivative: -

For each variable, the average of Predictor time Derivative and Corrector Time Derivative is Calculated

`((del(rho))/(delt))_i^(avg) = 1/2(((del(rho))/(delt))_i^p + ((del(rho))/(delt))_i^c)`

Similarly,

`((del(V))/(delt))_i^(avg),((del(e))/(delt))_i^(avg),((del(T))/(delt))_i^(avg) "and"((del(P))/(delt))_i^(avg)` are calculated

iv. Calculating Next time Step Variable: -

`rho^(t+trianglet) = rho_i ^t + ((del(rho))/(delt))_i^(avg) trianglet`

Similarly,

`V^(t+trianglet),e^(t+trianglet),T^(t+trianglet)"and "P^(t+trianglet)` are calculated

B. Adaptive Time Step Sizing: -

As in any Explicit Method, there is a stability Requirement, and for that Courant Number and for stability Courant Number must be less than or equal to 1. For the present study Courant Number 0.6 is used

`C = (a + V)/((trianglex)/(trianglet))`

`=> trianglet = C(trianglex)/(a + V)`

`=> trianglet_i = C(trianglex_i)/(a_i + V_i)`

Now,

For each time step, the minimum value of the `trianglet` should be chosen for stability

`=> trianglet = min(trianglet_1,trianglet_2,.......,trianglet_N)`

C. Boundary Conditions: -

At the inflow boundary, the Value of Velocity can't be taken as Zero(if Velocity will be zero then there will be no mass flow) but the values of Density and Temperature can be taken as the same as reservoir Values

Fig. 4 Diagrammatic Representation Characteristic Lines at Boundary Nodes(This image is taken from Computational Fluid Dynamics, author John D. Anderson, Jr., (1995) for study purpose only and with no intention to infringe the Copyright of the original Author and/or the Publication)

As on Node 1, the left running Characteristic is Subsonic and the Right Running Characteristic is Sonic/Supersonic. So the Velocity at Node 1 can be taken as Floating and determined using Node 2 and Node 3

`(V_2 - V_1)/(trianglex) = (V_3 - V_2)/(trianglex)`

`=> V_1 = 2V_2 - V_3`

As on Node N, the left and Right Running Characteristic is Supersonic. So the Density, Velocity and Temperature at Node N can be taken as Floating and determined using Node N-1 and N-2

`(rho_N -rho_(N-1))/(trianglex) = (rho_(N-1) -rho_(N-1=2))/(trianglex)`

`rho_N = 2rho_(N-1) - rho_(N-2)`

Similarly,

`V_N = 2V_(N-1) - V_(N-2)`

`T_N = 2T_(N-1) - T_(N-2)`

The Boundary Conditions are as follows: -

D. Using Solution Vectors, Flux Vectors, and Source Terms to Solve Conservational Form: -

In our Present Equations, the Conservational Forms of Governing Equations can be converted into `U`(Solution Vector), `F`(Flux Vector), and `J`(Source Terms). The Non-Dimensional Conservative Form of Governing Equations are converted into `U`, `F`, and `J` Forms

Non-Dimensional Conservative Forms Equations: -

Continuity Form `-> (del(rho^'A^'))/(delt^') + (del(rho^'A^'V^'))/(delx^') = 0`

Momentum Equation `-> (del(rho^'A^'V^'))/(delt^') + (del(rho^'A^'(V^')^2 + (1/gamma)P^'A^'))/(delx^') = (1/gamma)P^'(delA^')/(delx^')`

Energy Conservation Equation `-> (del(rho^'A^'((e^')/(gamma - 1) + (gamma/2)(V^')^2)))/(delt^') + (del(rho^'A^'V^'((e^')/(gamma - 1) + (gamma/2)(V^')^2) + P^'A^'V^'))/(delx^') = 0`

These Forms can be written as follows: -

Continuity Form `-> (delU_1)/(delt) = - (delF_1)/(delx)`

Momentum Equation `-> (delU_2)/(delt) = - (delF_2)/(delx) + J_2`

Energy Conservation Equation `-> (delU_3)/(delt) = - (delF_3)/(delx)`

From Comparision between both Forms we can write: -

`U_1 = rho^'A^'`

`U_2 = rho^'A^'V^'`

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

`F_1 = rho^'A^'V^'`

`F_2 = rho^'A^'(V^')^2 + (1/gamma)P^'A^'`

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

`J_2 = (1/gamma)P^'(delA^')/(delx^')`

5. Method and Explanation of Code: -

The Total solution has 3 code Blocks, two are functions and one is the main script. The two Functions are the Two forms(Conservative and Non-Conservative Forms) of Governing Equations.

A. Explanation of Conservative form Solution Function: -

This Function is Coded to Give the Output Matrix of `rho^'`, `V^'`, `T^'`, and `dot(m)^'` by inputting the Length of Nozzle, Number of Grid Points, Total time Steps, and Courant Number

i. Section A(Defining Variables for Computation): -

`->` Creating `dx`, x Vector, and Area Vector using the input Variables `L`, `N`

`->` Inputting Initial Conditions of Density, Velocity, and Temperature

`->` Converting the initial Conditions in U1, U2, and U3 forms as Explained in section 4.D.

i. Section B(Computation using MacCormack's Technique): -

`->` Creating a time loop from 1 to time Steps as per input.

`->` Storing Values of U1, U2, and U3 in U1_old, U2_old, and U3_old respectively

`->` Calculating the Values of `trianglet` using the method explained in section 4.B.

`->` Calculating Predictor Time Derivatives as explained in Section 4.A.i.a.

`->` Calculating Predictor Variable Value as explained in Section 4.A.i.b.

`->` Calculating Corrector Time Derivative as explained in Section 4.A.ii.a.

`->` Calculating Average Time Derivative as explained in Section 4.A.iii.

`->` Calculating Next time Step Variable as explained in Section 4.A.iv.

`->` Applying Boundary Conditions as explained in Section 4.C.

`->` Storing Value of `rho`,`V`,`T`,and `dot(m)` in OUTPUT Matrix

B. Explanation of Non-Conservative form Solution Function: -

This Function is Coded to Give the Output Matrix of `rho^'`, `V^'`, `T^'`, and `dot(m)^'` by inputting the Length of Nozzle, Number of Grid Points, Total time Steps, and Courant Number

i. Section A(Defining Variables for Computation): -

`->` Creating `dx`, x Vector, and Area Vector using the input Variables `L`, `N`

`->` Inputting Initial Conditions of Density, Velocity, and Temperature

ii. Section B(Computation using MacCormack's Technique): -

`->` Creating a time loop from 1 to time Steps as per input.

`->` Storing Values of `rho`, `V`, `T`, in `rho_(old)`, `V_(old)`, and `T_(old)` respectively

`->` Calculating the Values of `trianglet` using the method explained in section 4.B.

`->` Calculating Predictor Time Derivatives as explained in Section 4.A.i.a.

`->` Calculating Predictor Variable Value as explained in Section 4.A.i.b.

`->` Calculating Corrector Time Derivative as explained in Section 4.A.ii.a.

`->` Calculating Average Time Derivative as explained in Section 4.A.iii.

`->` Calculating Next time Step Variable as explained in Section 4.A.iv.

`->` Applying Boundary Conditions as explained in Section 4.C.

`->` Storing Value of `rho`,`V`,`T`,and `dot(m)` in OUTPUT Matrix

C. Explanation of Main Script: -

This Script is used to Slove the Quasi 1D Subsonic-Supersonic Nozzle flow by calling both the Functions

i. Section A(Defining Variables): -

`->` Defined Values of L, N, total time steps, `gamma`, and Courant Number

`->` Creating `dx`, x Vector, and Area Vector using the input Variables `L`, `N`

`->` Creating `m_("Analytical")`

function [rho_sol,V_sol,T_sol,m_sol] = W7_Assign1_Conser_fun(xbyL,N,time_steps,gam,C)

%% Defining Variables for Computations
 % Creating Variables
   dx = xbyL/(N-1);
   x = [0:dx:xbyL];
   A = 1 + 2.2*(x - 1.5).^2;

 % Initail Values
  % Density
    rho(1:floor((N-1)/6)) = 1;
    rho(floor((N-1)/6)+1:floor((N-1)/2)) = 1 - 0.366*(x(floor((N-1)/6)+1:floor((N-1)/2)) - 0.5);
    rho(floor((N-1)/2)+1:N) = 0.634 - 0.3879*(x(floor((N-1)/2)+1:N) - 1.5);
  % Velocity(Assuming mass flow rate = 0.59)
    V = 0.59./(rho.*A);
  % Temperature
    T(1:floor((N-1)/6)) = 1;
    T(floor((N-1)/6)+1:floor((N-1)/2)) = 1 - 0.167*(x(floor((N-1)/6)+1:floor((N-1)/2)) - 0.5);
    T(floor((N-1)/2)+1:N) = 0.833 - 0.3507*(x(floor((N-1)/2)+1:N) - 1.5);

 % Converting Variables in U form
   U1 = rho.*A;
   U2 = rho.*A.*V;
   U3 = rho.*(T/(gam -1) + (gam/2)*(V.^2)).*A;

%% Computation using MacCormack's Method
 % Time Loop
   for t = 1:time_steps
       % Storing Updated Variables
         U1_old = U1;
         U2_old = U2;
         U3_old = U3;

       % Calculating next dt
         for i = 1:N
             del_t(i) = C*(dx/(V(i) + sqrt(T(i))));
         end
         dt = min(del_t);

       % Predictor Method
        % Updating Variables F as per U
          F1 = U2;
          F2 = (U2.^2)./U1 + ((gam - 1)/gam)*(U3 - (gam/2)*((U2.^2)./U1));
          F3 = gam*(U2.*U3)./U1 - (gam*(gam - 1)/2)*(U2.^3)./(U1.^2);

        for i = 2:(N-1)
            % Calculating Space Derivatives
              dF1_dx = (F1(i+1) - F1(i))/dx;
              dF2_dx = (F2(i+1) - F2(i))/dx;
              dF3_dx = (F3(i+1) - F3(i))/dx;
              % Density and Temperature for Calculating J2
                rho_i = U1(i)/A(i);
                T_i = (gam - 1)*(U3(i)/U1(i) - (gam/2)*(U2(i)/U1(i))^2);
                J2 = (1/gam)*rho_i*T_i*((A(i+1) - A(i))/dx);

            % Calculating Time Derivatives
              dU1_dt_p(i) = -dF1_dx;
              dU2_dt_p(i) = -dF2_dx + J2;
              dU3_dt_p(i) = -dF3_dx;

            % Calculating Predictor Variables
              U1(i) = U1_old(i) + dU1_dt_p(i)*dt;
              U2(i) = U2_old(i) + dU2_dt_p(i)*dt;
              U3(i) = U3_old(i) + dU3_dt_p(i)*dt;
        end
       
       % Corrector Method
        % Updating Variables F as per U
          F1 = U2;
          F2 = (U2.^2)./U1 + ((gam - 1)/gam)*(U3 - (gam/2)*((U2.^2)./U1));
          F3 = gam*(U2.*U3)./U1 - (gam*(gam - 1)/2)*(U2.^3)./(U1.^2);

        for i = 2:(N-1)
            % Calculating Space Derivatives
              dF1_dx = (F1(i) - F1(i-1))/dx;
              dF2_dx = (F2(i) - F2(i-1))/dx;
              dF3_dx = (F3(i) - F3(i-1))/dx;
              % Density and Temperature for Calculating J2
                rho_i = U1(i)/A(i);
                T_i = (gam - 1)*(U3(i)/U1(i) - (gam/2)*(U2(i)/U1(i))^2);
                J2 = (1/gam)*rho_i*T_i*((A(i) - A(i-1))/dx);

            % Calculating Time Derivatives
              dU1_dt_c(i) = -dF1_dx;
              dU2_dt_c(i) = -dF2_dx + J2;
              dU3_dt_c(i) = -dF3_dx;
        end

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

       % Next Time Step Variables
         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
         U2(1) = 2*U2(2) - U2(3);
         V1 = U2(1)/U1(1);
         U3(1) = U1(1)*(T(1)/(gam-1) + gam/2*(V1^2));
         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);

       % Calculating Density, Velocity and Temperature from U
         rho = U1./A;
         V = U2./U1;
         T = (gam -1)*((U3./U1) - (gam/2)*(V.^2));

       % Storing Solutions in rho, V and T Solution Vectors
         rho_sol(t,:) = U1./A;
         V_sol(t,:) = U2./U1;
         T_sol(t,:) = (gam -1)*((U3./U1) - (gam/2)*(V.^2));
         m_sol(t,:) = rho.*V.*A;
   end
end

function [rho_sol,V_sol,T_sol,m_sol] = W7_Assign1_NonConser_fun(xbyL,N,time_steps,gam,C)

%% Defining Variables for Computations
 % Creating Variables
   dx = xbyL/(N-1);
   x = [0:dx:xbyL];
   A = 1 + 2.2*(x - 1.5).^2;

 % Initial Value
   rho = 1 - 0.3146*x;
   T = 1 - 0.2314*x;
   V = (0.1 + 1.09*x).*(T.^(0.5));
   lnA = log(A);

%% Computation using MacCormack's Method
 % Time Loop
   for t = 1:time_steps
       % Storing updated Variables
         rho_old = rho;
         V_old = V;
         T_old = T;

       % Calculating next dt
         for i = 1:N
             del_t(i) = C*(dx/(V(i) + sqrt(T(i))));
         end
         dt = min(del_t);

       % Predictor Method
         for i = 2:(N-1)
             % Calculating Space Derivative
               drho_dx = (rho_old(i+1) - rho_old(i))/dx;
               dV_dx = (V_old(i+1) - V_old(i))/dx;
               dT_dx = (T_old(i+1) - T_old(i))/dx;
               dlnA_dx = (lnA(i+1) - lnA(i))/dx;

             % Calculating Time Derivative
               drho_dt_p(i) = -rho_old(i)*dV_dx - rho_old(i)*V_old(i)*dlnA_dx - V_old(i)*drho_dx;
               dV_dt_p(i) = -V_old(i)*dV_dx - (1/gam)*(dT_dx + T_old(i)*drho_dx/rho_old(i));
               dT_dt_p(i) = -V_old(i)*dT_dx - T_old(i)*(gam - 1)*(dV_dx + V_old(i)*dlnA_dx);

             % Calculating Predictor Variables
               rho(i) = rho_old(i) + drho_dt_p(i)*dt;
               V(i) = V_old(i) + dV_dt_p(i)*dt;
               T(i) = T_old(i) + dT_dt_p(i)*dt;

         end

       % Corrector Method
         for i = 2:(N-1)
             % Calculating Space Derivative
               drho_dx = (rho(i) - rho(i-1))/dx;
               dV_dx = (V(i) - V(i-1))/dx;
               dT_dx = (T(i) - T(i-1))/dx;
               dlnA_dx = (lnA(i) - lnA(i-1))/dx;

             % Calculating Time Derivative
               drho_dt_c(i) = -rho(i)*dV_dx - rho(i)*V(i)*dlnA_dx - V(i)*drho_dx;
               dV_dt_c(i) = -V(i)*dV_dx - (1/gam)*(dT_dx + T(i)*drho_dx/rho(i));
               dT_dt_c(i) = -V(i)*dT_dx - T(i)*(gam - 1)*(dV_dx + V(i)*dlnA_dx);

         end

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

       % Next Time Step Variable
         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 Condition
         V(1) = 2*V(2) - V(3);
         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);

       % Storing Solutions in rho, V and T Solution Vectors
         rho_sol(t,:) = rho;
         V_sol(t,:) = V;
         T_sol(t,:) = T;
         m_sol(t,:) = rho.*V.*A;
   end
end

% To Solve Quasi 1D Sub-Supersonic Nozzle Flow using MacCormack's Technique
  clear all
  close all
  clc

%% Defining Variables
   xbyL = 3;
   N = 31;
   time_steps = 1400;
   gam = 1.4;
   C = 0.6;

 % Creating Variables
   dx = xbyL/(N-1);
   x = [0:dx:xbyL];
   m_Analytical = 0.579*ones(1,N);

%% Calling Functions for Solutions
   [rho_con,V_con,T_con,m_con] = W7_Assign1_Conser_fun(xbyL,N,time_steps,gam,C);
   [rho_noncon,V_noncon,T_noncon,m_noncon] = W7_Assign1_NonConser_fun(xbyL,N,time_steps,gam,C);

%% Calculating Variables and Plotting
  % Calculating Conervative Pressure and Mach Number
    P_con = rho_con.*T_con;
    mach_con = V_con./(T_con.^0.5);

  % Calculating Non-Conservative Pressure and Mach Number
    P_noncon = rho_noncon.*T_noncon;
    mach_noncon = V_noncon./(T_noncon.^0.5);

%% Animating Plot from 1st Time Step to 1400 Time Step
   figure(1)
    for i = 1:time_steps
        % Conservative Form Solution
             subplot(5,2,1)
             plot(x,rho_con(i,:))
             axis([0 3 0 1.4])
             title({['Conservative ','Form'];['Time Step = ',num2str(i)]})
             xlabel('Nozzle Length Ratio')
             ylabel('Density')
             subplot(5,2,3)
             plot(x,V_con(i,:))
             axis([0 3 0 2.2])
             xlabel('Nozzle Length Ratio')
             ylabel('Velocity')
             subplot(5,2,5)
             plot(x,T_con(i,:))
             axis([0 3 0 1.2])
             xlabel('Nozzle Length Ratio')
             ylabel('Temperature')
             subplot(5,2,7)
             plot(x,P_con(i,:))
             axis([0 3 0 1.2])
             xlabel('Nozzle Length Ratio')
             ylabel('Pressure')
             subplot(5,2,9)
             plot(x,mach_con(i,:))
             axis([0 3 0 3.5])
             xlabel('Nozzle Length Ratio')
             ylabel('Mach Number')

        % Non-Conservative Form Solution
             subplot(5,2,2)
             plot(x,rho_noncon(i,:))
             axis([0 3 0 1.4])
             title({['Non-Conservative ','Form'];['Time Step = ',num2str(i)]})
             xlabel('Nozzle Length Ratio')
             ylabel('Density')
             subplot(5,2,4)
             plot(x,V_noncon(i,:))
             axis([0 3 0 2.2])
             xlabel('Nozzle Length Ratio')
             ylabel('Velocity')
             subplot(5,2,6)
             plot(x,T_noncon(i,:))
             axis([0 3 0 1.2])
             xlabel('Nozzle Length Ratio')
             ylabel('Temperature')
             subplot(5,2,8)
             plot(x,P_noncon(i,:))
             axis([0 3 0 1.5])
             xlabel('Nozzle Length Ratio')
             ylabel('Pressure')
             subplot(5,2,10)
             plot(x,mach_noncon(i,:))
             axis([0 3 0 3.5])
             xlabel('Nozzle Length Ratio')
             ylabel('Mach Number')

           % Adding Time Delay
             pause(0.005)
    end

%% Plotting mass Flow rate for both forms
 % Conservative Form Mass Distribution at diffrent time Steps
   figure(2)
   plot(x,m_con(1,:),x,m_con(50,:),x,m_con(100,:),x,m_con(150,:),x,m_con(200,:),x,m_con(700,:),x,m_Analytical,'-bs')
   xlabel('Length of the Nozzle')
   ylabel('Non-Dimensional Mass Distribution')
   legend('0*dt','50*dt','100*dt','150*dt','200*dt','700*dt','Analytical Non-Dimensional Mass Flow Rate')
   title('Non-Dimensional Mass Distribution at Diffrent Time Step throughout the nozzle for Conservative Form')

 % Conservative Form Mass Distribution at diffrent time Steps
   figure(3)
   plot(x,m_noncon(1,:),x,m_noncon(50,:),x,m_noncon(100,:),x,m_noncon(150,:),x,m_noncon(200,:),x,m_noncon(700,:),x,m_Analytical,'-bs')
   xlabel('Length of the Nozzle')
   ylabel('Non-Dimensional Mass Distribution')
   legend('0*dt','50*dt','100*dt','150*dt','200*dt','700*dt','Analytical Non-Dimensional Mass Flow Rate')
   title('Non-Dimensional Mass Distribution at Diffrent Time Step throughout the nozzle for Non-Conservative Form')

 % Steady-State Mass Distribution for both Forms and Analytical Solution
   figure(4)
   plot(x,m_con(1400,:),x,m_noncon(1400,:),x,m_Analytical)
   xlabel('Length of the Nozzle')
   ylabel('Non-Dimensional Mass Distribution')
   legend('Conservative Mass Distribution','Non-Conservative Mass Distribution','Analytical Mass Distribution')
   title('Steady-State Mass Distribution')