Spray impingement modelling: Evaluation of the dissipative energy loss and in ﬂ uence of an enhanced near-wall treatment

The goal of the present research is to contribute to improve the knowledge about the spray impingement topic through a numerical study that is aimed at investigating the impact of using the dissipative energy terms that are available in the literature when they are embedded into a speci ﬁ c dispersion model. Comparing all the numerical approaches, a non-negligible disagreement is observed between the relationship proposed in the original model and the other ones drawn from the literature. This fact evidences the in ﬂ uence of the energy dissipated term on the secondary atomization outcome. The present work also provides a comprehensive study on the estimation of the energy dissipated during the splash event. This is a major contribution to the permanent literature since the few works found only addressed the spread regime. In addition, this paper gives an in-depth analysis on the in ﬂ uence that an enhanced treatment of the boundary layer in the region close to the wall may have in the simulation of such ﬂ ows. The work revealed that this near-wall droplets tracking method provides an alternative way to increase the accuracy of the dispersed phase and achieve more consistent results without the necessity of a direct mesh re ﬁ nement.


Introduction
Ever since the pioneering studies of Worthington [1,2] the impact of droplets on solid surfaces have captured the interest of investigators in numerous fields due to the host of applications where this phenomenon can be found.Therefore, a major scientific effort has been invested in the near decade in a comprehensive understanding of the mechanisms underlying the spray impingement process whether through experimental, numerical or theoretical analysis.
A number of semi-empirical models for use in numerical simulations have emerged in the literature with the purpose of describing the interaction between the impinging droplets and the wall.The first known attempt was carried out by Nabber and Reitz [3] employing the KIVA code [4].However, this model had some important limitations related to the conditions for the occurrence of each impingement regime which were specified using an analogy to an oblique liquid jet.Later, Senda et al. [5] developed an impingement model to predict the secondary atomization and liquid film formation resulting from the impact of incident droplets on a surface, as well as the heat transfer between the wall film and the heated wall.This work was mainly based on the experimental investigation on hot walls by Wachters and Westerling [6] and the predicted results, although close to those obtained experimentally, only one particular case of injection and impinging conditions was considered.Other models ensued, such as the one of Nagaoka et al. [7] for 3-D simulations of spray-valve interaction and transient film behaviour but they were only expected to be applicable to a wall whose temperature was above the fuel boiling point since the same experimental data were used.It was only in 1995 that a model to predict the outcomes of diesel spray droplets impacting on a wall with temperatures below the fuel boiling point appeared.The numerical simulation proposed by Bai and Gosman [8] was formulated using a combination of simple theoretical analysis and experimental data from a wide variety of sources and could now predict the splash effect, as well as the droplet dispersion more effectively.The work was enhanced later through the refinement of the quantitative criteria defining the boundaries between different regimes but also by improving the treatment of the post-impingement characteristics in order to allow the application of the model to a wider range of impingement conditions (Bai et al. [9]).The simulations produced satisfactory results for the test cases selected by the authors but no clear evidence of their general applicability could be inferred.Moreover, the splash event began to be regarded with great curiosity due to its relevance in the spray impingement process.In fact, recently [10][11][12], significant attention has been given to this regime either through defining transition criteria that better fit specific experimental configurations or by characterizing the behaviour of drops during all the stages of the regime (expansion of the lamella, crown formation and propagation, etc.).Related to this matter is the dissipative energy loss, which is a fundamental parameter to estimate the post-impingement characteristics or, more precisely, to evaluate the velocity of the secondary droplets.However, there is little literature available related to this particular parameter, as well as little agreement even for what it represents exactly in this regime.The majority of the dissipative energy loss relationships have been deduced for the spread regime, i.e., during the expansion of the lamella and until it reaches its maximum extentwithout splashing.In numerical simulations, this situation can be overcome through some simplifying assumptions but which obviously implies inaccuracy.
It is the goal of the present research to contribute to improve the knowledge about the post-impingement characteristics of the secondary droplets through a numerical study that is aimed at investigating the impact of using the dissipative energy terms that are available in the literature when they are embedded into a specific dispersion model for a simulation of practical interest.In fact, unlike in spreading, there is no general consensus about what exactly represents the energy dissipation in the splash regime.This critical quantity corresponds to an additional parcel that allows energy conservation between kinetic and surface energy of the incident and secondary droplets.However, most of the studies carried out until now that focus on the analysis and estimation of this parameter only considered the energy lossthrough viscous dissipationduring the spreading of the lamella, until it reaches its maximum extent.On the other hand, in the splash regime the energy dissipation may also be associated with the crown emergence and the detachment of secondary droplets from the lamella rim.Besides the original energy dissipation term proposed in the dispersion model used in this work (which is based on the model of Bai et al. [9]), also the equations for the energy dissipation based on the investigations of Chandra and Avedisian [13] and Pasandideh-Fard et al. [14] have been introduced into the simulation.The latter two equations have been adjusted and completed with theoretical principles and experimental data in order to allow their use in the splash event.The prediction results are tested against the experimental measurements of Arcoumanis et al. [15] for spray impact under two crossflow velocities (5 and 15 m/s).
This work follows from a series of previous studies seeking to refine a flexible dispersion model in some aspects that would allow converging towards the best computational solution with minimal time constraints and using modelling assumptions from the available literature.Therefore, beyond the evaluation of the energy dissipation terms, the transition criterion between the spread and splash regimes proposed by Cossali et al. [16], which was found to present the best results in the prediction of the size distributions of the upward moving droplets [17], has been used during the simulations.It is noteworthy that the first part of this paper represents an extension of a previous work [18] in which the same set of energy dissipation equations combined with the spread/splash transition criterion of Bai et al. [9] have been evaluated.
The accurate prediction of dilute two-phase flows requires the effective modelling of both continuous gas and dispersed phases.The complexity of this process lies exactly in the simultaneous existence of both phases.The present multiphase computational model incorporates the influence of the particle dispersion into the continuous phase by treating droplets as sources of mass, momentum and energy in the gas flow field.The interaction between the impinging drops and the crossflow is of great importance in the outcome resulting from the impact.Consequently, the treatment of the boundary layer has been one of the topics that has raised numerous questions in the scientific community.The velocity profile approximation and the refinement required for the effective capture of droplet trajectory in the near-wall zone are key aspects to take into account during the modelling of such flows.Hence, in the second part of this paper, the influence in the final results of a refined treatment of droplets tracking in the region very close to the wall is also analysed.

Mathematical model
The numerical results presented in this work are based on a hybrid approach, in which the fluid phase is treated with an Eulerian reference frame and the particle phase is treated with a Lagrangian reference frame.Only the main features are outlined in this section since a thorough description of the mathematical approach is available in Refs.[19] and [20].The fluid phase is treated as a continuum by solving the partial differential equations in a fixed reference frame which represents each particle and its properties of interest.The turbulence is modelled by mean of the well-known "k-ɛ" turbulence model, which was found to predict reasonably well the mean flow (e.g.Ref. [21]) and the QUICK scheme of Leonard [22] is used to evaluate the convection terms in the discretization process.On the other hand, the dispersed phase is treated using a Lagrangian reference frame in which the particle trajectories are obtained by solving the particle momentum equation through the Eulerian fluid velocity field.
Estimated from E D ¼ ∫ tc 0 ∫ V ϕdVdt≈ϕVtc where the dissipation function is given by: ϕ Improvement of the Chandra and Avedisian [13] theoretical model, by replacing h by δ bl ¼ 2d I = ffiffiffiffiffi ffi Re p and assuming t c = 8/3(d I /U I )

Table 2
The four relationships tested in this study.
Model A (Bai et al. [9]) The flow configuration is shown schematically in Fig. 1 and consists of a spray stream injected through the upper wall of a rectangular channel with an inclination of 20°(in relation to the vertical plane), in the direction of the stream flow.The cross section of the computational domain is 86 mm× 32 mm, whilst the channel length is 350 mm.The location of the injection point (Z in ) lies 50 mm downstream from the inlet plane (Z in /H = 1.563) and on the symmetry plane.Therefore, the six boundaries of the computational domain considered here are: an inlet and outlet plane, a plane of symmetry and three solid walls at the top, bottom and side of the channel.
Numerous assumptions must be made in order to allow the modelling of the droplet behaviour within an adequate time frame.In the dispersed phase, it is assumed that the particles are sufficiently dispersed so that the interaction between droplets is negligible.In addition, since no reliable atomization model is yet accessible, an empirical procedure is used to estimate the characteristics of the spray at the exit of the injector.However, due to the limited experimental data available, it is assumed that droplets are spherical, the spray is dispersed so that the droplet collision/coalescence can be neglected, and the droplet aerodynamic breakup and evaporation are ignored (the measurements were performed at ambient pressure and temperature).Thus, the measured probability density function (pdf), which was taken in a particular plane 80 mm downstream the nozzle, is then reproduced at the injector exit.
The model of Bai et al. [9] considers four impingement regimes: stick, rebound, spread and splash.The existence of these regimes depends on the properties of the impinging droplets and the conditions of the solid surface, including whether the latter is dry or wetted.Under both conditions, the spread-splash regime transition criteria were derived from the Stow and Hadfield [23] data, giving rise to a Critical Weber number dependent on the Laplace number.The boundaries stickrebound and rebound-spread (derived from Lee and Hanratty [24] data) for wetted walls were set with the critical Weber numbers of 2 and 20, respectively.
In a previous work [17], a qualitative comparison of several transition criteria between spread and splash regimes available in the literature has been performed.It was found that, for the particular case of the upward moving droplet size distributions, the transition criterion proposed by Cossali et al. [16] revealed a better agreement with measurements than the other correlations evaluated.However, the opposite was verified in the case of the velocity-size correlations of the upward moving droplets, where the velocity profiles were overestimated along the entire range of droplet diameters.In Ref. [16], the authors investigated the transition between spread and splash regimes by analysing a large number of pictures of droplets impacting on an aluminium plate with a non-dimensional film thickness (δ) between 0.08 and 1.2.They found that a correlation based on the Weber and Ohnesorge numbers could fit adequately their data: Therefore, Eq. ( 1) has been used instead of the one proposed in the original model for the boundary between the spread and splash regimes.The present dispersion model solves the conservation of energy equation for an impinging particle that can produce up to 6 splashing parcels in which the size follows characteristic distributions for secondary droplets.The kinetic and surface energy of the incident drop are considered to be equal to the kinetic and surface energy of the secondary droplets plus a term that corresponds to the dissipative energy loss.This latter parameter is perhaps the most critical quantity to determine accurately and is the main focus of the first part of this work.
In the original model, Bai and Gosman [8] deduced their own relationship for the dissipated energy in terms of the critical Weber number.Due to the under-estimation of the values in certain ranges of application, the relationship was later revised and limited by a postulated value of 80% of the kinetic incident energy based on the normal incident velocity of the droplets.However, the origin of the approach used to define the energy dissipation is somewhat unknown.
The principal study available related with the dissipative energy loss has been conducted by Chandra and Avedisian [13], which concluded that the dissipation of energy is directly proportional to the viscosity.In fact, this study was focused on estimating the maximum diameter of a liquid drop which spreads on a surface.The authors found that the viscosity largely controls the post-impact occurrences, which in turn rules the splash phenomena.The time scalecharacteristic of convectionis estimated by assuming it as the time taken for the droplet thickness, h, to go from its maximum value to zero, whereas the volume of liquid in the drop, after flattened out in a regular disc shape, is proportional to the thickness of the disc and the square number of the maximum extension diameter of the liquid film.
Since then, several investigators took the assumptions of Chandra and Avedisian [13] to estimate the energy dissipated and tried to improve the accuracy of the theoretical model.In that pursuit, Pasandideh-Fard et al. [14] replaced the splat film thickness in the dissipation function ϕ, which they found that it overestimated the maximum extend of the film d max (up to 40%), by the boundary layer thickness δ bl at the solid-liquid interface.In addition, the authors assumed a new time scale t c associated with the droplet spreading until its maximum extension is reached depending on the impact velocity of the incident droplet.
The dissipative energy loss relationships deduced by the investigators and their source equations are outlined in Table 1.
The thickness (h) of the liquid drop after flattened out in a regular disc shape when impact occurs, which is found in the relationship presented by Chandra and Avedisian [13], is estimated considering that the incident drop keeps the same volume from the instant immediately before impact until the last instant of the spreading phase.In addition, it is assumed that the splash regime only occurs at the moment of the crown emergence.Thus, from the work of Yarin and Weiss [25], d max is defined as being twice the incident drop diameter.Therefore, the equation defining the thickness of the disc may be written as follows: However, Roisman et al. [26] estimated the following disc thickness from a mass balance without giving much details about it: According to Pasandideh-Fard [14], Eq. (3) should not be used because it does not correspond to the equality of the volume of a drop of diameter d I to a cylinder of height h and diameter d max .However, Roisman et al. [26] did not clearly mention that this was the kind of mass balance used.So, admitting the previous value for maximum spreading diameter, the thickness of the disc becomes: Both Eqs. ( 2) and (4) for the disc thickness have been introduced into the relationship deduced by Chandra and Avedisian [13] (called Model C and Model B, respectively) which in addition to the correlations of Bai et al. [9] (Model A) and Pasandideh-Fard et al. [14] (Model D) make up the four energy dissipation relationships tested in this paper as highlighted in Table 2.
It is noteworthy that Models B and C are related by a factor that is the incident droplet diameter raised to the third power.As it will be seen later, this difference in both modelsand specifically this parameterdoes not have much importance on the results obtained.

Results and discussion
This chapter is divided into two sections: the first one addresses the study of the energy dissipation for two crossflow velocities; and the second analyses the influence of an enhanced droplets tracking in the near-wall region.The numerical predictions of the model, which computational method has been described in the previous section, are presented and discussed.Only the most meaningful results are presented here in order to ensure brevity and conciseness.

Energy dissipation
The numerical predictions presented in this section are compared and assessed with the experimental data of Arcoumanis et al. [15] for two cases of crossflow velocities: 5 and 15 m/s.Two different populations of droplets are distinguished from the results presented below: one set consists of droplets moving towards the wall whilst the other comprises the droplets moving in the opposite direction.In addition, two types of graphs are presented: the probability density function of droplet size; and the droplet mean velocities as function of their size.Fig. 2 shows the four different measurement locations, where the present results are compared with the experimental data of Arcoumanis et al. [15].The four different positions a, b, c, and d, are located in the centre of the wind tunnel and lie respectively at 12, 15, 20, and 25 mm downstream of the injector nozzle and within a horizontal plane 5 mm above the impingement wall.
Fig. 3 shows droplet mean velocities in the direction normal to the wall.A good consistency is observed between the four relationships studied but the comparison between the computational results and the experimental data shows that the former are still far from the results expected.It can be seen that at locations c and d, the maximum predicted size classes do not exceed 160 μm and 100 μm, respectively, whereas the measurements reach 400 μm at both locations.In fact, these discrepancies were also reported in Ref. [9] and may be connected to the uncertainties in the derived initial conditions for the construction of the spray, which produce erroneous incoming droplet estimations.
Fig. 4 illustrates the size distributions of the upward moving droplets through an air flow rate of 15 m/s.It is evident that the simulated results show close agreement with the measurements at all the locations.This agreement between both results is also enhanced by the use of the transition criterion proposed by Cossali et al. [16] as concluded in Ref. [17].The greater discrepancy is found at the location closer to the injector where a slight shift of the distribution to lower drop sizes is observed.This implies that the simulation exhibits a higher probability of occurrence of smaller droplets in that position.
The energy dissipation relationships introduced into the dispersion model affects the post-impingement characteristics of secondary droplets by balancing the energy conservation equation with a new value for the velocity.Hence, one may not expect that the general behaviour of the particles in their way down towards the wall, which are mainly the ones coming directly from the injector in their incident trajectory, can be effectively influenced by changing the value of the energy dissipation parameter.On the other hand, the same statement does not apply to the velocity-size correlations of the upwardmoving droplets, as will be discussed later in this analysis.
The droplet mean velocities in the direction normal to the wall for the ones moving upward are shown in Fig. 5.As it can be seen, the energy dissipation relationship established in Model A is the one that present the best results as far as the velocity profile is concerned for an air flow rate of 5 m/s.Despite that, it still results in larger velocity values relatively to experimental data, which in this case is largely due to the use of the transition criterion of Cossali et al. [16] as concluded in Ref. [17], and does not predict any droplets with size classes larger than about 180 μm, whilst the measurements cover a range up to 300 μm.The other newly defined models show reasonable agreement between them but deviate even further from the experimental data.Considering now the results with the 15 m/s crossflow velocity in Fig. 6, it is observed that Model D approaches the results of Model A, which lie close to those obtained experimentally.However, once more, the other relationships tested overestimate the normal velocity magnitude of droplets along the entire range of size classes.

Near-wall droplet tracking method
The dispersion of droplets is associated with the interaction between the droplets and the surrounding environment.The presence of a gas phase through a crossflow with a certain velocity enables the appearance of a boundary layer that influences the trajectory of droplets in the region near the impinging surface.In fact, the difference between the air velocity in contact with the top and bottom of the particle causes it to rotate [27], which gives rise to an additional significant force that leads to a different trajectory.Some authors [28] reported that in the near-wall region, where the shear rate is higher, the magnitude of this transverse lift force becomes as important as the drag force.Furthermore, upon entering the boundary layer, additional deforming forces are exerted on the drop which could decelerate and deform it slightly before impact.The modelling of such flows requires the accurate description of the boundary layer in the region near the wall, whether by approximating the velocity profile by a linear function or not.This decision is normally by the construction of the mesh which, although independent of numerical errors, may lead to an erroneous prediction of the droplet behaviour in this region.Therefore, the key question that remains is: how the refinement of the treatment of the droplet tracking in the region near the wall affects the results?
The comparison between the dimensionless vertical profiles calculated and the one measured by Arcoumanis et al. [15] for an air flow velocity of 5 m/s is presented in Fig. 7.The velocity profile has been taken in the vertical plane at X/H = 0.05 and upstream of the injector position (Z/H = 1.3).The horizontal velocity component and the vertical length scale have been made dimensionless by the maximum crossflow velocity and the height of the solution domain, respectively.Fig. 7 reveals that the numerical results are in acceptable agreement with the Laser Doppler measurements of the air flow.The critical discrepancy occurs in the last points of the experimental data in which the last numerical node lies at Y/H = 0.05 and thence to the wall the velocity profile is approximated using a linear function.A direct comparison is made between the velocity profiles in the boundary layer given by the turbulent theory, experimental data and numerical predictions, indicating an adequate agreement that suggested the use of the linear approximation for the sake of mathematical simplicity.Therefore, the procedures implemented are centred on the basis of a more refined tracking of droplets in the near wall region, through a decrease of the time step size which would lead to an increase not only in the accuracy but also in the computational cost.The process starts from the second node above the wall with a refinement established of up to five times the original model in each iteration all the way down until the drop reaches the surface.
The calculated and measured size distributions for the downwardmoving droplets with an air flow velocity of 15 m/s are depicted in Fig. 8.An improved performance of the modified model is evident particularly at locations c and d, where the old model over-estimates the mode (most frequent) of the predicted droplet diameters.In fact, the droplets tracked at the locations further away from the injection planewhich have larger angle of incidence and thus enabling them to go furtherseem to be more vulnerable to the modified model possibly due to the longer time period that they spend within the influence of the boundary layer.The enhanced treatment near the wall allows a more accurate prediction of the distribution of the droplet sizes in the instants before (the large majority of them) the impact in the case where the effect of the crossflow is more important.Also noticeable at location d is that the droplets with sizes greater than 120 μm are not predicted although they were detected in the measurements.Figs. 9 and 10 depict the outgoing droplet mean velocities in the direction normal to the wall considering the air flow velocity of 5 and 15 m/s, respectively.In the case of the lower crossflow velocity, the discrepancies found in the original simulation for droplet with size classes higher than about 200 μm are stabilized in the improved model.However, range of size classes covered by the numerical simulation is not as high as in the experimental data.The opposite is verified in Fig. 10, in which the maximum size classes predicted in the modified model are much higher than both measured and original model in all the locations.Nevertheless, in general, the behaviour of numerical results is in reasonable agreement with the measurements.

Conclusion
Simulating the complex dynamic phenomena that occur during the impact of the spray droplets is of great interest in numerous applications.In this sense, refine flexible dispersion models that can be adjusted through the use of adapted and more suitable empirical correlations in order to fit specific configurations with minimal time constraints would be a great asset both for industry and researchers.In this work, the dissipative energy terms for the splash event found in the literature are investigated and compared.Besides the relationship proposed by Bai et al. [9], other equations of energy dissipation based on investigations devoted to the study of this parameter are introduced in the simulation and the results are tested against experimental data.The results show that the dissipative energy term that is considered in the energy balance between incident and secondary droplets have an important effect on the secondary atomization outcome.This influence is more important in the velocity-size correlations since this parameter has a direct effect on the normal velocity of the droplets resulting from the splash regime.Therefore, the comparison of the four expressions presented in this work allows concluding that the original relationship of the base model is, by now, the best approximation of the energy dissipated during the splash event.In fact, the disagreement between the latter expression and the other ones drawn from the literature evidences the non-negligible effect of this parameter on the results.In addition, these discrepancies between the numerical predictions demonstrate that the expressions used to give rise to Models B, C and D do not represent accurately the energy dissipated during the splash regime since the velocity profile was over-estimated along the entire range of droplet diameters.This situation is connected to the underestimation of the dissipative term, which is consistent to the fact that the expressions were proposed for the spread regime.No other studies were found in the open literature reporting the estimation of the energy loss during both the expansion of the lamella, and crown emergence and detachment of the secondary droplets from the lamella rim.
There are still noticeable differences between the predicted results and the measurements in the velocity-size correlations that are related to the empirical procedure used to estimate the characteristics of the spray at the injector exit due to the lack of a reliable atomization model available.On the other hand, a good agreement has been verified between the calculated and measured size distributions for both droplets moving downward and upward.The latter case is also improved by the use of the transition criterion proposed by Cossali et al. [16] but the opposite is verified in the droplet mean velocity results (as reported in Ref. [17]).Another observation on the energy dissipation study is that the incident droplet diameter parameter does not seem to play an important role on the final outcome.In fact, despite two of the studied relationships differ from each other by a factor of proportionality that is the incident droplet diameter raised to the third power, the results obtained are very similar in all cases and conditions considered.
Concerning the results of the near-wall droplet tracking method, the new approach used in the present work provides an alternative to model more accurately the dispersed phase in this region close to the impinging surface instead of simply increase the mesh refinement, which can lead to difficulty resolving the law of the wall.In fact, this method enables the achievement of more consistent results without the necessity of altering the gas phase.Both the size distributions and the mean velocities of the downward-moving droplets of the modified model exhibited a better agreement with the measurement data than the original simulation, particularly for the higher crossflow velocity, which is when the influence of the boundary layer on the outcome is more important.It is worth mentioning that these improved results are also connected to a slight increase in the computational cost and time but, yet much lower than what we would expect with a direct mesh refinement.Bearing in mind the importance that the shearing air has on droplet dispersion [15,29], this conclusion becomes interesting and might be helpful in specific applications.

Fig. 4 .
Fig. 4. Size distributions of upward-moving droplets at four locations for a crossflow velocity of 15 m/s.

Fig. 5 .
Fig. 5. Velocity-size correlations of the upward-moving droplets at four locations for a crossflow velocity of 5 m/s.

Fig. 6 .Fig. 7 .
Fig. 6.Velocity-size correlations of the upward-moving droplets at four locations for a crossflow velocity of 15 m/s.

Fig. 8 .
Fig. 8. Size distributions of downward-moving droplets at four locations for a crossflow velocity of 15 m/s.

Fig. 9 .
Fig. 9. Velocity-size correlations of the upward-moving droplets at four locations for a crossflow velocity of 5 m/s.

Fig. 10 .
Fig. 10.Velocity-size correlations of the upward-moving droplets at four locations for a crossflow velocity of 15 m/s.

Table 1
Dissipative energy loss relationships and corresponding observations.