Courses by Software
Courses by Semester
Courses by Domain
Tool-focused Courses
Machine learning
POPULAR COURSES
Post Graduate Program in Hybrid Electric Vehicle Design and AnalysisPost Graduate Program in Computational Fluid DynamicsPost Graduate Program in CADPost Graduate Program in CAEPost Graduate Program in Manufacturing DesignPost Graduate Program in Computational Design and Pre-processingPost Graduate Program in Complete Passenger Car Design & Product DevelopmentSuccess Stories
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…
Aditya Aanand
updated on 16 Oct 2022
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 →ρ
ii. Velocity(Along x) →V
iii. Temperature →T
iv. Pressure →P
v. Internal Energy →e
vi. Mass Flow rate →.m
vii. Specific Heat Ratio →γ
viii. Specific Heat at Constant Volume →cv
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)=cvT
ii. Speed of Sound(a)=√γRT
iii. P = rhoRT`
iv. cv=Rγ-1
C. Convention Used to denote Different Regions and Type of Variables(Let any Property be ϕ)
i. To Denote the Property at Reservior →ϕo
ii. To Denote the Property at Throat →ϕ*
iii. To Denote the Property at Exit →ϕe
iv. To Denote Non-Dimensional Property →ϕ′
D. Subsonic-Supersonic Nozzle with AA*=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. ρ′=1-0.3146⋅x′{x′∈[0,3]}
b. T′=1-0.2314⋅x′{x′∈[0,3]}
c. V′=(0.1+1.09⋅x′)(T′)12{x′∈[0,3]}
ii. For Conservative Form: -
a. ρ′={1x∈[0,0.5)1-0.366(x′-0.5)x∈[0.5,1.5]0.634-0.3879(x′-1.5)x∈(1.5,3]
b. T′={1x∈[0,0.5)1-0.167(x′-0.5)x∈[0.5,1.5]0.833-0.3507(x′-1.5)x∈(1.5,3]
c. V′=.m′ρ′A′
d. .m′=0.590
F. The Analytical mass Flow Rate.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(.mLeft)=ρAV
Mass Leaving from the Right Side(.mRight)=ρAV+∂(ρAV)∂x
Hence, Net Mass decrease due to mass leaving(.mNet)=(.mRight)-(.mLeft)
⇒.mNet=ρAV+∂(ρAV)∂x-ρAV
⇒.mNet=∂(ρAV)∂x
b. The Net Change of mass change within the Control Volume: -
Net Mass decrease within the control Volume(.mtime)=-∂(ρA)∂t
As the Net decrease in mass within the Control Volume is happening due to the Net mass leaving
Hence,
.mNet=.mtime
⇒∂(ρAV)∂x=-∂(ρA)∂t
⇒∂(ρA)∂t+∂(ρAV)∂x=0
| Continuity Equation in Conservation Form | ∂(ρA)∂t+∂(ρAV)∂x=0 (eq.1) |
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
Fnet=PA+2(PdA2)-(P+dP)(A+dA)
⇒Fnet=-(AdP+dPdA)
As dPdA is a very small Quantity and so it can be neglected
⇒Fnet=-AdP (eq.2)
Now,
As d(PA)=PdA+AdP
⇒AdP=d(PA)-PdA (eq.3)
Now,
Putting (eq.3) Values in (eq.2)
⇒Fnet=-(d(PA)-PdA) (eq.4)
b. The Effect of Net Force Applied on the Control Volume in the x direction is given by max
Effect of Net Force = max
As the Control Volume is moving with the fluid,
⇒max=ρAdxDVDt
⇒max=(ρAdx)(∂V∂t+V∂V∂x)
⇒max=(ρA∂V∂t+ρAV∂V∂x)dx (eq.5)
Now,
As ∂(ρAV)∂t=V∂(ρA)∂t+ρA∂(V)∂t
⇒ρA∂(V)∂t=∂(ρAV)∂t-V∂(ρA)∂t (eq.6)
And
As ∂(ρAV2)∂x=V∂(ρAV)∂x+ρAV∂(V)∂x
⇒ρAV∂(V)∂x=∂(ρAV2)∂x-V∂(ρAV)∂x (eq.7)
Now,
Putting Values of (eq.6) and (eq.7) in (eq.5)
max=(∂(ρAV)∂t-V∂(ρA)∂t+∂(ρAV2)∂x-V∂(ρAV)∂x)dx
⇒max=(∂(ρAV)∂t+∂(ρAV2)∂x-V∂(ρA)∂t-V∂(ρAV)∂x)dx
⇒max=(∂(ρAV)∂t+∂(ρAV2)∂x-V(∂(ρA)∂t+∂(ρAV)∂x))dx
Now,
As from (eq.1) the value of ∂(ρA)∂t+∂(ρAV)∂x=0
⇒max=(∂(ρAV)∂t+∂(ρAV2)∂x)dx (eq.8)
As the Net Force acting on the System is causing the net Effect on the System
So,
max=Fnet
⇒(∂(ρAV)∂t+∂(ρAV2)∂x)dx=-(d(PA)-PdA)
Dividing by dx
⇒∂(ρAV)∂t+∂(ρAV2)∂x=-∂(PA)∂x+P∂A∂x
⇒∂(ρAV)∂t+∂(ρAV2+PA)∂x=P∂A∂x
| Momentum Equation in Conservative Form | ∂(ρAV)∂t+∂(ρAV2+PA)∂x=P∂A∂x (eq.9) |
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: -
.Eexternal=PAV+2(PVdA2)-(P+dP)(A+dA)(V+dV)
⇒.Eexternal=PAV+2(PVdA2)-(PAV+AVdP+PVdA+VdPdA+PAdV+AdPdV+PdAdV+dPdAdV)
As dPdAdV is very small quantities and so it can be neglected
⇒.Eexternal=-AVdP-VdPdA-PAdV-PdAdV-AdPdV
⇒.Eexternal=-AVdP -PAdV-(dA)VdP-(dA)PdV-AdPdV
⇒.Eexternal=-d(PAV) (eq.10)
b. The Net Rate of increase of energy of the control fluid moving with the flow: -
.Einternal=mDeDt⏟Internal Energy+mDV22Dt⏟Kinetic Energy
⇒.Einternal=ρAdx(D(e+V22)Dt)
⇒.Einternal=ρA(∂(e+V22)∂t+V∂(e+V22)∂x)dx
⇒.Einternal=(ρA∂(e+V22)∂t+ρAV∂(e+V22)∂x)dx (eq.11)
Now,
As ∂(ρA(e+V22))∂t=ρA∂(e+V22)∂t+(e+V22)∂(ρA)∂t
⇒ρA∂(e+V22)∂t=∂(ρA(e+V22))∂t-(e+V22)∂(ρA)∂t (eq.12)
And,
As ∂(ρAV(e+V22))∂x=ρAV∂(e+V22)∂x+(e+V22)∂(ρAV)∂x
⇒ρAV∂(e+V22)∂x=∂(ρAV(e+V22))∂x-(e+V22)∂(ρAV)∂x (eq.13)
Now,
Putting (eq.12) and (eq.13) in (eq.11)
.Einternal=(∂(ρA(e+V22))∂t-(e+V22)∂(ρA)∂t+∂(ρAV(e+V22))∂x-(e+V22)∂(ρAV)∂x)dx
⇒.Einternal=(∂(ρA(e+V22))∂t+∂(ρAV(e+V22))∂x-(e+V22)(∂(ρA)∂t+∂(ρAV)∂x))dx
Now,
As from (eq.1) the Value of ∂(ρA)∂t+∂(ρAV)∂x=0
⇒.Einternal=(∂(ρA(e+V22))∂t+∂(ρAV(e+V22))∂x)dx
As the Internal increase of total Energy of the Control Volume is happening due to the External Energy
So,
.Einternal=.Eexternal
⇒(∂(ρA(e+V22))∂t+∂(ρAV(e+V22))∂x)dx=-d(PAV)
Dividing by dx
⇒∂(ρA(e+V22))∂t+∂(ρAV(e+V22))∂x=-∂(PAV)∂x
⇒∂(ρA(e+V22))∂t+∂(ρAV(e+V22)+PAV)∂x=0
| Energy Conservation in Conservative Form | ∂(ρA(e+V22))∂t+∂(ρAV(e+V22)+PAV)∂x=0 (eq.14) |
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 ρ,V,T,P within the nozzle. The following are the Governing Equations for Subsonic-Supersonic Quasi 1D Nozzle in Conservative form: -
| Continuity Equation | ∂(ρA)∂t+∂(ρAV)∂x=0 |
| Momentum Equation | ∂(ρAV)∂t+∂(ρAV2+PA)∂x=P∂A∂x |
| Energy Conservation Equation | ∂(ρA(e+V22))∂t+∂(ρAV(e+V22)+PAV)∂x=0 |
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)
∂(ρA)∂t+∂(ρAV)∂x=0
⇒∂(ρA)∂t=-∂(ρAV)∂x
⇒∂(ρA)∂t=-[ρA∂(V)∂x+V∂(ρA)∂x]
⇒∂(ρA)∂t=-[ρA∂(V)∂x+ρV∂(A)∂x+AV∂(ρ)∂x]
As the Area of the Cross-Section of Nozzle doesn't change with time (⇒A=fn(x)only)
Hence,
⇒A∂(ρ)∂t=-[ρA∂(V)∂x+ρV∂(A)∂x+AV∂(ρ)∂x]
Dividing both sides with A
⇒∂(ρ)∂t=-[ρ∂(V)∂x+ρVA∂(A)∂x+V∂(ρ)∂x] (eq.15)
And,
As d(lnA)=dAA (eq.16)
Now,
Putting (eq.16) in (eq.15)
⇒∂(ρ)∂t=-[V∂(ρ)∂x+ρ∂(V)∂x+ρV∂(lnA)∂x]
| Continuity Equation in Non-Conservative Form | ∂(ρ)∂t=-[V∂(ρ)∂x+ρ∂(V)∂x+ρV∂(lnA)∂x] (eq.17) |
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)
∂(ρAV)∂t+∂(ρAV2+PA)∂x=P∂A∂x
⇒∂(ρAV)∂t+∂(ρAV2)∂x=-∂(PA)∂x+P∂A∂x (eq.18)
Now,
As ∂(ρAV)∂t=ρA∂(V)∂t+V∂(ρA)∂t (eq.19)
And,
As ∂(ρAV2)∂x=ρAV∂(V)∂x+V∂(ρAV)∂x (eq.20)
Now,
Putting (eq.19) and (eq.20) in (eq.18)
⇒ ρA∂(V)∂t+V∂(ρA)∂t+ ρAV∂(V)∂x+V∂(ρAV)∂x=-∂(PA)∂x+P∂A∂x
⇒ ρA∂(V)∂t+ ρAV∂(V)∂x+V[∂(ρA)∂t+∂(ρAV)∂x]=-[∂(PA)∂x-P∂A∂x]
As from (eq.1) the Value of ∂(ρA)∂t+∂(ρAV)∂x=0
⇒ ρA∂(V)∂t+ ρAV∂(V)∂x=-[∂(PA)∂x-P∂A∂x] (eq.21)
Now,
As ∂(PA)∂x=A∂P∂x+P∂A∂x
⇒∂(PA)∂x-P∂A∂x=A∂P∂x (eq.22)
Now,
Putting the Value of (eq.22) in (eq.21)
⇒ρA∂(V)∂t+ ρAV∂(V)∂x=-A∂P∂x (eq.23)
Dividing Both Sides with A
⇒ρ∂(V)∂t+ ρV∂(V)∂x=-∂P∂x
⇒ρ∂(V)∂t=- [ρV∂(V)∂x+∂P∂x]
| Momentum Equation in Non-Conservative Form | ρ∂(V)∂t=- [ρV∂(V)∂x+∂P∂x] (eq.24) |
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)
∂(ρA(e+V22))∂t+∂(ρAV(e+V22)+PAV)∂x=0
⇒∂(ρA(e+V22))∂t+∂(ρAV(e+V22))∂x=-∂(PAV)∂x
⇒∂(ρAe)∂t+∂(ρAV22)∂t+∂(ρAVe)∂x+∂(ρAV(V22))∂x=-∂(PAV)∂x
⇒∂(ρAe)∂t+∂(ρAVe)∂x+∂(ρA(V22))∂t+∂(ρAV(V22))∂x=-∂(PAV)∂x (eq.25)
Now,
As ∂(ρAe)∂t=ρA∂(e)∂t+e∂(ρA)∂t (eq.26)
And,
As ∂(ρAVe)∂x=ρAV∂(e)∂x+e∂(ρAV)∂x (eq.27)
Now,
Putting Values of (eq.26) and (eq.27) in (eq.25)
⇒ρA∂(e)∂t+e∂(ρA)∂t+ρAV∂(e)∂x+e∂(ρAV)∂x+∂(ρA(V22))∂t+∂(ρAV(V22))∂x=-∂(PAV)∂x
⇒ρA∂(e)∂t+ρAV∂(e)∂x+e[∂(ρA)∂t+∂(ρAV)∂x]+∂(ρA(V22))∂t+∂(ρAV(V22))∂x=-∂(PAV)∂x
Now,
As from (eq.1) the Value of ∂(ρA)∂t+∂(ρAV)∂x=0
⇒ρA∂(e)∂t+ρAV∂(e)∂x+∂(ρA(V22))∂t+∂(ρAV(V22))∂x=-∂(PAV)∂x (eq.28)
Now,
As ∂(ρA(V22))∂t=ρA∂(V22)∂t+(V22)∂(ρA)∂t (eq.29)
And,
As ∂(ρAV(V22))∂x=ρAV∂(V22)∂x+(V22)∂(ρAV)∂x (eq.30)
Now,
Putting the Values of (eq.29) and (eq.30) in (eq.28)
⇒ρA∂(e)∂t+ρAV∂(e)∂x+ρA∂(V22)∂t+(V22)∂(ρA)∂t+ρAV∂(V22)∂x+(V22)∂(ρAV)∂x=-∂(PAV)∂x
⇒ρA∂(e)∂t+ρAV∂(e)∂x+ρA∂(V22)∂t+ρAV∂(V22)∂x+(V22)[∂(ρA)∂t+∂(ρAV)∂x]=-∂(PAV)∂x
Now,
As from (eq.1) the Value of ∂(ρA)∂t+∂(ρAV)∂x=0
⇒ρA∂(e)∂t+ρAV∂(e)∂x+ρA∂(V22)∂t+ρAV∂(V22)∂x=-∂(PAV)∂x
⇒ρA∂(e)∂t+ρAV∂(e)∂x+ρAV∂(V)∂t+ρAV2∂(V)∂x=-∂(PAV)∂x
⇒ρA∂(e)∂t+ρAV∂(e)∂x+V[ρA∂(V)∂t+ρAV∂(V)∂x]=-∂(PAV)∂x
Now,
As from (eq.23) the Value of ρA∂(V)∂t+ ρAV∂(V)∂x=-A∂P∂x
⇒ρA∂(e)∂t+ρAV∂(e)∂x+V[-A∂P∂x]=-∂(PAV)∂x
⇒ρA∂(e)∂t+ρAV∂(e)∂x-AV∂P∂x=-∂(PAV)∂x
⇒ρA∂(e)∂t+ρAV∂(e)∂x=-[∂(PAV)∂x-AV∂P∂x] (eq.31)
Now,
As ∂(PAV)∂x=AV∂(P)∂x+PA∂(V)∂x+PV∂(A)∂x
⇒∂(PAV)∂x-AV∂(P)∂x=PA∂(V)∂x+PV∂(A)∂x (eq.32)
Now,
Putting Values of (eq.32) in (eq.31)
⇒ρA∂(e)∂t+ρAV∂(e)∂x=-[PA∂(V)∂x+PV∂(A)∂x]
Now,
Dividing A on both Sides
⇒ρ∂(e)∂t+ρV∂(e)∂x=-[P∂(V)∂x+PVA∂(A)∂x]
Now,
As from (eq.16) value of dAA=d(lnA)
⇒ρ∂(e)∂t+ρV∂(e)∂x=-P∂(V)∂x-PV∂(lnA)∂x
⇒ρ∂(e)∂t=-[ρV∂(e)∂x+P∂(V)∂x+PV∂(lnA)∂x]
| Energy Conservation Equation in Non-Conservation Form | ρ∂(e)∂t=-[ρV∂(e)∂x+P∂(V)∂x+PV∂(lnA)∂x] (eq.33) |
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 ρ,V,T,P within the nozzle. The following are the Governing Equations for Subsonic-Supersonic Quasi 1D Nozzle in Non-Conservative form: -
| Continuity Equation | ∂(ρ)∂t=-[V∂(ρ)∂x+ρ∂(V)∂x+ρV∂(lnA)∂x] |
| Momentum Equation | ρ∂(V)∂t=- [ρV∂(V)∂x+∂P∂x] |
| Energy Conservation Equation | ρ∂(e)∂t=-[ρV∂(e)∂x+P∂(V)∂x+PV∂(lnA)∂x] |
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. ρ′=ρρo
b. V′=Vao
c. T′=TTo
d. P′=PPo
e. A′=AA*
f. e′=eeo
g. m′=mmo
h. t′=tLao
i. x′=xL
i. Non-Dimensionalizing Conservative Form: -
a. Continuity Equation: -
∂(ρA)∂t+∂(ρAV)∂x=0
Applying the Non-Dimensional terms in the above Equation
⇒∂(ρ′A′)∂t′(ρoA*Lao)+∂(ρ′A′V′)∂x′(ρoA*aoL)=0
Now,
Dividing both sides by (ρoA*aoL)
⇒∂(ρ′A′)∂t′+∂(ρ′A′V′)∂x′=0
| Non-Dimensional Continuity Equation of Conservative Form | ∂(ρ′A′)∂t′+∂(ρ′A′V′)∂x′=0 (eq.34) |
b. Momentum Equation: -
∂(ρAV)∂t+∂(ρAV2+PA)∂x=P∂A∂x
Applying the Non-Dimensional terms in the above Equation
⇒∂(ρ′A′V′)∂t′(ρoA*aoLao)+∂(ρ′A′(V′)2(ρoA*a2o)+P′A′(PoA*))∂x′(L)=P′∂A′∂x′(PoA*L)
⇒∂(ρ′A′V′)∂t′(ρoA*a2oL)+∂(ρ′A′(V′)2(ρoA*a2o)+P′A′(PoA*))∂x′(L)=P′∂A′∂x′(PoA*L)
Now,
As Po=ρoRTo and ao=√γRTo
⇒∂(ρ′A′V′)∂t′(ρoA*(γRTo)L)+∂(ρ′A′(V′)2(ρoA*(γRTo))+P′A′((ρoRTo)A*))∂x′(L)=P′∂A′∂x′((ρoRTo)A*L)
⇒∂(ρ′A′V′)∂t′(γρoA*RToL)+∂(ρ′A′(V′)2γ+P′A′)∂x′(ρoA*RToL)=P′∂A′∂x′(ρoRToA*L)
Now,
Dividing both sides by (ρoRToA*γL)
⇒∂(ρ′A′V′)∂t′+∂(ρ′A′(V′)2+(1γ)P′A′)∂x′=(1γ)P′∂A′∂x′
| Non-Dimensional Momentum Equation of Conservative Form | ∂(ρ′A′V′)∂t′+∂(ρ′A′(V′)2+(1γ)P′A′)∂x′=(1γ)P′∂A′∂x′ (eq.35) |
c. Energy Conservation Equation: -
∂(ρA(e+V22))∂t+∂(ρAV(e+V22)+PAV)∂x=0
Applying the Non-Dimensional terms in the above Equation
⇒∂(ρ′A′(e′(eo)+(V′)22(a2o)))∂t′(ρoA*Lao)+∂(ρ′A′V′(e′(eo)+(V′)22(a2o))+P′A′V′(Po))∂x′(A*aoL)=0
Now,
Dividing both sides by A*aoL
⇒∂(ρ′A′(e′(eo)+(V′)22(a2o)))∂t′(ρo)+∂(ρ′A′V′(e′(eo)+(V′)22(a2o))(ρo)+P′A′V′(Po))∂x′=0
Now,
As Po=ρoRTo, ao=√γRTo and eo=cvTo
⇒∂(ρ′A′(e′(cvTo)+(V′)22(γRTo)))∂t′(ρo)+∂(ρ′A′V′(e′(cvTo)+(V′)22(γRTo))(ρo)+P′A′V′(ρoRTo))∂x′=0
And,
As cv=Rγ-1
⇒∂(ρ′A′(e′(RToγ-1)+(V′)22(γRTo)))∂t′(ρo)+∂(ρ′A′V′(e′(RToγ-1)+(V′)22(γRTo))(ρo)+P′A′V′(ρoRTo))∂x′=0
⇒∂(ρ′A′(e′γ-1+(γ)(V′)22))∂t′(ρoRTo)+∂(ρ′A′V′(e′γ-1+(γ)(V′)22)+P′A′V′)∂x′(ρoRTo)=0
Now,
Dividing both sides by ρoRTo
⇒∂(ρ′A′(e′γ-1+(γ2)(V′)2))∂t′+∂(ρ′A′V′(e′γ-1+(γ2)(V′)2)+P′A′V′)∂x′=0
| Non-Dimensional Energy Conservation Equation in Conservative Form | ∂(ρ′A′(e′γ-1+(γ2)(V′)2))∂t′+∂(ρ′A′V′(e′γ-1+(γ2)(V′)2)+P′A′V′)∂x′=0 (eq.36) |
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: -
| Continuity Equation | ∂(ρ′A′)∂t′+∂(ρ′A′V′)∂x′=0 |
| Momentum Equation | ∂(ρ′A′V′)∂t′+∂(ρ′A′(V′)2+(1γ)P′A′)∂x′=(1γ)P′∂A′∂x′ |
| Energy Conservation Equation | ∂(ρ′A′(e′γ-1+(γ2)(V′)2))∂t′+∂(ρ′A′V′(e′γ-1+(γ2)(V′)2)+P′A′V′)∂x′=0 |
ii. Non-Dimensionalizing Non-Conservative Form: -
a. Continuity Equation: -
∂(ρ)∂t=-[V∂(ρ)∂x+ρ∂(V)∂x+ρV∂(lnA)∂x]
Applying the Non-Dimensional terms in the above Equation
⇒∂(ρ′)∂t′(ρoLao)=-[V′∂(ρ′)∂x′+ρ′∂(V′)∂x′+ρ′V′∂(ln(A′A*))∂x′](ρoLao)
⇒∂(ρ′)∂t′(ρoLao)=-[V′∂(ρ′)∂x′+ρ′∂(V′)∂x′+ρ′V′∂(ln(A′)+ln(A*))∂x′](ρoLao)
⇒∂(ρ′)∂t′(ρoLao)=-[V′∂(ρ′)∂x′+ρ′∂(V′)∂x′+ρ′V′∂(lnA′)∂x′+ρ′V′∂(lnA*)∂x′](ρoLao)
Now,
As A*=constant⇒∂(lnA*)∂x′=0
⇒∂(ρ′)∂t′(ρoLao)=-[V′∂(ρ′)∂x′+ρ′∂(V′)∂x′+ρ′V′∂(lnA′)∂x′](ρoLao)
Now,
Dividing both sides with ρoaoL
⇒∂(ρ′)∂t′=-[V′∂(ρ′)∂x′+ρ′∂(V′)∂x′+ρ′V′∂(lnA′)∂x′]
| Non-Dimensional Continuity Equation of Non-Conservative Form | ∂(ρ′)∂t′=-[V′∂(ρ′)∂x′+ρ′∂(V′)∂x′+ρ′V′∂(lnA′)∂x′] (eq.37) |
b. Momentum Equation: -
ρ∂(V)∂t=- [ρV∂(V)∂x+∂P∂x]
Applying the Non-Dimensional terms in the above Equation
⇒ρ′∂(V′)∂t′(ρoaoLao)=- [ρ′V′∂(V′)∂x′(ρoa2oL)+∂P′∂x′(PoL)]
Now,
As Po=ρoRTo and ao=√γRTo
⇒ρ′∂(V′)∂t′(ρo(γRTo)L)=- [ρ′V′∂(V′)∂x′(ρo(γRTo)L)+∂P′∂x′(ρoRToL)]
⇒ρ′∂(V′)∂t′(γρoRToL)=- [ρ′V′∂(V′)∂x′(γ)+∂P′∂x′](ρoRToL)
Now,
Dividing both sides with γρoRToL
⇒ρ′∂(V′)∂t′=- [ρ′V′∂(V′)∂x′+(1γ)∂P′∂x′]
| Non-Dimensional Momentum Equation of Non-Conservative Form | ρ′∂(V′)∂t′=- [ρ′V′∂(V′)∂x′+(1γ)∂P′∂x′] (eq.38) |
c. Energy Conservation Equation: -
ρ∂(e)∂t=-[ρV∂(e)∂x+P∂(V)∂x+PV∂(lnA)∂x]
Applying the Non-Dimensional terms in the above Equation
⇒ρ′∂(e′)∂t′(ρoeoLao)=-[ρ′V′∂(e′)∂x′(ρoaoeoL)+P′∂(V′)∂x′(PoaoL)+P′V′∂(lnA′A*)∂x′(PoaoL)]
⇒ρ′∂(e′)∂t′(ρoeoaoL)=-[ρ′V′∂(e′)∂x′(ρoaoeoL)+P′∂(V′)∂x′(PoaoL)+P′V′∂(lnA′+lnA*)∂x′(PoaoL)]
⇒ρ′∂(e′)∂t′(ρoeoaoL)=-[ρ′V′∂(e′)∂x′(ρoaoeoL)+P′∂(V′)∂x′(PoaoL)+P′V′(∂(lnA′)∂x′+∂(lnA*)∂x′)(PoaoL)]
Now,
As A*=constant⇒∂(lnA*)∂x′=0
⇒ρ′∂(e′)∂t′(ρoeoaoL)=-[ρ′V′∂(e′)∂x′(ρoaoeoL)+P′∂(V′)∂x′(PoaoL)+P′V′∂(lnA′)∂x′(PoaoL)]
Now,
As Po=ρoRTo and eo=cvTo
⇒ρ′∂(e′)∂t′(ρoaocvToL)=-[ρ′V′∂(e′)∂x′(ρoaocvToL)+P′∂(V′)∂x′(ρoRToaoL)+P′V′∂(lnA′)∂x′(ρoRToaoL)]
Now,
As cv=Rγ-1
⇒ρ′∂(e′)∂t′(ρoaoRToL(γ-1))=-[ρ′V′∂(e′)∂x′(ρoaoRToL(γ-1))+P′∂(V′)∂x′(ρoRToaoL)+P′V′∂(lnA′)∂x′(ρoRToaoL)]
Now,
Dividing both sides with (ρoaoRToL(γ-1))
⇒ρ′∂(e′)∂t′=-[ρ′V′∂(e′)∂x′+(γ-1)P′∂(V′)∂x′+(γ-1)P′V′∂(lnA′)∂x′]
| Non-Dimensional Energy Conservation Equation in Conservative Form | ρ′∂(e′)∂t′=-[ρ′V′∂(e′)∂x′+(γ-1)P′∂(V′)∂x′+(γ-1)P′V′∂(lnA′)∂x′] (eq.39) |
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: -
| Continuity Equation | ∂(ρ′)∂t′=-[V′∂(ρ′)∂x′+ρ′∂(V′)∂x′+ρ′V′∂(lnA′)∂x′] |
|
Momentum Equation |
ρ′∂(V′)∂t′=- [ρ′V′∂(V′)∂x′+(1γ)∂P′∂x′] |
| Energy Conservation Equation | ρ′∂(e′)∂t′=-[ρ′V′∂(e′)∂x′+(γ-1)P′∂(V′)∂x′+(γ-1)P′V′∂(lnA′)∂x′] |
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: -
(∂(ρ)∂t)pi=-[Vtiρti+1-ρti△x+ρtiVti+1-Vti△x+ρtiVtilnAti+1-lnAti△x]
Similarly,
(∂(V)∂t)pi and (∂(e)∂t)pi are calculated
b. Calculating Predictor Variables: -
ˉρt+△t=ρti+(∂(ρ)∂t)pi△t
Similarly,
ˉVt+△t and ˉet+△t
And using Equations to calculate ˉPt+△t and ˉTt+△t
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: -
(∂(ρ)∂t)ci=-[ˉVt+△tiˉρt+△ti-ˉρt+△ti-1△x+ˉρt+△tiˉVt+△ti-ˉVt+△ti-1△x+ˉρt+△tiˉVt+△tilnˉAt+△ti-lnˉAt+△ti-1△x]
Similarly,
(∂(V)∂t)ci,(∂(e)∂t)ci,(∂(T)∂t)ciand(∂(P)∂t)ci is calculated
iii. Calculating Average Time Derivative: -
For each variable, the average of Predictor time Derivative and Corrector Time Derivative is Calculated
(∂(ρ)∂t)avgi=12((∂(ρ)∂t)pi+(∂(ρ)∂t)ci)
Similarly,
(∂(V)∂t)avgi,(∂(e)∂t)avgi,(∂(T)∂t)avgiand(∂(P)∂t)avgi are calculated
iv. Calculating Next time Step Variable: -
ρt+△t=ρti+(∂(ρ)∂t)avgi△t
Similarly,
Vt+△t,et+△t,Tt+△tand Pt+△t 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△x△t
⇒△t=C△xa+V
⇒△ti=C△xiai+Vi
Now,
For each time step, the minimum value of the △t should be chosen for stability
⇒△t=min(△t1,△t2,.......,△tN)
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
V2-V1△x=V3-V2△x
⇒V1=2V2-V3
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
ρN-ρN-1△x=ρN-1-ρN-1=2△x
ρN=2ρN-1-ρN-2
Similarly,
VN=2VN-1-VN-2
TN=2TN-1-TN-2
The Boundary Conditions are as follows: -
| ρ1 | ρo |
| V1 | 2V2-V3 |
| T1 | To |
| ρN | 2ρN-1-ρN-2 |
| VN | 2VN-1-VN-2 |
| TN | 2TN-1-TN-2 |
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 →∂(ρ′A′)∂t′+∂(ρ′A′V′)∂x′=0
Momentum Equation →∂(ρ′A′V′)∂t′+∂(ρ′A′(V′)2+(1γ)P′A′)∂x′=(1γ)P′∂A′∂x′
Energy Conservation Equation →∂(ρ′A′(e′γ-1+(γ2)(V′)2))∂t′+∂(ρ′A′V′(e′γ-1+(γ2)(V′)2)+P′A′V′)∂x′=0
These Forms can be written as follows: -
Continuity Form →∂U1∂t=-∂F1∂x
Momentum Equation →∂U2∂t=-∂F2∂x+J2
Energy Conservation Equation →∂U3∂t=-∂F3∂x
From Comparision between both Forms we can write: -
U1=ρ′A′
U2=ρ′A′V′
U3=ρ′A′(e′γ-1+(γ2)(V′)2)
F1=ρ′A′V′
F2=ρ′A′(V′)2+(1γ)P′A′
F3=ρ′A′V′(e′γ-1+(γ2)(V′)2)+P′A′V′
J2=(1γ)P′∂A′∂x′
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 ρ′, V′, T′, and .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 △t 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 ρ,V,T,and .m in OUTPUT Matrix
B. Explanation of Non-Conservative form Solution Function: -
This Function is Coded to Give the Output Matrix of ρ′, V′, T′, and .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 ρ, V, T, in ρold, Vold, and Told respectively
→ Calculating the Values of △t 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 ρ,V,T,and .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, γ, and Courant Number
→ Creating dx, x Vector, and Area Vector using the input Variables L, N
→ Creating mAnalytical
→ Called both the Functions and generated the Solution Matrix
→ Calculating Pressure and Mach Number at each time step and each node for both Forms of Governing Equations.
iv. Section D(Animating Plot from 1st Time Step to 1400 Time Step): -
→ Created Two Sections of Plot on the Left Side displayed Variables of Conservative Form and on the Right Side displayed Variables of Non-Conservative Form.
→ A Pause is added between each time step to create the Animation
v. Section E(Plotting mass Flow rate for both forms): -
→ Plotted mass Distribution throughout the Nozzle for different time steps for Conservative Form with Analytical mass flow rate.
→ Plotted mass Distribution throughout the Nozzle for different time steps for Non-Conservative Form with Analytical mass flow rate.
→ Plotted mass Distribution throughout the Nozzle at Steady-State from Conservative, Non-Conservative Forms with Analytical mass flow rate.
6. Code: -
A. W7_Assign1_Conser_fun: -
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
B. W7_Assign1_NonConser_fun: -
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
C. Main Solution Script using MacCormack's Technique on Quasi 1D Nozzle Flow: -
% 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')
7. OUTPUT Generated from Codes: -
A. All Variables at different Time Steps: -
i. 1st time step: -
ii. 50th time Step: -
iii. 200th time step: -
iv. 700th Time step: -
v. 1000th time step: -
vi. 1400th time step: -
B. Mass Distribution: -
i. Mass Distribution at different time step for Conservative Form: -
ii. Mass Distribution at different time step for Non-Conservative Form: -
iii. Mass Distribution at steady-State in both Forms with Analytical mass flow rate: -
8. Explanation of the OUTPUTS: -
A. Comparision between Conservative and Non-Conservative Form using Variables Observations: -
The Plots at different time steps for both the forms show(by only observation) that both COnservative and Non-Conservative forms of Governing Equations converge to the Steady-state at nearly the same time and the end result generated by both the forms are similar.
The main benefit of the Conservative Form is in Conditions of Shock Capture which was not the case above so both the forms don't show any deviation in time of convergence.
B. Comparision between Conservative and Non-Conservative Form using mass flow rate Observations: -
The major variation between both forms can be observed in the variation of mass distribution through the nozzle. The Logical answer of mass distribution through the nozzle at any given time step should be constant and in the case of the Conservative form, it can be seen that for most of the time steps the mass flow rate appears to be constant as the Solution variable U2 itself represented mass flow rate so realistic results are computed whereas in the case of Non-Conservative form the mass flow rate at any time step is indirectly generated using ρ, A, and V and as Density and Velocity have different Time derivatives so the deviation of the product of Density, Area, and Velocity at any given time step is significant to be observable in the mass distribution plot.
In the 1400th Time step(Steady-State) the mass distribution calculated by Conservative form throughout the nozzle can be seen as linear whereas because of indirect calculation of mass flow rate in Non-conservative form it is non-linear but still is in very low deviation based on the numerical value and so both the mass flow rate achieved form both the forms are acceptable.
Leave a comment
Thanks for choosing to leave a comment. Please keep in mind that all the comments are moderated as per our comment policy, and your email will not be published for privacy reasons. Please leave a personal & meaningful conversation.
Other comments...
Be the first to add a comment
Read more Projects by Aditya Aanand (17)
Week 10 - Heat Transfer Analysis on BTMS for 4X2 Battery Set using Ansys
AIM:- A. To Perform a Heat Transfer Analysis on A battery set by Inlet Velocities and Gap Between the Batteries. B. The Explain the Relevant Results Generated by the Analysis. Given:- A. The Battery Set Geometry: - The Following are the Parameters for the Battery Set: - i. The…
05 Nov 2022 03:19 PM IST
Week 9 - Parametric study on Gate valve.
AIM: - A. To Perform Simulation on the Gate Valve by Setting the Opening from 10% to 80% B. To Calculate the Mass Flow rate and the Flow Factor for each Case Given: - A. The Given CAD Models of Gate Valves: - …
31 Oct 2022 12:05 PM IST
Week 8 - Simulating Cyclone separator with Discrete Phase Modelling
AIM: - A. To Perform the Analysis of the Given Cyclone separator for different particle sizes and different inlet speeds. B. To Explain the relevant theories and to explain the generated results Given and Assumed: - A. The Geometry of the Given CAD Model: - B. Material and their Properties:…
30 Oct 2022 07:43 AM IST
Week 6 - CHT Analysis on a Graphics card
Problem Statement: - Perform a Steady State CHT analysis on the Given Graphic Card Model for 3 inlet Velocities 1 m/sec, 3 m/sec, and 5 m/sec for Course and Fine Mesh. Create an appropriate Mesh, define appropriate Materials for each Component and Perform a mesh-independent Study. Expected Results: - 1. Explain…
23 Oct 2022 08:49 AM IST
Related Courses
Skill-Lync offers industry relevant advanced engineering courses for engineering students by partnering with industry experts.
Our Company
4th Floor, BLOCK-B, Velachery - Tambaram Main Rd, Ram Nagar South, Madipakkam, Chennai, Tamil Nadu 600042.
Top Individual Courses
Top PG Programs
Skill-Lync Plus
Trending Blogs
© 2025 Skill-Lync Inc. All Rights Reserved.