CFD analysis of nozzle effect on jet formation

(1)

CFD analysis of nozzle effect on jet formation

Analysis of mixing performance, turbulence and flow instability

MUSTAFA GUDUCU

Master’s Thesis in Aeronautical Engineering Supervisor: Mireia Altimira

Examiner: Mihai Mihaescu

(2)

(3)

Abstract

This project is concerned with CFD simulations of jets is- sued from elliptical nozzles. The investigated jet flow in this project is turbulent flow emanating from microscopic nozzles into a combustion chamber. Jet flows are very common in engineering, medical and environmental applications and are for instance used in fuel injection systems, spray painting and drying. Jet flow devices are also very common in applications such as cutting, hydraulic drilling, cooling and heating. A better understanding of the flow phenomenons in jet flows are required in order to make these devices function and perform in a more efficient way.

The performance of diesel engines is strongly affected by the fuel spray, atomization and in turn the mixing process. This depends ultimately on the dimensions and geometry of the nozzle. The purpose of this project was therefore to investigate different elliptical nozzle geometries which also was compared to a conventional circular nozzle.

Three dimensional simulations have been performed to investigate flow quantities in the turbulent Reynold’s Aver- aged Navier Stokes and Large Eddy Simulation models in a single phase flow. Simulation of a two-phase flow with the Large Eddy Simulation model was also performed to investigate the inception and development of cavitation. The Volume of Fluid approach was used to describe the two- phase flow and Rayleigh-Plesset equation to solve bubble dynamics.

The mathematical models regarding those in single phase flow have been solved in the CFD software ANSYS FLU- ENT, while those in two-phase flow have been solved in the open source C++ toolbox OpenFOAM 2.0.0.

(4)

I would like to take this opportunity to thank my supervisor Mireia Altimira for given me the chance to work with this project. She believed in me though I did not have any earlier experience in applied Computational Fluid Dy- namics. I also want to thank her for all her encouragement and support during the time I was working on my master’s thesis. She has been an outstanding mentor and I have learned lots from her.

I also want to thank my examinator Mihai Mihaescu for arranging everything so I could get started with this master’s thesis.

I have also Alexander Nygård to thank, he spent a lot of time helping me when I was facing very difficult problems and I also want to thank him for all his advice, which improved my work with this master’s thesis.

Finally, I want to thank Johan Malmberg for all the interesting discussions we had in the working place, which was very rewarding.

(5)

Introduction & Background

A jet is a stream of fluid that is emanating from a nozzle or orifice into a surrounding medium containing fluid with lower velocity than the jet stream itself. Jet flows are very common in engineering, medical and environmental applications and are for instance used in fuel injection systems, spray painting and drying. Jet flow devices are also very common in applications such as cutting, hydraulic drilling, cooling and heating.

(a) Fuel injector (b) Cutting

(c) Painting (d) Cooling

Figure 1.1: Example of Jet application in industry. Pictures from [17]-[20]

Very small droplets are wanted in cooling applications to optimize heat transfer 1

(8)

and evaporation. Regarding painting applications, a uniform distribution of droplets is wanted. If the cooling layer is to thick then the costs increases. However if its to thin then it doesn’t fulfill its purpose.

One application of liquid jets that is particularly challenging is the injection of fuel in internal combustion engines. Some of these challenges are:

• The most important challenge of fuel injection in internal combustion engines is to decrease emissions while maintaining high engine efficiency. This depends essentially on the air-fuel mixture and in turn on the atomization of the fuel jets, cavitation, turbulence and flow inside the nozzle.

• The injection duration in a injection system is very short. The desired amount of fuel has to be delivered during this short time period from the fuel injection system at a correct time and rate [13]. Further, to ensure that all the fuel vaporizes which than burns from the compression in the combustion engines.

The droplets that does not evaporate have less chance to burn and this reduces the engine performance and increases emissions. This requires essentially a high level of atomization in a short period of time in order to obtain a very well mix of air-fuel mixture.

• Further, high injection pressure is required to obtain a high level of atomization. However, this can induce cavitation which may increase flow losses and lead to surface fatigue. Cavitation may also produce noise and vibration, which today is not a property that is desired. On the other hand, cavitation increases turbulence intensity which intensify the spray disintegration [13]. It also increases the spray cone angle which enhances the air-fuel mixture [2].

Modern diesel engines operate at very high pressures and uses injectors with very small diameters and length’s, of about 180 µm and 1 mm respectively [13]. The high pressure and small geometry give rise to very high jet velocities. With high pressure and velocity, the spray formation becomes highly influenced by the flow inside the jet nozzle, cavitation and turbulence[36]. The injection of fuel in combustion engines is a very delicate process. It requires a high level of atomization in order to obtain a very good mixture in a short time period. This depends ultimately on injector geometry, formation of spray and in turn of break-up of the liquid. The injection of fuel in combustion engines is therefore very interesting.

The investigated jet in this project is turbulent flow of diesel injection that issues from elliptical microscopic nozzles into the combustion chamber in diesel engines.

The purpose of this chapter is to introduce different flow properties that is relevant for this project, with some fundamental theory. Some previously conducted research on the topic that are treated in this master thesis will be discussed and this chapter ends with the objectives of this master’s thesis.

(9)

The regime in which an injector works determines the break-up mechanisms.

There are four different regimes according to [26];

1. Rayleigh regime. At very low velocities, the flow is dripping and there is no jet, just like the dripping of a faucet. A jet is formed with increase in velocity and the jet break-up is caused by the growth of axisymmetric oscillations of the jet surface and induced by surface tension. The average droplets are larger than the nozzle diameter in this regime.

2. The first wind-induced regime. With further increase in velocity, the surface tension between the jet and ambient fluid is amplified and the break-up process is accelerated, leading to a decrease of the break-up length. The average droplet size decreases in this regime being in the size of the nozzle diameter approximately.

3. The second wind-induced regime. The break-up mechanism in this regime is caused by unstable growth of short wavelength surface waves due to turbulence. The average droplets are now much smaller than the nozzle diameter.

4. Atomization regime. In this regime, the jet is completely disrupted at nozzle discharge forming a conical spray. The average droplets are much finer now.

a) b) c) d)

Figure 1.2: Break-up regimes. a) Rayleigh regime, b) first wind-induced regime, c) second wind-induced regime, d) atomization regime [14]

.

(10)

1.1 Cavitation

The combination of high injection pressures and the geometry of the injectors can induce cavitation. Cavitation occurs when the local liquid pressure is below the saturated vapour pressure, illustrated in Figure 1.3. At this point, the liquid will vaporize and give rise to the formation of bubbles of vapour. These bubbles will then collapse when entering a region of higher pressure and lead to the break-up of the liquid [24].

Figure 1.3: Cavitation [21]

.

This can have very devastating effects in a short time period to the injector walls due to shock waves generated by the implosion. This flow phenomena will create rougher surface on the injector walls which will increase the flow losses and may also lead to surface fatigue. It can also produce vibration and noise due to the implosions. But cavitation may not necessarily be negative, it can be very useful in many applications if it can be controlled. Some examples are high power ultrasonics used as alternatives to surgeries or destruction of kidney stones with shock waves due to the collapse of the cavitation bubbles. Cavitation has also a positive effect on the fuel spray quality, it was for instance explained by [13] that the collapse of cavitation bubbles inside the injector will increase the turbulence level and thereby intensify spray disintegration. It was also explained that the transition from turbulent to cavitating flow will increase the spray cone angle and a decrease in penetration length. It was also shown from numerical analysis made by [27] and [28] that cavitation leads to an increase of spray cone angle and shorter break-up lengths. Another numerical analysis regarding fuel injectors was made [21]

where they showed that the mean liquid velocity and r.m.s velocity increases due to turbulence and unsteadiness induced by cavitation. A numerical analysis regarding the influence of nozzle geometry on the flow inside the nozzle and cavitation was made by [4]. They concluded that elliptical nozzle with a vertically oriented major axis will have lower discharge coefficient and is less prone to cavitate.

However, it is uncertain if cavitation bubbles leads to jet break-up by increase of turbulence kinetic energy of the jet or by causing a direct local break-up of the jet. Either way, cavitation enhances the air-fuel mixture.

(11)

1.2. ATOMIZATION OF FUEL JETS

1.2 Atomization of fuel jets

In a diesel engine, air is allowed in to the cylinder through the intake valve which then closes at the end of the intake stroke. The piston will then start to compress the air and this will give rise to a very high temperature and pressure. Fuel is then added by the injection system which undergoes atomization and the air-fuel mixture is then ignited by the heat generated from the compression. Atomization is the process in which a maximum surface to volume ratio is achieved through transformation of a liquid stream into a spray of fine particles in a gaseous medium [14]. The atomization process is divided into two stages, namely the primary break- up and secondary break-up, illustrated in Figure 1.4.

Figure 1.4: Break-up of a full-cone diesel spray [13]

.

Primary break-up is the first disintegration of coherent liquid stream into drops and discrete structures. Possible break-up mechanisms are turbulence, cavitation, growth of surface waves due to aerodynamic forces and relaxation of the velocity profile [13]. With the latter meant that the velocity may change when entering the combustion chamber. With no wall boundary conditions, the viscous forces gives rise to acceleration of the outer jet region, causing instabilities and break-up of the other jet region. A further disintegration of these particles due to aerodynamic forces caused by the relative velocity between droplets and surrounding gas is the secondary break-up as explained by [13]. Further, as the jet moves downstream, it becomes more diluted due to entrainment. The velocities of the droplets are highest in the axial direction and decreases in the radial direction due to interaction with the entrained gas. It is in this dilute spray that evaporation occur and this regime is very important in the combustion process because finer particles gives better mixing as explained before. The droplets can also collide and change droplet properties.

Droplets can also bond and form larger droplets, this is called droplet coalescence.

(12)

1.3 Secondary flow

Turbulent flow through bent nozzles are very common in industry in e.g. combustion engines, chemical reactors and heat- and mass transfer systems. Secondary flow in bent nozzles are very common because it forces the flow to change the flow direction due to centrifugal forces, which generates an adverse pressure gradient. The pressure increases in the outer wall of the bend which also implies an decrease in velocity (Bernoullis theorem). The pressure on the inner wall of the bend however will show the opposite behaviour. This unbalanced condition giving rise to secondary flow is illustrated in Figure 1.5a showing a pair of two counter-rotating vortices. This was first discovered by Dean [16] and these vortices are therefore called Dean vortices.

Lower left side of bend.

High pressure region

Lower right side of bend.

High pressure region Upper left side of bend.

Low pressure region

Upper right side of bend.

Low pressure region

(a) Counter-rotating vortices [16]

(b) recirculation region [13]

Dean proposed a parameter that defines the secondary flow in curved nozzles, this is called the Dean number;

De = Re sR

R_c (1.1)

where Re is the Reynolds number, R the nozzle radius and R_c is the curvature radius. The secondary flow is strongly dependent on the Reynolds number and bend curvature since it is proportional to Reynolds number and inversely proportional to bend curvature. Further, this equation is also the ratio between the centrifugal and viscous forces. The bent nozzle would be a simplification of the actual injector and we want to reproduce the effects on the jet of this sudden direction change, keeping a rather simple configuration.

In the early 70’s, experimental analysis was made regarding mean quantities on the flow by [31]. It was concluded that the secondary motion is maximal at 30^◦ from bend inlet and decreases until it reaches 90^◦, where it reaches a steady value.

Later on, experimental analysis was made regarding flow in a 90^◦ bend by [32]

(13)

1.3. SECONDARY FLOW

where turbulence intensity and velocity field was investigated with Laser-Droplet- Velocimetry (LDV). It was shown that these flow properties was strongly dependent on the Reynolds number. Another numerical and experimental study was made by [10] regarding secondary flow in a 90^◦ turbulent nozzle elbow. They showed that the bend curvature has an utmost effect on the swirl intensity generated after the elbow. They showed that the swirl intensity increases exponentially after the elbow but dissipates quicker with a larger curvature radius. Further, a numerical analysis was made regarding flow in a U-bend where turbulence properties was investigated by [33]. It was shown that the r.m.s. values increases when the flow passes by the bend due to additional mean strain, which was related to the turning of the primary flow, which also lead to creation of secondary flow.

The lowest pressure are reached in the recirculation region due to reduction of cross-sectional area of the jet at the inlet edges, which also is called vena contracta, illustrated in Figure 1.5b. If the pressure in this region reaches the vapour pressure of the liquid, then the liquid is transformed into vapour. Further cavitation is enhanced by shear flow, which is caused due to large velocity gradients. The shear flow generates vortices and the static pressure in the center of these vortices is lower than in the surrounding liquid due to centrifugal forces. This may give rise to additional cavitation bubbles, as it was described by [13].

The axisymmetrical jet¹ is the most common configuration in industry but some research has been made on non-axisymmetrical jets, such as the elliptical jet. It was shown from experiments (discussed below in this section) that the use of elliptical jets enhances mass entrainment, mixing and spray characteristics [7],[3], [5] and [8].

However elliptical jets have a tendency to create instability in the flow [3],[5],[34]

and few research has been made on this subject.

It is worth to study the behaviour of the elliptical jets on the flow to gain a better understanding on the mechanisms driving mixing and entrainment to further improve an already beneficial configuration. Another study on single phase flow was made by [3] and they showed from experimental analysis that small aspect ratio elliptical jets increase the mass entrainment and that there is a significant difference compared to axisymmetric jets. They concluded that the difference in flow properties between the major- and minor axis was very large due to axis switching, e.g. a difference of 26% of initial momentum thickness. A similar experimental study was made by [8] and they showed that the axis-switching behaviour of the elliptical jets results in an unstable jet and shorter breakup lengths compared to the circular jet. But this depends also on the viscosity of the liquid, they showed that increasing the aspect ratio of the elliptical nozzle results in a instable flow if the viscosity is low.

A experimental analysis regarding a two-phase flow was made by [5] where they investigated the efficiency of mixing performance and mass entrainment for non-

1Symmetry around an axis, circular nozzles is typically axisymmetric while elliptical nozzles are non-axisymmetrical

(14)

circular jets. The simulations they made compared elliptical, rectangular, triangular and circular jets, all with the same equivalent diameters. They concluded that the non-axisymmetrical jet gives better mass entrainment as well as mixing performance than the circular jet due to axis-switching induced by mean streamwise vortices.

1.4 Objectives

The performance of diesel engines is strongly affected by the fuel spray, atomization and in turn the mixing process as mentioned earlier. These depends ultimately on the dimensions and geometry of the nozzle. The purpose of this project was therefore to investigate different elliptical nozzle geometries which also was compared to the conventional circular nozzle. Several different configurations have been simulated under conditions found in diesel injectors for these models regarding spray characteristics, cavitation and mixing performance. A whole load of research have been made on diesel injectors but only few on elliptical jets. The research found concerning the elliptical jets was either experimental investigations or simulations made in RANS (Reynolds Averaged Navier-Stokes). The turbulence RANS model gives the ensemble averaged flow fields and do not resolve the large unsteady flow structures. The simulations in this master thesis are made in RANS but also in LES (Large Eddy Simulation) in order to capture large unsteady flow structures that are responsible for jet mixing. The LES model is based on low-pass filtering of the Navier-Stokes equation and resolves the large scales of the turbulence flow and models the small scale. It is more expensive than RANS and is for most cases not practical to use because it requires a very fine grid and fine discretization both in time and space. The main purpose of this master thesis is thus to capture more information of the flow properties using LES with different nozzle geometries. The main objectives are listed below.

• The first objective was to investigate the influence of nozzle cross-section and curvature. This was done by modelling different geometries of straight and bent nozzles using RANS in a single phase flow. Further, the two configurations with higher TI was also simulated using LES in a single phase flow to study the large unsteady flow structures that are responsible for jet mixing.

• Then the simulations of single phase flow are compared with simulations of two phase flow, in which liquid is injected into vapour gas to investigate the effect of cavitation.

• Finally the impact of simplifications in the injector geometry and mathematical model was assessed in the final flow characteristics.

(15)

Chapter 2

Mathematical modeling

This chapter describes the mathematical model of the flow inside the injectors.

It begins with an overview of the governing equations and turbulence. Further, mathematical models of single phase flow is described for both RANS and LES.

Finally, mathematical models of multiphase flow and cavitation are described.

2.1 Governing equations

Computational fluid dynamics (CFD) is a tool to analyse fluid dynamics and provide numerical simulations by computational methods. The range of physical phenomena that can be solved with CFD is very broad. The governing equations of fluid dynamics and the basis for CFD are the Navier-Stokes equations. These equations are derived from the conservation of continuity and momentum equations. These equations are for incompressible flow expressed as:

ρ∂ui

∂t + ρ ∂

∂x_j(u_iuj) = −∂p

∂x_i + ∂

∂x_j

"

µ ∂ui

∂x_j +∂uj

∂x_i

!#

+ F_i (2.1)

∂ui

∂x_i = 0 (2.2)

where ρ is the density, µ the fluid viscosity and F_iis an additional source term of surface tension which implies a pressure jump across the interface. These differential equations are discretized in finite volume method in order to perform numerical simulations. The Navier-Stokes equations are used throughout the entire project but with some modifications depending on which model that is applied.

2.2 Turbulence modeling

The most accurate simulation available is the Direct numerical Simulation (DNS) because it solves the unsteady Navier-Stokes equation in the entire domain without any turbulence models, meaning that the whole domain of scales must be resolved

9

(16)

which is computationally very expensive. This model has therefore only been used for very simple cases. It is not practical to use this model for complex cases which requires large grids and short time steps in order to capture the flow properties, not even with high parallel computers. The LES model lies between RANS and DNS regarding computational cost. This method is less expensive than the DNS method but is still not practical for many applications since it requires a very fine grid near the walls in order to resolve the large scales of turbulence. This model is based on low-pass filtering of the Navier-Stokes equation and resolves the large scales of the turbulence flow and models the small scales.

Turbulence is the state where the flow becomes unstable, leading to random vorticity fluctuations and chaotic flow [15]. The structure of the rotational flow in turbulence are called eddies which varies from large scales to smallest Kolmogorov microscales. These are the smallest scales in which eddies are entirely dominated by viscosity. The largest eddies are inviscid because the turbulent scale, i.e. characteristic properties of the vortices are of the same order as the scales for the mean flow. This implies that the largest eddies are dominated by inertia effects. The energy from the mean flow will be extracted when transporting these eddies due to vortex stretching. Vortex stretching is the changes on vorticity by production of velocity gradients due to conservation of angular momentum. This will create smaller eddies which will be transported by the largest eddies due to the effect of vortex stretching and this process will continue until the eddies are entirely dominated by viscous effects. If the flow is laminar or turbulent is depending on what the Reynolds number is. The Reynolds number is defined as

Re = ρU L

µ (2.3)

where U and L are characteristic velocity and length respectively. The char- acteristic length is the hydraulic diameter which is D_h = 4A/P in which P is the perimeter of the jet nozzle. The turbulent flow in the jet will be simulated using the steady turbulent RANS model and LES. These are explained in the subsections below.

2.2.1 RANS turbulence models

The turbulent RANS model solves the ensemble or time-averaged Navier-Stokes equations and is thus based on averaged quantities. This model is most widely used because it is fast and computationally cheap but it is not the optimal model since it uses time-averaging quantities. This means that it can not capture all the flow properties, thus some information is lost. The RANS equations model are given as follows

ρ∂ ¯ui

∂t + ρ ∂

∂x_j( ¯uiu¯j) = −∂ ¯p

∂x_i + ∂

∂x_j

"

µ ∂ ¯ui

∂x_j +∂ ¯uj

∂x_i

! + ρτ_ij

#

(2.4)

(17)

2.2. TURBULENCE MODELING

∂ ¯ui

∂xi

= 0 (2.5)

where the last term in equation 2.4 is the Reynolds stress tensor, which is the internal stress acting on the mean turbulent flow. Further, τ_ij is the kinematic Reynolds stress tensor, defined as

τij = −u⁰_iu⁰_j (2.6)

The Reynolds stress tensor is related to the mean velocity gradients through the Boussinesq hypothesis, defined as

−ρu_iuj = µ_t ∂ ¯ui

∂x_j +∂ ¯uj

∂x_i

!

−2

3ρkδij (2.7)

here δ_ij is the Kronecker delta, k the turbulence kinetic energy and µ_t is the turbulent eddy viscosity.

Two-equation models

Further two-equation models is used, which calculates the turbulent viscosity and describes the turbulent properties. The turbulent k − and k − ω models are used in this project. These are described below.

• k − model

RANS simulations were made in FLUENT which provides three different k −  models: Standard, RNG and Realizable. The latter was chosen since the flow phenomenons that occur near the wall is important to consider and this model is suitable for resolving flows in the near-wall region. The two-equation model describes the turbulent flow with two differential transport equations for the kinetic energy k and dissipation rate .;

∂k

∂t + ¯u_j ∂k

∂xj

= ∂

∂xj

"

ν + ν_t

σk

∂k

∂xj

#

+ P_k− ρ (2.8)

∂

∂t + ¯uj

∂

∂x_j = ∂

∂x_j

" ν + νt

σ

∂

∂x_j

#

+ C₁S − C2

² k +√

ν (2.9)

where P_k describes the production of the turbulent kinetic energy due to mean velocity gradients;

Pk= τ_ij∂ ¯ui

∂xj

(2.10) and

(18)

C₁= max

0.43, η η + 5

; η = Sk

(2.11)

S = (2S_ijS_ij)^1/2= 1 2

∂ ¯u_j

∂xi

+ ∂ ¯u_i

∂xj

!

(2.12) The terms σ_k, σ, Ci are empirical constants and ν is the kinematic viscosity.

The dissipation rate in the realizable model is derived from the mean-square velocity fluctuation and the reason for this model to be realizable is due to the conditions that are enforced for the Reynold’s stresses. These conditions make the model more consistent with turbulent flow than the Standard- and RNG models. The turbulent eddy viscosity for k − turbulence models is defined as

µ_t= ρC_µk²

(2.13)

where is the dissipation of the turbulence kinetic energy and C_µis the empirical constant, defined as a function of the mean rotation tensor and the mean strain-rate tensor;

C_µ= 1

A0+ A_s^kU^∗ (2.14)

U^∗≡^qS_ijS_ij + ˜Ω_ijΩ˜_ij , Ω˜_ij = ¯Ω_ij − _ijkω_k− 2_ijkω_k (2.15)

A₀ = 4.04 , A_s=√

6cosφ (2.16)

φ = 1 3cos⁻¹

√

6SijS_jkS_ki S˜³

, S =˜ ^qS_ijS_ij (2.17)

• k − ω model

This model is much better than the k − model regarding boundary layer flows and the difference is that the k − ω model uses a transport equation where the turbulent dissipation rate is replaced with the specific dissipation rate ω. There are two different k − ω models in FLUENT, Standard [22] and SST (Shear Stress Transport model) [23]. The latter model was chosen since it predicts separation very accurately and it also uses a combination of both k − and k − ω models. The transport equations for this model are

∂k

∂t + u_j ∂k

∂x_j = P_k− β^∗kω + ∂

∂x_j

"

(ν + σ_kνt) ∂k

∂x_j

#

(2.18)

(19)

2.2. TURBULENCE MODELING

∂ω

∂t + u_j ∂ω

∂x_j = αS²− βω²+ ∂

∂x_j

"

(ν + σ_ων_t) ∂ω

∂x_j

#

+2(1 − F₁)σ_ω21 ω

∂k

∂x_i

∂ω

∂x_i (2.19)

The turbulent stress tensor and kinematic viscosity for this model is defined as

τij = µ_t

2S_ij −2 3

∂uk

∂x_k

−2

3ρkδij , µt= ρa1k

max(a₁ω, ΩF2) (2.20) in which

F₁ = tanh



 (

min

"

max

√ k

∂β^∗ωd,500ν d²ω

!

, 4ρσ_ω2k CD_kωd²

#)4

 (2.21)

F2= tanh





"

max 2

√ k

∂β^∗ωd,500ν d²ω

!#2

 (2.22)

The term CD_kω = max2ρσ_ω2_ω¹_∂x^k

j

ω

∂xj, 10⁻²⁰ and the closure constants are given in table 2.1.

Table 2.1: Closure constants

σ_ω σ_k β a₁

k −  0.856 1.00 0.0828 k − ω 0.65 0.85 0.075

SST 0.09 0.31

Near wall treatment

The turbulence in the nozzle is affected by the presence of the wall, where the flow has great velocity gradients and the mesh should be very fine in order to satisfy the physics. Further, the Boussinesq assumption that the Reynolds stress tensor is proportional to the strain rate tensor is not valid in flows with strong curvature.

The two-equation models have therefore problems to predict the flow that is affected by strong curvature and needs wall treatments. There are three region in the near wall layer; viscous sublayer, buffer layer and logarithmic layer as can be seen in Figure 2.1 The inner part of the near wall layer is the viscous sublayer and this region is dominated by viscous effects. The logarithmic layer is the region where the flow is fully turbulent and the buffer layer is the transition between these two.

The dimensionless viscous wall distance y⁺ should be less than 5 in the viscous

(20)

sublayer, between 5 and 60 in the buffer layer and greater than 60 in the logarithmic layer. The k − ω models applies throughout the boundary layer, however it needs a sufficient near wall resolution. The k − models on the other hand needs wall functions to solve the viscous sublayer.

Figure 2.1: Law of the wall .

The Enhanced Wall Treatment (EWT) in FLUENT was used as the near wall modeling method for the k − realizable model. This method combines a two- layer model with wall-functions. It resolves the viscous sublayer [12] when y⁺ ≤ 5 and uses wall functions when y⁺ ≥ 30 [12]. This requires a fine near-wall mesh in order to resolve the viscous sublayer and satisfy the physics. The dimensionless wall distance was calculated from

y⁺= yuτ

ν (2.23)

where the wall friction velocity u_τ is a function of the friction coefficient, which in turn is dependent on the Reynold’s number.

2.2.2 LES

The purpose of LES is to capture the large scales of turbulence and model the small scales. If the eddies are smaller than the grid sizes, then it is re-modeled by a subgrid scale (SGS). The filtering operation is defined as

Φ =¯

∞

Z

−∞

Φ(~x⁰)G(~x − ~x⁰)dx and

∞

Z

−∞

G(~x − ~x⁰)dx = 1 (2.24) where G is the spatial filter and the filtered Navier-Stokes equations is as follows

ρ ∂ ¯ui

∂t +∂( ¯uiu¯j)

∂x_j

!

= −∂ ¯p

∂x_i + ∂

∂x_j ν∂ ¯ui

∂x_j

!

−∂τij

∂x_j (2.25)

(21)

2.3. MULTIPHASE FLOW

where the turbulence stress tensor is defined as τ_ij = ρ(u_iu_j− ¯u_iu¯_j). The purpose of the SGS model is to transport energy from resolved grid scales to un-resolved grid scales. The most widely used SGS model in channel- and internal flows is the Smagorinsky model which is defined as

τij−1

3τijδij = −2µ_tS¯ij where µt= 2(C_s∆)²| ¯S|Sij (2.26) where the filter width is ∆ = (∆x∆y∆z)^1/3and C_sis the Smagorinsky constant which usually have values between 0.1 and 0.2.

2.3 Multiphase flow

Multiphase flows are flows in which there is a mixture of phases. The flow pattern of multiphase flow regimes are characterized by a dispersed flow or separated flow. In the case of dispersed flow, one phase is consisting of infinitesimal particles, bubbles or drops that are distributed in the other continuous phase. A separated flow on the other hand consists of two or more phases that are separated by interfaces.

Multiphase flows are widely encounted in industrial processes, some examples are anti-icing of aircrafts, food processing and fuel injection in combustion engines, the latter being one of the key issues in this project. The prediction of the fluid behaviour in a mixing process is a key issue to enhance the efficiency of these processes. The prediction of spray characteristics in combustion engines for example is utmost important for the mixing process in order for the engine to work efficiently but also reduce pollutants.

One approach regarding mathematical modeling flows is the "one-fluid" approach, which is used for two-phase flow simulations in this project. This approach describes the entire flow field with one set of governing equations, which are solved on regular fixed grid. This approach can either be used whit a Front-Tracking method or Front-Capturing method. The properties of each phase are generally different and tracking methods are used to locate the interface and jump conditions are modeled for interface exchanges. Front-Tracking methods use Lagrangian markers to identify the interface while Front-capturing methods use marker functions to locate the interface. The second method is applied in this project through the Volume of Fluid method (VOF), described in the next section.

2.3.1 Volume of fluid

The VOF method is a free surface technique used to capture phase interfaces. This model solves the fluid of the mixture by a single set of mass and momentum equations given below.

∂ρ

∂t + ∂

∂xi

(ρu_i) = S_i (2.27)

(22)

∂(ρu_i)

∂t + ∂

∂xj

(ρu_iuj) = −∂p

∂xi

+ ∂

∂xj

"

µ ∂(u_i)

∂xj

+∂(u_j)

∂xi

!#

+ F_i (2.28) However, the pressure jump conditions due to surface tension and interface curvature are modelled through an additional source term in the momentum equation.

Further, due to the mass transfer between phases, the continuity equation is no longer divergence free and a source term is included in the volume fraction equation.

The dynamics of the different phases are solved by a transport equation, in which a indicator function α is used to represent the volume fraction of one phase.

This transport equation is given as

∂α

∂t + ∂

∂xi

(αu_i) = S_i (2.29)

Mixture properties ρ and µ are averaged, using the volume fraction as

ρ = αρ_l+ (1 − α)ρ_g (2.30)

µ = αµl+ (1 − α)µ_g (2.31)

where l and g denotes liquid phase and gas phase respectively.

2.4 Cavitation

The simulations of cavitation in this project are based on the Rayleigh-Plesset [24]

equation for dynamics of bubble growth and collapse. This equation describes the growth of a single vapour bubble in liquid and is defined as

R_Bd²R_B dt² + 3

2

dR_B dt

2

= p_B− p_∞ pl

− 4v_l RB

dR_B

dt − 2σ ρlRB

(2.32) where R_Bis the bubble radius and σ is the liquid surface tension coefficient. The variables p_B is the bubble pressure, which is equal to the saturated vapour pressure p_v when the bubble only contains vapour. The parameter p_∞ is the local far-field pressure. By neglecting viscous and surface tension terms, the above equation can be integrated in time yielding the Rayleigh equation (2.33), which describes the growth of the bubble.

R˙B= sign(p_v− p_∞) s2

3

|p_vap− p∞|

ρ_l (2.33)

Further the volume fraction of vapour can be connected to the bubble radius through number of bubbles per volume of liquid [Schnaerr and Sauer]

(23)

2.4. CAVITATION

α =

4π 3 nBR³_B

1 +^4π₃ nBR³_B (2.34)

and the mass flow rate is expressed as m = C˙ ρ_vρ_l

ρ (1 − α)α3 ˙R_B

R_B (2.35)

where C is a fitting coefficient.

(24)

(25)

Chapter 3

Models & Setup

The purpose of this section is to introduce the geometry of the different models, discretization of the computational domain, boundary conditions, verification and validation of the simulations. This section ends with a comparison of flow quantities between different turbulent models.

3.1 Geometry

CFD simulations were carried out on microscopic nozzles of two different geometry configurations, these are straight- and bent nozzles. Simulations were made on three different cross-sections, circular and two different elliptical cross-sections.

Regarding the bent nozzles, two different orientations of the ellipses with respect to the bend axes has been considered. The geometry parameters of these are shown in Table 3.1. The cross-sections of these nozzles are taken from [4] and they correspond to actual diesel injectors. The ellipses H₁ and H₂ have horizontally oriented major axis while the ellipses V₁and V₂have vertically oriented major axis. The major axis is denoted as a and minor axis as b. The cross sectional areas are approximately the same for all nozzle types.

D_h

a b

Figure 3.1: Nozzle

Nozzle type a(µm) b(µm) Area(m²) D_h(µm)

Circular 85 85 2.270e⁻⁸ 85

H₁ 103 70 2.265e⁻⁸ 82.6

H2 120 60 2.262e⁻⁸ 77.8

V₁ 103 70 2.265e⁻⁸ 82.6

V₂ 120 60 2.262e⁻⁸ 77.8

Table 3.1: Nozzle types

Simulations have been made on quarter parts of the straight nozzles. This have been made to save computational time. However, the flow conditions are still valid since all normal gradients are zero at the symmetry planes and there is no flow

19

(26)

variations normal to the boundary planes. The simulations have been made on half parts of the bent nozzles since these models are only symmetric on one plane.

However, entire models have been used for the LES cases.

• Straight nozzles

Simulations of straight nozzles have been made for the three cross-sections; cir- cle, H1 and H2. The geometry of these are shown in Figure 3.2. All the dimension- less magnitudes presented from now on are calculated using the D_has characteristic length scale. The dimensions for the quarter parts as shown in the figure are; nozzle length has a size of 9.4D, chamber length are 17.6D and chamber width and height are both 7.6D.

Figure 3.2: Computational domain for RANS simulations of straight nozzles .

• Bent nozzles

Simulations of bent nozzles have been made for all the nozzle types, presented in Table 3.1. These models are illustrated in Figure 3.3.

The nozzle length before the bend is 3.5D and the nozzle length after the bend is of equal size as the nozzle length for the straight models. The length of the bend is 3.5D with a outer curvature radius of 2.4D. The geometry of the chamber is equal of that of the straight models. The numbers in Figure 3.3 (1-9) shown at different cross sections is used from now on to present results but with index S before. Different planes in the bend is also depicted with angles (0 − 60^◦) from bend inlet. These positions also apply to the straight nozzles starting from position 2.

3.2 Discretization of the computational domain

The element shapes in discretization of three-dimensional geometries are tetrahe- drons, hexahedrons, prism and pyramids. Further, the grids are characterized as either structured or unstructured mesh. Structured grids consist of regular connectivity, where all the vertices are topologically alike. Unstructured grids on the other

(27)

3.2. DISCRETIZATION OF THE COMPUTATIONAL DOMAIN

3.5D 2.4D

Ln=9.4D 3.5D

1 2 4 5 6

7 8 9

3

Cirkel

V2 V1 H1 H2

L=17.6D

30 60 0^o

o

0

Figure 3.3: Submerged jet of straight pipes .

hand consist of irregular connectivity, where the vertices have arbitrarily varying neighboring cells. Unstructured grids may therefore be computationally more expensive, however it provides more freedom and better fits to complicated domains.

The mesh for all the geometry configurations was constructed in Ansys work- bench meshing with pure hexahedral elements, in order to obtain structured mesh.

The mesh for all straight models have been constructed in the same way and the same goes for the bent models, this is necessary to make a good comparison. There- fore, only one mesh for each configuration is shown here. The mesh domain was divided into different volumes, this was necessary to obtain a pure structured hexahedral mesh.

The mesh for the straight model is illustrated in Figure 3.4, corresponding to the circular nozzle. The mesh on the nozzle inlet starts with a coarser mesh in the core and it becomes finer as it approaches the near wall region. The mesh on the nozzle wall is quite uniform until it approaches the nozzle inlet- and outlet. The mesh is much finer at these locations in order to obtain a smoother transition between the neighboring cells. The mesh on the chamber is finer at the core and becomes coarser as it goes further away. This results in a mesh of 1.2 million elements.

(28)

Figure 3.4: Mesh for case 1.1 .

• Bent nozzles

The mesh for the bent model is illustrated in Figure 3.5. The mesh for the bent model have been constructed in the same manner as for the straight model. Only difference is on the nozzle wall, where the mesh is finer at the bend. This mesh corresponds to H2 nozzle, which results in a mesh of 2.8 million elements.

3.3 Cases

Simulations have been made on different nozzle geometries with different modeling techniques. All cases that have been simulated are listed in Table 3.2 to simplify the reading.

In case 1, simulations have been made in FLUENT using RANS with both turbulence models. The grid sensitivity analysis was only made for case 1.1, discussed in section 3.5. In case 2, similar simulations have been made but now with the bent

(29)

3.3. CASES

Figure 3.5: Mesh for case 2.4 .

Table 3.2: Cases

1. RANS straight nozzles k k_ω Grid sen. Software

1.1 Circular X X X FLUENT

1.2 H1 X X FLUENT

1.3 H2 X X FLUENT

2. RANS bent nozzles

2.1 Circular X FLUENT

2.2 H1 X FLUENT

2.3 V1 X FLUENT

2.4 H2 X X X FLUENT

2.5 V2 X FLUENT

nozzles instead. The nozzles with vertically oriented major axis are simulated for this case as well. This is because the orientation of the jet have an major influence on the flow properties, as discussed earlier. One simulation have been made for the k − model in this case to compare the prediction of flow separation from each model.

(30)

3.4 Boundary conditions

The boundary conditions are illustrated in Figure 3.6. A total gauge pressure of 30 MPa and a initial gauge pressure of 10 MPa have been set as inlet condition. A gauge pressure of 6 MPa was imposed as an outlet condition. The no-slip condition is imposed at the walls and the normal velocity and normal gradients of all variables are zero at the symmetry planes. It should be noted that, in simulations of the bent nozzle, there is only one symmetry plane and that, in LES simulations, the entire geometry is considered. The diesel fuel that was used have viscosity of 0.0024 kg/m-s and density of 730 kg/m³.

Figure 3.6: Boundary conditions for straight pipe .

Simulations have been solved using the pressure-velocity coupling algorithm, SIMPLE [35]. Further a second order scheme is used for pressure, second order upwind scheme for momentum, turbulence kinetic energy and specific dissipation rate/turbulent dissipation rate for the RANS cases. Regarding LES, a bounded central differencing scheme is applied for the momentum equations. In the VOF method, a VanLeer scheme is used for volume fraction flux.

Regarding the simulation of two-phase flow, both liquid and vapour fuel has been considered incompressible. The vapour fuel have density and kinematic viscosity of 9.4 kg/m³ and 7.45e-07 m²/s respectively. The surface tension between the phases has been set as 0.029 kg/s² and the pressure of vapour fuel as 1329Pa. Further, parameters of the cavitation model has been set as C_c=1, C_v=1, n₀=1.6e13 m⁻³ and R₀=2µm.

3.5 Verification

Numerical modeling provides solutions which are only approximations of real physical properties. The formulation of mathematical models and simplifications of the models results in errors in numerical modeling, which has to be considered. One error is called round-off error, which exists due to limited number of computer digits available for storage, but this type of error is not major compared to other errors.

(31)

3.5. VERIFICATION

Another error is the iterative convergence error which exists because the simulations have to be stopped eventually. This error is minimized by letting the residuals go down to at least 10⁻⁶ except the continuity residual, because it depends on the initial solution. If the initial guess is good, the residual for continuity will not go down as much as the other residuals. Another error that is considered in this project is called discretization error, which arises when converting a continuous model to discrete form. The solution should become less sensitive and reach a continuum solution when the mesh is refined. Simulations were therefore made for two additional mesh resolutions, a coarser and finer mesh. The mesh resolution for the original mesh lies between these two. A sizing factor (s) of 1.25 was applied to the original mesh resolution, which resulted in a coarser mesh of 500k elements and a finer mesh of 2.3 million elements for the straight model. Regarding the bent model, a sizing factor of 1.25 was only applied at the most relevant regions, due to computational cost and time limitations. This resulted in a coarser mesh of 1.9 million elements and a finer mesh with 3.7 million elements.

The solution of the different mesh resolutions was compared by looking at the velocity profiles at mid chamber. Further a Grid Convergence Index (GCI) by [25]

was applied to compute the discretization error. The GCI is defined as GCI_{f ine} = F_S||

s^g− 1 (3.1)

This is the GCI for the finer mesh in which F_S is a safety factor and a value of 1.25 was used which correspond to a comparison between three or more grids. The parameter g is the order of convergence, which is second order in this case i.e. 2 and is the relative error which is defined as

 = ff ine− f_orig

f_orig (3.2)

where the quantity f is the massflow rate for the original- and finer mesh. The GCI for the coarser mesh can then be extrapolated through

GCIcourse= s^g· GCI_{f ine} (3.3)

This grid sensitivity study was made only for the circular and ellipse H2 cross- section of straight and bent models, respectively due to resemblances.

The velocity profiles for all three mesh resolutions regarding the straight model are shown in Figure 3.7. It is very clear from this figure that there are almost no differences at all between the velocity profiles, indicating that the solution is grid independent.

From the GCI calculations, values for the mesh index of 0.035% and 0.055% were obtained for the fine- and course mesh respectively. This gives a discretization error

(32)

Figure 3.7: Grid sensitivity analysis for cases 1.1 at mid chamber x/L = 0.5, where L is the chamber length.

which is below 2%. However since there are almost no differences, the intermediate mesh was chosen in order to save computational cost and time.

• Bent nozzles

The velocity profiles for all three mesh resolutions regarding the bent model are shown in Figure 3.8. The solution is also grid independent.

Figure 3.8: Grid sensitivity analysis for cases 2.1-2.3

From the GCI calculations, values for the mesh index of 0.015% and 0.023% were obtained for the fine- and coarse mesh respectively. This gives a discretization error

(33)

3.6. VALIDATION

which also is below 2%. The intermediate mesh, i.e. original mesh will therefore be chosen as the baseline model.

3.6 Validation

Simulations have been compared to experimental analysis and the validation of the simulations are discussed in this section. Further, the impact of simplifications in injector geometry and mathematical model will be assessed for the final flow characteristics.

Simulations were compared to experimental analysis [39] regarding discharge coefficient. They reported values betwen 0.77-0.79 for the same reynolds number around 10000. The experimental analysis were made on cylindrical nozzles with circular cross-sections. The discharge coefficient is defined as the actual mass flow divided by the theoretical mass flow:

Cd= m˙ A√

2ρ_l∆P (3.4)

where A is the cross-sectional area of the nozzle and ∆P is difference in pressure between inlet and outlet cross-sections. These parameters are given in Tabel 3.3:

Table 3.3: Flow properties for injection pressure of 30 MPa and backpressure of 6 MPa

∆P (MPa) m(g/s)˙ C_D(−)

Circular 23.8 3.27 0.773

H1 23.8 3.34 0.792

H2 23.8 3.35 0.795

V1 23.8 3.11 0.737

V2 23.8 2.93 0.694

V2 (Two-phase) 24 2.83 0.64

We can see in this table that ∆P is the same for all cases, except being slightly higher in the two-phase flow. All the cases in this table, except the last row are simulations in single phase flow in the bent cases. Further, the discharge coefficient decreases with decrease in mass flow rate, which is expected since the mass flow rate is directly proportional to the discharge coefficient. We can also see that the H2 nozzle gives highest value meaning that it has higher transformation of pressue energy into kinetic energy. The V2 nozzle shows lowest values because of higher flow losses. The discharge coefficient in the two-phase flow is also lower than in single phase flow due to cavitation that further weakens the flow capacity.

The discharge values being fearly close to the experimental values gives a reason to not simulate the entire injection domain with tens of millions elements. It would be enough to simulate a small part of the injection domain as in the bent cases, which would save a great deal of computational cost and time.

(34)

The second term in the equation of Rayleigh-plesset (2.33) is the contribution of liquid viscosity. It is seen that this term is proportional to bubble deformation rate and inversely proportional to bubbles radius. This term is therefore only significant for bubbles with small radii. The last term in the equation is the contribution of surface tension. This term is only significant for bubbles with small radii as well since the bubble radius appear in the denominator. These terms have therefore been neglected. According to [37] various characteristic time scales can be compared to assess the impact of the neglected terms on the dynamics of a cavitation bubble.

These time scales are as follows:

τν = a²

4ν_l (3.5)

τ_S = a rρ_la

2σ (3.6)

τp = a

s ρl

p_ref− p_ν (3.7)

where a is the characteristic length scale for the bubble radius, p_ref has been set as 1 bar and the indexes in the time scales correspond to viscosity, surface tension and pressure. The characteristic time scales are plotted as a function of the characteristic length scale in Figure 3.9

Figure 3.9: Characteristic time scales versus length scale.

It is seen in this figure that the dynamics of cavitation bubbles are controlled by pressure in a wide range of radii. Pressure is dominating already at 1.5 µm.

Viscosity becomes dominating for cavitation bubbles smaller than this value and surface tension is never dominating. The viscous and surface tension terms can therefore be neglected considering that the initial bubble radius is 2 µm. The

(35)

3.7. COMPARISON OF TURBULENCE MODELS

simplified Rayleigh-Plesset equation is therefore enough for solving the dynamics of cavitation bubbles.

3.7 Comparison of turbulence models

The two-equation turbulence models; SST k − ω and realizable k − were compared in this project, in order to choose a suitable model that fits best to the investigated cases. The comparison of these models for both straight nozzle and bent nozzle is described below. This study was only done for cases 1.3 and 2.4 from table 3.2.

The former is also compared to the LES model of case 3.1.

• Straight nozzle

The comparison of the two-equation models for the straight nozzles is shown in Figure 3.10. These comparisons are regarding mean velocity magnitude and turbulence intensity.

Mean Velocity magnitude Turbulence Intensity

k-

LES Realizable

SST

Figure 3.10: Contours of mean velocity magnitude and turbulence intensity for three different turbulent models for straight pipe cases.

In Figure 3.10, both turbulent two-equation models shows similar behaviour, however the realizable k − model shows higher turbulence levels at the corner of the nozzle discharge. Further, the velocity in the nozzle are quite similar for both models. However, the velocity in the core of the jet becomes lower for the k −  realizable model as it develops with increasing axial distance in the chamber.

Compared to the LES model, none of the turbulent two-equation models shows resemblances of the turbulence intensity. It is seen that the two-equation models shows high turbulence levels at the nozzle discharge while the LES model shows much lower value. However, the LES model shows much higher turbulence levels at mid chamber.

(36)

• Bent nozzle

The comparison of the two-equation models for the bent nozzles is shown in Figure 3.10.

Velocity Turbulence Intensity

k-

LES Realizable

SST

Figure 3.11: Contours of mean velocity magnitude and turbulence intensity for three different turbulent models for bent pipe cases.

The velocity contours in Figure 3.11 shows that the k − ω model is capable of capturing separation at the bend more accurately than the k − model. This is quite clear when comparing these two-equation models to the LES model.

However, the k − model shows higher resemblance to LES regarding the ve- locity at the lowerside of the nozzle. Regarding the turbulence intensity, all three models shows high turbulence levels at the bend which develops downwards to the lower corner at the nozzle discharge. However, the k − model shows much higher turbulence level at the bend.

The realizable k − model gives the most realizable solution since it uses a eddy viscosity formulation that is not constant, which prevents the normal stresses to become non-negative and the Schwarzs inequality for shear stresses are satisfied.

The SST k − ω model is a combination of the standard k − ω and k − models, the former used near the boundary while the latter is used in the free stream. Since µ_t is limited in equation (2.20), the production of shear stress and turbulent kinetic energy are not over-predicted for SST k − ω model. As soon as the jet enters the chamber, a shear layer is created between the jet and surrounding fluid. With increase in shear layer, an increase in turbulence is expected.

When comparing to LES, it seems that the SST k-omega model predicts separation better than the realizable model and since a great deal of separation will

(37)

3.7. COMPARISON OF TURBULENCE MODELS

be analysed in this project, the SST k − ω model have chosen to be the primary turbulence model.

(38)

CFD analysis of nozzle effect on jet formation