Numerical Models for Viscoelastic Liquid Atomization Spray

Atomization spray of non-Newtonian liquid plays a pivotal role in various engineering applications, especially for the energy utilization. To operate spray systems efficiently and well understand the effects of liquid rheological properties on the whole spray process, a comprehensive model using Euler-Lagrangian approaches was established to simulate the evolution of the atomization spray for viscoelastic liquid. Based on the Oldroyd model, the viscoelastic linear dispersion relation was introduced into the primary atomization; an extended viscoelastic version of Taylor analogy breakup (TAB) model was proposed; and the coalescence criteria was modified by rheological parameters, such as the relaxation time, the retardation time and the zero shear viscosity. The predicted results are validated with experimental data varying air-liquid mass flow ratio (ALR). Then, numerical calculations are conducted to investigate the characteristics of viscoelastic liquid atomization process. Results showed that the evolutionary trend of droplet mean diameter, Weber number and Ohnesorge number of viscoelastic liquids along with axial direction were qualitatively similar to that of Newtonian liquid. However, the mean size of polymer solution increased more gently than that of water at the downstream of the spray, which was beneficial to stable control of the desirable size in the applications. As concerned the effects of liquid physical properties, the surface tension played an important role in the primary atomization, which indicated the benefit of selecting the solvents with lower surface tension for finer atomization effects, while, for the evolution of atomization spray, larger relaxation time and zero shear viscosity increased droplet Sauter mean diameter (SMD) significantly. The zero shear viscosity was effective throughout the jet region, while the effect of relaxation time became weaken at the downstream of the spray field.


Introduction
Atomization of non-Newtonian liquid has become increasingly prevalent in various engineering applications, such as injection of gelled fuels for propulsion system [1], spray coating for painting [2,3] manufacture of pharmaceutical tablets [4], and materials processing for suspension plasma spray [5].The non-Newtonian droplets can provide attractive and variegated properties to meet the specific requirements.The constitution relation for non-Newtonian liquids are complicated since the relationship between shear stress and rate of strain exhibits non-linear behavior.For the widespread shear-thinning liquids, there are several rheological models to describe the variation of viscosity along with strain rate [6,7], such as power-law model, Carreau model and Oldroyd model [8].Among them, Oldroyd model with three constants can effectively describe the rheological characteristics for viscoelastic liquids.For the atomization of non-Newtonian liquid, the external two-phase flow outside of the nozzle orifice can be generally divided into two parts, the near-field region and the far-field domain, where specific physical phenomena govern the flow.
In the near-field region, the primary atomization in which the liquid stream disintegrates into ligaments and drops by interacting with the gas takes place.The studies on primary atomization can be classified into three approaches: experimental measurements, theoretical analysis and modeling calculations.
Based on the techniques of high speed camera and phase Doppler anemometry, experiments provided visualized approaches to reveal the morphology and detailed information of the primary atomization such as breakup length, geometry, spray angle and so on [9].Dumouchel [10] has reviewed the experimental investigate dedicated to the primary atomization step for Newtonian liquids.Aliseda et al. [11] extracted the images of liquid jet break-up which is closed to the nozzle orifice from high-speed visualizations for several Weber and Reynolds numbers for viscous and non-Newtonian liquids.They observed that at lower Weber numbers the Kelvin-Helmholtz instability grew slowly and several intact wavelengths were observed prior to break-up, In contrast at larger weber number, the Rayleigh-Taylor (R-T) instability was arresting.Acceleration of the interface resulted in the R-T instability creating ligaments of fluid which eventually break-up into droplets.This experimental finding could inspire us to use the R-T instability to describe the breakup of cylindrical ligaments.Harrison et al. [12] experimentally examined the spray cone angle with three types of polymer solutions and studied the cone angle as a function of type of polymer and reduced concentration.Chao et al. [13] evaluated the effect of polymeric additives on spray structure from photographs of the sprays.Their studies show the correlations between additive concentrations and the formation of droplets.These experimental investigations are very helpful for understanding the features of viscoelastic liquid breakup and give some guidance to build up models.
The theoretical analysis of primary atomization mainly focused on the derivation of dispersion relations between the surface wave growth rate and the wave number, which is based on two instability analysis: Rayleigh-Taylor (R-T) and the Kelvin-Helmholtz (K-H) instability.Various breakup modes resulted in different dispersion relations.Sirignano and Mehring [14] have reviewed the basic mechanism of distortion and disintegration of liquid streams.For the non-Newtonian liquid, especially viscoelastic, the linear instability dispersion relationship is more complex than the Newtonian liquid.Middleman [15] firstly performed a linearized stability analysis on viscoelastic jet to predict the stream disintegration.Goren and Gorttlieb [16] expanded the Weber's Newtonian liquid linear stability analysis by involving an equivalent liquid viscosity and including an unrelaxed axial tension into the dispersion equation which can be used to explain the stabilization of the viscoelastic jet.Joseph et al. [17] took the high relative velocity between liquid and the acceleration of droplets into consideration and obtained the dispersion relation based on R-T instability.Yang et al. [18] used linear stability analysis to investigate the instability of a three-dimensional viscoelastic liquid jet surrounded by a swirling air stream.
For the modeling of primary atomization, the direct numerical simulation (DNS) of multiphase flows which is based on the one-fluid formalism coupled with interface tracking algorithms [19][20][21][22] is a promising approach to thoroughly investigate the whole process.As reviewed by Villermaux [20] and Gorokhovski et al. [21], the DNS of turbulent primary atomization process is still in its infancy owing to the high computational cost and big numerical challenges.In addition, most of the numerical studies are focused on the Newtonian liquids, few investigations have devoted to viscoelastic liquids since the complicated rheological relations further enhanced the difficulties [22].Instead, based on some simplifications and assumptions of the gas-liquid interface, and involving the instability analysis, we can find a compromising way to establish a viable physical model and predict the size of ligaments/droplets after primary atomization directly and effectively.
At the far-field domain, the droplets produced by primary atomization are unstable in the turbulent spray, and will undergo a series of events such as secondary breakup, collision and coalescence.For non-Newtonian spray, the effects of viscoelasticity on the droplets evolution and distribution in the downstream field have attracted numerous industrial applications [23,24].For different applications, the requirements of droplet size and distribution are different, for instance, finer droplets are desirable for combustions and painting spray [25].On the contrary, in pesticide sprays [26] and anti-misting jet, too small drops are undesirable.The evolution of droplets at the downstream of the spray is critical to the final outcomes.However in the downstream of the spray, there is a scarcity in literature to investigate the evolution and distribution of atomized drops, especially for the atomization of viscoelastic liquids.
In the experimental aspect, most studies devoted to determine the features of non-Newtonian droplet secondary breakup.Rivera [27] used the high-speed photography to find that Non-Newtonian drop secondary breakup mode is qualitatively similar to those observed for Newtonian drops.The only major difference is that the bag stretches more before it ruptures and the resulting breakup fragments persist longer than their Newtonian counterparts.Dechelette et al. [28] experimentally investigated the combined effect of polymers and soluble surfactants on dynamics of jet breakup and formation of satellite drop.Ma et al. [29] found the particle size distribution was fit to the Rosin-Rammler function through 3-D phase Doppler methods.Hartranft and Settles [25] indicated that extensional viscosity had the dominant influence on sheet breakup features compared to Weber number, which means the traditional Weber number may not qualify to describe the breakup of non-Newtonian fluids.These experimental investigations can support us to extend the breakup mode of Newtonian drop to the viscoelastic version.
In the modeling aspects, for the secondary atomization, some investigations have devoted to build up physical models to illustrate droplet breakup.Xue et al. [30] used a two-dimensional computational approach to predict the effect of coupled surfactant and non-Newtonian mechanisms on the formation of satellite drops.Rivera [27] extended the Taylor analogy breakup (TAB) model to describe the inelastic non-Newtonian liquid with power-law model and validated the model with experimental data.It is noted that most of previous models are unable to account for the whole spray pattern and its evolution without considering the poly-dispersion of droplets as well as collision and coalescence [31].There is very few investigations devoted to establish the numerical models accounting for the evolution of viscoelastic liquid spray [32], due to the complexity of rheological behavior and statistical difficulties.
As such, the motivation of this paper is to logically expand the above mentioned research efforts and aim specially to shed light on the simulations of the downstream of the spray for viscoelastic liquid.The numerical model will consider the primary breakup, secondary breakup, droplet tracking and collision.For primary breakup, Goren and Joseph's dispersion relationship will be used to describe the linear instability of viscoelastic liquid.For secondary atomization, the TAB is extended to the viscoelastic version based on the rheological analysis.For droplet collision, the criterion of coalescence is modified by the rheological properties.The numerical calculations will try to give some detailed information about the evolutions of droplet distribution, size distribution and critical dimensionless numbers (such as Weber number and Ohnesorge number).Results will involve three parts.Firstly, the primary atomization model is validated and compared by experimental data.Then, the evolution of spray field is investigated in detailed by comparison of polymer solution and water.Finally, the effects of rheological properties on spray performance are discussed.

Mathematical Models
The mathematical model considers two parts.As shown in Figure 1, the first one is the primary breakup sub-model, which predicts gas velocity at the nozzle exit and the droplet mean size.The second sub-model calculates the evolutions of the gas flow and the droplets spray accounting for secondary breakup, collision and coalescence, momentum transfer.The initial droplets characteristics in the second sub-model are provided by the first sub-model calculation.In droplet-gas two-phase flow, the two-way coupling of gas and droplet phase is considered, while the effect of primary breakup jet on the environment gas is neglected.

Sub-Model of Primary Atomization
The primary atomization is a process that the liquid stream disintegrates into ligaments and drops by interacting with the gas.It is a complex gas-liquid flow.In order to simplify the process, we extracted the main features of the primary atomization from experimental observations.At present study, effervescent atomizer is chosen since it has so far been found effective for high viscosity liquids [33], in which the gas is introduced into the liquid at low pressure to form a bubbly two-phase mixture upstream of the orifice.As indicated by experimental observation, the flow pattern closed to the nozzle exit is similar to annular.Then, as illustrated in Figure 1, the primary breakup sub-model for effervescent atomizer can be established through three steps.
The first step is estimation of the annular sheets as well as the cylindrical ligaments.The momentum equation of the annular flow combined with the state equation can be written as: where ρ g is the gas density, g r is the radius of gas flow, l m  is the mass flow rate of liquid, and ALR is the air-liquid ratio by mass.The radius of gas flow can be written in terms of orifice radius orific r and void fraction ϑ , . The interface velocity slip ratio slip v under different flow rate is expressed as [34]: where C is the experimental coefficient of the mass flow rate scaling.The interface velocity slip ratio and void fraction are related by: , where R is the gas constant, T is the temperature, in P is the pressure inside the nozzle, and out P is the pressure outside the nozzle.

Sub-Model of Primary Atomization
The primary atomization is a process that the liquid stream disintegrates into ligaments and drops by interacting with the gas.It is a complex gas-liquid flow.In order to simplify the process, we extracted the main features of the primary atomization from experimental observations.At present study, effervescent atomizer is chosen since it has so far been found effective for high viscosity liquids [33], in which the gas is introduced into the liquid at low pressure to form a bubbly two-phase mixture upstream of the orifice.As indicated by experimental observation, the flow pattern closed to the nozzle exit is similar to annular.Then, as illustrated in Figure 1, the primary breakup sub-model for effervescent atomizer can be established through three steps.
The first step is estimation of the annular sheets as well as the cylindrical ligaments.The momentum equation of the annular flow combined with the state equation can be written as: where ρ g is the gas density, r g is the radius of gas flow, .m l is the mass flow rate of liquid, and ALR is the air-liquid ratio by mass.The radius of gas flow can be written in terms of orifice radius r orific and void fraction ϑ, r g = √ ϑr orifice .The interface velocity slip ratio v slip under different flow rate is expressed as [34]: where C is the experimental coefficient of the mass flow rate scaling.The interface velocity slip ratio and void fraction are related by: By solving Equations ( 1)-( 3), ρ g , ϑ and v slip can be calculated for different operating conditions.The thickness of annular liquid sheet is then calculated from δ lig = 4(r o − r g )/ √ π, which is also the diameter of the typical cylindrical ligament.Besides, the gas velocity inside the nozzle orifice can be expressed as V g1 = .m g /(A noz × ϑ × ρ g ), where A noz is the area of the nozzle orifice and .m g is the gas flow rate.According to the energy conservation equation, the gas velocity at the nozzle exit can be derived from V g2 = V g1 2 + 2RT log(P in /P out ), where R is the gas constant, T is the temperature, P in is the pressure inside the nozzle, and P out is the pressure outside the nozzle.Secondly, the formed ligaments are approximated by cylindrical jets.The ligaments will breakup into fragments at the wavelength of the most rapidly growing wave.The dispersion relations based Energies 2016, 9, 1079 5 of 17 on instability analysis are used to determine the critical wave number.Here, Joseph's dispersion relationship [17] is used to describe the linear instability of viscoelastic liquid.
Joseph's dispersion relationship is derived by considering the acceleration of a drop exposed to a high speed air stream.The equations of ligament motion coupled with boundary conditions are solved to deduce the expression of the disturbed interface: Equation ( 4) is an approximate analysis of Rayleigh-Taylor instability based on viscoelastic potential flow which gives the critical wavelength and growth rate to within less than 10% of the exact theory.In Equation ( 4), k is the magnitude of the wave number, α is the amplification rate, and σ is the liquid surface tension.µ ve is the effective shear viscosity of the liquid, can be assumed to be αρ l /2kµ ve .
. a is the acceleration of the liquid ligament: .
where c d is the drag coefficient.As the ligament Reynolds number Re lig = ρ g δ lig u g − u l /µ g is between 3 and 1000, here c d = 24(1 + Re lig 2/3 /6)/Re lig .Then, Equation ( 4) reduces to be a three-degree polynomial of α as follows and can be solved analytically: Based on the above dispersion relation, the critical breakup wave number that corresponds to the maximum growth rate can be used to calculate the wavelength of a typical ligament fragment.Finally, assuming that each fragment is stabilized to one spherical droplet under the influence of surface tension, the Sauter mean drop diameter SMD (defined as the ratio of volume-surface mean diameter ) can then be calculated from the conservation of mass:

Sub-Model of Droplet-Gas Two-Phase Flow
The second sub-model is based on hybrid Euler-Lagrangian coordinate system to describe the evolution of spray.The flow field is axisymmetric, and the parcels are tracked in three-dimensional coordinate.

Gas Phase
For the turbulent gas spray, the conservation equations of mass, momentum and energy are written as follows: Energies 2016, 9, 1079 where u g is the gas velocity vector, F is the momentum transfer between the gas and the discrete liquid phase, K is the turbulent kinetic energy per unit mass, σ s is the viscous stress tensor, q is the heat flux vector, ε is the viscous dissipation rate, and Q is the heat transfer between the gas and the liquid phase.Turbulence effects were modeled by the standard k − ε model.F and Q are defined as: , where α g is the fraction of gas in the computational cell, V cell is the cell volume, N p is the number of droplets in this parcel, and F p is the force acting on each droplet, Q p is the heat flux through surface of each droplet, and u l is the droplet velocity.
For the discrete droplets, the momentum equation can be given as: , where F p is the force exerted on droplet, m l is the mass of droplet, r l is the radius of droplet, g is the gravity force, and C D is the coefficient of the drag force and can be expressed as: Here the droplets are assumed to be spherical, and the effect of droplet deformation on drag coefficient is not taken into consideration.

Liquid Phase Secondary Atomization
In this section, we first analyze the constitutive relations of the Newtonian and viscoelastic liquids, introducing more physical properties to describe the viscoelastic contributions.Then, the TAB model which is based on Taylor's analogy between an oscillating droplet and a spring-mass system was extend to a non-Newtonian version to describe the deformation and distortion of viscoelastic liquids.
For the Newtonian fluid, the shear stress τ ij can be expressed as: where µ is the Newtonian viscosity, and d ij is the rate of deformation tensor.In cylindrical coordinates, d ij with (i, j) component can be taken as: For the viscoelastic fluid, according to Jeffreys' study [35], the rheological model can be expressed with three constants (µ, λ 1 , λ 2 ) as follows: Using the kinetic theory of transient mechanical response in small amplitude motions, the variables (u r , u z ) in Equations ( 9) and ( 10) can be written as f (r, z, t) = F(r)e ikz+αt , where k is the wave number, and α is the growth rate of a perturbation.As deduced by Middleman [15], Equations ( 9) and (10), respectively become the following: By comparing the above two dynamical equations, it concludes that they are identical, the Newtonian viscosity µ can be replaced by µ ve = µ(1 + αλ 2 )/(1 + αλ 1 ) for the viscoelastic model.Based on the above analysis, the original TAB model proposed by Amsden et al. [36] can be used as a foundation for the viscoelastic version.The derivation of viscoelastic version of TAB model include three aspects as below.
First is the modification of deformation equation.The TAB model is based on the Raleigh-Taylor instability analysis, describing the droplet distortion by a forced, damped, harmonic oscillator in which the forcing term is given by the aerodynamic interaction between the droplet and gas; the damping is due to the liquid viscosity, and the restoring force is supplied by the surface tension.Droplet distortion is denoted by the deformation parameter, ζ = 2τ/r, where τ denotes the change of radial cross-section from its equilibrium position and r is the initial drop radius.The deformation equation of ζ in terms of the normalized distortion parameter is: where ρ l is the density, µ l is the viscosity, σ is the surface tension, U rel is the relative drop-gas velocity, and the subscripts g and l denote the gas and liquid, respectively.It is assumed that the external aerodynamic force is normal to the drop and the turbulence fluctuations of the surrounding gas are ignored.In reality, there are various modes when droplets breakup as they are subjected to the stresses of the surrounding turbulence gas.However, it is impossible to describe all secondary breakup modes in a single model, therefore only the fundamental mode of oscillation is considered in TAB model.For a viscoelastic version of TAB model, the constant Newtonian shear viscosity µ l was replaced by a viscoelastic viscosity µ ve .In the expression of µ ve , the growth rate of a perturbation α can be taken by the rate of strain .γ: where λ 1 is the fluid relaxation time, λ 2 is the retardation time, and µ 0 is the zero shear viscosity.The rate of strain .γ is approximated as .γ = U rel /2r 0 .As indicated in Equation ( 14), on the one hand, µ ve will approach zero shear viscosity µ 0 when the rate of strain .γ is close to 0. On the other hand, µ ve will approach µ 0 λ 2 /λ 1 when .γ becomes larger.Therefore, when λ 2 < λ 1 and the rate of strain is increasing, µ ve < µ 0 , which means that the liquid is shear thinning.For dilute polymer solutions, the relationship of retardation time and relaxation time is λ 2 = λ 1 η solvent /η solution [37], where η solvent and η solution are the viscosity of solvent and solution, respectively.
Second, the breakup criterion is changed.In TAB model, drop breakup occurs when the normalized drop distortion ζ(t) exceeds the critical value of 1.For Newtonian droplet, the stability limit has been determined experimentally to be We g = ρ g U 2 rel r/σ = 6.For high viscous liquid, the effect of viscosity on the stability should be related to the breakup criterion.As proposed by Brodkey [38], the critical Weber number can be modified by introducing the Ohnesorge number: We g.crit = 6(1 + 1.077Oh 1.6 ), Oh = µ l / 2ρ l σr (15) when Oh ≤ 0.1, the value of this expression approximates to 6, which means the case of low viscosity for Newtonian liquid is recovered.When Oh is higher, the stability limit is increased and scales as µ 1.6 l in the asymptotic limit.Based on the Brodkey's empirical relation, an effective Weber number is introduced to instead of the regular one [38]: for the viscoelastic liquid, µ ve is used instead of µ l in Oh number.Third, the creation of the product droplets is also revised.The population dynamics is used to define the rate of droplet creation.For each breakup event, with the mass conservation principle, it is assumed that the number of product droplets is proportional to the number of critical parent drops, where the proportionality constant depends on the drop breakup regime: where m(t) denotes the mean mass of the product drop distribution, and the breakup frequency K bu depends on the drop breakup regimes.Breakup frequency K bu = 0.05ω as suggested by Amsden et al. [36] Energies 2016, 9, 1079 8 of 17 is used.The drop oscillation frequency ω is given by ω 2 = 8σ/(ρ l r 3 ) − 25µ 2 l /(4ρ 2 l r 4 ), where µ l is also replaced by µ ve for viscoelastic version.

Droplets Collision
Droplet collision is calculated by a statistical method, rather than a deterministic method since the real collisions process is highly stochastic.All the droplets in one computational parcel behave in the same manner, they either do or do not collide, which is depended on the probability of the collision whether larger than a random number.For collision result, an impact parameter, ψ = χ/(r 1 + r 2 ) is involved to judge the outcomes, where χ is the distance between the center of one drop and the relative velocity vector U rel , r 1 and r 2 are the radii of small and large droplets, respectively.The criteria, χ cr determines the transition boundary: drops coalescence when χ ≤ χ cr , and bounce when χ > χ cr .χ cr can be expressed as: where the collision Weber number for Newtonian liquid defined as: For viscoelastic liquid, We e f f l = We l /(1 + 1.077Oh 1.6 ) is used instead of We l .f (γ) is a function of the radius ratio γ: According to the conservation of momentum, the drop velocity after a bouncing collision is [36]: . The droplet velocity after coalescence is calculated as . m 1 and m 2 are the mass of the small and large droplets, respectively; and u 1 and u 2 are the velocity of the small and large droplets, respectively.The new radius after droplet coalescence is r new = (r 3 1 + r 3 2 ) 1/3 .

Numerical Setup
Figure 2 shows the computational mesh and boundary conditions of numerical model.The radius and axial lengths of the computational regime are 6 cm and 15 cm respectively, and 2π in the circular direction.The computational domain is in a cylindrical coordinate system with size of 56 (radial) × 65 (axial) × 32 (azimuthal).Finer meshes are used for the core region of the spray.The grid dependency has been validated in our previous studies [39,40] and the computational results are convergence under present mesh size.
Energies 2016, 9, 1079 9 of 17 the circular direction.The computational domain is in a cylindrical coordinate system with size of 56 (radial) × 65 (axial) × 32 (azimuthal).Finer meshes are used for the core region of the spray.The grid dependency has been validated in our previous studies [39,40] and the computational results are convergence under present mesh size.As presented in Table 1, there are five types of boundary conditions involved for the gas field.The parcels are introduced in the system at the nozzle exit, which has the original size calculated from primary breakup sub-model with Rosin-Rammler distribution: , where q co indicates the value of the width of the distribution.The Rosin-Rammler distribution is a prevalent size distribution used in the atomization field [41], and the experimental data also supported the Rosin-Rammler function for non-Newtonian fluid atomization spray [31].The number of parcels should be As presented in Table 1, there are five types of boundary conditions involved for the gas field.The parcels are introduced in the system at the nozzle exit, which has the original size calculated from Energies 2016, 9, 1079 9 of 17 primary breakup sub-model with Rosin-Rammler distribution: F R−R (D) = 1 − exp −(D/SMD 0 ) q co , where q co indicates the value of the width of the distribution.The Rosin-Rammler distribution is a prevalent size distribution used in the atomization field [41], and the experimental data also supported the Rosin-Rammler function for non-Newtonian fluid atomization spray [31].The number of parcels should be large enough to reduce statistical fluctuation.Usually, a total number of 10,000 computational parcels are injected into the flow for each case.
Table 1.Five types of boundary conditions involved in the model.

Validation of Primary Atomization
Figure 3 shows a comparison of experimental data [42] and model predictions for polymer solution atomization at P inj = 0.35 MPa, D noz = 1 mm and ALR ranging from 0.02 to 0.11.The physical and rheological properties are listed in Table 2.The polymer solution were formulated by adding poly(ethylene oxide) to solvent having a common composition of 60% glycerine-40% water solution by mass.Polyethylene oxide (PEO) molecular weights of 100,000 were used at mass concentrations of 0.001%, 0.01% and 0.15% [42].The dispersion relations proposed by Joseph are solved analytically.Results show that, with an increase in air-liquid mass flow ratio (ALR), droplet SMD will decline significantly, especially at lower ALR.The agreement between the predictions and measurements is qualitative achieved within 5%-20%.Joseph's model is effective to capture the influence of rheological properties on the breakup process of polymer solution and to quantitatively predict the resulting droplet sizes.However the shortcoming is that the divergence may occur when the liquid ligament is not similar to the flattened drop.Therefore, at lower ALR, the divergence is obviously since the liquid ligaments is slender.measurements is qualitative achieved within 5%-20%.Joseph's model is effective to capture the influence of rheological properties on the breakup process of polymer solution and to quantitatively predict the resulting droplet sizes.However the shortcoming is that the divergence may occur when the liquid ligament is not similar to the flattened drop.Therefore, at lower ALR, the divergence is obviously since the liquid ligaments is slender.

Evolution of Spray Field
Among the investigations of the far-field region, the particle size evolution along with the axial distance is the primary interest.Figure 4 firstly validates the numerical results of water's SMD with experimental data [43] and then compares the downstream variation of SMD for water and polymer solution.Both experiment and simulations were conducted based on the same operating conditions.The operating parameters and liquid physical properties are listed in Table 3.As shown in Figure 4, an agreement between the simulation results and experimental data has been obtained for the case of water, and the divergence of the experiment and simulation is between 2.5% and 25%.In Figure 4, the viscoelastic liquid added the difficulties in the breakup of ligaments

Evolution of Spray Field
Among the investigations of the far-field region, the particle size evolution along with the axial distance is the primary interest.Figure 4 firstly validates the numerical results of water's SMD with experimental data [43] and then compares the downstream variation of SMD for water and polymer solution.Both experiment and simulations were conducted based on the same operating conditions.The operating parameters and liquid physical properties are listed in Table 3. and droplets, impairing the atomization effect by increasing the drop size in the whole process, which is consistent with the results of other researchers [12,44].Further, in both the cases of water and polymer solution, the droplet diameter first decreases and then increases with the axial distance.This phenomenon can be explained by the competition between breakup and coalescence effects.Figure 5 shows the occurrence number of droplets breakup and collision.Close to the nozzle exit, with high kinetic energy and large velocity difference between the droplet and atomization gas, turbulence breakup dominates, causing the drop size to rapidly decrease, while, further downstream, as the kinetic energy decay and small velocity difference, turbulence breakup ends and droplet coalescence plays an important role, causing the drop size to gradually increase.Besides the common and well-known points, there are two noteworthy differences between  As shown in Figure 4, an agreement between the simulation results and experimental data has been obtained for the case of water, and the divergence of the experiment and simulation is between 2.5% and 25%.In Figure 4, the viscoelastic liquid added the difficulties in the breakup of ligaments and droplets, impairing the atomization effect by increasing the drop size in the whole process, which is consistent with the results of other researchers [12,44].Further, in both the cases of water and polymer solution, the droplet diameter first decreases and then increases with the axial distance.This phenomenon can be explained by the competition between breakup and coalescence effects.Figure 5 shows the occurrence number of droplets breakup and collision.Close to the nozzle exit, with high kinetic energy and large velocity difference between the droplet and atomization gas, turbulence breakup dominates, causing the drop size to rapidly decrease, while, further downstream, as the kinetic energy decay and small velocity difference, turbulence breakup ends and droplet coalescence plays an important role, causing the drop size to gradually increase.
breakup dominates, causing the drop size to rapidly decrease, while, further downstream, as the kinetic energy decay and small velocity difference, turbulence breakup ends and droplet coalescence plays an important role, causing the drop size to gradually increase.Besides the common and well-known points, there are two noteworthy differences between these two types of liquids.First, in Figure 4, the gap between the mean size of water and polymer solution is shortened along with the development of the jet.Before the nadir of the curve, the case of polymer solution decays more steeply than that of water.In contrast, after the nadir of the curve, the case of water rises more obviously than that of polymer solution.On the rising part of the curves, the increasing degree of water is about 123%, while the case of polymer solution is less than 57% (this value is calculated from ).This phenomenon indicates that the mean size of polymer solution increases more gently than that of water at the axial distance (5 100) / x noz L D -, which is benefit to control the desirable size in industrial applications.Secondly, in Figure 5, the occurrence numbers for water are higher than polymer solutions for both breakup and collision process in the upstream of the spray.This is because droplets breakup time and frequency depend on the competition between the aerodynamic force and restoring force.An increase in liquid Besides the common and well-known points, there are two noteworthy differences between these two types of liquids.First, in Figure 4, the gap between the mean size of water and polymer solution is shortened along with the development of the jet.Before the nadir of the curve, the case of polymer solution decays more steeply than that of water.In contrast, after the nadir of the curve, the case of water rises more obviously than that of polymer solution.On the rising part of the curves, the increasing degree of water is about 123%, while the case of polymer solution is less than 57% (this value is calculated from (SMD L x /D noz =140 − SMD L x /D noz =4 )/SMD L x /D noz =4 ).This phenomenon indicates that the mean size of polymer solution increases more gently than that of water at the axial distance (5-100)L x /D noz , which is benefit to control the desirable size in industrial applications.Secondly, in Figure 5, the occurrence numbers for water are higher than polymer solutions for both breakup and collision process in the upstream of the spray.This is because droplets breakup time and frequency depend on the competition between the aerodynamic force and restoring force.An increase in liquid viscosity will delay and weaken the breakup process, resulting in higher droplet size.The decrease of the frequency of droplet secondary breakup would also lead to the number of product droplets decreasing, therefore the collision occurrence number of polymer solutions is obviously less than that of water in the upstream region.However, when the axial distance beyond 40L x /D noz , the collision occurrence number for polymer solution exceeds that of water.Since the droplets collision process were influenced by the droplet size and concentration rate.Concentrated drop distribution causes the collision happening frequently.The outcome of collision was also affected by the liquid viscosity.For the polymer solution droplets, impact parameter increased due to higher viscosity, leading to an increase in possibility of droplet coalescence after collision.However, the change of collision occurrence number did not influence the trend of SMD in Figure 4 significantly, due to the occurrence number is too small compared to the total particle number.The most important factors influencing the droplet breakup can be grouped into two dimensionless numbers, the Weber number (We) and the Ohnesorge number (Oh).The average Weber number is defined as We = ρ g U rel 2 SMD/σ l .Here the modified weber number (Equation ( 16)) is used for the polymer solution.The average Ohnesorge number is expressed as Oh = µ l / √ ρ l σ l SMD.For viscoelastic liquid, the constant viscosity µ l was replaced by a viscoelastic viscosity µ ve .
The evolutions of Weber and Ohnesorge number along with the axial direction are illustrated in Figure 6.Weber number indicates the ratio of the inertial force to the surface tension force, in which the trend is basically determined by the gas-liquid relative velocity.At the nozzle exit, We is higher due to the large difference between the velocity of droplet and gas, causing the secondary breakup of droplets.In this period, We of polymer solution is larger than that of water due to the larger SMD at the beginning.Then at downstream of the spray, We dramatically declines as the kinetic energy dissipated and We of polymer solution is less than that of water since We for viscoelastic liquid is modified by the effective viscosity.The effective Weber number decreased as liquid viscosity increasing.The Ohnesorge number is the ratio of an internal viscosity force to an interfacial surface tension force, which significantly increases with the increase of liquid viscosity.For both water and polymer solution, the Oh curves have a peak at about 4 mm from the nozzle exit, which is corresponding to the nadir of the curve of droplet size (as shown in Figure 4).In short, the case of polymer solution has a greater magnitude of change for both We and Oh due to introducing the changes in viscosity.
For the polymer solution droplets, impact parameter increased due to higher viscosity, leading to an increase in possibility of droplet coalescence after collision.However, the change of collision occurrence number did not influence the trend of SMD in Figure 4 significantly, due to the occurrence number is too small compared to the total particle number.
The most important factors influencing the droplet breakup can be grouped into two dimensionless numbers, the Weber number (We ) and the Ohnesorge number ( O h ).The average Weber number is defined as . Here the modified weber number (Equation ( 16)) is used for the polymer solution.The average Ohnesorge number is expressed as μ / ρ σ l l l Oh SMD = .For viscoelastic liquid, the constant viscosity μ l was replaced by a viscoelastic viscosity μ ve .
The evolutions of Weber and Ohnesorge number along with the axial direction are illustrated in Figure 6.Weber number indicates the ratio of the inertial force to the surface tension force, in which the trend is basically determined by the gas-liquid relative velocity.At the nozzle exit, We is higher due to the large difference between the velocity of droplet and gas, causing the secondary breakup of droplets.In this period, We of polymer solution is larger than that of water due to the larger SMD at the beginning.Then at downstream of the spray, We dramatically declines as the kinetic energy dissipated and We of polymer solution is less than that of water since We for viscoelastic liquid is modified by the effective viscosity.The effective Weber number decreased as liquid viscosity increasing.The Ohnesorge number is the ratio of an internal viscosity force to an interfacial surface tension force, which significantly increases with the increase of liquid viscosity.For both water and polymer solution, the O h curves have a peak at about 4 mm from the nozzle exit, which is corresponding to the nadir of the curve of droplet size (as shown in Figure 4).In short, the case of polymer solution has a greater magnitude of change for both W e and O h due to introducing the changes in viscosity.In addition to drop size and critical numbers, liquid viscosity can also affect the spray pattern.Figure 7 shows the spatial distributions of droplet size and number for water and polymer solution, respectively.The cross-sectional plane is located at the axial distance of 50 mm.It can be seen from Figure 7 that the spatial distribution of droplets becomes wider for water.The distribution range of In addition to drop size and critical numbers, liquid viscosity can also affect the spray pattern.Figure 7 shows the spatial distributions of droplet size and number for water and polymer solution, respectively.The cross-sectional plane is located at the axial distance of 50 mm.It can be seen from Figure 7 that the spatial distribution of droplets becomes wider for water.The distribution range of water case is about 40% larger than the polymer solution.It can be explained as the ligaments and fragments are more likely to adhere due to viscoelastic properties.Furthermore, finer droplets are easier for the gas to accelerate and carry them outward of the centerline, resulting in the spray angle will expand.In turn, the distribution will also affect the droplet size.Concentrated drop distribution causes the collision happening frequently and the droplet size increasing, which is in accordance with that indicated in Figure 5.
water case is about 40% larger than the polymer solution.It can be explained as the ligaments and fragments are more likely to adhere due to viscoelastic properties.Furthermore, finer droplets are easier for the gas to accelerate and carry them outward of the centerline, resulting in the spray angle will expand.In turn, the distribution will also affect the droplet size.Concentrated drop distribution causes the collision happening frequently and the droplet size increasing, which is in accordance with that indicated in Figure 5.

Effects of Viscoelastic Parameters
Figure 8 shows the analysis of the effects of rheological properties on primary atomization.Figure 9 explores the evolutions of droplet size along with axial direction for various liquid properties.For polymer solutions, the rheology properties can be adjusted through the addition of polymer molecules.The liquid physical properties of base case are as the same as that of code 2 in Table 2.All cases in Figures 8 and 9 share the same operating conditions except the specified parameters in the legend of the figure.In order to compare the behavior of Newtonian and viscoelastic liquid, the calculated results of water spray have been added into Figures 8 and 9.
Obviously, from Figures 8 and 9, the increases in relaxation time and zero-shear viscosity cause the atomization quality deteriorates since the growth rate of the disturbance wave reduces and the liquid becomes stringier to resist the disruptive force, while small surface tension brings the droplet size to decrease significantly.The case with surface tension of

Effects of Viscoelastic Parameters
Figure 8 shows the analysis of the effects of rheological properties on primary atomization.Figure 9 explores the evolutions of droplet size along with axial direction for various liquid properties.For polymer solutions, the rheology properties can be adjusted through the addition of polymer molecules.The liquid physical properties of base case are as the same as that of code 2 in Table 2.All cases in Figures 8 and 9 share the same operating conditions except the specified parameters in the legend of the figure.In order to compare the behavior of Newtonian and viscoelastic liquid, the calculated results of water spray have been added into Figures 8 and 9.
brings the ligaments and droplet less stability.This finding shows the importance of selecting solvents with lower surface tension for the polymer solution when finer atomization effects are desirable.Moreover, the effects of relaxation time become weaken at the downstream of the spray field.The case with shorter relaxation time approaches to the base case since the influence of relaxation time depends on the gas-liquid relative velocity.When the relative velocity declined rapidly, the effective viscosity of viscoelastic liquid is mainly up to zero shear viscosity.

Conclusions
A comprehensive mathematical model has been proposed to simulate the outflow of the effervescent atomization for viscoelastic liquid, accounting for primary atomization droplet tracking, secondary breakup and droplet collision.The computational models have included two parts.The first one is the primary breakup sub-model.With the simplification and assumptions of the atomization of liquid stream, the primary breakup sub-model has been established based on the viscoelastic linear dispersion relation.The second sub-model can calculate the evolutions of the droplets spray, in which the gas-liquid two-phase flow has been simulated by the Euler-Lagrangian approach.Both the droplet secondary atomization model and collision model have been extended to the viscoelastic version by introducing rheological parameters.The main findings of the present study can be summarized in four points: solvents with lower surface tension for the polymer solution when finer atomization effects are desirable.Moreover, the effects of relaxation time become weaken at the downstream of the spray field.The case with shorter relaxation time approaches to the base case since the influence of relaxation time depends on the gas-liquid relative velocity.When the relative velocity declined rapidly, the effective viscosity of viscoelastic liquid is mainly up to zero shear viscosity.

Conclusions
A comprehensive mathematical model has been proposed to simulate the outflow of the effervescent atomization for viscoelastic liquid, accounting for primary atomization droplet tracking, secondary breakup and droplet collision.The computational models have included two parts.The first one is the primary breakup sub-model.With the simplification and assumptions of the atomization of liquid stream, the primary breakup sub-model has been established based on the viscoelastic linear dispersion relation.The second sub-model can calculate the evolutions of the droplets spray, in which the gas-liquid two-phase flow has been simulated by the Euler-Lagrangian approach.Both the droplet secondary atomization model and collision model have been extended to the viscoelastic version by introducing rheological parameters.The main findings of the present study can be summarized in four points: Obviously, from Figures 8 and 9, the increases in relaxation time and zero-shear viscosity cause the atomization quality deteriorates since the growth rate of the disturbance wave reduces and the liquid becomes stringier to resist the disruptive force, while small surface tension brings the droplet size to decrease significantly.The case with surface tension of 20 × 10 −3 J/m 2 is almost identical to the case of water.In both instability analysis and TAB model, the surface tension can resist the occurrence and development of instability and smooth the disturbance.That is, a decrease in surface tension brings the ligaments and droplet less stability.This finding shows the importance of selecting solvents with lower surface tension for the polymer solution when finer atomization effects are desirable.Moreover, the effects of relaxation time become weaken at the downstream of the spray field.The case with shorter relaxation time approaches to the base case since the influence of relaxation time depends on the gas-liquid relative velocity.When the relative velocity declined rapidly, the effective viscosity of viscoelastic liquid is mainly up to zero shear viscosity.

Conclusions
A comprehensive mathematical model has been proposed to simulate the outflow of the effervescent atomization for viscoelastic liquid, accounting for primary atomization droplet tracking, secondary breakup and droplet collision.The computational models have included two parts.The first one is the primary breakup sub-model.With the simplification and assumptions of the atomization of liquid stream, the primary breakup sub-model has been established based on the viscoelastic linear dispersion relation.The second sub-model can calculate the evolutions of the droplets spray, in which the gas-liquid two-phase flow has been simulated by the Euler-Lagrangian approach.Both the droplet secondary atomization model and collision model have been extended to the viscoelastic version by introducing rheological parameters.The main findings of the present study can be summarized in four points: (1) Based on the comparison with experimental data and the predicted results of other dispersion relation, the primary model using Joseph's relation effectively captured the influence of rheological properties on the breakup process of polymer solution and to quantitatively predict the resulting droplet sizes.(2) The evolutionary trend of droplet mean diameter, Weber number and Ohnesorge number of viscoelastic liquids along with axial direction were qualitatively similar to that of Newtonian liquids.The SMD for the case of polymer solutions were obviously larger than that of water throughout the spray field.While, the gap of SMD between water and polymer solution became shorten along with the development of the jet.This phenomenon indicated that the mean size of polymer solution increased more gently than that of water at the axial distance (5-100)L x /D noz , which benefited the stable control of the desirable size in the applications.(3) The changes of Weber and Ohnesorge number for polymer solution were more dramatic than that of water, since the viscosity of the viscoelastic liquid was influenced by the relative velocity of gas and droplets.Besides, the distribution range of water was about 40% larger than that of polymer solution in present study, which indicated the spray angle became narrow for the viscoelastic liquids.(4) In primary atomization, the surface tension was the dominating physical property for viscoelastic liquid rather than relaxation time and zero shear viscosity.When the finer atomization effects were desirable, it was important to select solvents with lower surface tension for the polymer solutions, while, concerning the evolution of atomization spray, larger relaxation time and zero shear viscosity increased droplet SMD significantly.The zero shear viscosity was effective throughout the jet region, while the effect of relaxation time became weaken at the downstream of the spray field.

Figure 2 .
Figure 2. Geometry and computational mesh: a middle section plane.(a) Middle section plane; and (b) Cross section plane.

Figure 2 .
Figure 2. Geometry and computational mesh: a middle section plane.(a) Middle section plane; and (b) Cross section plane.

Figure 3 .Table 2 .
Figure 3.Comparison of experimental data and model predictions for the Sauter mean diameter (SMD) versus air-liquid mass flow ratio (ALR) at 0.35 MPa inj P =

Table 3 .
Operating parameters and liquid physical properties of baseline case.

Figure 3 .
Figure 3.Comparison of experimental data and model predictions for the Sauter mean diameter (SMD) versus air-liquid mass flow ratio (ALR) at P inj = 0.35 MPa, D noz = 1 mm.

Figure 4 .
Figure 4. Droplet Sauter mean diameter along with the dimensionless axial distance for water and polymer solution.

Figure 5 .
Figure 5. Occurrence number of droplets breakup and collision.

Figure 4 .
Figure 4. Droplet Sauter mean diameter along with the dimensionless axial distance for water and polymer solution.

Figure 4 .
Figure 4. Droplet Sauter mean diameter along with the dimensionless axial distance for water and polymer solution.

Figure 5 .
Figure 5. Occurrence number of droplets breakup and collision.

Figure 5 .
Figure 5. Occurrence number of droplets breakup and collision.

Figure 6 .
Figure 6.Evolutions of Weber and Ohnesorge number along with axial direction for water and polymer solution.

Figure 6 .
Figure 6.Evolutions of Weber and Ohnesorge number along with axial direction for water and polymer solution.

Figure 7 .
Figure 7. Spatial distribution of the drop number and size at the axial distance of 50 mm.(a) Number distribution of water; (b) Size distribution of water ; (c) Number distribution of polymer solution; (d) Size distribution of polymer solution.
/m × is almost identical to the

Figure 7 .
Figure 7. Spatial distribution of the drop number and size at the axial distance of 50 mm.(a) Number distribution of water; (b) Size distribution of water ; (c) Number distribution of polymer solution; (d) Size distribution of polymer solution.

Figure 8 .
Figure 8. Droplet SMD versus ALR for various liquid properties.

Figure 9 .
Figure 9. Droplet SMD along the axial distance for various liquid properties.

Figure 8 .
Figure 8. Droplet SMD versus ALR for various liquid properties.

Figure 8 .
Figure 8. Droplet SMD versus ALR for various liquid properties.

Figure 9 .
Figure 9. Droplet SMD along the axial distance for various liquid properties.

Figure 9 .
Figure 9. Droplet SMD along the axial distance for various liquid properties.