A Computational Fluid Dynamics-based Correlation to Predict Droplet Sauter Mean Diameter for Σ − Y liq Atomization Model in Pressure Swirl Atomizer

DOI: http://dx.doi.org/10.24018/ejers.2021.6.7.2645 Vol 6 | Issue 7 | November 2021 69 Abstract — Liquid atomization is crucial to ensure efficient combustion as it is an inherent part of the injector system. The combustion of fuels relies on effective atomization to increase the surface area of the fuel and consequently achieve high rates of mixing and evaporation. Pressure swirl atomizers are inexpensive and reliable type of atomizer for fuel injection owing to its superior atomization characteristics and relatively simple geometry. The Sauter mean diameter (SMD) of atomizer contributes significantly to the combustion chamber performance. This paper presents a two-step strategy to predict droplet SMD for Ʃ − Yliq atomisation model in pressure swirl atomizer through the combination of experimentally validated Computation Fluid Dynamics (CFD) and Optimal Latin Hypercubes (OLHC) Design of Experiments (DoE) techniques. A three-dimensional Eulerian two-phase CFD model is developed to account for liquid and gas phases as a single continuum with high-density variation at large Reynolds and Weber numbers and validated against experimental measurements, before being employed to carry out a parametric study involving operating conditions and fluid properties of the pressure swirl atomizer. The atomizer is then represented in terms of four design variables, namely liquid viscosity, liquid velocity, surface tension and atomizer exit diameter. An 87-point OLHC DoE is constructed within the design variables space using a permutation genetic algorithm resulting in an accurate SMD prediction. Results show the newly developed SMD prediction is found to be superior compared with existing correlations and indicate significant improvement in the droplets SMD.


I. INTRODUCTION
Sauter mean diameter (e.g., SMD or D32) is one of the main atomizer characteristics and has been widely discussed as it plays an important role in internal combustion engines, gas turbines, furnaces, boilers, and agricultural sprays. Prediction of optimum SMD is not only essential in the design stage of these applications but also influences its performance [1], [2]. The focus is to keep SMD as low as possible due to higher efficiency and lower production of pollutions in combustion. Pressure swirl atomizer performs with various operating conditions, fluid properties, and geometric parameters in order to assess the best droplet SMD for a given application [3], [4]. This leads to a lot of analytical, empirical, semiempirical, and numerical SMD correlations being developed [5]. Fundamentally, the simplest correlation was formulated in the general form as [6]: where a, b, c, d are experimental constants determined from fitting data and , , 1 ∆ 1 are surface tension, kinematic viscosity, mass flow rate, and injection pressure of the atomizer, respectively. The first correlation for swirl atomizer was proposed by Ohnesorge [7] in the mid-thirties and takes the form where SMD is only related to injection pressure and mass flow. Radcliffe [8] and later on Jasuja [9] established further corrections namely:  (4) respectively which are similar and included liquid properties such as viscosity and surface tension. Simmons [10] correlated the Weber number with inlet pressure to define the following equation: for six swirl atomizers of small flow rates using water and kerosene. Correlation (5) where 0 orifice area, 1 inlet hole area and swirl chamber area. The effects of geometry were also proposed by Ballester et al. [12] by fitting experimental data, resulting with an expression of SMD as: = 0.43 0.55 −0.74 −0.05 −0.24 (7) It is noted that the above correlation equations (2)-(7) were obtained without theoretical considerations [5]. The mass flow rate, surface tension, and inlet pressure under the condition of the Weber number being larger than 10 were the inputs used by Kennedy [13] to derive a SMD correlation for six different atomizers with 25 types of liquids. His expression was found to be independent of liquid viscosity. In more convenient correlation by Lefebvre [14]:  (8) which included air density in addition to a range of surface tension from 27 to 73 mN/m, viscosity from 1 mm 2 /s to 18 mm 2 /s, mass flow rate and injection pressure in his correlation. He used liquids such as water, kerosene, gas oil, residual fuel oil, and diesel oil. It should be indicated that due to a disparity in the enumerated formulas above it is not possible to define one universal correlation based on the general form. Wang et al. [15], [16] even assert that within the same limits of fluid properties and nozzle operating conditions is not possible to define one universal SMD correlation based on this formulation.
Furthermore, in 1987 Wang and Lefebvre [16] defined the most complex semi-empirical correlation which was based on the theory. = 1 + 2 where 1 represents the first stage in the atomization process and 2 represents the final stage of atomization. This principle was also used by Wei and Yong [5] in 2014 to derive improved semi empirical correlation to predict Sauter Mean Diameter for pressureswirl atomizers Theoretical formula based on a hypothesis regarding the thickness of a plane disintegrating liquid sheet by Couto et al [17] in 1997 and define: They correlated liquid properties with atomizer dimensions and properties of the surrounding air. Their results were compared with other empirical correlations in their work. Further investigation was made by Musemic [18] in 2011 by testing 20 atomizers with water-glycerol mixture with a dynamic viscosity of up to 5 mPas and came out with the relation: He observed that a higher swirl ratio leads to smaller drop size and conclude that a higher swirl ratio means a higher ratio of the tangential velocity components to the axial velocity [19], [20]. Recent researchers in an attempt to have a comprehensive prediction of SMD used CFD tool with atomization model such as Stochastic Secondary Droplet (SSD) model to develop Sauter Mean Diameter (SMD correlations Dikshit et al. [21]. These results are often compared with empirical equations. However, the Σ− Yliq atomization model with the design of experiments has not been used in developing the SMD correlation. Chen et al. [22] in 2015 use historical data design under response surface methodology (RSM) to formulate the correlation of the SMD; however, it was for air-blast type atomizer and modified power-law relationship for SMD and proposed this equation to be:

II. GOVERNING EQUATION
An entirely Eulerian approach proposed by Vallet et al. [23] treats a two-phase medium as a single continuum where the dense phase is described similarly to a species in a multicomponent reactive mixture.
Let ̃ be the liquid mass fraction per unit mass of the twophase medium, then conservation equation for liquid mass fraction [24]: where ̅ is the Reynold-average density, ̃ is the Favreaveraged velocity of both phases. The mean density ̅ is related to the Favre averaged liquid mass fraction ̃ by: where and are the constant liquid and gas densities respectively. It is assumed that the pressure acting upon both phases is equal. ̇ is the mean rate of vaporization per unit surface of the liquid Ʃ is the mean surface area of the gasliquid interface per unit of two-surface media. Dispersion of the liquid by the turbulence is expressed in the equation by using the turbulent diffusivity and the turbulent Schmidt number as constant = 0.7.
Let Ʃ be the average surface area of the liquid-gas interface per unit mass of two-phase medium. The transport equation for the Ʃ can be written as [24]: where is the turbulent diffusivity, Ʃ is turbulent Schmidt number and is a constant Ʃ = = 0.7, is the rate of surface production and is proportional to turbulence time scale given by: Ʃ is equilibrium interface area and is related to equilibrium drop size [23], [24] by: where is a constant ŋ is the surface tension of the liquid. The atomization model (12) and (14) require a turbulent diffusivity and an integral scale and the standard kturbulence model was used to calculate these variables as well as providing closure of the fluid dynamics transport equations [25].
Once Ʃ is calculated, the Sauter mean diameter (SMD) can be found as: III. NUMERICAL METHOD The design cases were generated using a statistical Design of Experiments (DoE) technique known as Latin Hypercube Designs (LHD) in Matlab. This resulted in eighty-seven (87) design cases for seventy-three (73) different nozzle exit diameters. Felipe et al. [26] indicate that for four (4) independent variables for medium designs using LHD a minimum of seventy (70) cases are required. The design variables (DVs) used are liquid dynamic viscosity from 0.00031 to 0.2 Pa.s, surface tension, 0.02 to 0.075 mN/m, nozzle exit orifice diameter 0.0015 to 0.0035 m and liquid velocity from 1 to 6 m/s as shown in Table I. A commercial Computational Fluid Dynamics (CFD) software STAR-CD code was used to simulate the two-phase flow for the eightyseven (87) design cases for seventy-three (73) different nozzle exit diameters. The standard k-ε turbulence model was used to capture the turbulent behavior in the flow and the SMD was obtained for each designed case for the model. To establish the relationship between these operating conditions and fluid properties and the resulting droplet Sauter Mean Diameter (SMD) to obtain usable SMD correlation for the model, surface fitting of the data was carried out. The result was compared to two existing empirical SMD correlations such as Jasuja and Radcliffe. The where is nozzle exit diameter from the DoE, the inlet area from the geometry of the nozzle and , the diameter of the swirl chamber from the nozzle geometry.
Therefore, SMD correlation based on 4-factor DoE is formulated as shown below: =̇∆ where , a, b, c, and d are experimental constants determined from fitting the data. As indicated, the data for dependent variable SMD and independent variables , , ̇ and ∆ were generated from the DoE points. Using Matlab Isqcurvefit code in the optimization toolbox gives = 2.025, = 0.122, = 0.078, = 0.06 and = −0.42 with the normal of the residuals (resnorm) of the data representing the measure of the goodness of fit as 9.7261 which is less than the residual norm at the bad ending point. Hence, the fit is extremely good considering the variations in the liquid properties and operating conditions. Therefore, the SMD correlation for the model is: To generate the surface plots, two properties are held at constant levels to estimate the approximate SMDs or response variables for the correlations and also observe the relationship between these two predictor variables and SMD based on the new and empirical SMD correlations. In this case, the possible SMDs are known from any two properties selected. With mesh grid surface plot code in Matlab, the two continuous variables selected and the SMD correlations were defined with built-in functions and the response surface generated.

IV. RESULTS AND DISCUSSION
When the experimental constants in the new correlation were compared to the two existing empirical SMD correlations by Jasuja and Radcliffe, the result shows that the model under predicts all the constants except for the injection pressure with reference to the Jasuja SMD. However, over predicted the injection pressure in the Radcliffe case. The slight disparities in the experimental constants predicted may be attributed to the range of values used for the kinematic viscosity and the difficulty in adequately capturing surface tension. To further analyse the new SMD correlation, surface plots for various combinations of predictor variables against the SMD withhold values were carried out as shown and discussed below. Fig. 1 and 2 show the 3D surface plot of SMD against liquid properties with constant operating conditions for existing and new SMD correlations. The kinematic viscosity ranges from 1 mm 2 /s to 1000 mm 2 /s and the continuous variables for surface tension are between 20 mNm -1 and 75 mNm -1 to cover wide range liquids properties. The results show that the response surface or SMDs of the new model equation exhibit similar trend and shape with the existing correlations and the SMDs were smaller than those for the two existing correlations for the same conditions and constraints. This may be due to the smaller exponents or dependency of surface tension and the mass flow rate in the new SMD model which means the model could not adequately account for the surface tension.
The response surface shown in Fig. 2 is curved because the new SMD correlation contains quadratic terms that are statistically significant. The highest SMD approximately 110 µm for wrinkle resistance of the surface is seen in the upper left corner of the plot and corresponds with the combination of highest values of 1000 mm 2 /s and 75 mNm -1 for kinematic viscosity and surface tension respectively. The lowest values of SMDs for the surface are observed in the lower-left corner of the plot and this corresponds with low values of both kinematic viscosity and surface tension 1 mm 2 /s and 20 mNm -1 , respectively. The third and fourth predictor variables, mass flow rate and pressure differential, are not displayed in the surface plot since they are hold values. These results were obtained for the hold values for mass flow rate of 0.0074 kg/s with its experimental-power constant and pressure at 0.93MPa. The significance of this result is that liquids with lower fluid properties will give small droplet SMDs with other conditions remaining unchanged for a particular pressure swirl nozzle.  The continuous variables for pressure are generated between 0.5 and 10 bars based on the operating conditions of pressure swirl nozzles indicated by PNR Ltd in the user manual and the kinematic viscosity lies within 1 mm 2 /s and 1000 mm 2 /s. The average constant surface tension value is taken as 47.5 mNm -1 and the hold value for the mass flow rate is 0.0074 kg/s. From the three plots similar profiles were observed and the response surface for these predictors for the new fitted model was less than those in the existing SMD correlations. The results also show the peak SMD for this combination and conditions is in the upper right corner of the plot and occurs at the lower value of pressure at 0.5 bar and the highest value of 1000 mm 2 /s for the kinematic viscosity. With these results on these conditions, the local minimum SMDs are expected to emanate from pressure swirl nozzle with liquids with higher viscosity operating under lower pressure conditions. Fig. 5 and 6 demonstrate the surface plot of SMD against kinematic viscosity and mass flow rate with constant pressure operating conditions and surface tension for the new and existing SMD correlations. The minimum and maximum values for mass rate flow are 0.0017 kg/s and 0.1667 kg/s respectively and were taken from the user manual for nozzles by PNR Ltd at a constant pressure of 3 bars. The holding value for the liquid property is assumed as 47.5 mNm -1 for surface tension and the boundary space for kinematic viscosity ranges from 1 mm 2 /s to 1000 mm 2 /s. The plots show how the response variable, the SMD, relates to two of the four predictors or components based on the new fitted model SMD correlation. As indicated, the surface plot can only show three variables at a time, while holding the other two components at a constant level. The SMD observed on the plots are only valid for these fixed levels of these input variables. The plots show the same general trends and profiles for the existing and new SMD correlations and maximum droplet sizes are expected from a nozzle of higher flow rate with denser liquids based on these input conditions.

V. CONCLUSION
A usable SMD correlation for Ʃ − atomisation model was established which was consistent with existing correlations. It accommodated a wide range of fluid properties and nozzle operating conditions and easily computed the SMD design purposes. This would also serve as a reference or benchmark for researchers working to improve on this model within the fluid properties, operating conditions, and nozzle geometrical regimes considered. When the surface plots of various combinations of the predictor variables against SMD with typical values were performed, it was observed that the response variables for the new SMD model were smaller than the two existing SMD correlations which were desirable