Chapter one

Introduction and Literature Review

Nowadays, space has become an important scientific, economic, and strategic location. A
crowd of satellites orbits the Earth to monitor its atmosphere and oceans, predict the weather,
broadcast data such as global positioning to users all over the world, and observe closely war
zones. During the last century, space propulsion has been developed to access these orbits and
even space beyond, through chemical propulsion such as liquid and solid rocket engines and
electrical propulsion. As improvements of the existing systems are investigated, new solutions
like solar sails and nuclear propulsion are considered. Hybrid chemical propulsion, using both
liquid and solid propellants, has recently drawn back the attention on it with Spaceship One
successful flight.
Today’s rockets are remarkable collections of human ingenuity that have their roots in
science and technology of the past they are natural outgrowths of literally thousands of years of
experimentation and research on the rocket and rocket proportion. And all kinds of rockets are
work on one law and one principle, its Newton’s third law of motion.
Sir Isaac Newton first presented his three laws of motion in the "Principia Mathematical
Philosophise Naturalis" in 1686. His third law states that for every action (force) in nature there
is an equal and opposite reaction. In other words, if object A exerts a force on object B, then
object B also exerts an equal and opposite force on object A, notice that the forces are exerted on
different objects.
In aerospace engineering, the principle of action and reaction is very important. Newton's
third law explains the generation of thrust by a rocket engine. In a rocket engine, the hot exhaust
gas is produced through the combustion of a fuel with an oxidizer. The hot exhaust gas flows
through the rocket nozzle and is accelerated to the rear of the rocket. In reaction, a thrusting
force is produced on the engine mount. The thrust accelerates the rocket as described by
Newton's second law of motion.

1.1 Hybrid Rocket Motor

Hybrid rocket motor (HRMs) are a class of rocket propulsion systems that have garnered
considerable interest in recent years due to a unique combination of potential advantages as
compared to the more conventional all-solid and all-liquid propellant rockets. These advantages
include a high degree of environmental cleanliness, safety, controllability, simplicity, and ease of
manufacturing, all of which contribute to a low cost of manufacture and operation. And the word
"motor" is as common to hybrid rocket as the word "Engine".
As their name suggests, HRMs employ both a solid and a liquid propellant, most
commonly a liquid oxidizer and a solid fuel. In a “classical” HRMs, the fuel grain is cast in a
combustion chamber and typically contains one or more axial ports along which combustion
occurs. The oxidizer is injected into the fuel port(s) at the head end, and an exhaust nozzle
accelerates the gas or gas-particle combustion products rightward out of the aft end, producing
thrust. This basic configuration is illustrated in Fig. 1.1.
Hypersonic cruise vehicles must push the state of the art in many areas of technology to
achieve practicality and simply to survive the harsh aerothermal flight environment. Vehicles
must necessarily operate close to, but within, safe limits without large safety factors. As
hypersonic flight experience slowly accumulates, it is becoming increasingly clear just how
difficult it will be to achieve practical, safe, and routine hypersonic flight. Accordingly, the need
to achieve successful vehicle designs mandates highly accurate design information on
aerodynamics, aero heating, and propulsion system performance. The source of much of this
design information is the world’s complement of hypersonic speed rocket, which must satisfy
increasingly stringent flow quality requirements in order for designers to know accurately the
limits of hypersonic flight.

Fig. 1.1: Schematic of a hybrid rocket motor.

These engines potentially offer higher specific impulse ( I sp) and much greater
controllability than solid rocket motors, yet greater simplicity and thus potentially lower cost
than liquid rocket engines. They are also characterized by increased safety due to the fact that the
fuel and oxidizer are stored in different phases and are thus much less likely to form explosive
mixtures in the event of a malfunction. As such, they have generated considerable interest in
various rocket propulsion applications in recent years. However, unlike solid rockets, whose fuel
and oxidizer are pre-mixed, or liquid rockets, where propellant mixing and combustion occurs at
the droplet scale, hybrid combustion takes place in a macroscopic turbulent boundary layer down
the length of the fuel grain. This mechanism tends to lead to lower mixing and combustion
efficiencies. The nature of this turbulent diffusion flame combustion mechanism also makes its
internal ballistics somewhat more difficult to model accurately than either a traditional solid- or
liquid propellant rocket.
Finally, the fuel surface regression rate of an HRMs also tends to be significantly lower
than that of a solid-propellant rocket motor, often necessitating either complicated fuel grain
geometries or an otherwise more elaborate engine design in order to provide a sufficient burning
surface area to achieve target thrust levels. This becomes progressively more problematic as
engines are scaled to larger sizes, and the added complexity can negate many of the advantages
of a classical HRMs.
1.2 Type of nozzle:
The nozzle is a duct of varying cross-sectional area in which the velocity of fluid increases
with the corresponding drop in pressure, there are three types [1] in general are mentioned
below:
1.2.1 Convergent nozzle:
In the convergent nozzle, the cross-sectional area decreases from the Inlet section to the
outlet section.

Fig. 1.2 convergent nozzle

1.2.2 Divergent nozzle:

In the divergent nozzle, the cross-sectional increases from the inlet section to the outlet
section.

Fig. 1.3 divergent nozzle

1.2.3 Convergent-Divergent nozzle (C.D. nozzle):

In the convergent-divergent nozzle or De Laval nozzle, the cross-sectional area first
decreases from the inlet section to the throat and then increases from its throat to the outlet
section.
 Fig. 1.4 convergent-divergent nozzle
1.3 Different rocket nozzle configurations:
1.3.1 Conical nozzle:
The conical nozzle was often used in early rocket applications because of its simplicity
and ease of construction. The cone gets its name from the fact that the walls diverge at a constant
angle.

Fig. 1.5 Conical nozzle.

1.3.2 Contoured nozzle:

A contoured nozzle or (Bell nozzle) is the most commonly used nozzle shape, offering
significant advantages over the conical nozzle, both in size and performance. Near the throat, the
nozzle diverges at a relatively large angle but the degree of divergence tapers off further
downstream of the nozzle. Near the nozzle exit, the divergence angle is very small. In this way, it
minimizes weight while maximizing performance. The most important design issue is to contour
the nozzle to avoid oblique shocks and to maximize performance.

Fig. 1.6 Contoured nozzle.

1.3.3 Annular nozzle:

The throat section of such nozzle is an annulus formed between a central plug and an
outer tailpipe. Downstream of the throat, the exhaust gases are made to expand in a diverging
annulus formed between the diverging tailpipe and the converging central plug. The thrust loss
due to flow divergence in an annular nozzle, with diverging walls formed by conical surfaces,
would be less than that of a conical nozzle of the same area ratio and length [13].
 Fig. 1.7 Annular nozzle.

1.3.4 Plug nozzle:

The throat of such a nozzle, as shown in Fig. 1.8, is located as an annulus at the outer
diameter. At the outer edge of the annulus, the exhaust gases expand abruptly to the ambient
pressure. The plug surface may be designed in order to produce uniform exit flow parallel to the
nozzle axis. The cowl-lip diameter would be the same as the exit diameter of a uniform flow of a
C-D nozzle expanding the gases to the same ambient pressure. However, the plug-type would be
much shorter than the equivalent C-D one.

Fig. 1.8 Pplug nozzle.

1.3.5 Expansion-Deflection nozzle:

E-D type nozzles are constructed as a compact combustion chamber with the throat
section annulus located close to the nozzle centerline. M e exhaust gases are issue from the throat
in an outward direction and expand around the shoulder of the central plug, the nozzle wall
contour would then turn them into a nearly axial direction.
 Fig. 1.9 Expansion-Deflection nozzle

1.4 Historical Background & Literature review:

One of the major parameters in the design of rocket motors is the determination of the
optimum contour of the exhaust nozzle. Because the performance of rocket motors is highly
dependent on the aerodynamic design of the expansion nozzle, the main design parameters being
the contour shape and the area ratio. The optimal design of traditional bell-type nozzles for given
operating conditions (i.e. chamber and ambient pressures) is already supported by accurate and
validated tools. This has a great effect on the determination of the thrust developed by the rocket.
A convergent-divergent De Laval nozzle is typically used to convert the heat liberated in the
combustion chamber into kinetic energy for propulsion. Thrust is mainly produced by the
momentum imparted to the products of combustion when discharging through the exhaust
nozzle.
During their passage, the gases are continuously accelerated from low subsonic to high
supersonic velocities. The C-D nozzle may be divided into three parts:
1. The convergent subsonic section whose design influences the mass flow of the exhaust
gases and, to some extent, the combustion efficiency achieved in the chamber.
2. The throat section which determines, with the operating conditions in the combustion
chamber, the mass flow rate through the nozzle.
3. The divergent supersonic section whose wall configuration, together with the exit area,
determines the additional velocity imparted by the expansion of gases through this
portion of the nozzle.
This division is most appropriate because of the different effects each part has in
determining the thrust developed, and because of the different methods of analysis which have to
be used for computing the flow field in each one of the regions

1.4.1 Existing survey for nozzle:

In the early years of rocket and missile nozzle design, conical nozzles were designed
because the methodology for designing efficient contour nozzles did not exist. Although they
were necessary at that time, conical nozzles were an inadequate solution. Such nozzles resulted
in performance losses because they did not produce axially directed exhaust. The angularity in
the conical nozzle leads to incorrect nozzle design for wind tunnel applications. In some cases,
these nozzles were excessively heavy because of their length, which leads to enormous drag.
This is one of the main reasons for supersonic flights were not possible in the early days.
Therefore, the design of an efficient nozzle and parallel flow exhaust was a challenging task.
From the report G. Satyanarayan et al [2] it is studied and analyzed that for some inlets and
boundary conditions the rectangular nozzle gives an increased velocity of 23.93% compared to
the square nozzle and 24.47% compared to the circular nozzle. This is due to the pressure as well
as the temperature drop into the rectangular nozzle s greater than the other two nozzles. The
rectangular nozzle gives an increased pressure and temperature drop of about 22.93% and
42.56% compared to the square nozzle. Also, the rectangular nozzle gives increased temperature
and pressure drop of about 23.97% and 43.68% compared to the circular nozzle, and hence the
rectangular cross-section nozzle the best cross-section compared to circular and square for the
same input and boundary condition.
Venkatesh V. et al [3] his paper observed that the result obtained by using Fluent After the
simulating the transient gas flow by coupled explicit solver it gives a 2-D result which shows that
the contoured nozzle gives a greater Mach number at exit compared to the conical nozzle. As the
contoured nozzle gives a maximum expansion ratio compared to the conical nozzle. Hence the
conical nozzle has to be used at sea level whereas the contoured nozzle to be used at a higher
altitude since a greater expansion ratio is required at higher altitude for a given length of the
nozzle. But contour nozzle has complex geometry that is the profile may be parabolic or cubic
therefore it is difficult to fabricate.
From the report Arjun Kundu et al [4] it is observed that the result obtained after the CFD
analysis shows that a smaller exit diameter gives a greater Mach number compared to the larger
diameter for the same inlet and boundary condition. The first nozzle with an exit diameter 60 cm
increase in velocity from inlet to outlet is 137.37 m/s but in the second nozzle with an exit
diameter 53 cm the increase in velocity from inlet to exit is 152.6 m/s. This represents that
velocity at exit increases with a decrease in diameter. The nozzle with an exit diameter 53 cm
gives a maximum Mach number at exit compared to the nozzle with 60 cm diameter.
Karna S. Patel [5] his paper provides a discussion on the study of variation of Mach
number at the exit of the nozzle with variation in the divergent angle with the help of CFD
analysis and the objective of the research is to investigate the best suited divergent angle. As well
as the phenomenon of oblique shock is visualized to find the suitable divergent angle to prevent
shock.
The report of Nikhil D. Deshpande et al [6]. The theoretical and CFD analysis of the De-
Laval nozzle have been carried out which was aiming to provide theoretical formulae to calculate
velocity, pressure, the temperature at every section of the nozzle. The validation of these is
carried out using CFD software ANSYS FLUENT. Results obtained theoretically are verified
with the results obtained by CFD and it was found that Divergence Angle (Degree) Mach
number results are identical. There are some variations in the results because the theoretical
study does not account for the boundary layer effect, shock wave, radial velocity component,
wall friction, and so on. But the variation is quite insignificant.
Mohan Kumar G et al [7] worked on “design and optimization of De-Laval nozzle” and his
paper had been carried out to prevent shock-induced flow separation. The aim was to achieve
maximum thrust without flow separation due to flow separation. For maximum thrust and
efficiency, the direction of flow of the stream through the nozzle should be axial. Flow
conditions were selected based on pressure, temperature, and gases that are available at the exit
of the combustion chamber. To prevent shock-induced flow separation, overexpansion
conditions, and optimum expansion conditions are studied. And it was found that optimum
expansion is critical for the efficient operation of the nozzle.
From the report of Sourabh Shriwas et al [8] reduction of jet noise in the exhaust
nozzle of aircraft engine has been studied. It was found that various methods could be adopted
for noise reduction. Such as the use of the Chevron nozzle which consists of a sharp saw tooth at
the aft which causes smooth mixing of hot steam with the bypass steam hence turbulence is
reduced. Also making use of Nanomaterials with a higher density is capable of absorbing noise,
increasing the bypass ratio of turbofan, acoustic liners can also be implemented for noise
reduction.
Vignesh. Metal [9] from his work on “conceptual Design of Short Take-Off Supersonic
Aircraft with cold flow Nozzle “is carried out in 2013. The rectangular C-D nozzle is an altered
type DE-Laval nozzle having a reduced frontal area. The ultimate idea is to reduce takeoff
distance. By placing the C-D nozzle under the wing of military aircraft flow can be accelerated at
the choked flow condition resulting in a sudden increase in thrust and short takeoff. A numerical
study was carried out which highlighted that for the design shock waves are not formed in the
divergent section and also takeoff distance reduces by 50%.

1.5 Research problem:

1.6 Objectives:
The broad objectives of this study are:
1. Develop a two-dimensional method to predict shock and expansion free nozzle
shapes for the divergent section of the supersonic nozzle of the hybrid rocket.
2. Include numerical investigation of flow-through and after the nozzle of a hybrid
rocket motor.
3. Use computational fluid dynamics to determine whether the calculated profiles
produce uniform flow and free of shock in the supersonic section.
4. Study the comparison of conical nozzle and bell nozzle.
5. Study the effects of length and area change on the pattern.

1.7 Scope of the work:

Engineering Gas Dynamic is one of the basic sciences in mechanical engineering which
deals with the study of the motion of gases. The gas-dynamic equation was used here to study
and design all types of nozzles used in hybrid rocket engine and the operation they do no it. The
computational fluid dynamic was used here as an effective tool for a numerical investigation of
flow.
Method of characteristics (MOC) is used here for designing a supersonic nozzle. The
minimum length supersonic contour nozzle for rocket application is designed and the length of
the nozzle is optimized using MatLab 2016 software and OriginPro 9.0 is statictcs program used
for analyzing the datd. The nozzle contour is designed for different Mach number and altitude
conditions. A solver is to be chosen for meshing is Gambit it is very helpful which was chosen.
The analyzed will be in ANSYS 17.5 CFD full packeg. And finally, we have compared the
results we have got with published results for validation.
 Chapter two
Theory of Rocket

2.1 Fundamentals:
A convergent-divergent nozzle is typically used to convert the heat liberated in the
combustion chamber into kinetic energy for propulsion. Thrust is mainly produced by the
momentum imparted to the products of combustion when discharging through the exhaust
nozzle.
During their passage, the gases are continuously accelerated from low subsonic to high
supersonic velocities. The convergent-divergent nozzle may be divided into three parts:
1. The convergent subsonic section whose design influences the mass flow of the exhaust
gases and, to some extent, the combustion efficiency achieved in the chamber.
2. The throat section which determines, with the operating conditions in the combustion
chamber, the mass flow rate through the nozzle.
3. The divergent supersonic section whose wall configuration, together with the exit area,
determines the additional velocity imparted by the expansion of gases through this
portion of the nozzle as shown in the Fig.2.1 below.

Fig. 2.1 showing approximate flow velocity (v), together with

the effect on temperature (T) and pressure (p).
 As Anderson demonstrates in his book Modern Compressible Flow with Historical
Perspective [16], the following equation can be derived from the quasi-one-dimensional flow.
ⅆA du
=(M ¿¿ 2−1) ¿
A u (2.1)

The equation states that, for a compressible isentropic flow through a duct, an increase in
cross-sectional area A results in a decrease in velocity u if the Mach number is lower than unity.
This is very familiar to us in everyday life; who hasn’t squeezed the end of a garden hose to
increase the velocity of the water jet coming out of it? What is not so familiar to most of us is the
other truth this equation reveals: that an increase in area A causes an increase in velocity if the
Mach number is greater than one.
Gustav de Laval experienced this first hand in 1888 when he significantly improved the
performance of his impulse steam turbine by adding a divergent section to the convergent nozzle
that propelled the turbine blades, creating the first convergent-divergent nozzle. However, he
was not familiar with the equation above and, in fact, initially, it was disputed whether the flow
was indeed supersonic within the divergent section. The existence of supersonic flow inside the
divergent section was not confirmed until Aurel Stodola published the first scientific experiments
on convergent-divergent nozzles in 1903, confirming the shock theory of Rankine and Hugoniot
from a few decades earlier [12].
Although De Laval’s and Stodola’s work with nozzles was intended to improve the
performance of impulse steam turbines, the primary application that motivates the current work
in rocketry. A rocket is a flight vehicle that propels itself by using the reacting forces created by
expelling part of its mass at high speeds. The mass expelled is the rocket propellant, composed of
a fuel and an oxidizer that is burned within a combustion chamber to create the hot and highly
pressurized exhaust gas which is expelled from the nozzle exit at high speeds to achieve
maximum thrust. The thrust obeys Newton’s second law, which states that the rate of change of
momentum of an object is equal to the force acting on it:

ⅆm
F= ν
ⅆt e (2.2)

This method of propulsion has been known to man since at least the 13th century when
the first solid propellant rockets were used by the Chinese as a part of their” fire” arrows [1] In
1813, the English mathematician William Moore published equations describing the motion of
rockets, including the famous delta-v equation. Moore’s efforts were intended for military
applications, but the first man to envision the use of rockets for space exploration was the
Russian schoolteacher and visionary Konstantin Tsiolkovsky, who also independently derived
the delta-v equation and published it in 1903.
The delta-v equation describes how the change in a space vehicle velocity (∆V) in the
absence of gravity or any other body forces, is connected to its initial mass (m 0) including the
propellant, its final mass (mf) after the propellant has been expelled [16], and the effective
propellant exhaust
velocity.
m0 (2.3)
ΔV =v e ln
mf
 It is important to hold the initial mass (m 0) of a rocket to a minimum to reduce the energy
needed to reach the orbit. The force equation Eq. (2.2) and the delta-v equation (Eq.2.3) both tell
us that m0 can be lowered while keeping F or ∆V constant, by increasing the exhaust velocity, v e.
It was the American rocket pioneer Robert H. Goddard who first realized that the de Laval
nozzle could be used to achieve high rocket exhaust velocities, and he proceeded to build the
world’s first liquid rocket engine in 1926 [1].

The thrust, F, created by a nozzle is given as:

F=ṁ ν e +( p ¿ ¿ e− pa ) A ¿
(2.4)
Whereṁ is the mass flow rate through the nozzle and v e is the exhaust velocity, pe is the
static pressure at the nozzle exit, pa is the ambient pressure and Ae is the nozzle’s exit cross-
sectional area.
Note that Eq (2.4) in Eq (2.2) plus an added term accounting for the pressure mismatch
between the ambient air and the exhaust. For an ideally expanded flow, the pressure at the nozzle
exit (pe) is equal to the ambient pressure, (pa) making the second term in Eq (2.4) zero.
Two fundamental parameters that are frequently used in rocket engine design and theory
are the specific impulse (Isp) and the expansion ratio. The specific impulse is defined as:

F
I sp =
ṁ (2.5)

It is a simple measure of rocket engine performance that compares the engine thrust to its
mass flow rate. For an ideally expanded nozzle flow I sp = v e and in that case, we can replace v e
with Isp in Eq (2.3) and truly appreciate its importance to rocket efficiency. In the form above, its
unit is m/s but it is commonly normalized with Earth’s gravitational constant g = 9.81m/s 2,
giving it a unit of seconds [1]. A typical first stage engine has a specific impulse between 300s
and 450s.
The expansion area is defined as the ratio between the exit cross-sectional area and throat
cross-sectional area:

Ae
ε=
At (2.6)

“A high expansion ratio is desired as it maximizes exhaust velocities. However, at sea

level, the ambient air pressure limits this due to potential flow separation, as will be covered in
the next section. A practical limitation exists as well, as a high expansion ratio means a high exit
diameter with associated packaging challenges." [15]

2.2 Nozzle operation:

The purpose of any engine nozzle is to direct the exhaust of a rocket engine in one
direction, generating thrust in the opposite direction. The exhaust, a high-temperature mix of
gases, has an effectively random momentum distribution (i.e., the exhaust pushes in any
direction it can). If the exhaust is allowed to escape in this form, only a small part of the flow
will be moving in the correct direction and thus contribute to forward thrust. The nozzle redirects
the exhaust moving in the wrong direction so that it generates thrust in the correct direction.
As shown in Fig. 2.2, state (1) the flow through the nozzle when it is completely
subsonic. The flow accelerates out of the chamber through the converging section, reaching its
maximum (subsonic) speed at the throat [16]. The flow then decelerates through the diverging
section and exhausts into the ambient as a subsonic jet. Lowering the backpressure in this state
increases the flow speed everywhere in the nozzle.

Fig. 2.2 Flow in a Converging-Diverging nozzle.

In-state (2) the flow is the same as the state (1) but a slightly higher led to a decrease in
the pressure and an increase in the velocity.
Further lowering p e results in the state (3). The flow pattern is exactly the same as in
subsonic flow, except that the flow speed at the throat has just reached Mach 1. Flow through the
nozzle is now choked since further reductions in the backpressure can't move the point of M=1
away from the throat. However, the flow pattern in the diverging section does change as the
backpressure is lowered further.
In the state (4) as p e is lowered below that needed to just choke the flow a region of
supersonic flow forms just downstream of the throat. Unlike a subsonic flow, the supersonic flow
accelerates as the area gets bigger. This region of supersonic acceleration is terminated by a
normal shock wave. The shock wave produces a near-instantaneous deceleration of the flow to
subsonic speed. This subsonic flow then decelerates through the remainder of the diverging
section and exhausts as a subsonic jet. In this regime, if the backpressure is lowered or raised the
length of supersonic flow in the diverging section before the shock wave increases or decreases,
respectively.
 State (5) show that if p e is lowered enough the supersonic region may be extended all the
way down the nozzle until the shock is sitting at the nozzle exit. Because of the very long region
of acceleration (the entire nozzle length), the flow speed just before the shock will be very large.
However, after the shock, the flow in the jet will still be subsonic.
Lowering the backpressure further causes the shock to bend out into the jet. and a
complex pattern of shocks and reflections is set up in the jet which will now involve a mixture of
subsonic and supersonic flow, or (if the backpressure is low enough) just supersonic flow.
Because the shock is no longer perpendicular to the flow near the nozzle walls, it deflects it
inward as it leaves the exit producing an initially contracting jet. We refer to this as an over-
expanded flow because in this case, the pressure at the nozzle exit is lower than that in the
ambient (the back pressure)- i.e. the flow has been expanded by the nozzle too much, as shows
state (5).
A further lowering of the back-pressure changes and weakens the wave pattern in the jet.
Eventually, the backpressure will be lowered enough so that it is now equal to the pressure at the
nozzle exit. In this case, the waves in the jet disappear altogether. and the jet will be uniformly
supersonic. This situation, since it is often desirable, is referred to as the 'design condition ‘,
Pe=Pa.
Finally, in-state (7) if the backpressure is lowered even further, we will create a new
imbalance between the exit and back pressures (exit pressure greater than backpressure). In this
situation, called under-expanded, expansion waves that produce gradual turning and acceleration
in the jet form at the nozzle exit, initially turning the flow at the jet edge outward in a plume and
setting up a different type of complex wave pattern.

2.3 Nozzle flow separation:

Flow separation in supersonic nozzles is a basic fluid dynamic phenomenon that occurs at a
certain pressure, resulting in shock formatting and shock turbulent boundary layer interaction
inside the nozzle.
When the pressure ratio is pushed further from them the designed one, either at highly
overexpanded nozzles or in a transitional phase when full-flowing is not still established, flow
violently separates from the nozzle wall upstream due to an adverse pressure gradient as shown
in Fig. 2.3. The separated boundary-layer flow continues as a free jet and exits the nozzle. This
separation mode named free shock separation (FSS) was observed early as reported in the
number of empirical models in the late ‘50s and ‘60s and it is well understood till now. With the
further development of the high-power nozzles, namely TOC and compressed parabolic as for
example of J-2S engines asymmetric FSS has been reported and investigated. In the ‘70s during
the cold flow tests of the J-2S, Nave and Coffey first reported the phenomenon of asymmetric
flow separation that yields tilted separation surface and causes the highest value of the side loads.
Particularly, pressure downstream of the separation showed unsteady behavior with the strong
oscillations and values quite above the ambient one which they attributed to the reattachment of
the separated flow to the nozzle. Since the separation extension is limited by this, they named it
restricted shock separation (RSS).
The same behavior is noted during the development of the contemporary used TOC
launcher nozzles as Vulcain 2 nozzle and reported in works as of Hagemann and Frey [16], Deck
[10], Reijasse, Schimizu et al., Roquefort. Hagemann and Frey [16] related the RSS problem
with the internal shock and cap shock pattern which, as mentioned, are characteristic for the
TOC and compressed parabolic nozzles. The strong side loads due to RSS are one of the major
concerns of the modern launcher rocket engine design. In the current study, the TOC nozzle is
investigated for pressure ratios under which FSS, RSS, and shock patterns occur at the hysteresis
regime.

Fig. 2.3: Schematic drawing and nomenclature of free shock separation [23].

2.4 Method of characteristic:

The Method of Characteristics (MOC) is a numerical procedure appropriate for solving
two-dimensional compressible flow problems. By using this technique flow properties such as
direction and velocity, can be calculated at distinct points throughout a flow field. The method of
characteristics, implemented in computer algorithms, is an important element of supersonic
computational fluid dynamics software. These calculations can be executed manually, with the
aid of spreadsheet programming or technical computing software is the number of characteristic
lines increase, so do the is a point, and the manual calculations can become exceedingly tedious
Omer Musa [25].
 The method of Characteristics was developed by the mathematicians Jaques Saloman
Hadamard in 1903 and by Tullio Levi-Civita in 1932 [16]. The method of characteristics uses a
technique of following propagation paths in order to find a solution to partial differential
equations.
The physical condition of a two-dimensional, steady, isentropic, irrotational flow can be
expressed mathematically by the nonlinear differential equation of the velocity potential. The
method of characteristics is a mathematical formulation that can be used to find solutions to the
aforementioned velocity potential, satisfying given boundary conditions for which the governing
partial differential equations (PDEs) become ordinary differential equations (ODEs).
The name comes from a method used to solve hyperbolic partial differential equation,
finding “characteristic lines" (combinations of the independent variables) along which the partial
differential equation reduces to a set of ordinary differential equations, or even, in some cases, to
algebraic equations which are easier to solve. The applications of the method of characteristics
for nozzle flows are not limited to the design of contours. The method may also be used to
analyze the flow field inside a known contour as well. The method is also not limited to the flow
within the nozzle. The approach can be extended to analyze the exhaust plume for both under
expanded and over-expanded nozzle flow using the free pressure boundary of the exhaust plume
[7].

Fig. 2.4 analysis the flow inside the nozzle with

MOC.

2.5 MOC Applied to Contour Design:

The MOC technique meshes very cleanly with the needs of a contour design procedure.
Further, MOC is easier to apply to the direct-design problem than to the analysis of a known
contour because the boundary conditions are easier to apply in the design problem. With the help
of Fig. 4, the direct-design procedure may be summarized as follows. First, the nozzle exit
conditions are specified in terms of the exit Mach number and the nozzle exit radius. With the
assumption of uniform and parallel exit flow along with the characteristic nozzle in Fig. 4, the
entire solution (ν, θ, μ, and M) along the nozzle is known. Grid points may be distributed along
with the nozzle. Two additional boundary conditions are needed: (1) a centerline Mach number
distribution, and (2) an upstream boundary condition near the throat. A typical throat boundary
condition is the right-running throat characteristic TI in Fig. 4, which may be constructed if a
flow-field solution is available for a specified throat contour up to point T. Note that for the
design application, a characteristic can serve as a boundary condition. Other starting line
alternatives with varying degrees of suitability include a straight, vertical, slightly supersonic
(uniform flow) line at the geometric throat.
 a straight, vertical, slightly supersonic line at proper quasi-1-D area ratio and x-station
with flow angle linearly varying from zero on the centerline to wall angle.
 a spherical source-flow start line at a specified station.
 a transonic solution for circular-arc throat along a vertical line at a specified station.
 a transonic solution for circular-arc throat along a straight line with a specified slope
angle.

Fig. 2.5 MOC Applied to Direct Design Problem

With the inlet and outlet boundary conditions defined, the centerline Mach number
distribution may be chosen, preferably for compatibility with the inlet and outlet. The centerline
Mach number distribution should be monotonic and should match some derivatives of the inlet
and exit boundary conditions (typically derivatives of Mach number with respect to x up to third
order, but frequently less).
With the boundary values computed at selected grid points along with the nozzle in Fig. 2.5
the MOC design procedure can begin Fig. 2.5 Starting at point C on the centerline, interior
characteristic intersection points can be constructed, first marching up characteristic CD to
construct the next upstream characteristic. A numerical procedure similar to that for the analysis
problem can be applied with modifications, but this procedure will not be detailed here. As each
characteristic is constructed, mass flux can be integrated along with the characteristic, the left
characteristic being treated as a surface through which mass is flowing until a mass flow rate
equal to that crossing CD is reached. This defines a point on the nozzle wall. When the entire
characteristic network has been so erected, the complete nozzle contour is fully defined. Note
that the MOC-based contour design does not proceed upstream of the throat. Consider also that
the inviscid and perfect-gas contour so designed does not depend on tunnel stagnation
conditions.

Fig. 2.6 Construction of Characteristic Solution

The MOC was used here to design a Contour nozzle, where a MatLab program was used
to draw a profile shape of the nozzle by using a gas dynamic equation.

2.6 Computational Fluid Dynamics:

It is a branch of fluid mechanics using a numerical method to solve fluid problems using
computer simulation to visualize how a gas or liquid flows as well as how the gas or liquid
affects objects as it flows. Computational fluid dynamics is based on the Naiver-Stokes
equations. These equations describe how the velocity, pressure, temperature, and density of a
moving fluid are related.
Computational fluid dynamics has been around since the early 20th century and many
people are familiar with it as a tool for analyzing airflow around cars and aircraft. As the cooling
infrastructure of server rooms has increased in complexity, CFD has also become a useful tool in
the data center for analyzing thermal properties and modeling airflow. CFD software requires
information about the size, content, and layout of the data center. It uses this information to
create a 3D mathematical model on a grid that can be rotated and viewed from different angles.
CFD modeling can help an administrator identify hot spots and learn where cold air is being
wasted or air is mixing.
 Simply by changing variables, the administrator can visualize how cold air will flow
through the data center under a number of different circumstances. This knowledge can help the
administrator optimize the efficiency of the existing cooling infrastructure and predict the
effectiveness of a particular layout of IT equipment. For example, if an administrator wanted to
take one rack of hard drive storage and split the hard drives over two racks, a CFD program
could simulate the change and help the administrator understand what adjustments would need to
be made to deal with the additional heat load before any time or money has been spent.

2.6.1 Methodology:
In all of these approaches, the same basic procedure is followed.
 During preprocessing
o The geometry and physical bounds of the problem can be defined using computer-
aided design (CAD). From there, data can be suitably processed (cleaned-up) and
the fluid volume (or fluid domain) is extracted.
o The volume occupied by the fluid is divided into discrete cells (the mesh). The
mesh may be uniform or non-uniform, structured, or unstructured, consisting of a
combination of hexahedral, tetrahedral, prismatic, pyramidal, or polyhedral
elements.
o The physical modeling is defined – for example, the equations of fluid motion +
enthalpy + radiation + species conservation
o Boundary conditions are defined. This involves specifying the fluid behavior and
properties at all bounding surfaces of the fluid domain. For transient problems, the
initial conditions are also defined.
 The simulation is started and the equations are solved iteratively as a steady-state or
transient.
 Finally, a postprocessor is used for the analysis and visualization of the resulting solution.

2.6.2 Governing Equations:

For a Newtonian fluid, the governing equations of fluid flow describing the conservation
of mass, momentum, and energy in Cartesian coordinate systems can be written as follows [29]
Conservation of Mass
∂ ρ ∂( ρu) ∂( ρν) ∂( ρw)
+ + + =0
∂t ∂x ∂y ∂z (2.7)

Conservation of Momentum
We choose they-coordinate to be in the vertical direction opposite to the direction in gravity.
The buoyancy forces act in the y-direction. The conservation of momentum equations in the
Cartesian coordinates can be written as
x-Momentum Equation
 ∂ ∂( ρuu) ∂ (ρvu) ∂( ρuw ) −∂ P ∂ 2 u ∂2 u ∂2 u
∂t
( ρu)+
∂x
+
∂y
+
∂z
=
∂x
+μ +( +
∂ x2 ∂ y2 ∂ z2 ) (2.8)

y-Momentum Equation

∂ ∂ ∂ ∂ −∂ P ∂2 v ∂2 v ∂2 v
∂t
ρ∨+
∂x
ρuV + ρvv + ρwv =
y ∂z ∂y
+μ + (
+
∂ x2 ∂ y2 ∂ z2 )
− ρgβ(T ¿¿ ∞−T )¿
(2.9)

z-Momentum Equation

∂ ∂ ( ρuw) ∂(ρvw ) ∂( ρww) −∂ P ∂2 w ∂2 w ∂2 w

∂t
( ρw )+
∂x
+
∂y
+
∂z
=
∂x
+μ + (
+
∂ x2 ∂ y2 ∂ z2 ) (2.10)

Conservation of Energy

∂ ∂ ∂ ∂ ∂2 T ∂2 T ∂2 T
∂t
( ρ c p T )+ ( ρuc p T )+
∂x ∂y ∂z (
( ρv c p T )+ (ρw c p T )=k + +
∂ x2 ∂ y2 ∂ z2 )
+q
(2.11)

2.6.3 Naiver-Stokes Equations:

The Naiver-Stokes equations describe the full range of scale present in fluid flows,
including the finest scales of turbulence [17]. However, with the exception of cases where the
equations can be greatly simplified, there exists no known analytical solution to the equations.
Computational Fluid Dynamics (CFD) is therefore an important research field and an
engineering tool that deals with solving the Naiver-Stokes equations numerically. Direct
Numerical Simulations, commonly known as DNS, attempt to solve the flow equations presented
in the preceding section, and, in order to achieve accurate solutions, the DNS needs to resolve
the motion of the finest turbulent eddies.
This requires the computational grid used to obtain the numerical solution to be fine
enough to capture the geometrical features of the smallest turbulent eddies, and the time step the
solver takes needs to be small enough to capture their dynamics. For high speed, wall-bounded
flows of engineering interest, the range of scales is so wide that the computational cost of
resolving the whole turbulence spectra is too high to be economically feasible or even possible to
perform. This has led to the field of turbulence modeling, where the objective is to either model
the effect of the whole turbulence spectrum on the mean flow or to model the effects of the parts
of the spectrum that contain the finer scales of turbulence on the resolved spectrum.

2.6.4 Turbulence Model:

2.6.4.1 The k-ω Turbulence Model:
The k-ω model is based on the transport equations for the turbulence kinetic energy k
and the specific dissipation rate ω. The k-ω model is supposed “more accurate than k-ε in the
near-wall layers, and has therefore been successful for flows with more moderate adverse
pressure gradients, but fails for flows with pressure-induced separation” [29] .Shear-Stress
Transport modification, k-ω is widely used in industrial, commercial, and research codes;
however, because of its free-stream sensitivity to the values of ω, it has yet to overtake k-ε in
popularity. The transport equations for k and ω are given as:
∂ ρk ∂ ρu i k ∂
+ = Γ
∂t ∂ xi ∂ xJ k (2.12)

∂ ( ρk ) ∂ ( ρu i k ) ∂ ∂k ~
∂t
+
∂ xi
=
(
∂xj
Γk
∂ xj)+ Gk −Y k + S k
(2.13)

∂(ρω) ∂( ρui ω) ∂ ∂ω ~
∂t
+
∂ xi
= Γ
∂ xj k ∂ xj ( )
+ Gω −Y ω+ S ω + Dω
(2.14)

~
where Γ k and Γ w represents the effective diffusivities, S is the strain rate magnitude, Gk
~
represents the generation of turbulence kinetic energy due to mean velocity gradients, Gω
represents the generation of ω, Y represents the dissipation of its respective variable due to
turbulence, Sk and Sω are user-defined source terms, and Dω is the representative cross-
diffusion term :
Dω=2 ¿
(2.15)

where F1 is the blending function and σ w ,2 = 1.168.

2.6.4.2 The SST Turbulent Model:
The SST k-ω model was proposed by Menter in [17]. The SST formulation combines
the best of the two models i.e. the use of the k -ω formulation in the inner parts of the boundary
layer: which makes the model directly usable all the way down to the wall through the viscous
sub-layer, hence the SST k - ω model can be used as a ’Low-Re turbulence model’ without any
extra damping functions. The SST formulation also switches to a k - ε behavior in the free-
stream and thereby avoids the common k -ω problem that the model is too sensitive to the inlet
free-stream turbulence properties. This model shows good behavior in adverse pressure
gradient and separated flows.
The transport equations for SST k - ω reads:
∂(ρω) ∂ γ ∂ ∂k
∂t
+
∂ xj
( ρωU j )= ρ P−β ¿ ρωω+
μt ∂ xj
( μ+σ ω μ ) −(∂ xj )
+2 ¿
(2.16)

∂(ρω) ∂ γ ∂ ∂k
∂t
+
∂ xj
( ρωU j )= ρ P−β ¿ ρωω+
μt ∂ xj
( (
μ+σ ω μ ) −
∂ xj )
+2 ¿ (2.17)
With,
 ρk /ω
μt =
Ω F2
(
ma x 1 ,
a1 ω ) (2.18)

Modeling the effects of turbulence on the main flow is usually done by decomposing the
flow variables according to Reynolds decomposition,

ϕ=ϕ́+ ϕ' (2.19)

Where the flow variable is φ decomposed into its mean component φ and its fluctuating
component ϕ '. This decomposition, or separation of scales, can be applied to the Naiver-Stokes
equations to derive the Reynolds-averaged Naiver-Stokes equations (RANS). when dealing with
compressible flows, it is convenient to use Favre-filtering to avoid having to model extra terms
in the Reynolds-averaged continuity equation [17]. A Favre-filtering flow quantity φ~ is density-
weighted such that.
ρϕ
´
ϕ́=
ρ́ (2.20)

and the following decomposition applies:

ϕ=ϕ́+ ϕ ¿
(2.21)
Where ϕ ¿ can be thought of as the high-frequency component and φ~ as the low
frequency component of the flow variable φ
Using this the unsteady Favre-filtered Reynolds-averaged Naiver-Stokes (URANS) equations
can be formulated:
∂ ρ́ ∂( ρ́ ~
u j)
+ =0
∂t ∂ xj (2.22)

∂( ρ́ ~
ui ) ∂( ρ́ ~
ui ~u j ) ∂ Ṕ ∂ σ´ij ∂ τ ij
+ x
=¿ + +
∂t ∂ j ∂ xi ∂ xi ∂ x j (2.23)

∂ ρ́ ~e0 ∂ ρ́ ~e 0 ~
u j −∂ ṕ ~
uj ∂ μ ∂~ T ∂ 1 ∂
∂t
+
∂ xj
=
∂xj
+
∂ xj
Cp(pr ∂ x j
+q tj +)∂x j
( ui´σ ij )−
2 ∂x j
ui ui u j´−~
ρ(¿~ u i ui ~
uj ) ¿
(2.24)

The Favre-filtered viscous stress tensor reads

2 (2.25)
σ´ij =μ 2 ~
(S ij− ~S δ )
3 mm ij
and the Favre-filtered strain rate tensor is
~ ~
~ 1 ∂ui ∂ uj
Sij = ( +
2 ∂ x j ∂ xi ) (2.26)

Two additional terms have now appeared, the turbulent stress tensor, τij, and the turbulent
heat flux term, q tj The turbulent stress tensor is
τ ij =−ρ́ (~
ui u j − ~
ui ~
u j)

~
~u ~ ~ ~ ~ ( u} {widetilde {u}} rsub {j}} widetilde {+ {widetilde {u}} rsub {i} {u} rsub {j} rsup { ~ } {u} rsub {j} rsup { (2.27)
)+ u⏟
¿−ρ́ (⏟
(i u j −u i u j ) +⏟
I
i
II
i
III
)
Terms I-IIIcalled the Leonard stresses, cross stresses and Reynolds stresses, respectively.
Similar to τij the turbulent heat flux is given as
~ ~ ~
q tj=−C p ρ́−~T u j −~
T~u j=−C p ρ́( ~
T~u j −~
T~u j+ T } {widetilde {u}} rsub {j}} + widetilde {widetilde {T} {u} rsub {j} rsup { + T } {u} rsub {j} rsup { )
(2.28)

Two terms remain to be addressed in the energy equation, i.e. the product ui σ ij can be
replaced with ui ¿ + τ ij ) and the last term in the energy equation is considered negligible and is
dropped. The Naiver-Stokes equations form a closed set of equations but in the derivations of
the URANS equations, two new terms in the momentum and the energy equations appear that
need to be modeled. Those are the turbulent stress terms, which include the unresolved high-
frequency components of velocity and temperature, respectively, and are therefore related to
unresolved turbulence.
These terms need to be modeled, and the following section presents the turbulent
model’s formulation that was used in this thesis.
2.6.4.3 The" k - Ɛ " Turbulence model:
The "k - Ɛ " model by Launder and Spalding [17] is an eddy viscosity model and, like
other eddy viscosity models, it is based on the Boussinesq assumption that proposes a linear
relationship between the turbulent stresses and the mean flow gradients through an eddy
viscosity, µt, also known as the turbulent viscosity. Based on this assumption, the turbulent
stress tensor can be modeled as
~
~ 2 ∂ uk 2 ~
(
τ ij =μ t 2 S ij −
3 ∂ xk )
δ ij − ρ́ k δ ij
3 (2.29)

and the turbulent heat flux tensor is model is based on a simple diffusion term:
 μt ∂ ~Tt
q j=C p (2.30)
p rt ∂ x j
where Pt is the turbulent Prandtl number. The k - Ɛ turbulence model gives the turbulent
viscosity as
~2
k
μt =C μ ρ́ ~
ε (2.31)

~
The modeled turbulent kinetic energy, k , and dissipation ~ε, are obtained by solving their
respective modeled transport equations:

∂ ρ́ ~k ∂ μ t ∂ ~k
+
∂t ∂ x i (~
ρ́ k u j− μ+
( ) )
σk ∂ x j
=Pk − ρ́ ~ε
(2.32)

And,

∂ ρ́ ~ε ∂ μ ∂ ~ε
+
∂t ∂ xj ( ( ) )
~
ρ́ ~ε u j− μ+ t
σε ∂x j
=¿
(2.33)

Pk is the turbulence production term and Cµ, σ k, σ ε, C ε 1 and C ε 2are model constants. Regarding
boundary conditions, some relatively low values are usually prescribed for k and ~ε at inflow
boundaries at walls, k = 0 and ∂ ε /∂ y = 0 where y is the direction normal to the wall. The
production of turbulence is modeled according to

2 ∂~uk 2 ∂ ~ui
((
P k = μt 2 ~
S ij −
3 ∂ xk )
δ ij − ρ́ ~k δ ij
3 )∂ xi (2.34)

and, in the current work, a realizable constraint is applied to the turbulent viscosity:
~2 C ρ́ ~k
k
μt =min(¿ C μ ρ́ ~ , ~r ~ ) ¿ (2.35)
ε √ sij s ij
where Cr is constant.
A standard log-law wall function was used in all simulations for two reasons. Firstly, the
linearized solver proved to be unstable when a low-Reynolds number turbulence model was
used, most likely because of high gradients near the wall. Secondly, the separation point was
accurately predicted using the wall function so its use was continued in the DES and URANS
simulations to save computational time. A low-Re turbulence model requires a fine enough grid
resolution in the wall-normal direction to be able to resolve the flow gradients in the boundary
layer all the way to the wall, requiring the first grid point to be at y+ < 1. Wall functions serve
as a compromise between accuracy and computational cost and only require the first grid point
outside of the wall to be located in the logarithmic layer [17]; the first grid point can be
between y+ > 30 and y+ < 100. When the log-law wall function is used, it is assumed that the
velocity profile follows the log-law, as in a fully developed boundary layer10:
~
ut 1
= ln ¿¿ ¿
ur k (2.36)

where κ is the Von Kárman constant, B is a constant. uτ = p ¿) is the friction velocity, τ w is the
wall shear stress and
ρ́ y ur
+¿= ¿
μ
y (2.37)

is the normalized wall distance. The log-law wall function is activated if y+ & 10:7 and, when
activated, the wall shear stress, Dw , is adjusted to fulfill Equation 2.21, see Davidson [17] for
more detail. When the turbulence model is activated, pressure p and temperature T~ are
computed according to:
~
u ~u ~
(
ṕ=(γ −1) ρ́ ~e 0− k k −k
2 ) (2.38)

And,

1 ~ ~ u ~
u ~
~
T=
Cv 2 (
e0 − k k − k ) (2.39)

Where
γ =C p ∕ C ν
(2.40)

2.6.5 Turbulent combustion model:

In order to model combustion, the Eddy Dissipation Model (EDM) has been selected as
the best compromise between accuracy and computational effort. Being the EDM essentially
based on the hypothesis of fast chemistry, the influence of the chemical system used by the
model is particularly related to the number of chemical species used. Following the analysis
described in [17] the following model has been selected to simulate the paraffin/N2O
combustion:

151N2O+H2+25C2H4
147N2+41H2O+28CO+22CO2+10OH+10O2+8NO+H2
as reported in [4], if a lower number of chemical products is used for the combustion reaction,
the flame temperature can become higher than the flame adiabatic temperature, which results in
unphysical solutions. This formulation of the chemical reaction is used to model a one-phase
combustion process, where all the reactants are in a gaseous state.
2.6.6 Discretization methods:
The stability of the selected discretization is generally established numerically rather than
analytically as with simple linear problems. Special care must also be taken to ensure that the
discretization handles discontinuous solutions gracefully. The Euler equations and Naiver–
Stokes equations both admit shocks and contact surfaces.
Some of the discretization methods being used are:
Finite volume method:
The finite volume method (FVM) is a common approach used in CFD codes, as it has an
advantage in memory usage and solution speed, especially for large problems, high Reynolds
number turbulent flows, and source term dominated flows (like combustion).[17]
In the finite volume method, the governing partial differential equations (typically the
Naiver-Stokes equations, the mass, and energy conservation equations, and the turbulence
equations) are recast in a conservative form, and then solved over discrete control volumes. This
discretization guarantees the conservation of fluxes through a particular control volume. The
finite volume equation yields governing equations in the form,
∂
∭ Q ⅆ V + ∬ F ⅆ A=0
∂t (2.41)

where Q is the vector of conserved variables F is the vector of fluxes V is the volume of the
control volume element, and A is the surface area of the control volume element.
Finite element method:
The finite element method (FEM) is used in structural analysis of solids but is also
applicable to fluids. However, the FEM formulation requires special care to ensure a
conservative solution. The FEM formulation has been adapted for use with fluid dynamics
governing equations. Although FEM must be carefully formulated to be conservative, it is much
more stable than the finite volume approach. [17] However, FEM can require more memory and
has slower solution times than FVM.
In this method, a weighted residual equation is formed:
Ri= ∭ W i Q ⅆ V e
(2.42)
where R is the equation residual at an element vertex I, Q is the conservation equation expressed
on an element basis, W is the weight factor, and V is the volume of the element.
Finite difference method:
The finite difference method (FDM) has historical importance and is simple to program.
It is currently only used in a few specialized codes, which handle complex geometry with high
accuracy and efficiency by using embedded boundaries or overlapping grids (with the solution
interpolated across each grid).
∂Q ∂ F ∂G ∂ H
+ + + =0
∂t ∂ x ∂ y ∂z (2.43)

where Q is the vector of conserved variables, and F, G, and H are the fluxes in the x, y, and z
directions respectively.
The use of computational flow dynamics (CFD) is a numerical methodology for solving
the governing equations of fluid flow. The governing equations of fluid flow are partial
differential equations; when discretized on a mesh, they transform into algebraic equations which
can be solved by a finite-difference/finite-volume algorithm. the following sections will briefly
describe the governing equations, turbulence model, flow conditions, and properties employed in
this work.
2.6.7 Modelling species:
ANSYS Fluent can model the mixing and transport of chemical species by solving
conservation equations describing convection, diffusion, and reaction sources for each
component species. Multiple simultaneous chemical reactions can be modeled, with reactions
occurring in the fluid phase (volumetric reactions) and/or on wall or particle surfaces, and in the
porous region. And the Eddy Dissipation model had been used.
Species Transport Equations:
When you choose to solve conservation equations for chemical species, ANSYS Fluent
predicts the local mass fraction of each species, through the solution of a convection-diffusion
equation for the i thspecies. This conservation equation takes the following general form:

∂(ρ Y i)
+∇ ¿ (2.44)
∂t

where Ri is the net rate of production of species i by chemical reaction and Siis the rate of
creation by addition from the dispersed phase plus any user-defined sources.

2.6.8 Eddy Dissipation Model:

Under some combustion conditions, fuels burn quickly and the overall rate of reaction is
controlled by turbulent mixing. In high-temperature non-premixed flames, for example,
turbulence slowly convects/mixes fuel and oxidizer into the reaction zones where they burn
quickly. In certain premixed flames, the turbulence slowly convects/mixes cold reactants and hot
products into the reaction zones, where the reaction occurs rapidly. In such cases, one
approximation is to assume the combustion is mixing limited, allowing neglect of the complex
chemical kinetic rates and instead of assuming instantaneous burn upon mixing. For mixed-is-
burned approximation, ANSYS Fluent provides a turbulence-chemistry interaction model, based
on the work of Magnussen and Hjertager called the eddy-dissipation model. With this model, the
net rate of production of species due to reaction is given by the smaller (that is, limiting value) of
the two expressions below:

ε Y
Ri , r=v 'i ,r M w ,i Aρ mⅈ n R ' R
k (
v R ,r M w , R ) (2.45)

ε ∑Y p
Ri , r=v 'i ,r M w ,i ABρ
k
( N
'
∑ M w , j v ' j ,r
l
) (2.46)

where,
Y p= the mass fraction of any product species, P
Y R=the mass fraction of a particular reactant, R
A= an empirical constant equal to 4.0.
B= an empirical constant equal to 0.5.

2.7 Numerical investigation Method:

The RANS equations and turbulence models (as well as the radiative transfer model)
create a system of seven equations that need to be solved numerically. An analytical solution for
these equations is hard, therefore, an iterative numerical solution method is used on a mesh to
approximate the partial differential equations into approximate algebraic equations. The
linearized algebraic equations iteratively converge to the nonlinear solutions by employing a
suitable algorithm built-in ANSYS Fluent. A convergence criterion is specified to achieve
acceptable accuracy. When all the flow properties in all cells of the mesh reach the convergence
criteria, the solution is considered “converged” and the iterative process ends. And this is
required to be less than 0.001 which is ANSYS Fluent default parameter to achieve accuracy.
 Chapter three
Analysis And Calculation
3.1 Introduction:
In order to design a rocket nozzle with efficient thrust characteristics, supersonic flows
must be achieved. This requires the design of a converging-diverging nozzle. The goal of a
converging-diverging nozzle is to accelerate the combustion gases produced by burning
propellants to Mach 1 at the throat. This is the physical velocity limit in a converging nozzle
which is known as "choked" flow. Once supersonic velocities are achieved at the throat, the fluid
will actually accelerate through the diverging section of the nozzle.
 A favorable pressure gradient must be exerted across the nozzle. At design conditions, the
backpressure is similar to nozzle exit pressure and flow accelerates throughout the nozzle. But
backpressure is varying with changing altitude as the rocket is flight.
This section is dealing with the design two-type rocket nozzle Conical nozzles are self-
explanatory in that they are a cone shape defined by their half angle. And Bell shape has a
parabolic shape the method of characteristics (MOC) is used to design bell nozzles. The
streamlines of the expanding gas are calculated to produce a bell-shaped nozzle.
In the nozzle design, the most important part is the divergent section because it plays a
very important role in converting presser energy into kinetic energy. Some consideration, as
follows, should be satisfied in nozzle design.
o shortest possible nozzle length for minimum space envelope, weight, wall friction losses,
and Colling requirements.
o minimum separation and turbulence losses within the nozzle.
o uniform, parallel, axial gas flow at the nozzle exit for maximum momentum vector.
o Easy to manufacturing.
Assumption:
1. For simplicity, the combustion gas is assumed to be an ideal gas.
2. The gas flow is isentropic (constant entropy), frictionless, and adiabatic (it is little or no
heat gained or lost).
3. The gas flow is constant (steady) during the period of the propellant burn.
4. The gas flow is along a straight line from the gas inlet to the exhaust gas exit (along the
nozzle’s axis of symmetry.
5. The gas flow behavior is compressible since the flow is at very high velocities.

3.2 Design of Conical nozzle:

Conical nozzle contour is the simplest contour which provides satisfactory in most
respects, was used almost in space application. a conical nozzle allows ease of manufacture and
flexibility of converting an existing design to a higher or lower expansion area ratio without
major redesign.
The design typical conical nozzle is based on Hunter [18] who designed a 2D convergent-
divergent axisymmetric conical nozzle shown in Fig. 4.12,throat section has the contour of a
circular arc with a radius R. Throat nominal area of A t = 2.785 × 10-3 m2 and the area ratio of
Ae/At = 1.79 The half-angle of the nozzle convergent cone section can range from 20 to 45 deg.
The divergent cone half-angle α varies from approximately 10 to 18 deg. The following equation
can express the length of the conical nozzle:

R t ( √ ϵ−1 )+ R (sec α −1)

Ln=
tan a (3.1)
 The conical nozzle with an (11o) divergent half-angle has become almost a standard because
it is a good compromise based on weight, length, and performance.

Fig. 3.1 Conical nozzle design points.

The design property is shown in the table

Table 3.1 Design properties of Conical Nozzle.

Parameter value unit
Inlet diameter 1.386 m
Throat diameter 0.55 m
Outlet diameter 0.972 m
Area ratio 1.79 -
Total length 4.55 m
Divergent length 1.89 m
Divergent angle 11.01 o

3.3 Design of Bell nozzle:

The bell-shaped or contour nozzle is probably the most commonly used shaped rocket
engine nozzle. The Method of Characteristics (MOC) is a numerical procedure appropriate for
solving two-dimensional compressible flow problems [16]. The nozzle designed with MOC is
able to achieve uniform exit flow conditions [3]. The bell-shaped or contour nozzle had a high
angle expansion section (20 to 50°) right behind the nozzle throat; this was followed by a gradual
reversal of the nozzle contour slope so that the divergence angle was smaller at the nozzle exit
(usually less than a 10° half-angle) [19].
To carry out numerical simulation of a bell type nozzle, computer code is developed in the
MATLAB program. The program is present in the appendix (1). The code uses MOC and the
stream function to define the higher efficiency nozzle contours for isentropic, irrotational
supersonic flows of any working fluid for any user-defined exit Mach number. The program
starts with reading data from the user, such as enter the values of exit Mach number, throat
radius, the number of divisions required at the initial expansion curve, and the value of specific
heat (γ). In the next step, the program evaluates the angle of divergence for the initial expansion
contour downstream of the throat based on the exit Mach number that the user entered. After
evaluating the values of θ, x, y, ν, M, and µ for the points chosen on the initial curvature; using
these values, the program evaluates all the properties on the intersection points of the
characteristic lines based on theoretical formulas. At last, one can get the x, y coordinates of the
divergent section of the nozzle based on obtained properties θ, ν, and M. as shown in the Fig..
The present design is carried out for an outlet diameter of 0.972 m and throat diameter 0.55 m.
The details of design parameters, obtained from the program, for flow (Fluid: Combustion gases
at 400K).
The Fig. 3.2 shows the design of the 2D axisymmetric (only upper part) divergent section
in the countered or bell-shaped nozzle. The design was done using MATLAB vertex at end of
the blue line present the nozzle wall this vertex takes as X, Y point.

Fig. 3.2 Wall contour nozzle in MATLAB.

the design property is shown in the table below:

Table 3.2 Design properties of Bell Nozzle

Parameter value unit
Inlet diameter 1.386 m
Throat diameter 0.55 m
Outlet diameter 0.972 m
Area ratio 1.79 -
Total length 4.55 m
Divergent length 1.89 m
Divergent angle 11.01 o

3.4 The methodology of comparing:

It is known that the analytical solution for the nozzle gives the same value for all types of
converging-diverging nozzles and that is not correct, also assuming that the working fluid is air
and this assumption has a degree of error.
In this chapter, we will design and study two types of converging-diverging nozzles
Conical nozzle and Bell-shaped or contour nozzle Fig. 3.3 all nozzle at the same area ratio 1.79
(as highlighted) and same total length.

Table 3.3 Designs parameters for the two nozzles.

Comparing is done here with the same design condition see table and same operating condition.
Parameter Conical nozzle Bell nozzle
Inlet diameter 1.386 m 1.386 m
Throat diameter 0.55 m 0.55 m
Outlet diameter 0.972 m 0.972 m
Area ratio 1.79 1.79
Total length 4.55 m 4.55 m
Divergent length 1.89 m 1.89 m
The most important parameter for study and analysis hybrid rocket nozzle is nozzle
pressure ratio (NPR) it is the ratio between the pressure at the nozzle inlet (P o) and the pressure
at the nozzle outlet (Po) it controls the formulation of shocks, their strength, the shape and type of
flow (over and under expansion), and little influence on Mach number, so the parameter (NPR)
will be mentioned a lot in the research.
From Table 3.4 below. obviously, that area ratio plays an important role in changing
Mach number for a nozzle has a 1.79 area ratio as highlighted Mach number =2.07.
 Table 3.4 Properties against Mach number

Computational Fluid Dynamics (CFD) was used for analysis by using ANSYS Fluent
17.2 for modeling and simulation and to calculate the actual value of the Mach number for each
type of nozzle.

Fig. 3.3 The two types of nozzles (Conical & Bell).

The physical space is transformed into a computational space in order to promote the
numerical efficiency, and the physical boundary conditions are applied conveniently using the
following coordinates transformation
η = η (x, y), ξ = ξ (x, y ) (3.2)
After the coordinate transformation, equation (3.2) becomes:
^ ∂
∂Q ∂ ^ ^
+ F−^
^ F v + G− G v =0
∂t ∂ξ ∂η (3.3)

^
WhereF G

^ −1 Q , F =J −1 ξ x F+ ξ y G ,
^
Q=J
G
^ () ( η x F+ η y G) (3.4)

F
^v −1 ξ x F v + ξ y G v

( ) (
^
Gv
=J
η x F v+ η y G v ) (3.5)

3.5 The effect of Length, Temperature, and Area ratio changes on Mach
number:
It is known that the design of the conical nozzle is related to design equations and the
length of the nozzle is calculated by the equation (3.1), as a result, the part affecting the flow
characteristics is the divergent section and the convergent section does not have any effect, so the
length of divergent section will only change without changing the convergent section all other
design condition will remain constant like inlet, throat, and outlet diameter, area ratio = 1.79 for
all designs.
In this study the length divergent section will be changed four times, twice smaller 3.55,4.0
m, and twice greater than the original length 5.0, 5.55 m, note that the original length of the
divergent section is 4.55 m.

3.5 4.0 4.55 5.0 5.55

5

In this section two, temperature effects will be discussed, to see how it effects on Mach
number and other properties, the design temperature (Combustion gases temperature) on which
the nozzle is designed is = 400 K also it will be increased to know the amount of change that the
temperature change causes.
In section three the effect of changing the nozzle Area ratio on the other flow parameter.
This investigation happens with a constant of the other nozzle design parameter the area ratios
are taken in this study is 1.79 this is the design area ratio 2.3 and 1.55 was used here to see the
effect of this change.

3.6 Boundary condition and flow domain:

In order to run a simulation of the flow in with supersonic nozzles a suitable boundary
condition most be used as shown in the Fig. 3.4. And for more accurate and realistic simulation
and to find out the effects that occur behind the nozzle an addition was made, nozzle flow
domain is a computational domain added after nozzle for more accurate and realistic simulation.
The nozzles must be built virtually so that a mesh can be generated in the fluid region.
By entering the points which define the nozzle’s contour we imported these points into Gambit,
to account for the effects of supersonic nozzle exhaust and to illustrate the behavior of the flow
respect to the ambient conditions a large open area is modeled both upstream and downstream of
the nozzle. The nozzle inlet section is extended upstream of the nozzle to properly model the
stagnation conditions at the inlet section. Gambit is a mesh generating program used to mesh the
whole fluid domain of the simulation. All points are connected to produce a 2D axisymmetric
virtual geometry Fig. 3.2 shows the typical geometry.

Fig. 3.4 Computational domain and boundary conditions.

Table 3.5 Boundary condition.

The choice of boundary condition should be made with high precision because the
accuracy of the solution depends on the choice of boundary condition. The boundary condition is
used here is shown in the table.
Boundary Type Condition
Nozzle inlet Pressure inlet P0 = NPR × Pa, T0
Nozzle wall Wall No-slip and adiabatic
Ambient inflow Pressure inlet P0 = KPa, T0 =300 K
Symmetry symmetry 2D Planar
Outflow Pressure outlet Pout = KPa
 First of all, the boundary condition is used here is "Pressure inlet" at the nozzle inlet
which is define an inflow condition based on the known pressure value (P 0) at the inlet, No-slip
and adiabatic wall condition is imposed at the "Nozzle wall" of the nozzle, Ambient inflow is
Pressure inlet to but the value of pressure is equal to ambient pressure, the symmetry boundary
condition is used to reduces the computational calculations which lead to reducing the solution
time, Pressure outlet is define an outflow condition based on the flow pressure (P) at the outlet.
3.7 Numerical Mesh:
The used mesh is 2D, structured, and axisymmetric. The region near the axis has smaller
cells to capture details of the plume pattern. The throat region is the densest part of the mesh see
Fig. 3.5. In the region near the wall, the gradient of quantities is considerably high and requires
fine grids close to the wall to capture the changes in quantities. The boundary layer consists of 20
cells with a total thickness of ~0.6mm.fig illustrates mesh generation.

Fig. 3.5 a closer look at the meshed nozzle.

3.7.1 Grid independence study:

Having determined the geometry of the nozzle from the preliminary and detailed design
phases, the Gambit Mesher was employed to model this geometry, treating it as an axisymmetric
problem. Initially, a relatively coarse structured mesh was generated consisting of quadrilateral
elements, but having sufficient density to resolve the important geometry of the nozzle. To
facilitate mesh sensitivity analysis, four meshes were considered: standard, coarse, fine, and very
fine. All four were structured grids consisting of the quadrilateral element. These are illustrated
in Fig. 3.6a,3.6b, 3.6c,3.6d respectively.
 Fig. 3.6a Standard mesh.

Fig. 3.6b Coarse mesh.

Fig. 3.6c Fine mesh.

Fig. 3.6d Very fine mesh.

The standard mesh consisted of 16925cells, with a maximum equisize skew of 0.125655
and Minimum Orthogonal Quality of 0.874345. The coarse mesh included 4740cells, with a
maximum equisize skew of 0.141282 and Minimum Orthogonal Quality of 0.858718. The fine
mesh included 39000 cells, with a maximum equisize skew of 0.127948 and Minimum
Orthogonal Quality of 0.872052. The very fine mesh included 64000 cells, with a maximum
equisize skew of 0.127948 and Minimum Orthogonal Quality of 0.858718.

For the purposes of verifying the nozzle mesh the plot of static pressure in all adapted
meshes at the nozzle wall was plotted as in Fig. 3.7. This Fig. explain different type of mesh
from coarser to a very fine mesh.

Fig. 3.7 Mesh Comparation.

From Fig. 3.7 all mesh seemed to act as one but as soon as the different at the position 3.5
appears, the coarse mesh continues to decrease, the standard mesh increases first but the rest of
the meshes increase at the same time and they act as one, and we have found that the coarser
mesh does not work for this type of application with the comparison of another type and another
type acts as same. For this reason, we well used the fine mesh in this research. Because it has a
reasonable number of cells.
Now that the geometry has meshed, it can be imported into ANSYS FLUENT, the fluid
flow simulation program. Once imported, the solver type, material and properties, operating
conditions, and boundary conditions must all be defined. In many cases, the time step, controlled
by the Courant number, must be reduced or modulated throughout the simulation to facilitate
convergence.

3.8 CFD Setup:

The ANSYS Fluent 17.5 CFD package was used, employing the control volume
technique to solve the governing Naiver-Stokes equations of continuity, momentum, and energy.
the steady-state solver was used, and the axial symmetry of the problem was invoked to remove
the need to solve the full 3-dimensional equations and greatly reduce the computational
requirements. Due to the fact that the flow of greatest interest was supersonic, FLUENT’s
“Density-based” (or coupled) solver was used, employing the energy equation as the transport
equation for temperature, the continuity equation as the transport equation for density, and the
ideal gas law for the computation of the pressure.
Fluent offers the choice of either a Density-based implicit solver or a Density-based
explicit solver, each relating to how the governing equations are linearized. In the implicit
scheme, the unknown value of a given variable in each cell is computed from a relation that
includes both existing and unknown values from the neighboring cells. Since this means that
each unknown will appear in multiple equations the system must be solved simultaneously. This
approach offers broader stability characteristics and tends to allow a converged steady-state
solution to obtain much faster. The scheme is also unconditionally stable regardless of the
courant number chosen, though lower courant numbers can be required during startup when the
change in the solution is highly nonlinear.
Although being comparatively memory-intensive, the method’s stability characteristics
motivated its selection for this thesis. As the nozzle flow introduces the added complexity of
supersonic flow and the shock waves this implies, a second-order accurate solution was desired,
since first-order solutions tend not to adequately resolve these regions. Accordingly, second-
order up winding was used for the flow, and the QUICK scheme was employed for the modified
turbulent viscosity. As the nozzle flow introduces the added complexity of supersonic flow and
the shock waves this implies, a second-order accurate solution was desired, since first-order
solutions tend not to adequately resolve these regions. Accordingly, second-order up winding
was used for the flow, and the QUICK scheme was employed for the modified turbulent
viscosity.

Table 3.6 CFD setup.

Set-up
Solver Type Density-Based
Space 2D
Velocity Formation: Absolute
Time steady
2D space Axisymmetric
Models
 Energy checked
Viscous k-epsilon (2 Equ)
k-epsilon Model realizable
Near wall treatment Enhanced wall
treatment
Materials
Name air
properties density Ideal gas
cp 1006.43 J/kg*K
Molecular Weight: 28.966 kg/kmol
viscosity Sutherland law
Operating
Conditions:
Pressure: Operating Pressure: 0 Pa
Gravity: Not Checked
Reference Pressure X(m): 0
Location:
Y(m): 0
Boundary
condition
Pressure Inlet: Gauge Total Pressure: 1000000 Pa
Supersonic/Initial 1000000 Pa
Gauge
Pressure:
Direction Normal to boundary
Specification Method:
Total Temperature: 400K
Turbulent intensity 5
(%)
Turbulent viscosity 10
ratio
Pressure Outlet: Gauge Pressure:
Backflow Direction Normal to Boundary
Specification
Method:
Average pressure Not Checked
specification
Target Mass- flow Not Checked
Rate:
Backflow Total 300 k
Temperature:
Backflow Turbulent 5
intensity (%)
Backflow Turbulent 10
viscosity ratio
 Ambient Pressure Gauge Total Pressure:
Inlet:
Supersonic/Initial
Gauge
Pressure:
Direction Normal to boundary
Specification Method:
Total Temperature: 300
Turbulent intensity 5
(%)
Turbulent viscosity 10
ratio
Reference values
Compute from Ambient pressure inlet
Solution methods
formulation Implicit
Flux type Roe-FDS
Spatial Gradient Least squares cell-
discretization based
Flow Second-order upwind
Turbulent kinetic Second-order upwind
energy
Turbulent dissipation Second-order upwind
rate
Solution control
Courant number 5
Under relaxation
factors
Turbulent kinetic 0.8
energy
Turbulent dissipation 0.8
rate
Turbulent viscosity 1
solid 1
Monitors
Residuals
Continuity 1x10^-6
X- velocity 1x10^-6
Y- velocity 1x10^-6
Energy 1x10^-6
k 1x10^-6
Epsilon 1x10^-6
Solution
initialization
Initialization method Hybrid initialization
3.9 Validation:
To ensure the validity and accuracy of the employed numerical investigation, the
validation step is very important, a rigorous comparison between the present computation and
Hunter’s [18] experimental data obtained for a 2D convergent-divergent nozzle at the nozzle
pressure ratio (NPR) of 2.4 is implemented. The test nozzle used in Hunter’s investigation is a
subscale, non-axisymmetric, 2D convergent-divergent nozzle with the throat nominal area of A t
= 2.785 × 10-3 m2 and the area ratio of Ae/At = 1.79. The geometric details of the 2D convergent
divergent nozzle are shown in Fig. 3.8.
The schlieren flow visualization at NPR = 2.4 for the studied 2D convergent-divergent
nozzle is illustrated in Fig. 4. It comprises a well-defined lambda shock structure and a free
shock flow that begins from the leading branch of the lambda shock structure and extends
downstream. The abscissa is normalized by the throat position (x t) measured away from the
nozzle inlet. Fig. 5 includes the diameter of the Mach disk and two angles, one between the
nozzle wall and the oblique shock and the other between the oblique and reflected shocks, that
all have been measured experimentally by Hunter [18].
 Fig. 3.8 Schlieren flow visualization at NPR = 2.4 for the 2D-CD nozzle. Photographed by
Hunter [18], NASA Langley Research Center. Compared with our CFD result.
As shown in Fig. 3.8 the comparison between our result Fig. 3.4A colored nozzle and hunter
experimental result Fig. 3.4B black and white the nozzle is shown axisymmetric (only upper part
is showed) because in engineering problem its make the solution easy and reduce the desired
time for solving. In the Hunter result, the flow Plame does not exist because the device using in
visualization is limited to nozzle only.
By plot, the result to be more clear as to show in Fig. 3.9 the experimental data form hunter
[18] in the red line and the CFD data in the black line clearly show that the experimental and
CFD data is similar.

Fig. 3.9 Comparison of Experimental and Computational Pressure Data.

 Chapter four
Results and Discussion
4.1 Air Nozzle Result:
The first experiment was performed on air without any combustion and that was to see
whether the real manufacturing nozzle would yield to 1 Mach or not. In the convergent-divergent
nozzle, the Mach number must be greater than unity (M>1) to achieve the desired speed
(supersonic speed). And the major difference between the nozzle working on-air and working on
combustion gases is the temperature difference and density difference the result of both will be
presented.
The CFD analysis of a rocket engine nozzle has been conducted to understand the
phenomena of subsonic flow through it at various divergent lengths and areas. A two
-dimensional axisymmetric model is used for the analysis and the governing equations were
solved using the finite-volume method in ANSYS Fluent ® software. The variations in the
parameters like the Mach number, static pressure, and static temperate are being analyzed. The
phenomena of oblique shock are visualized and the normal shock is also illustrated.
Computational Fluid Dynamics (CFD) is an engineering tool that assists experimentation. Its
scope is not limited to fluid dynamics; CFD could be applied to any process which involves
transport phenomena with it. To solve an engineering problem, we can make use of various
methods like the analytical method, experimental methods using prototypes. The analytical
method is very complicated and difficult. The experimental methods are very costly. If any errors
in the design were detected during the prototype testing, another prototype is to be made
clarifying all the errors and again tested. This is a time-consuming as well as a cost consuming
process. Flow instabilities might be created inside the nozzle due to the formation of shocks
which reduce the exit Mach number as well as the thrust of the engine. This could be eliminated
by choosing the right design.

4.2 Result and Discussion:

The result we have got using ANSYS Fluent it’s an analysis result of the conical nozzle,
bell nozzle, the effect of changing length, and temperature on Mach number.
4.2.1 Conical Nozzle Result:
The design typical conical nozzle is based on Hunter [18] who designed 2D convergent-
divergent axisymmetric conical nozzle. The result that was Hunter [18] has come to are a
practical (experimental) result and we here doing the practice using CFD simulation. The
boundary condition using for simulation is:
 Inlet pressure = 1 MPa
 Inlet temperature= 400 K
  Outlet temperature = 300 K

The Fig. (4.1) show the contour of Mach number the results show, for NPR = 2.412, a
nozzle shock with a pronounced lambda foot structure and fully detached separation layer
extended from the leading lambda shock in downstream. The separation region formed
downstream can be considered as a consequence of the adverse pressure gradient through the
shock, which forces the incoming boundary layer to separate. The oblique shock structures are of
the weak type resulting in low supersonic flow downstream while the flow immediately past the
Mach stem is subsonic.

Free Get Exhaust

Mach Disc

Oblique
shocks

Fig. 4.1a Mach number contour for a complete nozzle.

Fig. 4.1b Mach number contour for an axisymmetric conical nozzle.

 The Fig. 4.2 show the contour of static pressure, Fluent defines the total pressure as a gauge
pressure with respect to operating pressure. From Fig. 4.2, the static pressure sho wn at the
nozzle varies from 1000 kPa to 416 kPa. From the inlet throughout to the exit of the nozzle,
pressure drops significantly. Until the shock pressure becomes very high. And the Mach number
on the other side falls down the cause of shock apparent. Shock wave in rockets applications is
not desirable and preferably avid it.

Fig. 4.2 Contou r of static pressure.

The Fig. 4.3 show the contour of density, the density of gases com from combustion is
higher and continually decrees until the shock wave see that a slight increase.

Fig. 4.3 Contour of Density.

As mentioned previously, the exhaust flow pattern is dependent on whether the flow is over-
expanded or under-expanded [20]. In over-expanded flow co nditions, the ambient pressure is
higher than the nozzle exit pressure, while for under-expanded flow conditions the ambient
pressure is less than the nozzle exit pressure. In the case of over-expanded conditions,
the exhaust flow adapts to the ambient through a system of oblique shocks and expansion waves.
However, in the case of under expanded conditions, the exhaust gas continues its expansion
behind the nozzle exit. Fig. 4.4 shows a good qualitative comparison between the present
numerically predicted and the observed flow pattern downstream the nozzle exits for both flow
conditions.

a- Under expanded flow condition b- Over expanded flow condition

Fig. 4.4 Qualitative Comparison between predicted and observed flow
pattern downstream nozzle exit.

The above results indicate that the shock location inside the nozzle and the separation point
is affected by the NPR. Therefore, a wide range of NPRs has been simulated using the K-W
turbulence model. The predictions of the shock position as well as the separation point are
plotted and compared with the experimental measurements as shown in Fig. 3.8.The results
indicate that by increasing the NPR, the shock position as well as the separation point move
downstream. By increasing the NPR to be about five, it is found that the shock location does not
change nearby the nozzle exit.
The dimensionless shear stress (τ/τ in) distribution for different NPRs is plotted in Fig
(4.5), where sin is the shear stress at the nozzle inlet. For separated flow at low NPR, the results
indicate that the flow did not attach the nozzle surface and the free shear layer started at the
trailing lambda foot is completely detaching past the shock. As the separation point becomes
nearly near the effective nozzle exit, the lambda shock system adjusted to satisfy the continuity
of pressure and flow direction. This trend is indicated by a slight increase in the shear stress near
the nozzle exit. These results can also be seen in the Mach number graph, as shown in Fig. 4.1,
For NPR = 2.412, the shear stress plotted in Fig. 4.5, a positive velocity distribution is obtained
nearby the effective nozzle exit. The Mach contours, shown in Fig. 4.2, illustrate the existence of
a lambda shock system which consists of a normal shock wave at the nozzle core and an oblique
shock at the nozzle walls as a result of the boundary layer shock interaction.

Fig. 4.5 The computed shear stress distribution for different NPR.

4.2.2 Bell Nozzle Result:

In this type, we designed the nozzle by using MatLab through MOC the program designs
the convergent part of the nozzle that comes parabolic shape in bell nozzle as opposed to the
straight one in the conical nozzle. By using the same boundary condition used in the previous
type.
 Inlet pressure = 1 MPa
 Inlet temperature= 400 K
 Outlet temperature = 300 K
 The Fig. 4.6 show the contour of Mach number for NPR=2.412, The Mach number is
subsonic at the inlet After crossing the throat portion the flow becomes supersonic. The exit
Mach number has some variations. In the boundary layer, the flow becomes supersonic but at the
midpoint of the exit, the flow becomes less supersonic. Mach velocity vectors show the change
of Mach numbers with the velocity vectors. Velocity vectors show the change of Mach numbers.

Shock Diamond

Fig. 4.6a Contour of Mach number for a complete Bell nozzle type.

Fig. 4.6b Contour of Mach number for an axisymmetric Bell nozzle type.
 Noctic that from Fig. 4.6 shock diamond was appearance at the end of the nozzle shock
diamond or Mach diamond are formation of shock wave patterns that apper when the rocker
oprateed in an atomspher. The “diamonds” are actually a complex folw field made visible by
abrupt changes in local density and pressuer. Shock diamond form when the supersonic exaust
form a nozzle is slightly over expanded, meaning that the ststic pressuer of gasses exiting the
nozzle is less than the ambint pressuer, Fig. 4.6c show the actuall skock dimond comparing with
simulation one.

Fig. 4.6c Show the real shock dimoand [24].

The contour Fig. 4.7 show the velocity magnitude for a Bell nozzle type as shown clearly, it’s an
over-expanded flow condition.

Fig. 4.7 Contour of Velocity magnitude for an axisymmetric Bell nozzle type.
 The Static temperature Fig. 4.8a is 400 K at the nozzle and varies gradually through the
axial distance. From the inlet to out temperature decreases. But for the total temperature Fig.
4.8b, it remains the same in the whole nozzle. But at the exit portion, the total temperature drops
slightly.

Fig. 4.8a Contour of Static temperature for an axisymmetric Bell nozzle type.

Fig. 4.8b Contour of the Total temperature for an axisymmetric Bell nozzle type.
For comparing to the type of nozzle at the same design parameter and same boundary
condition bell nozzle type is shows superior because the probability of shock appearance is less
than the conical nozzle. And for the same design parameter bell nozzle give a higher Mach
number than the conical nozzle.

4.2.3 The effect of Length, Temperature, and Area ratio changes on Mach
number:

Table 4.1 Design properties of Conical Nozzle.

 the study was done on conical nozzle type the original length of this nozzle is 4.55 m and
the other study length is 3.55, 4.0, 5.0 and 5.55 m the only length of the divergent part was
changed and the other design parameter remains constant as see in the table below:

Parameter value unit

Inlet diameter 1.386 m
Throat diameter 0.55 m
Outlet diameter 0.972 m
Area ratio 1.79 -
Total length 4.55 m
Divergent angle 11.01 o

The boundary condition is used for this study well remain constant as used above. we
changing the length and meshed the geometry and simulate in ANSYS Fluent and to make the
result easier to understand we used OriginPro to plot the result and make a comparison. As
shown in Fig. 4.9, Fig. 10 the result of Density, Static pressure, Total pressure, Mach number,
Static temperature, and Total temperature.

Fig. 4.9 Result of Density, Static pressure plots with different lengths.
 Fig. 4.10 Result of Total pressure and Static temperature, plots with
different lengths.

Fig. 4.11 show the plots of different nozzles parameter compared with different lengths. As
shown in the plot Fig.s any line represents the specific length of nozzle comparing with others, in
the density plot we find that the 5.55 nozzle length has a maximum value at point 0.6 of x/l x-
axis. The other important thing shows in plot Fig. 4.11 Mach number as show the 5.55 nozzle
length has a maximum value of Mach number and the 3.55 nozzle length have a minimum value
of Mach number. From that Fig. out the increase of length lead to an increase in Mach number,
and a decrease the length lead to a slight decrease in Mach number.

Fig. 4.11 Result of Mach number and Total temperature compared with different nozzle lengths.

Combustion temperature has a significant effect on flow patterns and nozzle efficiency
because it needs to be carefully controlled to give optimum performance for the nozzle. And to
what is the effect four temperature degree was compared the first 400K (the design temperature)
and 1000K, 1500K, 2000K and 3000K.
The plot Fig. 4.12 show the comparisons result between different temperature degree. Any
line represents a specific temperature with a specific parameter, Density, Static pressure, Total
pressure, Mach number, Static temperature, and Total temperature.

Fig. 4.12 Result of Mach number and Density compared with different nozzle temperature.

The result obtained from the plot Fig. 4.13 is show that the maximum value of density is
presented in the design temperature at T=400K and with continuous increase temperature leads
to a decrease the density and this is a bad Indicator because with a decrease in density the
volume of gases will increase. Fig. 4.13 shows the result of Static pressure and Total temperature
compared with different nozzle temperature.
 Fig. 4.13 Result of Static pressure and Total temperature compared with different nozzle
temperature.

The static pressure plot shown obviously the changing temperature does not affect in
pressure as seeing clearly all temperatures are identical. And with increasing temperature, the
Mach number slightly decrease.
In the end, the increase in temperature does not desirable because it needs special materials
to deal with this temperature.

Finally, the Fig. 4.14 represents the effect of changing nozzle Area ratios on its flow pattern
and flow characteristic. Three area ratios are chosen first the area ratio 1.79 which is the design
nozzle area ratio and 2.3, 1.55 to the effect of increasing and decrease. We notice from the plots
above the static pressure in beginning take different paths in area ratio values but soon as the
unified path takes place which happened in position 0.8.
 Fig. 4.14 Effect of nozzle Area ratios on flow properatic.

However, in changing area ratio values the same happens, by looking at the Fig. 4.14 the
area ratio 1.79 is the designed value and we see if we get lower than the designed value the static
pressure still decreasing to a value below all the values then it suddenly increases that sudden
decrease means a weak shock wave appears and it gets stronger whenever we get lower than the
designed value.
In both cases, we saw whenever we change the designed point or value shocks appear in
the nozzle and lead to change in the Mach number and lowering the efficiency of the nozzle.

4.2.4 Species Reaction Nozzle:

On the basis of the encouraging results obtained in cold gas conditions and taking into
account the good agreement between numerical and experimental results, two-dimensional CFD
computations have been used to verify the performances of the proposed nozzle concept in hot
firing conditions. paraffin wax as fuel and nitrox as an oxidizer was used with seven species (H2,
O2, OH, H, O, H2O, N2) have been used to. The same mesh, previously used for cold gas
computations, has been used for the hot gas cases. This choice is justified on the basis of
previous works, carried out with the same CFD tool, and by the fact that the original mesh is
already well refined, in particular in the area around the nozzle tip for the specific purpose of this
study. In addition, a sensitivity study has been performed on hot gas cases. Boundary conditions
have been set in terms of feeding conditions for the combustion chamber (propellant pressure
and temperature level, mixture ratio) and for the secondary flow (pressure and temperature
levels, pure Hydrogen flow). Ambient air, at sea level conditions, is taken into account as
surrounding fluid (pressure and temperature levels, chemical composition).

Steady-state and transient (start-up and shut-down) computations have been carried out. For
both types of computations, gaseous propellant inlet conditions have been chosen. Indeed,
despite the capability of the CPS-C code to handle variable thermodynamic properties over a
very large temperature and pressure ranges (including supercritical and transcortical
thermodynamic conditions), the computational complexity related to phase change mechanism
has been omitted, as these physical phenomena are not directly linked to the flow separation
issue. Indeed, the combustion process taking place in the combustion chamber is mostly
completed before the flow enters the sonic throat.
 Previous studies demonstrated that macroscopic combustion heterogeneities (due to injector
path or to mixture ratio distribution across the injection plate) can have a minor effect on the
local position of the flow separation line along the nozzle circumference. Nevertheless, this
aspect is covered in this study by the transient computation: the combustion pressure ranges from
few bars to the full-flowing conditions and the obtained results can be intended as the local worst
case for the effectiveness of the proposed concept.
Fig. 4.15 show the plots of static pressuer and Mach number compar between using the
compress air as working fluid and combustion gases as a working fluid we find that the
combustion gases gave us slitly law Mach number resultant form slity decrice in pressure in
point 5.5 at x/t axis. Combustion gasess nozzle flow give us the same result as contur because of
the the result shown as plot.

Fig. 4.15 Static pressuer and Mach number compare with two deffernt gasses.

Chapter Five
Conclusion And Recommendation
5.1 Conclusion:
The numerical simulations of compressible flow passing through a 2D convergent-divergent
nozzle with different geometries and different nozzle pressure ratios are presented. The nozzle is
assumed to have impermeable and adiabatic walls with a flow straightener in the upstream side
and is connected to a plenum surrounding the nozzle geometry and extended in the downstream
direction (adding Fairfield after nozzle) for more accuracy and realistic simulation. The predicted
results are obtained by solving the RANS equations for compressible flow in its conservative
form coupled with both the energy equation and the equation of a state. Serval lengths and
temperate degrees were applied.
We saw the pressure ratio plays an important role in Mach number, shocks waves, the value
of Mach number, and the efficiency in general. Using the paraffin wax as fuel and nitrox as an
oxidizer release a huge amount of heat and the temperature of the burnt gas was very high
because of that the material of the nozzle surface must be taken into account because the high
temperature of burnt gases can cause melting and leads to a big failure. We also concluded from
previous chapters if we increse the designed length it well lead to increase the Mach number but
cose to shock appernt mack that increasing uesless and decreasing the design length and area
ratios of the nozzle cause to shock waves appear and changinging the design temperuter lead to
decrease Mach number and decrease the efficiency and lowering the thrust which is not
preferred.
Based on the analysis of the results of numerical simulations the following concluding
remarks can be drawn:
In this geometry maximum Mach number is in the same place of minimum pressure.
Different types of an isentropic parameter such as temperature, pressure, Mach number
has been calculated. The pressure of fluid decreased and Mach is increased.
The converging-diverging geometry does not let Mach number to reach more than
unity(M=1) at the throat. Mach number continues rising after throat.
CFD calculations might have different outcomes with different schemes and this CFD
result of one case may not be the spitting image of the same case through a different type
of turbulent model.
For smaller diameters there are low chances of flow separation hence thrust exerted on
the body is larger in the case of a C-D nozzle with a smaller diameter than the larger
diameter.

5.2 Recommendation:
1- after designing and calculating the dimensions of the nozzle, engineers should be
accurate in manufacturing, because any change in this dimension can cause a serious
problem and lead to failure.
2- Pressure differences must take into account because with different altitudes pressure
changes, and the designers must take this consideration.
3- Depending on the fuel and combustion temperature the material of the wall surface is
chosen avoiding to melting of the nozzle surface.
4- Three-dimensional analyses should be considered in further steps.
References:
[1] Senthil, S. "Gas Dynamics and Jet Propulsion." (2009).
[2] G. SATYANARAYANA, 2. Ch. VARUN, 3. S.S. NAIDU, “CFD analysis C-D nozzle”, ACTA
TECHNICA CORVINIENSIS- Bulletin of engineering, Tome VI, ISSN 2067-3809, 2013
[3] Arjun Kundu, Devyanshu Prasad, Sarfraj Ahmed, “Effect of Exit Diameter on the performance of
Converging-Diverging Annular Nozzle using CFD”, International Journal of Innoavatine Research in
Science, Engineering and technology, Vol. 5, Issue 6, June 2016.
[4] Venkatesh.V, C Jaya pal Reddy, “Modeling and Simulation of Supersonic nozzle using
Computational Fluid Dynamics”, International Journal of Novel Research in Interdisciplinary Studies,
Vol. 2, Issue 6, pp: (16-27), 2015.
[5] Karna S. Patel, “Flow Analysis and Optimization of Supersonic Rocket Engine Nozzle at Various
Divergent Angle using Computational Fluid Dynamics (CFD)”, IOSR Journal of Mechanical and Civil
Engineering (IOSR-JMCE) e-ISSN: 2278-1684,p-ISSN: 2320-334X, Volume 11, Issue 6 Ver. IV
(NovDec. 2014), PP 01-10
[6] Nikhil D. Deshpande, Suyash S. Vidwans, Pratik R. Mahale, Rutuja S. Joshi, K. R. Jagtap,
“Theoretical and CFD analysis of De Laval Nozzle”, International Journal of Mechanical And Production
Engineering, ISSN: 2320-2092, Volume- 2, Issue-4 April, 2014
[7] Mohan Kumar G, Dominic Xavier Fernando, R. Muthu Kumar, “Design and Optimization of De
Laval Nozzle to prevent shock induced Flow Separation”, Advances in Aerospace Science and
Applications, ISSN 2277-3223 Volume 3, Number 2 (2013), pp. 119-124
[8] Sourabh Shriwas, Sachin Shah, Prashant Shah, Kalpit P. Kaurase, “Reduction Of Jet Noise In the
Aircraft Nozzle”, International Journal Of Research In Aeronautical And Mechanical Engineering ISSN
(ONLINE): 2321-3051 Vol.3 Issue.3, March 2015
[9] Vignesh. M , Sathish. k , Adrin Issai Arasu. L ,Dharun Lingam. K , and Sanal Kumar. V. R,
“Conceptual Design of Short Takeoff Supersonic Aircraft with Cold Flow Nozzles”, 3 rd International
Conference on Trends in Mechanical and Industrial Engineering (ICTMIE'2013) , Kuala Lumpur
(Malaysia), January 8-9, 2013
[10] R. Lárusson, N. Andersson, L.-E. Eriksson, and J. Östlund. Comparison of Eigenmode Extraction
Techniques for Separated Nozzle Flows. Proc. of the 50 th AIAA/ASME/SAE/ASEE Joint Propulsion
Conference (2014). July 28-30, Cleveland, Ohio, USA
[11] R. Lárusson, N. Andersson, L.-E. Eriksson, and J. Östlund. Linear Stability Analysis Using the
Arnoldi Eigenmode Extraction Technique Applied to Separated Nozzle Flow. Proc. of the 49th
AIAA/ASME/SAE/ASEE Joint Propulsion Conference (2013). July 14-17, San Jose, CA, USA
[12] Saleh B, Koglbauer G, Wendland M, Fischer J. Working fluids for low-temperature organic Rankine
cycles. Energy. 2007 Jul 1;32(7):1210-21.
[13] G. V. R. Rao. Exhaust Nozzle Contours for Optimum Thrust. Journal of Jet Propulsion 28 (1958)
[14] R. Lárusson, N. Andersson, and J. Östlund. Hybrid RANS-LES Simulation of Separated Nozzle
Flow. Proc. of the 52nd AIAA/ASME/SAE/ASEE Joint Propulsion Conference (2016). July 25-27, Salt
Lake City, UT, USA
[15] G. Hagemann, J. Alting, and D. Preclik. Scalability for Rocket Nozzle Flows Based on Subscale and
Full-Scale Testing. Journal of Propulsion and Power 19 (2003)
[16] R. Lárusson, H. E. Hafsteinsson, N. Andersson, and L.-E. Eriksson. Investigation of Supersonic Jet
Flow Using Modal Decomposition. Proc. of the 20th AIAA/CEAS Aeroacoustics Conference (2014).
June 16-20, Atlanta, GA, USA
[16] J. D. Anderson. Modern Compressible Flow with Historical Perspective. 3rd edition. McGraw-Hill,
2003
[17] Mukherjee S, Biswas G. Computational fluid dynamics. Alpha Science International Limited; 2014.
[18] C. A. Hunter. Experimental Investigation of Separated Nozzle Flows. Journal of Propulsion and
Power 20 (2004).
[20] K.M. Pandey and A.P. Singh., CFD Analysis of Conical Nozzle for Mach 3 at Various Angles of
Divergence with Fluent Software, International Journal of Chemical Engineering and Applications, Vol.
1, No. 2, August 2010, pp.179-185.
[21] The Space Report 2016. Space Foundation. 2016. url: http://www.spacefoundation.
org/sites/default/files/downloads/The_Space_Report_2014_Overview_ TOC_Exhibits.pdf.
[22] M. J. Zucrow. Aircraft and Missile Propulsion. John Wiley & Sons, 1958.
[23] S. J. Dick, ed. Historical studies in the societal impact of spaceflight. National Aeronautics and
Space Administration, 2015. url: https://www.nasa.gov/ connect/ebooks/index.html.
[24] Guentert, E. C. and Neumann, H. E., \Design of Axisymmetric Exhaust Nozzles By Method of
Characteristics Incorporating A Variable Isentropic Exponent," NASA TR R-33, 1959.
[25] Omer Musa. Design of supersonic nozzle using Method Of Characteristic, Master thesis 2013 ,
Sudan university of science and technology.
[26] F. R. Menter and M. Knutz. Adaptation of Eddy-Viscosity Turbulence Models to Unsteady
Separated Flow Behind Vehicles. Proc. of the Symposium on the aerodynamics of heavy vehicles: trucks,
buses and trains. (2002).
[27] C. Loh and K. Zaman. Numerical investigation of transonic resonance with a convergent-divergent
nozzle. AIAA Journal 40.12 (2002), pp. 2393–2401.
[28] L. Davidson. Fluid mechanics, turbulent flow and turbulence modeling. Chalmers University of
Technology. 2017. url: http://www.tfd.chalmers.se/~lada/postscript_files/solids-and-fluids_turbulent-
flow_turbulence-modelling.pdf.
[29] cfd
 APPENDIX A
MATLAB Supersonic Bell Nozzle Design Code

 THE MAIN PROGRAM:

clc; close all; clear all;

%INPUT VALUES
p_1 = xlsread('PARAMS.xlsx','PARAMETERS','B2'); %CHAMBER PRESSURE
T_1 = xlsread('PARAMS.xlsx','PARAMETERS','B3'); %CHAMBER TEMP
FT = xlsread('PARAMS.xlsx','PARAMETERS','B4'); %DESIRED THRUST OR....
m_dot = xlsread('PARAMS.xlsx','PARAMETERS','B5'); %DESIRED MASS FLOW
RATE....
ALT = xlsread('PARAMS.xlsx','PARAMETERS','B6'); %ALTITUDE
g = xlsread('PARAMS.xlsx','PARAMETERS','B7'); %GAMMA
R = xlsread('PARAMS.xlsx','PARAMETERS','B8'); %GAS CONSTANT

%% exit pressure
if (11000>ALT) && (ALT<25000)
T = -56.46; %C
p_o = 1000*(22.65*exp(1.73-0.000157*ALT));
elseif ALT>=25000
T = -131.21 + 0.00299*ALT ;
p_o = 1000*(2.488*((T+273.1)/216.6)^-11.388);
else
T = 15.04 - 0.00649*ALT;
p_o = 1000*(101.29*((T+273.1)/288.08)^5.256);
end

%% begin calculation
PR = p_o/p_1;
PR2 = (p_o/p_1)^((g-1)/g);
TT = (2*g*R*T_1)/(g-1);
p_t = ((2/(g+1))^(g/(g-1)))*2.068;
v_t = sqrt((2*g*R*T_1)/(g+1));
v_e = sqrt(TT*(1-PR2));

if m_dot==0
m_dot=FT/v_e;
elseif FT==0
 FT = m_dot/v_e;
else
fprintf('You can either set desired thrust OR mass flow rate')
end

T_e = T_1*(p_o/p_1)^((g-1)/g);
a_e = sqrt(g*R*T_e);

Me = v_e/a_e;

% MOC
TR = 0.559; %throat radius (mm)
RTOD = 180/pi;
DTOR = pi/180;
P = []; %x axis points
%% PM FUNCTION
A = sqrt((g+1)/(g-1));
B = (g-1)/(g+1);
v_PM = @(x) A*atan(sqrt(B*(x^2-1))) - atan(sqrt(x^2-1));

%% CALCULATE T_MAX, BREAK UP INTO DIVISIONS

T_max = 0.5*v_PM(Me)*RTOD;
DT = (90-T_max) - fix(90-T_max);
T(1) = DT*DTOR;
n = T_max*2;

for m = 2:n+1
T(m) = (DT + (m-1))*DTOR;
%Mach from T(i) using T(i) = v_PM (FALSE POSITION)
x_int = [1 1.01*Me];
func = @(x) T(m) - v_PM(x);
M(m) = fzero(func,x_int);
P(m) = 0 + TR*tan(T(m)); %X-AXIS POINTS
%RRSLOPES
RR(m) = -TR/P(m);
%LR slopes
LR(m) = tan(T(m)+asin(1/M(m)));
SL(m) = -RR(m);
end
%% PLOTTING
P(1) = [];
l = length(P);

for j = 1:l
P1 = [0 TR];
P2 = [P(j) 0];
plot(P2,P1,'k')
hold on
xlabel('CENTERLINE')
ylabel('RADIUS')
end
hold on;
LR(1) = []; RR(1) = [];
SL(1) = [];
F = RR(m-1);

for c = 1:length(P)-1
x(c) = (TR+SL(c)*P(c))/(SL(c)-F);
y(c) = F*x(c)+TR;
X_P = [P(c) x(c)];
Y_P = [0 y(c)];
plot(X_P,Y_P,'b');
end

hold on

%% FIRST WALL SECTION

TM = T_max*DTOR;
xw(1) = (TR+SL(1)*P(1))/(SL(1)-tan(TM));
yw(1) = tan(TM)*xw(1)+TR;
X_P2 = [P(1) xw];
Y_P2 = [P(2) yw];
plot(X_P2,Y_P2,'g');
%DIVIDE (delta slopes)
DTW = tan(TM)/(length(P)-1);
s(1) = tan(TM);
b(1) = TR;
v_PM(Me)

for k = 2:length(P)-1
s(k) = tan(TM)-(k-1)*DTW; %slope
b(k) = yw(k-1)-s(k)*xw(k-1); %y-int
xw(k) = (b(k)+SL(k)*P(k))/(SL(k)-s(k));
yw(k) = s(k)*xw(k)+b(k);
X_P3 = [x(k) xw(k)];
Y_P3 = [y(k) yw(k)];
plot(X_P3,Y_P3,'b');
fprintf ('%.4f %.4f 0\n',xw(k),yw(k))
end
hold on
LAST (PTS)
xf = (b(length(b))+SL(length(SL))*P(length(P)))/SL(length(SL))
yf = b(length(b))
X_F = [P(length(P)) xf]
Y_F = [0 yf]
plot(X_F,Y_F,'k');

xw = [0 xw];
yw = [TR yw];

xlswrite('PARAMS.xlsx',transpose(xw),'PTS','A1:A62');
xlswrite('PARAMS.xlsx',transpose(yw),'PTS','B1:B62');