Numerical study of bubbly flow in a swirl atomizer Numerical study of bubbly (cid:13)ow in a swirl atomizer

In the present work, we extend our previous research on swirl nozzles by introducing bubbles at the nozzle inlet. A large-scale hollow cone pressure-swirl atomizer is studied using scale-resolving simulations. The present (cid:13)ow conditions target a Reynolds number range of 600 (cid:20) Re (cid:20) 910 and gas-to-total volumetric (cid:13)ow rate ratios between 0.07 (cid:20) (cid:12) (cid:20) 0.33 with (cid:12) = 0 as an experimental and computational reference. The computational setup has relevance to high-viscosity bio-fuel injection processes. Flow rate ratio and bubble diameter sweeps are carried out to study their e(cid:11)ect on the inner-nozzle (cid:13)ow and the liquid (cid:12)lm characteristics outside the nozzle. The present (cid:13)ow system is shown to pose highly versatile physics including bubble coalescence, bubble-vortex interaction, and faster liquid (cid:12)lm destabilization relative to (cid:12) = 0 case. The main results are as follows. (1) (cid:12) is found to have a signi(cid:12)cant e(cid:11)ect on the bimodal bubble volume probability density function (PDF) inside the swirl chamber. Also, the total resolved interfacial area of the near-ori(cid:12)ce liquid (cid:12)lm increases with (cid:12) . (2) At representative value of (cid:12) = 0.2, the exact bubble size at the inlet is observed to have only a minor e(cid:11)ect on the swirl chamber (cid:13)ow and liquid (cid:12)lm characteristics. (3) The bubble-free ( (cid:12)


I. INTRODUCTION
Here, two-phase gas-liquid flow is studied inside a large-scale pressure-swirl atomizer used for injection of high-viscosity, bio-mass based liquid fuels in boiler applications. Studying in-nozzle two-phase flows has relevance for various scenarios related to injector operation.
For example, the injected liquid can be aerated with gas in order to enhance atomization of the subsequent liquid spray (effervescent atomization) 1,2 . In addition, the gas (vapor) may be generated by boiling process inside the nozzle. For example, modern recovery boilers are often operated in conditions susceptible for flash boiling 3 . In addition to having a strong impact on the spray quality 4 , flashing affects the in-nozzle flow as significant vapor generation may occur already inside the nozzle (approximate void fraction as high as 50-70% at nozzle exit) 3 . Despite the strong influence on the spray characteristics, the effect of the gas presence on the nozzle flow has thus far received little attention.
Pressure-swirl atomizers utilize swirling motion to enhance atomization. In a swirl nozzle, the flow enters a swirl chamber through tangential or axial (vaned) inlet ports which creates a strong vortex characterized by the formation of a gaseous core in typical operating conditions. Customarily, swirl nozzles are operated with purely single-phase liquid feeds. The spray undergoes different transitional modes before reaching a fully developed hollow cone flow regime when the inlet velocity, and therefore the swirl strength, is increased 5 . The hollow cone mode leads to increased spray angle, discharge velocity and amount of air-to-liquid interaction which, in turn, enhance atomization 6 . For review of research on swirl atomizers, see Vijay, Moorthi, and Manivannan 6 and Kang et al. 7 .
Correspondingly, Lin, Kennedy, and Jackson 11 concluded that the aerated-liquid spray is strongly related to the structure of the internal two-phase flow inside a small-scale atomizer (orifice diameter d o ∼ 2 mm) by testing out various discharge passage lengths and converging angles. Unsteady effects in effervescent atomizers have been studied for example by Sen et al. 12  In addition, Eulerian-Eulerian two-fluid models have been used to study flows in effervescent atomizers 21,22 .
In the context of pressure-swirl atomizers, the interaction of the main vortex with the two-phase flow is important. Bubbles in a swirling flow start to migrate towards the center of rotation due to centrifugal force which is larger for liquid 23, 24 . Maneshian, Javadi, and Taeibi Rahni 23 studied numerically the evolution of single bubbles in 3d under the combined effect of rotating flow and gravitational acceleration for a range of Morton and Bond numbers. Similarly, Hsiao, Ma, and Chahine 25 studied a swirl separation process in different gravity environments. Migration towards center of rotation was also observed by Van Nierop et al. 24 , who studied small air bubbles in rotating flow in order to determine the drag and lift coefficients of the bubbles. Effervescent atomization in swirl nozzle context has been investigated by Ochowiak 26 . They studied the discharge coefficient (C d ) of a small-scale (d o ∼ 2 mm) effervescent-swirl atomizer with various fluids and orifice shapes. The nozzle geometry (orifice shape and type) was noted to greatly affect C d . They also demonstrated diminishing C d with GLR, while noting that at higher GLR the shape of the orifice did not play an important role.
Recently, Laurila et al. 27 investigated a large-scale pressure-swirl atomizer (d o ∼ 2 cm) with significant asymmetric features in contrast to the more commonly studied small-scale symmetric nozzle types. Numerical simulations were conducted at a Reynolds number range 3 of 420 ≤ Re ≤ 5300 and several flow modes were identified. S-shaped, transitional and fully developed hollow cone liquid films were observed at Re = 420, 830 and ≥1660, respectively.
Later, the work was extended by conducting both experimental and numerical investigations at a focused Reynolds number range (600 ≤ Re ≤ 910) 28 . At this range, the flow was found to be in a transitional state without a fully developed air core consistent with the previous findings.
In the present work, we extend our previous studies 27,28 by introducing bubbly flow at the inlet of the nozzle. Understanding the characteristics of the in-nozzle two-phase flow has relevance for injector operation and subsequent spray formation. The literature highlights the importance of nozzle geometry in the case of effervescent atomization. However, there are only a few studies on internal flow in effervescent-swirl atomizers and especially on largescale designs such as the current geometry. The presently studied scenario features complex flow physics including developing two-phase pipe flow with bubble coalescence, interaction of gas structures with swirling flow and destabilization of the ejected liquid film by bubbles.
The main goals of this study are to: 1. Relate and compare the present bubbly flow scenarios to the corresponding cases without added gas.
2. Investigate the effect of gas-to-total volumetric flow rate ratio on the general flow characteristics.
3. Study the effect of inlet bubble diameter on the flow dynamics including two-phase flow structures.

A. Governing equations and numerical methods
Here, the same numerical method as in our previous publications 27, 28 is utilized. The volume-of-fluid method 29 is used for the time dependent solution of the two-phase flow. In VOF, the phases are identified by the indicator field, α, which obtains a value of one in the liquid and zero in the gas phase, i.e. α = α l . Based on the incompressible Navier-Stokes 4 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS
and the momentum equation which incorporates the effects of the gravitational, viscous and surface tension forces Above, u, ρ, p, µ, σ and κ are the velocity, density, pressure, viscosity, surface tension and curvature, respectively. The density and viscosity are mixture quantities and they are calculated as ρ = αρ l + (1 − α)ρ g and µ = αµ l + (1 − α)µ g based on the phase-specific values (l for liquid, g for gas). In Eq. 2, the surface tension force is modeled using the Continuum

B. Nozzle geometry
The studied large-scale asymmetric swirl nozzle is shown in Fig. 1. The flow enters the nozzle through an inlet pipe (I) which is aligned 45 • relative to the swirl chamber axis. The enlargement region (II) connects the pipe to the swirl chamber (III) through a single inlet port. As the entry is tangential to the chamber wall, the flow is forced into a swirling motion.

5
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.  This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

C. Computational mesh
The used computational domain and mesh are shown in Fig. 3. Here, the same meshing strategy is used as in our previous work 27, 28 . We utilize an unstructured mesh consisting of

D. Boundary conditions
The resolved bubbly inlet condition is implemented with a uniform Dirichlet condition for the velocity (u = U b , where U b is the mean bulk velocity) and a mapping procedure for and thus the slip ratio is S = u g /u l = 1, the flow rate ratio β is equivalent to the volumeaveraged void fraction inside the reference cylinder. Therefore, the number of bubbles inside the reference volume is determined such that the volume-averaged void fraction is equal to β. The details of the bubbly inlet conditions in the different simulation cases are given in This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI:10.1063/5.0028963
At the outlet boundary, a mixed boundary condition is applied for the variables depending on the flow direction. If the flow direction is towards the domain, Dirichlet conditions are used for both the velocity (u extrapolated from inside the domain) and the indicator field (α = 0). If the flow direction is outwards from the domain, zero gradient Neumann conditions are used for both fields. In both cases, a Dirichlet condition (p = 0 Pa) is used for the total pressure. At the walls, no-slip condition is imposed for the velocity, while the pressure and indicator fields have Neumann conditions.

E. Cases and parameters
In this paper, flow rate ratio and bubble diameter sweeps are carried out. The main flow parameters of the studied cases are shown in Table I, while the corresponding bubbly inlet conditions are illustrated in Fig. 4. To characterize the flows, we use the gas-to-total volumetric flow rate ratio, β = Q g /(Q l + Q g ) 39,40 . In the table, the total mass flow rate (ṁ), the total volumetric flow rate (Q = Q l + Q g ), and β are related byṁ The Reynolds number, Re = ρ l U b d p /µ l , is referred to the liquid properties, the mean bulk velocity, and the inlet pipe diameter (d p ).
In the flow rate ratio sweep (cases V1-V3), the total mass flow is kept constant, while This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.  The particular choice of the mass/volumetric flow rates of the current cases with twophase inlet flow is based on the previous cases of Laurila et al. 28 . The previous work provides experimental references with single-phase inlet flow for the current setup. In addition to the experimental data, we repeat the simulations of the previous work with inlet conditions consistent with the present approach of uniform inlet velocity (cases S1-S4). In this paper, we use the term single-phase to refer to the cases S1-S4 with single-phase inlet flow.

PLEASE CITE THIS ARTICLE AS
In addition to the flow parameters given in Table I, the other presently used parameters are also based on the previous experiments/simulations of Laurila et al. 28 . The liquid/gas densities and viscosities are ρ l = 1230 kg/m 3 , ρ g = 1 kg/m 3 , µ l = 192 mPas and µ g = 0.0148 mPas, respectively, while the surface tension is σ = 0.064 N/m.

10
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

A. Inlet pipe
First, bubble dynamics is investigated inside the inlet pipe. A developing two-phase flow is observed inside the inlet pipe as is illustrated by the instantaneous snapshots of the phase interface (α = 0.5 iso-surface) shown in Fig. 5. As imposed by the inlet conditions, the gas phase initially consists of spherical bubbles which start to rapidly deform in the flow field. In the vicinity of the wall, the developing viscous boundary layer stretches and breaks up the bubbles (I in Fig. 5). In addition, the bubbles are observed to migrate towards the center of the pipe thus forming a nearly gas free wall layer. Such behavior is consistent with previous research on downflow in a vertical channel [41][42][43] . Near the pipe center, the bubbles start to deform and agglomerate. The agglomeration is affected by wake effects between the bubbles. Typically, the trailing bubble is entrained and merged in the wake of the preceding bubble (II). The merging process is similar to a collision between two rising bubbles 44,45 . In such a process, customarily the trailing bubble undergoes more deformation than the leading one 46 . The agglomeration leads to formation of large gas structures which exhibit slug-like behavior (III). The time evolution of the volume-averaged void fraction and total resolved interfacial area inside the inlet pipe are shown in Fig. 6. Here, the simulations with two-phase inlet flow are initialized from a flow field where the nozzle is filled with liquid. Hence, there is 11 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.  to quantify the probability of finding gas volume in gas structures of size d i . More precisely, time-averaged PDFs are calculated based on the instantaneous bubble counts N i,j collected at time instances j. As the total gas volume in a size class i at time j is V i,j ∝ N i,j d 3 i , the normalized probability density function is calculated as:

PLEASE CITE THIS ARTICLE AS
where ∆d = d i − d i−1 is the width of the diameter size class which is constant for all classes.
The PDFs in Fig. 8 (a)  14 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

B. Enlargement region and swirl chamber
The enlargement region redirects the flow from the inlet pipe and forces the swirling motion inside the swirl chamber. Mean velocity fields together with exemplary instantaneous phase interface iso-surfaces are shown in Fig. 9. In the enlargement, two distinct regions are identified. In the high velocity region (I), the flow attaches near the bottom of the

15
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. is visible from the tangential mean velocity profiles shown in Fig. 10. The plots show the two-phase (β > 0) velocity profiles together with the single-phase (β = 0) references of corresponding mass and volumetric flows at different distances from the chamber roof.

PLEASE CITE THIS ARTICLE AS
In addition, previous experimental results 28 are shown when available. Overall, the mean velocity profiles inside the chamber are similar in both the single and two-phase cases. This can be attributed to the fact that the same total volumetric flow rate is applied in both cases. Upon closer inspection of the profiles, it is noted that the local deviations from the single-phase references become more pronounced with increasing β. For example at β = 0.33 (30 mm), larger deviations at y > 0 can be observed due to the higher local presence of gas.
However, it can be concluded that the main velocity characteristics inside the chamber are only moderately sensitive to β.
The bubble volume PDFs in the swirl chamber are depicted in Fig. 11. The PDFs are calculated with the same procedure as described in Sec. III A. In all cases, the distributions are bimodal. Contrary to the inlet pipe (see Fig. 8), a significant portion of the gas volume in the swirl chamber is contained in bubbles with small effective diameter (first peak at This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. migration towards the center of rotation is consistent with previous literature 23, 24 and caused by the density difference and centrifugal acceleration. In Fig. 11 (a), there is a correlation between β and the size of the largest gas structures. The second, higher peak, which corresponds to the large structures, moves towards higher d i as β is increased. In addition, the relative portion between small and large structures shifts towards the large structures with increasing β. In contrast, the PDFs in Fig. 11 (b) show very close correspondence to each other. Interestingly, the inlet bubble diameter (5-15 mm) does not have a significant impact on the PDF. Consistent with this observation, also the polydispersed case behaves similarly.

C. Discharge orifice
Time-averaged velocity and void fraction profiles at the discharge orifice are shown in  To further quantify the effective gaseous core diameter (d ef f ), and therefore, the film 18 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. rate. In the two-phase cases (β > 0), the increase in the total flow rate Q is due to increase in β as the quantities are related by Q = Q l /(1 − β) where Q l is approximately constant.

PLEASE CITE THIS ARTICLE AS
The effective diameter is estimated to be between 4.4-6.1, 5.9-7.5 and 7.5-8.9 for β = 0.07, 0.20 and 0.33, respectively. In addition, d ef f in the single-phase (β = 0) cases is shown in the figure together with previous results of Laurila et al. 28 (same flow parameters except for the parabolic velocity inlet condition). As expected, the inclusion of gas to the inlet flow increases the core diameter relative to the cases with single-phase inlet flow. However, the core size in the two-phase cases does not increase significantly faster as a function of Q than in the single-phase cases. This may be partly attributed to the increasing axial velocity at orifice center where predominantly gas is present (see Fig. 12 (a) and (c)), i.e. the increased gas volumetric flow can be accounted for either by an increase in the gas velocity or gas cross-sectional area. In the present case, the results suggest the former.
The gaseous core existence factor, F ave , is shown in Fig. 13 (b). Here, F ave ∈ [0,1] is defined as the fraction of time when a gaseous core of any size is present at the orifice plane. F ave is unity in the two-phase cases since the plane is, in practice, always occupied by some gas structure. In cases with single-phase inlet flow, the increasing core existence factor indicates a core in a developing state. The observed values are in good agreement 20 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI:10.1063/5.0028963
with previous results 28 . In the single-phase cases, the larger variability of d ef f at low Q is partly attributed to the low existence factor which decreases the effective averaging time.

D. Liquid film
The effect of the gas presence on the ejected liquid film is highly important with regard to subsequent droplet formation and spray characteristics. As illustrated in Fig. 14

21
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

22
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI:10.1063/5.0028963
volume PDFs in the swirl chamber were noted to be bimodal and strongly affected by β. Both the size of the largest gas structures and the scale-separation between the small and large bubbles increased with β. At the discharge orifice (d o ), the effective gas core diameter (d ef f ) increased with β and was found to be in the range 0.2 < d ef f /d o < 0.4. In addition, the core was always present at β > 0, whereas at β = 0 the core was intermittent.
2. It was found that β significantly affected the liquid film outside the nozzle. In summary, we have provided numerical evidence that bubble inclusion at the inlet affects global liquid film characteristics with relevance to atomization. We highlight the rich physics inside effervescent-swirl nozzles including bubble coalescence, bubble accumulation inside the main vortex, and faster destabilization of the liquid film with increasing β.

ACKNOWLEDGMENTS
The authors acknowledge ANDRITZ Oy. for financial support. The computational resources for this study were provided by CSC -Finnish IT Center for Science and the Aalto Science-IT project.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request. Restrictions apply to the availability of the nozzle geom-

24
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.
etry, which was used under license for this study. Data is available from the authors upon reasonable request and with the permission of ANDRITZ Oy.

Appendix A: Validation
The numerical method is validated against the widely used 54-57 benchmark cases of Hysing et al. 38 . The tests consist of two rising bubble cases with different surface tension, density and viscosity ratios. The main flow parameters are given in Table II

25
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.  Table III. The selected cell sizes are based on our previous work (single-phase inlet flow), where resolution similar to M3 was concluded to adequately resolve the developing air core in the orifice vicinity 28 . In the mesh assessment, the time step is varied to maintain Courant number of Co = 0.5, while all other parameters are kept constant between the runs. For M1-M3, the initial conditions are interpolated from the densest mesh.

PLEASE CITE THIS ARTICLE AS
Mean velocity profiles inside the swirl chamber and at the discharge orifice are shown in Fig. 19 (a)-(c). It is noted, that the mean profiles are relatively well captured already 26 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.   The bubble volume PDFs in the swirl chamber are shown in Fig. 19 (f). The relative 27

PLEASE CITE THIS ARTICLE AS
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI:10.1063/5.0028963
weight of the small diameter bubbles of the total gas volume is increased with increasing mesh resolution. In general, the cell size is directly connected to the minimum bubble size that can be resolved. Hence, the smallest resolved bubble is limited by ∆x. We acknowledge that the smallest bubbles present in the current setup cannot be resolved with reasonable grid resolution.