- Email: [email protected]
Detailed numerical simulation of swirling primary atomization using a mass conservative level set method
Accepted Manuscript
Detailed Numerical Simulation of Swirling Primary Atomization Using a Mass Conservative Level Set Method Changxiao Shao , Kun Luo , Yue Yang , Jianren Fan PII: DOI: Reference:
S0301-9322(16)30620-6 10.1016/j.ijmultiphaseflow.2016.10.010 IJMF 2489
To appear in:
International Journal of Multiphase Flow
Received date: Revised date: Accepted date:
13 May 2015 28 June 2016 22 October 2016
Please cite this article as: Changxiao Shao , Kun Luo , Yue Yang , Jianren Fan , Detailed Numerical Simulation of Swirling Primary Atomization Using a Mass Conservative Level Set Method , International Journal of Multiphase Flow (2016), doi: 10.1016/j.ijmultiphaseflow.2016.10.010
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
Highlights
We report detailed numerical simulations of swirling liquid atomization by using a recently developed mass conservative level set method.
Through comparing the sheet thickness, the breakup length and the cone angle, the numerical
CR IP T
convergence of the global characteristics of the swirling two phase flow has been obtained. The numerical results show that turbulent inflow can induce liquid sheet breakup near the nozzle exit, reduce the stiffness of the liquid sheet, and lead to the statistically homogeneous
AN US
distribution of small-scale liquid structures in the radial direction.
Compared with the single-phase jet, the two-phase jet exhibits the chaotic velocity filed downstream that can enhance the mixing of droplets and ambient gas, and the precessing vortex core (PVC) is not observed in the center of the two-phase jet. The preferential alignment of
i
M
with the intermediate strain rate indicates that the
AC
CE
PT
ED
fluctuating velocity in the recirculation zone is statistically similar to isotropic turbulence.
1
ACCEPTED MANUSCRIPT
Detailed Numerical Simulation of Swirling Primary Atomization Using
AN US
Method
CR IP T
a Mass Conservative Level Set
By
Changxiao Shao1, Kun Luo*, Yue Yang2, Jianren Fan1 1
M
State Key Laboratory of Clean Energy Utilization
Zhejiang University, Hangzhou 310027, P. R. China State Key Laboratory of Turbulence and Complex Systems
ED
2
AC
CE
PT
Peking University, Beijing 100871, P. R. China
Submitted to
International Journal of Multiphase Flow
*Corresponding author, E-mail: [email protected], Tel: 86-0571-87951764 2
ACCEPTED MANUSCRIPT
Abstract: We report detailed numerical simulations of swirling liquid atomization. A recently developed mass conservative level set method is employed to capture the gas-liquid interface and a ghost fluid method is utilized to deal with the jump conditions across the interface. The swirl and atomization characteristics of two-phase
CR IP T
annular swirling jets with the influence of turbulent inflow are investigated. Through comparing the sheet thickness, the breakup length and the cone angle, the numerical convergence of the global characteristics of the swirling two phase flow has been
AN US
obtained. The numerical results show that turbulent inflow can induce liquid sheet breakup near the nozzle exit, reduce the stiffness of the liquid sheet, and lead to the statistically homogeneous distribution of small-scale liquid structures in the radial
M
direction. Compared with the single-phase jet, the two-phase jet exhibits the chaotic velocity filed downstream that can enhance the mixing of droplets and ambient gas,
ED
and the precessing vortex core (PVC) is not observed in the center of the two-phase
PT
jet. In addition, the recirculation zone is smaller and farther from the nozzle exit for the turbulent inflow case than that from the laminar inflow case, and the preferential
i
CE
alignment of
with the intermediate strain rate indicates that the fluctuating
AC
velocity in the recirculation zone is statistically similar to isotropic turbulence. The interaction of the liquid-gas interface and vortices shows the preferential normal alignment of the vorticity and the normal of the interface, and the liquid sheet can generate high shear layers to produce anisotropic small-scale fluctuations. Keywords: primary atomization; swirl jet; level set method; turbulence; interface
3
ACCEPTED MANUSCRIPT
1 Introduction Liquid atomization is a common phenomenon in a variety of scientific and engineering applications, such as mixing, spraying, printing, food processing, agriculture, pharmaceutical process and combustion devices in gas turbine. The
CR IP T
atomization process is complex, involving highly turbulent and convoluted interfaces as well as breakup and coalescence of liquid masses, so that its mechanism has not yet been well understood [1]. It is of importance to explore the complex details of the
AN US
liquid atomization.
The atomization process of liquid can be divided into primary atomization and secondary atomization. Detailed experimental analysis of primary atomization remains a challenge due to the difficulty posed by dense spray measurements [2].
M
Alternatively, numerical simulations with interface capturing methods have been
ED
developed to study the primary atomization. The major interface capturing methods include the volume of fluid (VOF) method and the level set (LS) method. The VOF
PT
method is an Euler-type method with good mass conservation, but the implementation
CE
of this method is complicated owing to the reconstruction of the interface [3].
AC
Meanwhile, the LS method can easily describe topological changes under complex motions [4] but most of the existing LS methods are lack of the mass conservation property.
Both VOF and LS methods have been utilized to investigate primary atomization in straight or cross jet flows, such as in Kim et al. [5]-[6], Herrmann et al. [7], Xiao et al. [16], Shinjo and Umemura [18], Desjardins et al. [21], Ménard et al. [25]-Error! 4
ACCEPTED MANUSCRIPT
Reference source not found., Li et al. Error! Reference source not found.-Error! Reference source not found., Fuster et al. [30], M. Arienti and M. Sussman [30], Chen et al. [31].
However, there are very few of studies on primary atomization in swirling jet flow except Li et al. [28] that investigated a realistic swirling flow and qualitatively
CR IP T
demonstrated the validity and feasibility of the Coupled LS and VOF approach with experimental results. The swirling flow is significantly important when it is employed in gas turbine combustors as an approach to enhance fuel-air mixing and flame
AN US
stabilization. As demonstrated by Siamas et al. [46], many of the existing studies on annular swirling liquid jets focused on experimental visualizations and simplified mathematical models, which are difficult or insufficient to reveal the complex details
M
of liquid breakup and atomization in a two-phase environment.
For single-phase swirling flows, the swirling characteristics [42], such as vortex
ED
breakdown [47] and the precessing vortex core (PVC) [48], are extensively
PT
investigated. The formation of the recirculation zone, a form of vortex breakdown, can serve to stabilize flame. Many theories of vortex breakdown have been proposed,
CE
e.g., the generation of negative azimuthal vorticity in Brown et al. [54], the axial
AC
pressure exceeding the streamwise momentum flux in Mahesh [55] and inertial wave propagation in Squire [56]. The PVC is known as the rotating vortical structures in the center of the swirling jet, and its presence is significant for the instability and the combustion system. The observed PVC is generally considered as a secondary instability related to the post vortex breakdown [48] and can be suppressed in the presence of density gradient. On the other hand, the corresponding phenomena or the 5
ACCEPTED MANUSCRIPT
results of swirling characteristics can be different in a two-phase annular swirling jet, because the liquid sheet may significantly alter the swirling characteristics. In order to investigate the primary atomization in a swirl combustor, the atomization characteristics of an annular swirling jet are simulated using the
CR IP T
high-fidelity detailed numerical simulation with a recently developed mass conservative level set method. In addition, the influence of inflow turbulence on the global atomization characteristics is discussed. The outline of the present study is
AN US
organized as follows. In section 2, we provide the governing equations and describe the level set method. In section 3, we present numerical configurations, resolution considerations and the inflow conditions. In section 4, we discuss the atomization
ED
vorticity and interface.
M
characteristics, the influence of inflow turbulence and the correlation between
2 Numerical Methods
PT
2.1 Governing equations
CE
The Navier-Stokes equations for gas and liquid phases read
u u u p u ut g , t u 0 , t
AC
(1) (2)
where u is the velocity, ρ is the density, p is the pressure, g is the gravitational acceleration and µ is the dynamic viscosity. The material properties of gas and liquid are constant, i.e., ρ = ρl, µ = µl in liquid phase and ρ = ρg, µ = µg in gas phase and they are subjected to jump conditions at the interface, namely, l g and 6
ACCEPTED MANUSCRIPT
l g . The velocity at the interface is continuous, i.e., u ul ug 0 . The pressure across the interface can be expressed as
p 2 nt u n
,
(3)
where σ is the surface tension, κ is the interface curvature and n is the interface
CR IP T
normal. 2.2 Interface capturing method
The level set method is used to capture the interface, which is implicitly given by
AN US
the iso-surface of the smooth level set function G x, t 0 . Generally, the level set function G is imposed to be the signed distance function to the interface, i.e., G x, t x x ,
(4)
M
where x is the location at the interface that is closest to x. The level set function is defined to be positive for the liquid phase and negative for the gas phase.
ED
However, the level set method suffers the issue of mass loss. Instead of the
PT
signed distance function, the hyperbolic tangent function proposed by Olsson and Kreiss [34][35] is employed here, G x, t 1 , 2
x, t 1 tanh
CE
2
is the thickness of the profile. The evolution of the interface is implicitly
AC
where
(5)
captured by the level set equation,
u 0 . t
(6)
The level set equation is solved by the fifth-order upstream central scheme for spatial discretization and the second-order semi-implicit Crank-Nicolson scheme for time integration. It is noted that will no longer remain the hyperbolic tangent 7
ACCEPTED MANUSCRIPT
profile even at the well-resolved velocity, which can contaminate the interface and lead to mass loss. Consequently, an additional re-initialization equation [35] needs to be introduced to overcome this drawback by solving the following equation to a steady state,
where
(7)
CR IP T
1 n n n ,
is pseudo time. Although, this re-initialization equation can significantly
reduce the mass loss compared with the original re-initialization equation, this
AN US
approach still suffers from mass loss in under-resolved regions. Therefore, we use a mass conserving LS method [32] to further improve the mass conservation of the liquid phase by correcting the volume fraction based on the local curvature at the
M
interface. More information about the mass remedy procedure can be found in Luo et al. [32] in detail.
ED
2.3 Ghost fluid method for the jump conditions
PT
The jump conditions for pressure gradient in the Poisson equation as well as jump condition for the viscous terms are taken into account by the ghost fluid method.
CE
The details of this method can be found in Desjardins et al. [36]. Consider an
AC
interface located at x between the two grid locations xi and xi+1, where xi+1 is within the liquid phase. The pressure jump in the pressure Poisson equation is then written as
1 p p 1 p p l ,i 1 g ,i g ,i 1 * g g ,i p p 1 * 2 , 2 x x g ,i x x * where g l 1 and x xi / x .
8
(8)
ACCEPTED MANUSCRIPT
The pressure jump in Eq. (3) only includes the surface tension force. The continuum surface force (CSF) method is used for the viscous term and the spatial discretization is the second-order finite difference scheme. 2.4 Numerical schemes
CR IP T
The discretization of the Navier-Stokes equations is based on the staggered uniform grid, in which the pressure p and the LS function G are stored at the cell centers, while velocity is stored at the face centers. The spatial discretization of the
AN US
Navier-Stokes equations is performed using the second-order finite central difference schemes. The second-order semi-implicit iterative procedure [33] for time integration is utilized, which is economical, stable and accurate. The iteration can be expressed as (9)
M
f ukn11 u n tf 1 u n ukn 1 1 t ukn11 ukn 1 , 2 2 u
ED
where f is the right hand side of Navier-Stokes equations, and f / u is the Jacobian. The low and upper indices k and
n denote the kth Newton-Raphson sub-iterative
PT
step and the time step, respectively.
CE
The time step t at time t n is determined by the CFL conditions
AC
CFL dy CFL dz CFL dx 2 CFL dy 2 CFL dz 2 t min CFL dx , , , , , , (10) u v w 4 4 4 where CFL is the CFL number that is set to 1. The density and viscosity fields are computed as
g l g ,
(11)
g l g .
(12)
A projection method is used for temporal integration of the Navier-Stokes 9
ACCEPTED MANUSCRIPT
equations. The details are listed as below. (1) Compute an intermediate velocity field u* :
un u* un t un un
(13)
CR IP T
Note that the viscosity is computed following Eq. (7) when the interface crosses a mesh cell.
(2) Compute the pressure field by solving the Poisson equation with GFM: * p n 1 n 1 u . t
AN US
(14)
(3) Compute the real divergence-free velocity field by using the pressure gradient:
u t *
p n 1
n1
.
ED
3 Numerical Simulations
(15)
M
u
n 1
3.1 Flow configuration and computational cases
PT
To mimic fuel atomization in gas turbine combustor, we simulate the swirling jet
CE
injected into surrounding still air using the mass conserving LS method. Fig. 1 shows a schematic of the flow configuration. The exit shape of the nozzle is annular and the
AC
diameters of the outer and inner circle are Dout = 0.4 mm and Din = 0.2 mm, respectively. The initial condition of the liquid phase is a semi-sphere membrane. The liquid is swirling at the exit of the nozzle and the swirl number is defined as S U / U x , where U and U x are the azimuthal velocity and axial velocity,
respectively. The bulk azimuthal velocity and axial velocity are both set to 10 m/s, hence the swirl number is S = 1. The computational domain is 10Dout 10Dout 10Dout. 10
ACCEPTED MANUSCRIPT
The jet direction aligns with the x-direction. The inflow condition at the annular nozzle is velocity inlet and the outflow condition is convective outlet. The boundary conditions at the y and z directions are periodic conditions. Three important forces, including the inertia force l L2V 2 , the surface tension
CR IP T
force L , and the viscous force l LV act on the liquid in atomization (excluding gravity force), where L and V are the characteristic length and velocity, respectively. From these forces, two independent non-dimensional groupings are the Reynolds
AN US
number Re l LV / l and the Weber number We l LV 2 / . Another important number is the Ohnesorge number Oh We / Re l / l L .
The major parameters are listed in Table 1. The characteristic length and velocity
M
are set to be the outer diameter and axial velocity, respectively, resulting in the Re, We and Oh of 2000, 222 and 0.015. The density ratio, viscosity ratio, and the swirl
ED
number are chosen to mimic the realistic conditions in gas turbine combustor, but the
PT
Reynolds number and the Weber number are limited to 103 and 102, respectively, due to our available computer resources.
CE
The computational domain is a rectangular region. The Cartesian grid system is
AC
used for all cases. The cases are listed in Table 2. Cases 1 and 3 are performed to investigate the grid resolution effect on the results. Cases 2 and 3 are performed to investigate the effect of turbulence on the swirling atomization. Case 3 is also conducted to investigate the two-phase characteristics compared to the corresponding single-phase jet. The computation is carried out on the National Supercomputer in Tianjin. Cases 1 and 2 used 448 cores for about 50 hours, and Case 3 used 2048 cores 11
ACCEPTED MANUSCRIPT
for about 160 hours. 3.2 Resolution considerations In order to ensure the adequate resolution of the turbulent atomization, the smallest turbulent scale and liquid structures need to be resolved. The smallest
CR IP T
turbulent scale that needs to be resolved is the Kolmogorov length scale. The liquid structures in the present simulation include thin liquid sheet, ligaments and droplets, and the phase interface thickness. These two issues lead to a very strict constraint on
AN US
the mesh size. As demonstrated in [37][22], such resolution requirements make simulations by using the current state-of-the-art methods unfeasible in the context of the detailed numerical simulation. Therefore, we only focus on the global features of
computational resources.
M
turbulent swirling atomization with the highest resolution depending on our available
ED
Similar to that in [22], the resolution of turbulence is related to the number of
PT
grid points per jet diameter. For the coarse case, the grid resolution x is 7.8µm resulting in 50 mesh points per Dout. For the fine case, the grid resolution x is
CE
3.9µm resulting in 100 mesh points per Dout. The resolution of the liquid structures
AC
follows the consideration of Ménard et al. [25]. The grid-based Weber number is defined as Wex lU 2 x / and it is suggested that no further breakup occurs if
Wex is smaller than about 10. In the present simulations, the grid Weber number Wex is 17.36 for the coarse case and 8.68 for the fine case.
3.3 Turbulent inflow condition Non-turbulent uniform profile inflow is first specified in the simulation for the 12
ACCEPTED MANUSCRIPT
liquid jet. This enables the simulation to exclude the effect of inflow turbulence on the swirling atomization. Following [15]-[17], the influence of turbulent inflow condition is then investigated. In order to generate the realistic turbulent inflow condition, a precursor simulation of a periodic swirling pipe flow is performed using the liquid
CR IP T
properties. The cylindrical mesh is used and the computational domain is [0, 10D], [0.25D, 0.5D], [0, 2 ] in x, r, directions, respectively. The velocities in the azimuthal and axial directions are both set to 10 m/s. The domain is discretized on a
AN US
192 64 48 mesh. Fig. 2 shows the contours of instantaneous velocity magnitude with random fluctuations in the longitudinal direction of the pipe flow and the coupling with the external flow. Once the flow reaches the statistically stationary state, and L are
M
the velocity field is recorded within the time period U / L , where U
streamwise velocity magnitude and 10D, respectively. The Cartesian velocities in the
ED
x, y, z directions then read
PT
U U x , V U cos arctan
y , W U sin arctan z
y . z
(16)
CE
For the laminar inflow case, Eq. (16) is directly applied to the velocity boundary condition. For the turbulent inflow case, a swirling turbulent pipe flow simulation is
AC
first conducted and then stored as the time-dependent inlet velocity of the external flow in the combustor for swirling atomization. Since the scale of the pipe is equal to the outer diameter of the inflow of the combustor, the stored pipe inflow date file can be directly read for the external flow calculations. Moreover, the mean velocity and root mean square (rms) profiles obtained at the nozzle exit are illustrated in Fig. 3. The mean velocity increases sharply from the wall. 13
ACCEPTED MANUSCRIPT
The turbulent intensity is about 20% of the mean velocity and the rms velocity reaches a maximum value near the wall. Although it is not validated by the experimental results, the present velocity and rms profiles are similar to that in [15]-[17], so we will investigate the effect of turbulent inflow on the atomization by
CR IP T
utilizing this inflow condition.
4 Results and Discussions 4.1 Global atomization characteristics
AN US
In the later stage of the sheet breakup for the laminar case, some finger-like structures appear at the rim of the sheet, which can be seen in the simulation results in Fig. 4(a) with the coarse mesh and in Fig. 4(b) with the fine mesh. The experimental
M
result in Taylor [38] is also shown in Fig. 4(c) for comparison. In [38], a liquid sheet of kerosene issuing from the orifice with the diameter 0.3 mm of a swirl atomizer
ED
under pressure of 1.36 atm into air at pressure 0.32 atm. The visualizations of the
PT
simulation and experimental results are qualitative similar despite of the different conditions. However, the finger-like structures in the numerical simulation results are
CE
generally longer than those in the experiment. This difference seems to be attributed
AC
to the liquid property or the insufficient grid resolution that leads to unphysical breakup. It can be seen from Figs. 4(a) and 4(b) that ligament length is significantly reduced in the case by refining grid resolution, which implies that the grid resolution has an impact on the breakup of swirling sheet in the numerical results. On the other hand, continuously increasing the grid resolution seems unfeasible considering that the sheet thickness may be down to tens of nanometers. Therefore, the present study 14
ACCEPTED MANUSCRIPT
focuses on the large-scale characteristics of spray formation and the effects of grid resolution on the overall characteristics, including the sheet thickness, breakup length, recirculation zone and spray cone angle. 4.1.1 Liquid sheet thickness
CR IP T
The thickness of the liquid sheet at the nozzle exit has an impact on the drop sizes in the primary atomization. Taylor [38] speculated that the spiral variation of the sheet thickness at the nozzle exit and the oscillating air core could lead to the
AN US
variation of the liquid thickness in the conical sheet. The multi-dimensional spray modeling [39] predicts that the sheet thickness at the breakup length is related to the cone angle, the maximum growth rate and the liquid velocity.
M
Fig. 5 shows the overall interface spatial and temporal evolution of the laminar swirling atomization in Case 1. Here, the dimensionless time T = tU/D is used where
ED
U is the liquid axial velocity. At the first stage of the injection, the liquid structure
PT
forms the mushroom-shaped tip at T = 1.25 in Fig.5 due to the liquid penetration by impingement against the still air, which is similar to that shown in Shinjo et al. [18].
CE
Some ligaments are first generated at this rolled-up structures and then breaks into
AC
droplets at the outer recirculation zone (ORZ) surrounding the conical sheet and at T = 2.50 in Fig.5. The liquid film spreads downstream smoothly until T = 3.75. In the present simulation, the sheet variation shows a limited dependence on the
radial disturbances for the laminar inflow condition, which can be evident in the smooth interface sheet in Fig. 5. The thickness h of the film can be calculated from the mass flow rate equation [1] 15
ACCEPTED MANUSCRIPT
m lVave h 2 x tan h ,
(17)
where Vave is the average axial velocity along the liquid sheet, x and are the stream-wise position and the half cone angle, respectively. Eq. (17) is based on the assumption that liquid sheet does not rupture before the breakup length. Although the
CR IP T
simulation result shows a few of droplets at the ORZ pinched off from the tip of the conical sheet in Fig. 5, it is still reasonable to assume the corresponding volume is much smaller than that of the liquid sheet.
AN US
The spatial variations of the sheet thickness for laminar inflow conditions at different grid resolutions at T = 12.5 are shown in Fig. 6. The thickness of the sheet decreases sharply near the nozzle exit with the increasing stream-wise position and
M
then decreases slowly. Before the location x/Dout = 1.5, the thicknesses of the sheet for different grid resolutions are collapsed to each other. It is noted that the thickness of
x tan
PT
h x tan
ED
the sheet for the fine case is in accordance with the physically reasonable solution 2
m lVave
(18) in the
CE
of Eq. (17), where m 7.536 104 kg / m3 , Vave 8.15m / s and 45
present simulation. The difference mainly locates at the film breakup zone, comparing
AC
the coarse and fine meshes of the laminar cases. For the numerical simulations, the liquid sheet may occur non-physical breakup at the scale of the grid size. Although 10243 grids have been adopted in the present simulation, there are about 2.5 cells across this thickness where the sheet begins to breakup, which still doesn’t guarantee the sufficient resolution of the sheet breakup. However, the results prove to capture the major characteristics of the thickness of the sheet. Moreover, the liquid sheet 16
ACCEPTED MANUSCRIPT
breaks at about x/Dout = 2 for the coarse mesh case as shown in Fig. 6. The sheet thickness at this position is h 2.5x , which should be carefully investigated. For the fine mesh case, the sheet thickness at the breakup position is about h 5x . 4.1.2 Breakup length and recirculation zone
CR IP T
The breakup length of the liquid sheet is defined as the distance from the nozzle exit to the position where the sheet begins to break. The recirculation zone is defined as a flow transition with a free stagnation point/region on the axis followed by a
AN US
reverse flow. These two features can be observed in Fig. 7, which shows the overlaps of the interface among different plane-cuts for the coarse mesh (a) and the fine mesh (b) at the statistically stationary state T = 12.5. The plane-cuts along the azimuthal
M
direction are projected on the x-y plane in Fig. 7. It illustrates that the breakup of the liquid length occurs around x/Dout = 2~3 for Case 1 with the coarse mesh and x/Dout =
ED
3~4 for Case 3 with the fine mesh. The recirculation zones are formed around x/Dout =
PT
2~5 for both cases. Therefore, the grid resolution has an impact on the breakup length, but refining grid resolution does not change the general swirling-flow characteristics
CE
significantly.
AC
4.1.3 Spray cone angle The spray cone angle is an important spray characteristic for evaluating atomizer
performance. The cone angle is defined as the angle between the tangents to the envelope of the spray. Qualitatively, the spray cone angle increases with time and reaches to a certain value. In the first stage, the axisymmetric vortex pair is generated, entraining the surrounding gas to the swirl core, and bend the liquid film towards the 17
ACCEPTED MANUSCRIPT
axis. In the later stage, the vortex breaks down and the spray cone angle increases. Quantitatively, the temporal evolutions of the spray angle with different grid resolutions are shown in Fig. 8. It shows that the evolutions of the spray angle are almost the same for different grid resolutions.
CR IP T
In summary, the grid resolution has effects on the breakup zone of the liquid sheet, but both Case 1 and Case 3 can capture the major characteristics of the swirling flow.
AN US
4.2 Turbulent Inflow Effect 4.2.1 Flow structures
We compare qualitative differences of the swirl atomization between laminar
M
inflow and turbulent inflow in Figs. 5 and 9. In the turbulent inflow condition, the helical wave (see Fig. 9, T = 1.25), which is not observed in the laminar inflow
ED
condition, appears due to the inlet turbulence disturbance. The helical wave and
PT
turbulent liquid inflow make the conical sheet significantly wrinkled. The conical sheet develops for a short distance and then immediately ruptures into ligaments and
CE
droplets (see Fig. 9, T = 3.75 – 6.25). Note that the breakup length of the conical sheet
AC
is much shorter than that in the laminar inflow condition. Moreover, the mean length of the ligaments is significantly reduced due to the turbulent straining motion that enhances the rupture of the ligaments. From Figs. 5 and 9, we conclude that turbulent inflow enhances the flow instability, which results in the droplets disperse uniformly in the radial direction after atomization. In order to further investigate the turbulent inflow effect, we overlap different 18
ACCEPTED MANUSCRIPT
plane-cuts of the interface with the turbulent inflow condition at T = 12.5 in Fig. 10. At x/Dout = 0.5-1.5 near the nozzle exit, the liquid sheets do not overlap with each other, which implies that the sheet thickness fluctuates along the azimuthal direction and it is different from the laminar inflow case. In the turbulent inflow case, the
CR IP T
small-scale liquid structures after breakup are gathered around x/Dout = 2 and y/Dout = 3. While for the laminar inflow case, those structures are gathered in the downstream of x/Dout = 4. This shows that the turbulent inflow can promote the liquid sheet
AN US
breakup near the nozzle exit and enhance the radial movement of the liquid structures. Additionally, the instantaneous interfaces at T = 12.5 for different inflow conditions viewed from upstream are shown in Fig. 11. The liquid sheet atomization
M
in the laminar inflow condition displays a two-step cascade, including the process from the sheet to ligaments and from ligaments to droplets, which is different from
ED
that in [40][41]. The ligaments hardly occur in the turbulent inflow condition in
PT
Fig.11(b). Therefore, the turbulence inflow significantly enhances the liquid atomization, which agrees with the result in the crossflow simulation [15].
CE
The different inflow conditions lead to different characteristics of liquid
AC
dispersion. In the laminar inflow case, the liquid structures experience the two-step cascade from the sheet through ligaments to droplets. The liquid dispersion for the turbulent inflow case is quite different. The liquid sheet starts to undulate along the azimuthal direction in the upstream, and then liquid drops formed undergo the radial dispersion, both toward and outward the centerline. The liquid structures in the downstream shows approximately a statistically uniform distribution in the radial 19
ACCEPTED MANUSCRIPT
direction. It is demonstrated that the turbulence induced in the liquid inflow can weaken the stiffness of the sheet and leads to more even distribution of the liquid structures in the radial direction. 4.2.2 Recirculation characteristics
CR IP T
In order to understand the structure characteristics of the gas-liquid swirling recirculation zone generated by vortex breakdown, the mean velocity profiles on x-y plane-cuts in Fig. 12. Although there are extensive studies on vortex breakdown in
AN US
single phase flow [42][43][44], there are few studies on vortex breakdown in gas-liquid swirling flow. In [42] for the single phase flow, vortex breakdown has been defined as an abrupt flow transition with a free stagnation point/region on the axis
M
followed by a reverse flow and a fully turbulent region. Fig. 12 shows the instantaneous axial velocity field of the swirling gas-liquid flow in the laminar inflow
ED
and in the turbulent inflow condition at T = 15.0. Compared to the single phase swirl
PT
flow, the gas-liquid swirl flow has intense fluctuations. Two distinct features can be observed in Fig. 12. First, the distance of the
CE
recirculation zone location from the nozzle is longer for the laminar inflow condition.
AC
From Fig. 12(a), the axial velocity is approximately zero when x/Dout = 1.5. The central recirculation zone occurs at x/Dout = 1.5 - 4.5 for laminar inflow, while the central recirculation zone is at x/Dout = 1.0 - 3.0 for turbulent inflow. It is interesting to note that the maximum negative axial velocity occurs at the liquid sheet breakup into filaments and droplets, x/Dout = 2.5 for laminar inflow and x/Dout = 1.5 for turbulent inflow. Second, the size of the breakdown in the laminar inflow case is nearly twice as 20
ACCEPTED MANUSCRIPT
large as that in the turbulent inflow case. The possible reason for the differences is that the liquid sheet can maintain its momentum in the radial and axial directions in the laminar inflow case, so the momentum exchange from liquid phase to gas phase is in further downstream. A primary axisymmetric vortex pair can be observed at the
CR IP T
location where liquid sheet breaks into droplets for the laminar inflow condition. However, the recirculation zone is suppressed in the turbulent inflow case due to the turbulent sheet breakup at a short distance from the nozzle and the subsequent vortex
AN US
breakdown.
Fig. 13 shows the axial and radial velocity profiles at different axial positions for laminar and turbulent inflows at T = 15.00. In Fig. 13(a), the axial velocity decreases
M
in the axial direction between the inlet and the stagnation point in both laminar and turbulent inflow cases owing to the adverse pressure gradient [45]. Moreover, the
ED
radial distance of the recirculation zone in the turbulent inflow case is shorter than
PT
that in the laminar inflow case, which agrees with the observation in Fig. 12. 4.3 Characteristics of vortices and interfaces
CE
4.3.1 Precessing vortex core (PVC)
AC
The precessing vortex core (PVC) is known as the rotating vortical structures in
the center of the swirling jet. Its presence is significant for the instability and the combustion system. Previous studies found that PVC occurs when the swirl number S > 0.6-0.7 [49] in single phase flow. However, it is an open question if the PVC exists in two-phase annular swirling jets. To identity PVC, we display the velocity vectors for the single-phase jet, the laminar and turbulent inflow cases of the two-phase jet on the 21
ACCEPTED MANUSCRIPT
plane-cuts at x/Dout = 3. For the single phase case, the PVC is evident in the center of the jet in Fig. 14(a) with the occurrence of the maximum axial velocity spinning around the origin, which agrees with the simulation result in [46] and the experimental observations in [50]. On the other hand, there is no obvious PVC in both
CR IP T
two-phase jet cases with laminar and turbulent inflow conditions, which was also observed in [46] where the swirl number is 0.4. It appears that the two-phase cases have strong radial momentum to weaken the vortical flow in the center of the jet.
AN US
In the two-phase cases, the complex vortical structures appear in the downstream. For the laminar inflow case, the clockwise regular velocity vectors occur at near nozzle exit (not shown), where the liquid sheet has not broken up. In the downstream,
M
the velocity vectors along the large circular rim are oriented to the positive radial direction, which indicates the liquid drops disperse to the outward of the sheet. In the
ED
inner zone of the liquid sheet, the velocity field is chaotic with numerous vortical
PT
structures in Fig. 14(b), which can enhance mixing of the drops and the ambient gas. For the turbulent inflow case, the velocity vectors at the downstream in Fig. 14(c) are
CE
more chaotic that those in the laminar inflow case.
AC
In contrast, the single phase case exhibits very regular flow structures. The
clockwise rotation of the velocity vectors appears at the upstream and the symmetric triple-vortices with the center are formed in Fig. 14(a). This can be attributed to the helical wave development along the swirling jet. In Fig. 14(a) the radial zone induced by swirling flow at z/Dout = [-2, 2], so the radial development in the single phase case is compact compared to the two-phase cases. With the same inlet velocity, the liquid 22
ACCEPTED MANUSCRIPT
phase has much larger momentum than that in the single gas-phase case to penetrate to streamwise and radial directions, leading to the wide spread of liquid structures. 4.3.2 Correlation of the interface and vortices The gas-liquid interfaces in Case 3 at T = 6.25 are shown in Fig. 15, along with
CR IP T
the vortical structures characterized as the iso-surfaces of the Q-criterion. Here, Q Aij Aij Sij Sij / 2 is the second invariant of the velocity gradient tensor, where Aij and Sij are the asymmetric and symmetric parts of the velocity gradient tensor,
AN US
respectively. It can be seen that high vorticities concentrate on the breakup region. The preferred normal alignment of the LS scalar gradient with the vorticity field near the interface can be quantified by the cosine of the angle between vorticity vector and
,
M
LS scalar gradient vector as
(19)
ED
where verticity vector u . Fig. 16 shows the PDF of at T = 6.25 that has
PT
a strong peak at 0 . Here, 0 means strong perpendicularity and 1
CE
no perpendicularity. The strong perpendicularity of the interface and the vortices is also in accord with the result in [53] for homogeneous isotropic turbulence.
AC
In order to examine the dynamics in small-scale turbulence, the enstrophy
evolution equation for two phase flow is utilized as
p 2 S D 1 2 , S Dt 2
(2 0)
t where D / Dt / t u , S u u / 2 , and the terms on the right hand side
23
ACCEPTED MANUSCRIPT
are the enstrophy production due to vortex stretching, the baroclinic term, and viscous dissipation term, respectively. The LS function G can be seen as a passive scalar, and the small-scale passive scalar structure can be analyzed in terms of the evolution equation for the scalar
CR IP T
gradient as D 1 2 Gi Gi Sij G j , Dt 2
(21)
where Gi G / xi , and the term on the right-hand side is the production.
AN US
The turbulent production terms for ii and Gi Gi are i Sij j and Gi Sij G j , respectively. The alignments among the vorticity i , the scalar gradient Gi and the principal rate of Sij are significant to determine the turbulent production terms
M
[57][58]. Fig. 17 shows the PDFs for the cosine of the angle between the vorticity and
ED
the principal strain rate directions in an area of the recirculation zone (encircled by the solid lines in Fig. 7b and Fig. 12a), where ai, bi, ci are the eigenvectors corresponding
are
real
and
CE
c
PT
to the directions of the eigenvalues a, b, c, respectively. Since Sij is symmetric, a, b, satisfy
abc 0
(solenoidal
condition)
and
a b c (a 0 c) . Here, cos 1 means strong alignment and cos 0
AC
no alignment. It shows in Fig.17 that
i is aligned with the intermediate (b) strain
rate, which is consistent with the findings of Ashurst et al. (1987) for isotropic turbulence and Abe et al. (2009) of far-wall region for turbulent channel flow. In the recirculation zone, the turbulent eddies are generated by the wakes formed due to droplets dynamics. The gas phase flow in the recirculation region is modulated and exhibits similar characteristics to the isotropic turbulence as shown in Fig. 17. This 24
ACCEPTED MANUSCRIPT
indicates that the present recirculation zone of swirling two-phase jet also contain the similar small-scale structures as those in isotropic turbulence. In the near-sheet regions (encircled by the dashed lines in Fig. 7b and Fig. 12a),
i is also found to be aligned with the intermediate (b) strain rate as shown in Fig.
CR IP T
18, which agrees with the findings of Abe et al. (2009) for the near-wall region for turbulent channel flow and Shinjo et al. (2015) [59] for the near-droplet regions for dense spray. Moreover, the alignments of LS gradient Gi with the extensional (a)
AN US
and the compressive (c) strain rates (Fig. 19) tend to an orientation angle of 45° ( cos 0.7 ). Both trends of the alignments among
i , Gi and the principal strain
rate are qualitatively similar to those in the near-wall zone for the turbulent channel
M
flow[58], which indicates that the liquid sheet can generate high shear layers to
5 Conclusions
ED
produce anisotropic small-scale fluctuations.
PT
We report the spray characteristics of the swirling liquid primary atomization. The primary atomization is simulated using the detailed numerical simulation with the
CE
recently developed mass conservative LS method and the ghost fluid method handling
AC
jump conditions. The interface is represented by an iso-surface of a hyperbolic tangent function which is kept with a re-initialization process. Through comparing the sheet thickness, the breakup length and the cone angle, we demonstrate that the grid resolution is satisfactory from the convergence of the major characteristics of the swirling two-phase flow. The liquid sheet thickness for the fine case agrees well with theoretical analysis. 25
ACCEPTED MANUSCRIPT
Compared with the single phase jet, the two-phase jet exhibits the chaotic velocity field downstream. The momentum carried by the liquid-phase fluctuates intensely, which can suppress the generation of PVC in the center of the jet. Compared to the laminar inflow case, the turbulent inflow reduce the stiffness of the
CR IP T
liquid sheet leading to statistically homogeneous distribution of the liquid structures in the radial direction. In addition, the recirculation zone is smaller and farther from the nozzle exit for the turbulent inflow case. The interaction of the liquid-gas interface
AN US
and vortices is characterized by the preferential normal alignment between the gradient of the LS scalar and vorticity. The alignments among the vorticity, the level-set scalar gradient and the principal rate of the symmetric part of the velocity
M
gradient tensor agree with the experimental results, and indicate that small-scale dynamics in the recirculation zone is similar to that in isotropic turbulence and is
ED
shows anisotropic statistics near the liquid sheet.
PT
Acknowledgement
CE
This work is financially supported by the National Natural Science Foundation of China (Nos. 91541202, 51222602, and 51276163). Yue Yang acknowledges the
AC
support by the National Natural Science Foundation of China (Nos. 11522215 and 91541204) and the Young Thousand Talent Program of China.
References [1] N. Ashgriz, Handbook of atomization and sprays: Theory and Applications, (2011) [2] X. Li, M. C. Soteriou, W. Kim, J. M. Cohen, High fidelity simulation of the 26
ACCEPTED MANUSCRIPT
spray generated by a realistic swirling flow injector, Journal of Engineering for Gas Turbines and Power, 136(2014)071503-1-10 [3] G. Tryggvason, R. Scardovelli, S. Zaleski, Direct numerical simulations of gas-liquid multiphase flows, Cambridge University Press, 2011
CR IP T
[4] S. Osher, J. A. Sethian, Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulation. J. Comp. Phys., 79(1988)12 [5] D. Kim, O. Desjardins, M. Herrmann, P. Moin. Toward two-phase simulation of
AN US
the primary breakup of a round liquid jet by a coaxial flow of gas, Center for Turbulence Research, Annual Research Briefs, 2006
[6] D. Kim, P. Moin, Numerical simulation of the breakup of a round liquid jet by a
M
coaxial flow of gas with a subgrid Lagrangian breakup model, Center for Turbulence Research, Annual Research Briefs, 2011
ED
[7] M. Herrmann. A balanced force refined level set grid method for two-phase
PT
flows on unstructured flow solver grids. J. Comput. Phys., 227,2674-2706(2008) [8] M. Herrmann. On simulating primary atomization using the refined level set grid
CE
method. Atom. Sprays, 21(2011)283-301
AC
[9] M. Herrmann, M. GoroK-Hovski, An outline of an LES subgrid model for liquid/gas phase interface dynamics, Center for Turbulence Research, Proceedings of the Summer Program, 2008
[10] M. Herrmann. Detailed numerical simulations of the primary atomization of a turbulent liquid jet in crossflow. Journal of Engineering for Gas Turbines and Power, 132(2009)061506-1-10 27
ACCEPTED MANUSCRIPT
[11] M. Herrmann. Detailed simulations of the breakup processes of turbulent liquid jets in subsonic crossflows. In 11th Triennial International Annual Conference on Liquid Atomization and Spray Systems, ICLASS 2009 [12] M. Herrmann. A parallel Eulerian nterface tracking/Lagrangian point particle
CR IP T
multi-scale coupling procedure. J. Comput. Phys., 229(2010)745-759 [13] M. Herrmann, M. Arienti, M. Soteriou. The impact of density ratio on the liquid core dynamics of a turbulent liquid jet injected into a crossflow, Journal of
AN US
Engineering for Gas Turbines and Power, 133(2011)061501-1-9
[14] M. Herrmann. The influence of density ratio on the primary atomization of a turbulent liquid jet in crossflow. Proceedings of the Combustion Institute, 33
M
(2011) 2079-2088
[15] F. Xiao, M. Dianat, J. J. McGuirk, Large eddy simulation of liquid-jet primary
ED
breakup in air crossflow. AIAA Journal, 51(2013)2878-2893
PT
[16] F. Xiao, M. Dianat, J. J. McGuirk, Large eddy simulation of single droplet and liquid jet primary breakup using a coupled level set/volume of fluid method.
CE
Atomization and Sprays, 24(2014)281-302
AC
[17] F. Xiao, M. Dianat, J. J. McGuirk, LES of turbulent liquid jet primary breakup in turbulent coaxial air flow, International Journal of Multiphase Flow, 60(2014) 103-118
[18] J. Shinjo, A. Umemura, Simulation of liquid jet primary breakup: dynamics of ligament and droplet formation, International Journal of Multiphase Flow, 36 (2010)513-532 28
ACCEPTED MANUSCRIPT
[19] J. Shinjo, A. Umemura, Detailed simulation of primary atomization mechanisms in Diesel jet sprays (isolated identification of liquid jet tip effects), Proceedings of the Combustion Institute, 33(2011)2089-2097 [20] J. Shinjo, A. Umemura, Surface instability and primary atomization
CR IP T
characteristics of straight liquid jet sprays, International Journal of Multiphase Flow, 37(2011)1294-1304
[21] O. Desjardins, V. Moureau, H. Pitsch, An accurate conservative level set/ghost
AN US
fluid method for simulation turbulent atomization, J. Comput. Phys., 227(2008) 8395-8416
[22] O. Desjardins, H. Pitsch, A spectrally refined interface approach for simulating
M
multiphase flows, J. Comput. Phys., 228(2009)1658-1677
[23] O. Desjardins, H. Pitsch, Detailed numerical investigation of turbulent
ED
atomization of liquid jets, Atomization and sprays, 20(2010)311-336
PT
[24] O. Desjardins, J. O. McCaslin, M. Owkes, P. Brady, Direct numerical and large-eddy simulation of primary atomization in complex geometries,
CE
Atomization and sprays, 23(2013)1001-1048
AC
[25] T. Ménard, S. Tanguy, A. Berlemont, Coupling level set/VOF/ghost fluid methods: Validation and application to 3D simulation of the primary break-up of a liquid jet, International Journal of Multiphase Flow, 33(2007)510-524
[26] R. Lebas, T. Ménard, P. A. Beau, A. Berlemont, F. X. Demoulin, Numerical simulation of primary break-up and atomization: DNS and modeling study, International Journal of Multiphase Flow, 35(2009)247-260 29
ACCEPTED MANUSCRIPT
[27] X. Li, M. C. Soteriou, M. Arienti, M. M. Sussman, High-fidelity simulation of atomization and evaporation in a liquid jet in cross-flow, 49th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition, 2011
CR IP T
[28] X. Li and M. C. Soteriou, High-fidelity simulation of fuel atomization in a realistic swirling flow injector, Atomization and Sprays, 23(2013)1049-1078
[29] D. Fuster, A. Bagué, T. Boeck, L. L. Moyne, A. Leboissetier, S. Popinet, P. Ray,
AN US
R. Scardovelli, S. Zaleski, Simulation of primary atomization with an octree adaptive mesh refinement and VOF method, International Journal of Multiphase Flow, 35(2009)550-565
M
[30] M. Arienti and M. Sussman, An embedded level set method for sharp-interface multiphase simulations of diesel injectors, International Journal of Multiphase
ED
Flow, 59(2014)1-14
PT
[31] X. Chen, D. Ma, V. Yang, S. Popinet, High-fidelity simulations of impinging jet atomization, Atomization and Sprays, 23(2013)1079-1101
CE
[32] K. Luo, C. X. Shao, Y. Yang, J. R. Fan. A mass conserving level set method for
AC
direct numerical simulation of liquid atomization. J. Comput. Phys., 298 (2015) 495-519
[33] C. D. Pierce, P. Moin. Progress-variable approach for large eddy simulation of turbulent combustion. Technical Report TF80, Flow Physics and Computation Division, Dept. Mech. Eng., Standford University, 2001 [34] E. Olsson and G. Kreiss. A conservative level set method for two phase flow. J. 30
ACCEPTED MANUSCRIPT
Comput. Phys., 210,225-246(2005) [35] E. Olsson, G. Kreiss, and S. Zahedi. A conservative level set method for two phase flow ii. J. Comput. Phys., 225,785-807(2007) [36] O. Desjardins, V. Moureau, H. Pitsch, An accurate conservative level set/ghost
CR IP T
fluid method for simulating turbulent atomization, J. Comput. Phys., 227(2008) 8395-8416
[37] M. GoroK-Hovski and M. Herrmann, Modeling primary atomization, Annu. Rev.
AN US
Fluid Mech., 40(2008)343-366
[38] G. Taylor, The dynamics of thin sheets of fluid. III. Disintegration of fluid sheets, Proc. R. Soc. Lond. A, 253(1959)313-321
M
[39] P. K. Senecal, D. P. Schmidt, I. Nouar, C. J. Rutland, R. D. Reitz, M. L. Corradini, Modeling high-speed viscous liquid sheet atomization, International
ED
Journal of Multiphase Flow, 25(1999)1073-1097
PT
[40] Y. Wang, K. S. Im, K. Fezzaa, Similarity between the primary and secondary air-assisted liquid jet breakup mechanisms, Physical Review Letters, 100(2008)
CE
154502
AC
[41] E. A. Novikov and D. G. Dommermuth, Distribution of droplets in a turbulent spray, Physical Review E, 56(1997)5479-5482
[42] H. Liang and T. Maxworthy, An experimental investigation of swirling jets, J. Fluid Mech., 525(2005)115-159 [43] O. L. Negro and T. O’Doherty, Vortex breakdown: a review, Progress in Energy and Combustion Science, 27(2001)431-481 31
ACCEPTED MANUSCRIPT
[44] M. R. Ruith, P. Chen, E. Meiburg, T. Maxworthy, Three-dimensional vortex breakdown in swirling jets and wakes: direct numerical simulation, J. Fluid Mech., 486(2003)331-378 [45] E. Krause, The solution to the problem of vortex breakdown, In: Lecture notes
CR IP T
in physics, vol.371, Berlin: Springer, (1990)35-50 [46] G. A. Siamas, X. Jiang, L. Wrobel, Direct numerical simulation of the near-field dynamics of annular gas-liquid two-phase jets, Physics of Fluids, 21 (2009)
AN US
042103
[47] W. Kollmann, A. S. H. Ooi, M. S. Chong, J. Soria, Direct numerical simulations of vortex breakdown in swirling jets, J. Turbul., 2 (2001) N5
M
[48] K. Manoharan, S. Hansford, J. O’Connor, S. Hemchandra, Instability mechanism in a swirl flow combustor: Precession of vortex core and influence
ED
of density gradient, Proceedings of ASME Turbo Expo, 2015
PT
[49] A. K. Gupta, D. J. Lilley, N. Syred, Swirl flows, Tunbridge Wells, UK: Abacus Press, 1984
CE
[50] N. Syred, A review of oscillation mechanisms and the role the precessing vortex
AC
core (PVC) in swirl combustion systems, Prog. Energy Combust. Sci., 32 (2006) 93-161
[51] K. Luo, H. Pitsch, M. G. Pai, O. Desjardins, Direct numerical simulations and analysis of three-dimensional n-heptane spray flames in a model swirl combustor, Proceedings of the Combustion Institute, 33, 2143-2152 (2011) [52] P. Moin and S. V. Apte, Large-eddy simulation of realistic gas turbine 32
ACCEPTED MANUSCRIPT
combustors, AIAA Journal, 44, 698-708 (2006) [53] Y. Yang, D. I. Pullin, I. B. Moreno. Multi-scale geometric analysis of Lagrangian structures in isotropic turbulence. J. Fluid Mech.,654,233-270,(2010) [54] Brown, G. & Lopez, J. Axisymmetric vortex breakdown Part 2. Physical
CR IP T
mechanisms. J. Fluid Mech. 221, 553 (1990) [55] Mahesh, K. A model for the onset of breakdown in an axisymmetric compressible vortex. Phys. Fluids 8, 3338 (1996)
AN US
[56] Squire, H. B. Analysis of the vortex breakdown phenomenon. In Miszellaneen der Angewandten Mechanik, p. 306. Akademie-Verlag, Berlin (1962) [57] W. T. Ashurst, A. R. Kerstein, R. M. Kerr, C. H. Gibson, Alignment of vorticity
M
and scalar gradient with strain rate in simulated Navier-Stokes turbulence, Physics of Fluids, 30, 2343 (1987)
ED
[58] H. Abe, R. A. Antonia, H. Kawamura, Correlation between small-scale velocity
(2009)
PT
and scalar fluctuations in a turbulent channel flow, J. Fluid Mech., 627, 1-32
CE
[59] J. Shinjo, J. Xia, A. Umemura, Droplet/ligament modulation of local small-scale
AC
turbulence and scalar mixing in a dense fuel spray, Proceedings of the Combustion Institute, 35 (2015) 1595-1602
33
CR IP T
ACCEPTED MANUSCRIPT
PT
ED
M
AN US
Fig. 1: Schematic of flow configuration used for swirling film atomization.
AC
CE
Fig. 2: Contours of the instantaneous velocity magnitude inthe longitudinal directions of the pipe flow and the coupling with the external flow.
34
ED
M
AN US
Fig.3: Mean velocity (a) and rms (b) profiles at the liquid inlet.
CR IP T
ACCEPTED MANUSCRIPT
AC
CE
PT
Fig. 4. The breakup of the gas/liquid interface, (a) Simulation result from Case 1 with a coarse mesh; (b) Simulation result from Case 3 with a fine mesh; (c) Experimental result in Taylor [38]
Fig. 5: The interface evolution of the laminar swirling atomization for Case 3(T= 0, 1.25, 2.50, 3.75, 5.00, and 6.25 from left to right).
35
CR IP T
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
AN US
Fig. 6: The spatial variations of the sheet thickness for the laminar inflow condition at different grid resolutions at T= 12.5.
Fig. 7: The overlap of the interfaces on different plane-cuts for coarse mesh (a) and fine mesh (b) at T= 12.5
36
CR IP T
ACCEPTED MANUSCRIPT
M
AN US
Fig. 8: The temporal evolution of the spray angle for different grid resolutions.
AC
CE
PT
ED
Fig. 9: The interface evolution of the turbulent swirling atomization for Case 2(T= 0, 1.25, 2.50, 3.75, 5.00, and 6.25 from left to right).
Fig. 10: The overlap of the interfaces on different plane-cuts for the turbulent inflow condition at T= 12.5 37
CR IP T
ACCEPTED MANUSCRIPT
PT
ED
M
AN US
Fig. 11: The instantaneous interfaces at T= 12.5 for different inflow conditions at a clockwise rotating motion viewed from upstream. (a) laminar inflow; (b) turbulent inflow.
AC
CE
Fig. 12: The instantaneous axial velocity field of the swirling gas-liquid flow (a) in the laminar inflow case and (b) in the turbulent inflow condition at T= 15.00. The black lines denote the stream lines.
38
AN US
CR IP T
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
Fig. 13: The velocity profiles at different axial positions at T= 15.00. (a) axial velocity, (b) radial velocity. The solid and dash lines denote for laminar inflow and turbulent inflow, respectively.
39
AC
CE
PT
ED
M
AN US
CR IP T
ACCEPTED MANUSCRIPT
Fig. 14: Instantaneous velocity vector and axial velocity contour at x/Dout= 3 at T= 12.5 (a) single phase; (b) laminar inflow; (c) turbulent inflow.
40
AN US
CR IP T
ACCEPTED MANUSCRIPT
CE
PT
ED
M
Fig. 15: Interface-vortices interaction for fine case at T= 6.25 (white color for interface and green color for iso-surfaces of the Q criterion).
AC
Fig. 16: PDF of for T= 6.25.
41
CR IP T
ACCEPTED MANUSCRIPT
ED
M
AN US
Fig. 17: PDFs of the alignment between the principal strain rate and i in the recirculation zone.
Fig. 18: PDFs of the alignment between the principal strain rate and i in near regions of the
AC
CE
PT
liquid sheet structures.
Fig. 19: PDFs of the alignment between the principal strain rate and the LS gradient Gi in near regions of the liquid sheet structures.
42
ACCEPTED MANUSCRIPT
Table 1: Comparison of major parameters between realistic gas turbine combustor or the model combustor [28][51][52] and the parameters in the present detailed numerical simulation Parameters Turbine combustor Present simulation 90o
Swirl number
1.0
Density ratio
O(10-1000)
Dynamic viscosity ratio
O(10-100)
Surface tension
0.02-0.04
Reynolds number
O(106)
20 10
0.036
AN US
2000 222
M
O(10 )
Domain size
Number of grids
1 2 3
10D 10D 10D 10D 10D 10D 10D 10D 10D
5123 5123 10243
ED
Case
Grid resolution (µm) 7.8 7.8 3.9
AC
CE
PT
1.0
6
Weber number
Table 2: Simulation cases
90o
CR IP T
Spray cone angle
43
Inflow Laminar Turbulent Laminar
Detailed Numerical Simulation of Swirling Primary Atomization Using a Mass Conservative Level Set Method Changxiao Shao , Kun Luo , Yue Yang , Jianren Fan PII: DOI: Reference:
S0301-9322(16)30620-6 10.1016/j.ijmultiphaseflow.2016.10.010 IJMF 2489
To appear in:
International Journal of Multiphase Flow
Received date: Revised date: Accepted date:
13 May 2015 28 June 2016 22 October 2016
Please cite this article as: Changxiao Shao , Kun Luo , Yue Yang , Jianren Fan , Detailed Numerical Simulation of Swirling Primary Atomization Using a Mass Conservative Level Set Method , International Journal of Multiphase Flow (2016), doi: 10.1016/j.ijmultiphaseflow.2016.10.010
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
Highlights
We report detailed numerical simulations of swirling liquid atomization by using a recently developed mass conservative level set method.
Through comparing the sheet thickness, the breakup length and the cone angle, the numerical
CR IP T
convergence of the global characteristics of the swirling two phase flow has been obtained. The numerical results show that turbulent inflow can induce liquid sheet breakup near the nozzle exit, reduce the stiffness of the liquid sheet, and lead to the statistically homogeneous
AN US
distribution of small-scale liquid structures in the radial direction.
Compared with the single-phase jet, the two-phase jet exhibits the chaotic velocity filed downstream that can enhance the mixing of droplets and ambient gas, and the precessing vortex core (PVC) is not observed in the center of the two-phase jet. The preferential alignment of
i
M
with the intermediate strain rate indicates that the
AC
CE
PT
ED
fluctuating velocity in the recirculation zone is statistically similar to isotropic turbulence.
1
ACCEPTED MANUSCRIPT
Detailed Numerical Simulation of Swirling Primary Atomization Using
AN US
Method
CR IP T
a Mass Conservative Level Set
By
Changxiao Shao1, Kun Luo*, Yue Yang2, Jianren Fan1 1
M
State Key Laboratory of Clean Energy Utilization
Zhejiang University, Hangzhou 310027, P. R. China State Key Laboratory of Turbulence and Complex Systems
ED
2
AC
CE
PT
Peking University, Beijing 100871, P. R. China
Submitted to
International Journal of Multiphase Flow
*Corresponding author, E-mail: [email protected], Tel: 86-0571-87951764 2
ACCEPTED MANUSCRIPT
Abstract: We report detailed numerical simulations of swirling liquid atomization. A recently developed mass conservative level set method is employed to capture the gas-liquid interface and a ghost fluid method is utilized to deal with the jump conditions across the interface. The swirl and atomization characteristics of two-phase
CR IP T
annular swirling jets with the influence of turbulent inflow are investigated. Through comparing the sheet thickness, the breakup length and the cone angle, the numerical convergence of the global characteristics of the swirling two phase flow has been
AN US
obtained. The numerical results show that turbulent inflow can induce liquid sheet breakup near the nozzle exit, reduce the stiffness of the liquid sheet, and lead to the statistically homogeneous distribution of small-scale liquid structures in the radial
M
direction. Compared with the single-phase jet, the two-phase jet exhibits the chaotic velocity filed downstream that can enhance the mixing of droplets and ambient gas,
ED
and the precessing vortex core (PVC) is not observed in the center of the two-phase
PT
jet. In addition, the recirculation zone is smaller and farther from the nozzle exit for the turbulent inflow case than that from the laminar inflow case, and the preferential
i
CE
alignment of
with the intermediate strain rate indicates that the fluctuating
AC
velocity in the recirculation zone is statistically similar to isotropic turbulence. The interaction of the liquid-gas interface and vortices shows the preferential normal alignment of the vorticity and the normal of the interface, and the liquid sheet can generate high shear layers to produce anisotropic small-scale fluctuations. Keywords: primary atomization; swirl jet; level set method; turbulence; interface
3
ACCEPTED MANUSCRIPT
1 Introduction Liquid atomization is a common phenomenon in a variety of scientific and engineering applications, such as mixing, spraying, printing, food processing, agriculture, pharmaceutical process and combustion devices in gas turbine. The
CR IP T
atomization process is complex, involving highly turbulent and convoluted interfaces as well as breakup and coalescence of liquid masses, so that its mechanism has not yet been well understood [1]. It is of importance to explore the complex details of the
AN US
liquid atomization.
The atomization process of liquid can be divided into primary atomization and secondary atomization. Detailed experimental analysis of primary atomization remains a challenge due to the difficulty posed by dense spray measurements [2].
M
Alternatively, numerical simulations with interface capturing methods have been
ED
developed to study the primary atomization. The major interface capturing methods include the volume of fluid (VOF) method and the level set (LS) method. The VOF
PT
method is an Euler-type method with good mass conservation, but the implementation
CE
of this method is complicated owing to the reconstruction of the interface [3].
AC
Meanwhile, the LS method can easily describe topological changes under complex motions [4] but most of the existing LS methods are lack of the mass conservation property.
Both VOF and LS methods have been utilized to investigate primary atomization in straight or cross jet flows, such as in Kim et al. [5]-[6], Herrmann et al. [7], Xiao et al. [16], Shinjo and Umemura [18], Desjardins et al. [21], Ménard et al. [25]-Error! 4
ACCEPTED MANUSCRIPT
Reference source not found., Li et al. Error! Reference source not found.-Error! Reference source not found., Fuster et al. [30], M. Arienti and M. Sussman [30], Chen et al. [31].
However, there are very few of studies on primary atomization in swirling jet flow except Li et al. [28] that investigated a realistic swirling flow and qualitatively
CR IP T
demonstrated the validity and feasibility of the Coupled LS and VOF approach with experimental results. The swirling flow is significantly important when it is employed in gas turbine combustors as an approach to enhance fuel-air mixing and flame
AN US
stabilization. As demonstrated by Siamas et al. [46], many of the existing studies on annular swirling liquid jets focused on experimental visualizations and simplified mathematical models, which are difficult or insufficient to reveal the complex details
M
of liquid breakup and atomization in a two-phase environment.
For single-phase swirling flows, the swirling characteristics [42], such as vortex
ED
breakdown [47] and the precessing vortex core (PVC) [48], are extensively
PT
investigated. The formation of the recirculation zone, a form of vortex breakdown, can serve to stabilize flame. Many theories of vortex breakdown have been proposed,
CE
e.g., the generation of negative azimuthal vorticity in Brown et al. [54], the axial
AC
pressure exceeding the streamwise momentum flux in Mahesh [55] and inertial wave propagation in Squire [56]. The PVC is known as the rotating vortical structures in the center of the swirling jet, and its presence is significant for the instability and the combustion system. The observed PVC is generally considered as a secondary instability related to the post vortex breakdown [48] and can be suppressed in the presence of density gradient. On the other hand, the corresponding phenomena or the 5
ACCEPTED MANUSCRIPT
results of swirling characteristics can be different in a two-phase annular swirling jet, because the liquid sheet may significantly alter the swirling characteristics. In order to investigate the primary atomization in a swirl combustor, the atomization characteristics of an annular swirling jet are simulated using the
CR IP T
high-fidelity detailed numerical simulation with a recently developed mass conservative level set method. In addition, the influence of inflow turbulence on the global atomization characteristics is discussed. The outline of the present study is
AN US
organized as follows. In section 2, we provide the governing equations and describe the level set method. In section 3, we present numerical configurations, resolution considerations and the inflow conditions. In section 4, we discuss the atomization
ED
vorticity and interface.
M
characteristics, the influence of inflow turbulence and the correlation between
2 Numerical Methods
PT
2.1 Governing equations
CE
The Navier-Stokes equations for gas and liquid phases read
u u u p u ut g , t u 0 , t
AC
(1) (2)
where u is the velocity, ρ is the density, p is the pressure, g is the gravitational acceleration and µ is the dynamic viscosity. The material properties of gas and liquid are constant, i.e., ρ = ρl, µ = µl in liquid phase and ρ = ρg, µ = µg in gas phase and they are subjected to jump conditions at the interface, namely, l g and 6
ACCEPTED MANUSCRIPT
l g . The velocity at the interface is continuous, i.e., u ul ug 0 . The pressure across the interface can be expressed as
p 2 nt u n
,
(3)
where σ is the surface tension, κ is the interface curvature and n is the interface
CR IP T
normal. 2.2 Interface capturing method
The level set method is used to capture the interface, which is implicitly given by
AN US
the iso-surface of the smooth level set function G x, t 0 . Generally, the level set function G is imposed to be the signed distance function to the interface, i.e., G x, t x x ,
(4)
M
where x is the location at the interface that is closest to x. The level set function is defined to be positive for the liquid phase and negative for the gas phase.
ED
However, the level set method suffers the issue of mass loss. Instead of the
PT
signed distance function, the hyperbolic tangent function proposed by Olsson and Kreiss [34][35] is employed here, G x, t 1 , 2
x, t 1 tanh
CE
2
is the thickness of the profile. The evolution of the interface is implicitly
AC
where
(5)
captured by the level set equation,
u 0 . t
(6)
The level set equation is solved by the fifth-order upstream central scheme for spatial discretization and the second-order semi-implicit Crank-Nicolson scheme for time integration. It is noted that will no longer remain the hyperbolic tangent 7
ACCEPTED MANUSCRIPT
profile even at the well-resolved velocity, which can contaminate the interface and lead to mass loss. Consequently, an additional re-initialization equation [35] needs to be introduced to overcome this drawback by solving the following equation to a steady state,
where
(7)
CR IP T
1 n n n ,
is pseudo time. Although, this re-initialization equation can significantly
reduce the mass loss compared with the original re-initialization equation, this
AN US
approach still suffers from mass loss in under-resolved regions. Therefore, we use a mass conserving LS method [32] to further improve the mass conservation of the liquid phase by correcting the volume fraction based on the local curvature at the
M
interface. More information about the mass remedy procedure can be found in Luo et al. [32] in detail.
ED
2.3 Ghost fluid method for the jump conditions
PT
The jump conditions for pressure gradient in the Poisson equation as well as jump condition for the viscous terms are taken into account by the ghost fluid method.
CE
The details of this method can be found in Desjardins et al. [36]. Consider an
AC
interface located at x between the two grid locations xi and xi+1, where xi+1 is within the liquid phase. The pressure jump in the pressure Poisson equation is then written as
1 p p 1 p p l ,i 1 g ,i g ,i 1 * g g ,i p p 1 * 2 , 2 x x g ,i x x * where g l 1 and x xi / x .
8
(8)
ACCEPTED MANUSCRIPT
The pressure jump in Eq. (3) only includes the surface tension force. The continuum surface force (CSF) method is used for the viscous term and the spatial discretization is the second-order finite difference scheme. 2.4 Numerical schemes
CR IP T
The discretization of the Navier-Stokes equations is based on the staggered uniform grid, in which the pressure p and the LS function G are stored at the cell centers, while velocity is stored at the face centers. The spatial discretization of the
AN US
Navier-Stokes equations is performed using the second-order finite central difference schemes. The second-order semi-implicit iterative procedure [33] for time integration is utilized, which is economical, stable and accurate. The iteration can be expressed as (9)
M
f ukn11 u n tf 1 u n ukn 1 1 t ukn11 ukn 1 , 2 2 u
ED
where f is the right hand side of Navier-Stokes equations, and f / u is the Jacobian. The low and upper indices k and
n denote the kth Newton-Raphson sub-iterative
PT
step and the time step, respectively.
CE
The time step t at time t n is determined by the CFL conditions
AC
CFL dy CFL dz CFL dx 2 CFL dy 2 CFL dz 2 t min CFL dx , , , , , , (10) u v w 4 4 4 where CFL is the CFL number that is set to 1. The density and viscosity fields are computed as
g l g ,
(11)
g l g .
(12)
A projection method is used for temporal integration of the Navier-Stokes 9
ACCEPTED MANUSCRIPT
equations. The details are listed as below. (1) Compute an intermediate velocity field u* :
un u* un t un un
(13)
CR IP T
Note that the viscosity is computed following Eq. (7) when the interface crosses a mesh cell.
(2) Compute the pressure field by solving the Poisson equation with GFM: * p n 1 n 1 u . t
AN US
(14)
(3) Compute the real divergence-free velocity field by using the pressure gradient:
u t *
p n 1
n1
.
ED
3 Numerical Simulations
(15)
M
u
n 1
3.1 Flow configuration and computational cases
PT
To mimic fuel atomization in gas turbine combustor, we simulate the swirling jet
CE
injected into surrounding still air using the mass conserving LS method. Fig. 1 shows a schematic of the flow configuration. The exit shape of the nozzle is annular and the
AC
diameters of the outer and inner circle are Dout = 0.4 mm and Din = 0.2 mm, respectively. The initial condition of the liquid phase is a semi-sphere membrane. The liquid is swirling at the exit of the nozzle and the swirl number is defined as S U / U x , where U and U x are the azimuthal velocity and axial velocity,
respectively. The bulk azimuthal velocity and axial velocity are both set to 10 m/s, hence the swirl number is S = 1. The computational domain is 10Dout 10Dout 10Dout. 10
ACCEPTED MANUSCRIPT
The jet direction aligns with the x-direction. The inflow condition at the annular nozzle is velocity inlet and the outflow condition is convective outlet. The boundary conditions at the y and z directions are periodic conditions. Three important forces, including the inertia force l L2V 2 , the surface tension
CR IP T
force L , and the viscous force l LV act on the liquid in atomization (excluding gravity force), where L and V are the characteristic length and velocity, respectively. From these forces, two independent non-dimensional groupings are the Reynolds
AN US
number Re l LV / l and the Weber number We l LV 2 / . Another important number is the Ohnesorge number Oh We / Re l / l L .
The major parameters are listed in Table 1. The characteristic length and velocity
M
are set to be the outer diameter and axial velocity, respectively, resulting in the Re, We and Oh of 2000, 222 and 0.015. The density ratio, viscosity ratio, and the swirl
ED
number are chosen to mimic the realistic conditions in gas turbine combustor, but the
PT
Reynolds number and the Weber number are limited to 103 and 102, respectively, due to our available computer resources.
CE
The computational domain is a rectangular region. The Cartesian grid system is
AC
used for all cases. The cases are listed in Table 2. Cases 1 and 3 are performed to investigate the grid resolution effect on the results. Cases 2 and 3 are performed to investigate the effect of turbulence on the swirling atomization. Case 3 is also conducted to investigate the two-phase characteristics compared to the corresponding single-phase jet. The computation is carried out on the National Supercomputer in Tianjin. Cases 1 and 2 used 448 cores for about 50 hours, and Case 3 used 2048 cores 11
ACCEPTED MANUSCRIPT
for about 160 hours. 3.2 Resolution considerations In order to ensure the adequate resolution of the turbulent atomization, the smallest turbulent scale and liquid structures need to be resolved. The smallest
CR IP T
turbulent scale that needs to be resolved is the Kolmogorov length scale. The liquid structures in the present simulation include thin liquid sheet, ligaments and droplets, and the phase interface thickness. These two issues lead to a very strict constraint on
AN US
the mesh size. As demonstrated in [37][22], such resolution requirements make simulations by using the current state-of-the-art methods unfeasible in the context of the detailed numerical simulation. Therefore, we only focus on the global features of
computational resources.
M
turbulent swirling atomization with the highest resolution depending on our available
ED
Similar to that in [22], the resolution of turbulence is related to the number of
PT
grid points per jet diameter. For the coarse case, the grid resolution x is 7.8µm resulting in 50 mesh points per Dout. For the fine case, the grid resolution x is
CE
3.9µm resulting in 100 mesh points per Dout. The resolution of the liquid structures
AC
follows the consideration of Ménard et al. [25]. The grid-based Weber number is defined as Wex lU 2 x / and it is suggested that no further breakup occurs if
Wex is smaller than about 10. In the present simulations, the grid Weber number Wex is 17.36 for the coarse case and 8.68 for the fine case.
3.3 Turbulent inflow condition Non-turbulent uniform profile inflow is first specified in the simulation for the 12
ACCEPTED MANUSCRIPT
liquid jet. This enables the simulation to exclude the effect of inflow turbulence on the swirling atomization. Following [15]-[17], the influence of turbulent inflow condition is then investigated. In order to generate the realistic turbulent inflow condition, a precursor simulation of a periodic swirling pipe flow is performed using the liquid
CR IP T
properties. The cylindrical mesh is used and the computational domain is [0, 10D], [0.25D, 0.5D], [0, 2 ] in x, r, directions, respectively. The velocities in the azimuthal and axial directions are both set to 10 m/s. The domain is discretized on a
AN US
192 64 48 mesh. Fig. 2 shows the contours of instantaneous velocity magnitude with random fluctuations in the longitudinal direction of the pipe flow and the coupling with the external flow. Once the flow reaches the statistically stationary state, and L are
M
the velocity field is recorded within the time period U / L , where U
streamwise velocity magnitude and 10D, respectively. The Cartesian velocities in the
ED
x, y, z directions then read
PT
U U x , V U cos arctan
y , W U sin arctan z
y . z
(16)
CE
For the laminar inflow case, Eq. (16) is directly applied to the velocity boundary condition. For the turbulent inflow case, a swirling turbulent pipe flow simulation is
AC
first conducted and then stored as the time-dependent inlet velocity of the external flow in the combustor for swirling atomization. Since the scale of the pipe is equal to the outer diameter of the inflow of the combustor, the stored pipe inflow date file can be directly read for the external flow calculations. Moreover, the mean velocity and root mean square (rms) profiles obtained at the nozzle exit are illustrated in Fig. 3. The mean velocity increases sharply from the wall. 13
ACCEPTED MANUSCRIPT
The turbulent intensity is about 20% of the mean velocity and the rms velocity reaches a maximum value near the wall. Although it is not validated by the experimental results, the present velocity and rms profiles are similar to that in [15]-[17], so we will investigate the effect of turbulent inflow on the atomization by
CR IP T
utilizing this inflow condition.
4 Results and Discussions 4.1 Global atomization characteristics
AN US
In the later stage of the sheet breakup for the laminar case, some finger-like structures appear at the rim of the sheet, which can be seen in the simulation results in Fig. 4(a) with the coarse mesh and in Fig. 4(b) with the fine mesh. The experimental
M
result in Taylor [38] is also shown in Fig. 4(c) for comparison. In [38], a liquid sheet of kerosene issuing from the orifice with the diameter 0.3 mm of a swirl atomizer
ED
under pressure of 1.36 atm into air at pressure 0.32 atm. The visualizations of the
PT
simulation and experimental results are qualitative similar despite of the different conditions. However, the finger-like structures in the numerical simulation results are
CE
generally longer than those in the experiment. This difference seems to be attributed
AC
to the liquid property or the insufficient grid resolution that leads to unphysical breakup. It can be seen from Figs. 4(a) and 4(b) that ligament length is significantly reduced in the case by refining grid resolution, which implies that the grid resolution has an impact on the breakup of swirling sheet in the numerical results. On the other hand, continuously increasing the grid resolution seems unfeasible considering that the sheet thickness may be down to tens of nanometers. Therefore, the present study 14
ACCEPTED MANUSCRIPT
focuses on the large-scale characteristics of spray formation and the effects of grid resolution on the overall characteristics, including the sheet thickness, breakup length, recirculation zone and spray cone angle. 4.1.1 Liquid sheet thickness
CR IP T
The thickness of the liquid sheet at the nozzle exit has an impact on the drop sizes in the primary atomization. Taylor [38] speculated that the spiral variation of the sheet thickness at the nozzle exit and the oscillating air core could lead to the
AN US
variation of the liquid thickness in the conical sheet. The multi-dimensional spray modeling [39] predicts that the sheet thickness at the breakup length is related to the cone angle, the maximum growth rate and the liquid velocity.
M
Fig. 5 shows the overall interface spatial and temporal evolution of the laminar swirling atomization in Case 1. Here, the dimensionless time T = tU/D is used where
ED
U is the liquid axial velocity. At the first stage of the injection, the liquid structure
PT
forms the mushroom-shaped tip at T = 1.25 in Fig.5 due to the liquid penetration by impingement against the still air, which is similar to that shown in Shinjo et al. [18].
CE
Some ligaments are first generated at this rolled-up structures and then breaks into
AC
droplets at the outer recirculation zone (ORZ) surrounding the conical sheet and at T = 2.50 in Fig.5. The liquid film spreads downstream smoothly until T = 3.75. In the present simulation, the sheet variation shows a limited dependence on the
radial disturbances for the laminar inflow condition, which can be evident in the smooth interface sheet in Fig. 5. The thickness h of the film can be calculated from the mass flow rate equation [1] 15
ACCEPTED MANUSCRIPT
m lVave h 2 x tan h ,
(17)
where Vave is the average axial velocity along the liquid sheet, x and are the stream-wise position and the half cone angle, respectively. Eq. (17) is based on the assumption that liquid sheet does not rupture before the breakup length. Although the
CR IP T
simulation result shows a few of droplets at the ORZ pinched off from the tip of the conical sheet in Fig. 5, it is still reasonable to assume the corresponding volume is much smaller than that of the liquid sheet.
AN US
The spatial variations of the sheet thickness for laminar inflow conditions at different grid resolutions at T = 12.5 are shown in Fig. 6. The thickness of the sheet decreases sharply near the nozzle exit with the increasing stream-wise position and
M
then decreases slowly. Before the location x/Dout = 1.5, the thicknesses of the sheet for different grid resolutions are collapsed to each other. It is noted that the thickness of
x tan
PT
h x tan
ED
the sheet for the fine case is in accordance with the physically reasonable solution 2
m lVave
(18) in the
CE
of Eq. (17), where m 7.536 104 kg / m3 , Vave 8.15m / s and 45
present simulation. The difference mainly locates at the film breakup zone, comparing
AC
the coarse and fine meshes of the laminar cases. For the numerical simulations, the liquid sheet may occur non-physical breakup at the scale of the grid size. Although 10243 grids have been adopted in the present simulation, there are about 2.5 cells across this thickness where the sheet begins to breakup, which still doesn’t guarantee the sufficient resolution of the sheet breakup. However, the results prove to capture the major characteristics of the thickness of the sheet. Moreover, the liquid sheet 16
ACCEPTED MANUSCRIPT
breaks at about x/Dout = 2 for the coarse mesh case as shown in Fig. 6. The sheet thickness at this position is h 2.5x , which should be carefully investigated. For the fine mesh case, the sheet thickness at the breakup position is about h 5x . 4.1.2 Breakup length and recirculation zone
CR IP T
The breakup length of the liquid sheet is defined as the distance from the nozzle exit to the position where the sheet begins to break. The recirculation zone is defined as a flow transition with a free stagnation point/region on the axis followed by a
AN US
reverse flow. These two features can be observed in Fig. 7, which shows the overlaps of the interface among different plane-cuts for the coarse mesh (a) and the fine mesh (b) at the statistically stationary state T = 12.5. The plane-cuts along the azimuthal
M
direction are projected on the x-y plane in Fig. 7. It illustrates that the breakup of the liquid length occurs around x/Dout = 2~3 for Case 1 with the coarse mesh and x/Dout =
ED
3~4 for Case 3 with the fine mesh. The recirculation zones are formed around x/Dout =
PT
2~5 for both cases. Therefore, the grid resolution has an impact on the breakup length, but refining grid resolution does not change the general swirling-flow characteristics
CE
significantly.
AC
4.1.3 Spray cone angle The spray cone angle is an important spray characteristic for evaluating atomizer
performance. The cone angle is defined as the angle between the tangents to the envelope of the spray. Qualitatively, the spray cone angle increases with time and reaches to a certain value. In the first stage, the axisymmetric vortex pair is generated, entraining the surrounding gas to the swirl core, and bend the liquid film towards the 17
ACCEPTED MANUSCRIPT
axis. In the later stage, the vortex breaks down and the spray cone angle increases. Quantitatively, the temporal evolutions of the spray angle with different grid resolutions are shown in Fig. 8. It shows that the evolutions of the spray angle are almost the same for different grid resolutions.
CR IP T
In summary, the grid resolution has effects on the breakup zone of the liquid sheet, but both Case 1 and Case 3 can capture the major characteristics of the swirling flow.
AN US
4.2 Turbulent Inflow Effect 4.2.1 Flow structures
We compare qualitative differences of the swirl atomization between laminar
M
inflow and turbulent inflow in Figs. 5 and 9. In the turbulent inflow condition, the helical wave (see Fig. 9, T = 1.25), which is not observed in the laminar inflow
ED
condition, appears due to the inlet turbulence disturbance. The helical wave and
PT
turbulent liquid inflow make the conical sheet significantly wrinkled. The conical sheet develops for a short distance and then immediately ruptures into ligaments and
CE
droplets (see Fig. 9, T = 3.75 – 6.25). Note that the breakup length of the conical sheet
AC
is much shorter than that in the laminar inflow condition. Moreover, the mean length of the ligaments is significantly reduced due to the turbulent straining motion that enhances the rupture of the ligaments. From Figs. 5 and 9, we conclude that turbulent inflow enhances the flow instability, which results in the droplets disperse uniformly in the radial direction after atomization. In order to further investigate the turbulent inflow effect, we overlap different 18
ACCEPTED MANUSCRIPT
plane-cuts of the interface with the turbulent inflow condition at T = 12.5 in Fig. 10. At x/Dout = 0.5-1.5 near the nozzle exit, the liquid sheets do not overlap with each other, which implies that the sheet thickness fluctuates along the azimuthal direction and it is different from the laminar inflow case. In the turbulent inflow case, the
CR IP T
small-scale liquid structures after breakup are gathered around x/Dout = 2 and y/Dout = 3. While for the laminar inflow case, those structures are gathered in the downstream of x/Dout = 4. This shows that the turbulent inflow can promote the liquid sheet
AN US
breakup near the nozzle exit and enhance the radial movement of the liquid structures. Additionally, the instantaneous interfaces at T = 12.5 for different inflow conditions viewed from upstream are shown in Fig. 11. The liquid sheet atomization
M
in the laminar inflow condition displays a two-step cascade, including the process from the sheet to ligaments and from ligaments to droplets, which is different from
ED
that in [40][41]. The ligaments hardly occur in the turbulent inflow condition in
PT
Fig.11(b). Therefore, the turbulence inflow significantly enhances the liquid atomization, which agrees with the result in the crossflow simulation [15].
CE
The different inflow conditions lead to different characteristics of liquid
AC
dispersion. In the laminar inflow case, the liquid structures experience the two-step cascade from the sheet through ligaments to droplets. The liquid dispersion for the turbulent inflow case is quite different. The liquid sheet starts to undulate along the azimuthal direction in the upstream, and then liquid drops formed undergo the radial dispersion, both toward and outward the centerline. The liquid structures in the downstream shows approximately a statistically uniform distribution in the radial 19
ACCEPTED MANUSCRIPT
direction. It is demonstrated that the turbulence induced in the liquid inflow can weaken the stiffness of the sheet and leads to more even distribution of the liquid structures in the radial direction. 4.2.2 Recirculation characteristics
CR IP T
In order to understand the structure characteristics of the gas-liquid swirling recirculation zone generated by vortex breakdown, the mean velocity profiles on x-y plane-cuts in Fig. 12. Although there are extensive studies on vortex breakdown in
AN US
single phase flow [42][43][44], there are few studies on vortex breakdown in gas-liquid swirling flow. In [42] for the single phase flow, vortex breakdown has been defined as an abrupt flow transition with a free stagnation point/region on the axis
M
followed by a reverse flow and a fully turbulent region. Fig. 12 shows the instantaneous axial velocity field of the swirling gas-liquid flow in the laminar inflow
ED
and in the turbulent inflow condition at T = 15.0. Compared to the single phase swirl
PT
flow, the gas-liquid swirl flow has intense fluctuations. Two distinct features can be observed in Fig. 12. First, the distance of the
CE
recirculation zone location from the nozzle is longer for the laminar inflow condition.
AC
From Fig. 12(a), the axial velocity is approximately zero when x/Dout = 1.5. The central recirculation zone occurs at x/Dout = 1.5 - 4.5 for laminar inflow, while the central recirculation zone is at x/Dout = 1.0 - 3.0 for turbulent inflow. It is interesting to note that the maximum negative axial velocity occurs at the liquid sheet breakup into filaments and droplets, x/Dout = 2.5 for laminar inflow and x/Dout = 1.5 for turbulent inflow. Second, the size of the breakdown in the laminar inflow case is nearly twice as 20
ACCEPTED MANUSCRIPT
large as that in the turbulent inflow case. The possible reason for the differences is that the liquid sheet can maintain its momentum in the radial and axial directions in the laminar inflow case, so the momentum exchange from liquid phase to gas phase is in further downstream. A primary axisymmetric vortex pair can be observed at the
CR IP T
location where liquid sheet breaks into droplets for the laminar inflow condition. However, the recirculation zone is suppressed in the turbulent inflow case due to the turbulent sheet breakup at a short distance from the nozzle and the subsequent vortex
AN US
breakdown.
Fig. 13 shows the axial and radial velocity profiles at different axial positions for laminar and turbulent inflows at T = 15.00. In Fig. 13(a), the axial velocity decreases
M
in the axial direction between the inlet and the stagnation point in both laminar and turbulent inflow cases owing to the adverse pressure gradient [45]. Moreover, the
ED
radial distance of the recirculation zone in the turbulent inflow case is shorter than
PT
that in the laminar inflow case, which agrees with the observation in Fig. 12. 4.3 Characteristics of vortices and interfaces
CE
4.3.1 Precessing vortex core (PVC)
AC
The precessing vortex core (PVC) is known as the rotating vortical structures in
the center of the swirling jet. Its presence is significant for the instability and the combustion system. Previous studies found that PVC occurs when the swirl number S > 0.6-0.7 [49] in single phase flow. However, it is an open question if the PVC exists in two-phase annular swirling jets. To identity PVC, we display the velocity vectors for the single-phase jet, the laminar and turbulent inflow cases of the two-phase jet on the 21
ACCEPTED MANUSCRIPT
plane-cuts at x/Dout = 3. For the single phase case, the PVC is evident in the center of the jet in Fig. 14(a) with the occurrence of the maximum axial velocity spinning around the origin, which agrees with the simulation result in [46] and the experimental observations in [50]. On the other hand, there is no obvious PVC in both
CR IP T
two-phase jet cases with laminar and turbulent inflow conditions, which was also observed in [46] where the swirl number is 0.4. It appears that the two-phase cases have strong radial momentum to weaken the vortical flow in the center of the jet.
AN US
In the two-phase cases, the complex vortical structures appear in the downstream. For the laminar inflow case, the clockwise regular velocity vectors occur at near nozzle exit (not shown), where the liquid sheet has not broken up. In the downstream,
M
the velocity vectors along the large circular rim are oriented to the positive radial direction, which indicates the liquid drops disperse to the outward of the sheet. In the
ED
inner zone of the liquid sheet, the velocity field is chaotic with numerous vortical
PT
structures in Fig. 14(b), which can enhance mixing of the drops and the ambient gas. For the turbulent inflow case, the velocity vectors at the downstream in Fig. 14(c) are
CE
more chaotic that those in the laminar inflow case.
AC
In contrast, the single phase case exhibits very regular flow structures. The
clockwise rotation of the velocity vectors appears at the upstream and the symmetric triple-vortices with the center are formed in Fig. 14(a). This can be attributed to the helical wave development along the swirling jet. In Fig. 14(a) the radial zone induced by swirling flow at z/Dout = [-2, 2], so the radial development in the single phase case is compact compared to the two-phase cases. With the same inlet velocity, the liquid 22
ACCEPTED MANUSCRIPT
phase has much larger momentum than that in the single gas-phase case to penetrate to streamwise and radial directions, leading to the wide spread of liquid structures. 4.3.2 Correlation of the interface and vortices The gas-liquid interfaces in Case 3 at T = 6.25 are shown in Fig. 15, along with
CR IP T
the vortical structures characterized as the iso-surfaces of the Q-criterion. Here, Q Aij Aij Sij Sij / 2 is the second invariant of the velocity gradient tensor, where Aij and Sij are the asymmetric and symmetric parts of the velocity gradient tensor,
AN US
respectively. It can be seen that high vorticities concentrate on the breakup region. The preferred normal alignment of the LS scalar gradient with the vorticity field near the interface can be quantified by the cosine of the angle between vorticity vector and
,
M
LS scalar gradient vector as
(19)
ED
where verticity vector u . Fig. 16 shows the PDF of at T = 6.25 that has
PT
a strong peak at 0 . Here, 0 means strong perpendicularity and 1
CE
no perpendicularity. The strong perpendicularity of the interface and the vortices is also in accord with the result in [53] for homogeneous isotropic turbulence.
AC
In order to examine the dynamics in small-scale turbulence, the enstrophy
evolution equation for two phase flow is utilized as
p 2 S D 1 2 , S Dt 2
(2 0)
t where D / Dt / t u , S u u / 2 , and the terms on the right hand side
23
ACCEPTED MANUSCRIPT
are the enstrophy production due to vortex stretching, the baroclinic term, and viscous dissipation term, respectively. The LS function G can be seen as a passive scalar, and the small-scale passive scalar structure can be analyzed in terms of the evolution equation for the scalar
CR IP T
gradient as D 1 2 Gi Gi Sij G j , Dt 2
(21)
where Gi G / xi , and the term on the right-hand side is the production.
AN US
The turbulent production terms for ii and Gi Gi are i Sij j and Gi Sij G j , respectively. The alignments among the vorticity i , the scalar gradient Gi and the principal rate of Sij are significant to determine the turbulent production terms
M
[57][58]. Fig. 17 shows the PDFs for the cosine of the angle between the vorticity and
ED
the principal strain rate directions in an area of the recirculation zone (encircled by the solid lines in Fig. 7b and Fig. 12a), where ai, bi, ci are the eigenvectors corresponding
are
real
and
CE
c
PT
to the directions of the eigenvalues a, b, c, respectively. Since Sij is symmetric, a, b, satisfy
abc 0
(solenoidal
condition)
and
a b c (a 0 c) . Here, cos 1 means strong alignment and cos 0
AC
no alignment. It shows in Fig.17 that
i is aligned with the intermediate (b) strain
rate, which is consistent with the findings of Ashurst et al. (1987) for isotropic turbulence and Abe et al. (2009) of far-wall region for turbulent channel flow. In the recirculation zone, the turbulent eddies are generated by the wakes formed due to droplets dynamics. The gas phase flow in the recirculation region is modulated and exhibits similar characteristics to the isotropic turbulence as shown in Fig. 17. This 24
ACCEPTED MANUSCRIPT
indicates that the present recirculation zone of swirling two-phase jet also contain the similar small-scale structures as those in isotropic turbulence. In the near-sheet regions (encircled by the dashed lines in Fig. 7b and Fig. 12a),
i is also found to be aligned with the intermediate (b) strain rate as shown in Fig.
CR IP T
18, which agrees with the findings of Abe et al. (2009) for the near-wall region for turbulent channel flow and Shinjo et al. (2015) [59] for the near-droplet regions for dense spray. Moreover, the alignments of LS gradient Gi with the extensional (a)
AN US
and the compressive (c) strain rates (Fig. 19) tend to an orientation angle of 45° ( cos 0.7 ). Both trends of the alignments among
i , Gi and the principal strain
rate are qualitatively similar to those in the near-wall zone for the turbulent channel
M
flow[58], which indicates that the liquid sheet can generate high shear layers to
5 Conclusions
ED
produce anisotropic small-scale fluctuations.
PT
We report the spray characteristics of the swirling liquid primary atomization. The primary atomization is simulated using the detailed numerical simulation with the
CE
recently developed mass conservative LS method and the ghost fluid method handling
AC
jump conditions. The interface is represented by an iso-surface of a hyperbolic tangent function which is kept with a re-initialization process. Through comparing the sheet thickness, the breakup length and the cone angle, we demonstrate that the grid resolution is satisfactory from the convergence of the major characteristics of the swirling two-phase flow. The liquid sheet thickness for the fine case agrees well with theoretical analysis. 25
ACCEPTED MANUSCRIPT
Compared with the single phase jet, the two-phase jet exhibits the chaotic velocity field downstream. The momentum carried by the liquid-phase fluctuates intensely, which can suppress the generation of PVC in the center of the jet. Compared to the laminar inflow case, the turbulent inflow reduce the stiffness of the
CR IP T
liquid sheet leading to statistically homogeneous distribution of the liquid structures in the radial direction. In addition, the recirculation zone is smaller and farther from the nozzle exit for the turbulent inflow case. The interaction of the liquid-gas interface
AN US
and vortices is characterized by the preferential normal alignment between the gradient of the LS scalar and vorticity. The alignments among the vorticity, the level-set scalar gradient and the principal rate of the symmetric part of the velocity
M
gradient tensor agree with the experimental results, and indicate that small-scale dynamics in the recirculation zone is similar to that in isotropic turbulence and is
ED
shows anisotropic statistics near the liquid sheet.
PT
Acknowledgement
CE
This work is financially supported by the National Natural Science Foundation of China (Nos. 91541202, 51222602, and 51276163). Yue Yang acknowledges the
AC
support by the National Natural Science Foundation of China (Nos. 11522215 and 91541204) and the Young Thousand Talent Program of China.
References [1] N. Ashgriz, Handbook of atomization and sprays: Theory and Applications, (2011) [2] X. Li, M. C. Soteriou, W. Kim, J. M. Cohen, High fidelity simulation of the 26
ACCEPTED MANUSCRIPT
spray generated by a realistic swirling flow injector, Journal of Engineering for Gas Turbines and Power, 136(2014)071503-1-10 [3] G. Tryggvason, R. Scardovelli, S. Zaleski, Direct numerical simulations of gas-liquid multiphase flows, Cambridge University Press, 2011
CR IP T
[4] S. Osher, J. A. Sethian, Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulation. J. Comp. Phys., 79(1988)12 [5] D. Kim, O. Desjardins, M. Herrmann, P. Moin. Toward two-phase simulation of
AN US
the primary breakup of a round liquid jet by a coaxial flow of gas, Center for Turbulence Research, Annual Research Briefs, 2006
[6] D. Kim, P. Moin, Numerical simulation of the breakup of a round liquid jet by a
M
coaxial flow of gas with a subgrid Lagrangian breakup model, Center for Turbulence Research, Annual Research Briefs, 2011
ED
[7] M. Herrmann. A balanced force refined level set grid method for two-phase
PT
flows on unstructured flow solver grids. J. Comput. Phys., 227,2674-2706(2008) [8] M. Herrmann. On simulating primary atomization using the refined level set grid
CE
method. Atom. Sprays, 21(2011)283-301
AC
[9] M. Herrmann, M. GoroK-Hovski, An outline of an LES subgrid model for liquid/gas phase interface dynamics, Center for Turbulence Research, Proceedings of the Summer Program, 2008
[10] M. Herrmann. Detailed numerical simulations of the primary atomization of a turbulent liquid jet in crossflow. Journal of Engineering for Gas Turbines and Power, 132(2009)061506-1-10 27
ACCEPTED MANUSCRIPT
[11] M. Herrmann. Detailed simulations of the breakup processes of turbulent liquid jets in subsonic crossflows. In 11th Triennial International Annual Conference on Liquid Atomization and Spray Systems, ICLASS 2009 [12] M. Herrmann. A parallel Eulerian nterface tracking/Lagrangian point particle
CR IP T
multi-scale coupling procedure. J. Comput. Phys., 229(2010)745-759 [13] M. Herrmann, M. Arienti, M. Soteriou. The impact of density ratio on the liquid core dynamics of a turbulent liquid jet injected into a crossflow, Journal of
AN US
Engineering for Gas Turbines and Power, 133(2011)061501-1-9
[14] M. Herrmann. The influence of density ratio on the primary atomization of a turbulent liquid jet in crossflow. Proceedings of the Combustion Institute, 33
M
(2011) 2079-2088
[15] F. Xiao, M. Dianat, J. J. McGuirk, Large eddy simulation of liquid-jet primary
ED
breakup in air crossflow. AIAA Journal, 51(2013)2878-2893
PT
[16] F. Xiao, M. Dianat, J. J. McGuirk, Large eddy simulation of single droplet and liquid jet primary breakup using a coupled level set/volume of fluid method.
CE
Atomization and Sprays, 24(2014)281-302
AC
[17] F. Xiao, M. Dianat, J. J. McGuirk, LES of turbulent liquid jet primary breakup in turbulent coaxial air flow, International Journal of Multiphase Flow, 60(2014) 103-118
[18] J. Shinjo, A. Umemura, Simulation of liquid jet primary breakup: dynamics of ligament and droplet formation, International Journal of Multiphase Flow, 36 (2010)513-532 28
ACCEPTED MANUSCRIPT
[19] J. Shinjo, A. Umemura, Detailed simulation of primary atomization mechanisms in Diesel jet sprays (isolated identification of liquid jet tip effects), Proceedings of the Combustion Institute, 33(2011)2089-2097 [20] J. Shinjo, A. Umemura, Surface instability and primary atomization
CR IP T
characteristics of straight liquid jet sprays, International Journal of Multiphase Flow, 37(2011)1294-1304
[21] O. Desjardins, V. Moureau, H. Pitsch, An accurate conservative level set/ghost
AN US
fluid method for simulation turbulent atomization, J. Comput. Phys., 227(2008) 8395-8416
[22] O. Desjardins, H. Pitsch, A spectrally refined interface approach for simulating
M
multiphase flows, J. Comput. Phys., 228(2009)1658-1677
[23] O. Desjardins, H. Pitsch, Detailed numerical investigation of turbulent
ED
atomization of liquid jets, Atomization and sprays, 20(2010)311-336
PT
[24] O. Desjardins, J. O. McCaslin, M. Owkes, P. Brady, Direct numerical and large-eddy simulation of primary atomization in complex geometries,
CE
Atomization and sprays, 23(2013)1001-1048
AC
[25] T. Ménard, S. Tanguy, A. Berlemont, Coupling level set/VOF/ghost fluid methods: Validation and application to 3D simulation of the primary break-up of a liquid jet, International Journal of Multiphase Flow, 33(2007)510-524
[26] R. Lebas, T. Ménard, P. A. Beau, A. Berlemont, F. X. Demoulin, Numerical simulation of primary break-up and atomization: DNS and modeling study, International Journal of Multiphase Flow, 35(2009)247-260 29
ACCEPTED MANUSCRIPT
[27] X. Li, M. C. Soteriou, M. Arienti, M. M. Sussman, High-fidelity simulation of atomization and evaporation in a liquid jet in cross-flow, 49th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition, 2011
CR IP T
[28] X. Li and M. C. Soteriou, High-fidelity simulation of fuel atomization in a realistic swirling flow injector, Atomization and Sprays, 23(2013)1049-1078
[29] D. Fuster, A. Bagué, T. Boeck, L. L. Moyne, A. Leboissetier, S. Popinet, P. Ray,
AN US
R. Scardovelli, S. Zaleski, Simulation of primary atomization with an octree adaptive mesh refinement and VOF method, International Journal of Multiphase Flow, 35(2009)550-565
M
[30] M. Arienti and M. Sussman, An embedded level set method for sharp-interface multiphase simulations of diesel injectors, International Journal of Multiphase
ED
Flow, 59(2014)1-14
PT
[31] X. Chen, D. Ma, V. Yang, S. Popinet, High-fidelity simulations of impinging jet atomization, Atomization and Sprays, 23(2013)1079-1101
CE
[32] K. Luo, C. X. Shao, Y. Yang, J. R. Fan. A mass conserving level set method for
AC
direct numerical simulation of liquid atomization. J. Comput. Phys., 298 (2015) 495-519
[33] C. D. Pierce, P. Moin. Progress-variable approach for large eddy simulation of turbulent combustion. Technical Report TF80, Flow Physics and Computation Division, Dept. Mech. Eng., Standford University, 2001 [34] E. Olsson and G. Kreiss. A conservative level set method for two phase flow. J. 30
ACCEPTED MANUSCRIPT
Comput. Phys., 210,225-246(2005) [35] E. Olsson, G. Kreiss, and S. Zahedi. A conservative level set method for two phase flow ii. J. Comput. Phys., 225,785-807(2007) [36] O. Desjardins, V. Moureau, H. Pitsch, An accurate conservative level set/ghost
CR IP T
fluid method for simulating turbulent atomization, J. Comput. Phys., 227(2008) 8395-8416
[37] M. GoroK-Hovski and M. Herrmann, Modeling primary atomization, Annu. Rev.
AN US
Fluid Mech., 40(2008)343-366
[38] G. Taylor, The dynamics of thin sheets of fluid. III. Disintegration of fluid sheets, Proc. R. Soc. Lond. A, 253(1959)313-321
M
[39] P. K. Senecal, D. P. Schmidt, I. Nouar, C. J. Rutland, R. D. Reitz, M. L. Corradini, Modeling high-speed viscous liquid sheet atomization, International
ED
Journal of Multiphase Flow, 25(1999)1073-1097
PT
[40] Y. Wang, K. S. Im, K. Fezzaa, Similarity between the primary and secondary air-assisted liquid jet breakup mechanisms, Physical Review Letters, 100(2008)
CE
154502
AC
[41] E. A. Novikov and D. G. Dommermuth, Distribution of droplets in a turbulent spray, Physical Review E, 56(1997)5479-5482
[42] H. Liang and T. Maxworthy, An experimental investigation of swirling jets, J. Fluid Mech., 525(2005)115-159 [43] O. L. Negro and T. O’Doherty, Vortex breakdown: a review, Progress in Energy and Combustion Science, 27(2001)431-481 31
ACCEPTED MANUSCRIPT
[44] M. R. Ruith, P. Chen, E. Meiburg, T. Maxworthy, Three-dimensional vortex breakdown in swirling jets and wakes: direct numerical simulation, J. Fluid Mech., 486(2003)331-378 [45] E. Krause, The solution to the problem of vortex breakdown, In: Lecture notes
CR IP T
in physics, vol.371, Berlin: Springer, (1990)35-50 [46] G. A. Siamas, X. Jiang, L. Wrobel, Direct numerical simulation of the near-field dynamics of annular gas-liquid two-phase jets, Physics of Fluids, 21 (2009)
AN US
042103
[47] W. Kollmann, A. S. H. Ooi, M. S. Chong, J. Soria, Direct numerical simulations of vortex breakdown in swirling jets, J. Turbul., 2 (2001) N5
M
[48] K. Manoharan, S. Hansford, J. O’Connor, S. Hemchandra, Instability mechanism in a swirl flow combustor: Precession of vortex core and influence
ED
of density gradient, Proceedings of ASME Turbo Expo, 2015
PT
[49] A. K. Gupta, D. J. Lilley, N. Syred, Swirl flows, Tunbridge Wells, UK: Abacus Press, 1984
CE
[50] N. Syred, A review of oscillation mechanisms and the role the precessing vortex
AC
core (PVC) in swirl combustion systems, Prog. Energy Combust. Sci., 32 (2006) 93-161
[51] K. Luo, H. Pitsch, M. G. Pai, O. Desjardins, Direct numerical simulations and analysis of three-dimensional n-heptane spray flames in a model swirl combustor, Proceedings of the Combustion Institute, 33, 2143-2152 (2011) [52] P. Moin and S. V. Apte, Large-eddy simulation of realistic gas turbine 32
ACCEPTED MANUSCRIPT
combustors, AIAA Journal, 44, 698-708 (2006) [53] Y. Yang, D. I. Pullin, I. B. Moreno. Multi-scale geometric analysis of Lagrangian structures in isotropic turbulence. J. Fluid Mech.,654,233-270,(2010) [54] Brown, G. & Lopez, J. Axisymmetric vortex breakdown Part 2. Physical
CR IP T
mechanisms. J. Fluid Mech. 221, 553 (1990) [55] Mahesh, K. A model for the onset of breakdown in an axisymmetric compressible vortex. Phys. Fluids 8, 3338 (1996)
AN US
[56] Squire, H. B. Analysis of the vortex breakdown phenomenon. In Miszellaneen der Angewandten Mechanik, p. 306. Akademie-Verlag, Berlin (1962) [57] W. T. Ashurst, A. R. Kerstein, R. M. Kerr, C. H. Gibson, Alignment of vorticity
M
and scalar gradient with strain rate in simulated Navier-Stokes turbulence, Physics of Fluids, 30, 2343 (1987)
ED
[58] H. Abe, R. A. Antonia, H. Kawamura, Correlation between small-scale velocity
(2009)
PT
and scalar fluctuations in a turbulent channel flow, J. Fluid Mech., 627, 1-32
CE
[59] J. Shinjo, J. Xia, A. Umemura, Droplet/ligament modulation of local small-scale
AC
turbulence and scalar mixing in a dense fuel spray, Proceedings of the Combustion Institute, 35 (2015) 1595-1602
33
CR IP T
ACCEPTED MANUSCRIPT
PT
ED
M
AN US
Fig. 1: Schematic of flow configuration used for swirling film atomization.
AC
CE
Fig. 2: Contours of the instantaneous velocity magnitude inthe longitudinal directions of the pipe flow and the coupling with the external flow.
34
ED
M
AN US
Fig.3: Mean velocity (a) and rms (b) profiles at the liquid inlet.
CR IP T
ACCEPTED MANUSCRIPT
AC
CE
PT
Fig. 4. The breakup of the gas/liquid interface, (a) Simulation result from Case 1 with a coarse mesh; (b) Simulation result from Case 3 with a fine mesh; (c) Experimental result in Taylor [38]
Fig. 5: The interface evolution of the laminar swirling atomization for Case 3(T= 0, 1.25, 2.50, 3.75, 5.00, and 6.25 from left to right).
35
CR IP T
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
AN US
Fig. 6: The spatial variations of the sheet thickness for the laminar inflow condition at different grid resolutions at T= 12.5.
Fig. 7: The overlap of the interfaces on different plane-cuts for coarse mesh (a) and fine mesh (b) at T= 12.5
36
CR IP T
ACCEPTED MANUSCRIPT
M
AN US
Fig. 8: The temporal evolution of the spray angle for different grid resolutions.
AC
CE
PT
ED
Fig. 9: The interface evolution of the turbulent swirling atomization for Case 2(T= 0, 1.25, 2.50, 3.75, 5.00, and 6.25 from left to right).
Fig. 10: The overlap of the interfaces on different plane-cuts for the turbulent inflow condition at T= 12.5 37
CR IP T
ACCEPTED MANUSCRIPT
PT
ED
M
AN US
Fig. 11: The instantaneous interfaces at T= 12.5 for different inflow conditions at a clockwise rotating motion viewed from upstream. (a) laminar inflow; (b) turbulent inflow.
AC
CE
Fig. 12: The instantaneous axial velocity field of the swirling gas-liquid flow (a) in the laminar inflow case and (b) in the turbulent inflow condition at T= 15.00. The black lines denote the stream lines.
38
AN US
CR IP T
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
Fig. 13: The velocity profiles at different axial positions at T= 15.00. (a) axial velocity, (b) radial velocity. The solid and dash lines denote for laminar inflow and turbulent inflow, respectively.
39
AC
CE
PT
ED
M
AN US
CR IP T
ACCEPTED MANUSCRIPT
Fig. 14: Instantaneous velocity vector and axial velocity contour at x/Dout= 3 at T= 12.5 (a) single phase; (b) laminar inflow; (c) turbulent inflow.
40
AN US
CR IP T
ACCEPTED MANUSCRIPT
CE
PT
ED
M
Fig. 15: Interface-vortices interaction for fine case at T= 6.25 (white color for interface and green color for iso-surfaces of the Q criterion).
AC
Fig. 16: PDF of for T= 6.25.
41
CR IP T
ACCEPTED MANUSCRIPT
ED
M
AN US
Fig. 17: PDFs of the alignment between the principal strain rate and i in the recirculation zone.
Fig. 18: PDFs of the alignment between the principal strain rate and i in near regions of the
AC
CE
PT
liquid sheet structures.
Fig. 19: PDFs of the alignment between the principal strain rate and the LS gradient Gi in near regions of the liquid sheet structures.
42
ACCEPTED MANUSCRIPT
Table 1: Comparison of major parameters between realistic gas turbine combustor or the model combustor [28][51][52] and the parameters in the present detailed numerical simulation Parameters Turbine combustor Present simulation 90o
Swirl number
1.0
Density ratio
O(10-1000)
Dynamic viscosity ratio
O(10-100)
Surface tension
0.02-0.04
Reynolds number
O(106)
20 10
0.036
AN US
2000 222
M
O(10 )
Domain size
Number of grids
1 2 3
10D 10D 10D 10D 10D 10D 10D 10D 10D
5123 5123 10243
ED
Case
Grid resolution (µm) 7.8 7.8 3.9
AC
CE
PT
1.0
6
Weber number
Table 2: Simulation cases
90o
CR IP T
Spray cone angle
43
Inflow Laminar Turbulent Laminar