Probabilistic derivation of droplet velocity using quadrature method of moments

Nomenclature

D

drop diameter (m)

D n

normalization diameter

k

parameter of gamma distrbiution

M

spray moments

n ( D )

the number size distribution function

n v

the number of droplets per unit volume

r 32

Sauter mean radius

r crit

critical radius

Re

Reynolds number

S

the source term

t

time

u

the droplets velocity

V

the droplet volume

x c

any real number

Greek symbols

a

modeling velocity parameter

b

modeling velocity constant

Γ

gamma function

τ

Lattice relaxation time

1 Introduction

The main difficulty encountered in modeling spray flows, and more generally of two-phase discharges, lies in the jumps of discontinuity experienced by the physical quantities at the passage of the different interfaces. To correctly treat these discontinuities, two levels of modeling can be envisaged. The first way to model a two-phase medium is direct numerical simulation (DNS), which consists of indirectly solving the local instantaneous equations and the different jump conditions of the two-phase medium.

All the flow details, both in space and time, are thus solved numerically. Unfortunately, the counterpart of such precision in describing the two-phase medium is a very high numerical cost. Despite the recent improvement in computing resources and the high degree of maturity attained by current numerical methods, it is still not feasible to simulate complex industrial configurations using the DNS. A second way of modeling a two-phase medium is to describe it in average magnitudes. Using an averaged model makes it possible to reduce the computation time considerably. Although only large scales are solved, the degree of precision of such a model is generally sufficient for many industrial applications. The difficulty is no longer in a numerical order. It is moved to the closures. Many models are indeed needed to retrieve the information lost during the averaging process, especially concerning the geometry of the interfaces and the interactions between phases.

The averaged modeling of two-phase flows has been the subject of many developments over the last 50 years. The most commonly used model is undoubtedly the two-fluid model [1] because of its potential to describe a wide range of two-phase flows, from dispersed to phase-separated flows. The resolution scale chosen for averaged scale via the two-fluid model is the ultimate goal of obtaining a simple model to implement and digitally efficiently all kinds of configurations and complex geometries of spray droplets [2]. To describe a technical spray, polydispersity has to be considered. The two-fluid method can only reduce the polydispersity and the mapping of the polydispersity with the multi-fluid approach [3]. Although this is too computationally intensive, the application of these methods for simulating technical sprays is unimaginable due to the economic feasibility.

Another way to take advantage of the parallelizability of the Euler description use and polydispersity is to consider form-the-moment methods. It describes the evolution of droplet size distribution in space and time. The moment methods use the possibility that distributions are also characterized by their moments. Instead of the drop size distribution and describing and calculating the distribution function descriptively, the distribution moments are considered. Based on the kinetic spray equation, the moment methods were developed. The method of moment (MOM) was one of the first methods [4] used in chemistry.

McGraw extends this quadratic approximation method to the quadrature method of moments (QMOMs) [5,6]. The distribution function is approximated by abscissa and weighting, also called quadrature variables. Using these abscissas and weights, the moments of the distribution function can be reached and transported. Marchisio applies this method [7–9] and put it based on direct quadrature method of moments (DQMOMs) [10]. Instead of transporting the moments as in QMOM and using an algorithm, they reconstructed the moments in each time step; the quadratic approximation variables are directly transported in the DQMOM. Extended areas of application range from simulations of soot formation [11], solution of the Fokker Planck equation [12], and fluidized beds [13] to polydisperse sprays. The latter is from [14–16]. Fox et al. [17] introduced the so-called ratio constraints for DQMOM to represent complete vaporization. The new method performs verification calculations and compares them with results from an Euler-Lagrange analysis. From the works [17–19], it becomes clear that the computational effort for the DQMOM is lower than for the Euler-Lagrange method. Pigou et al. [20] extended QMOM to solve the publishing balance equations. Askari et al. [21] extended the quadratic moment method to simulate the heat and mass transfer in two-phase flow in bubbly flow.

The standard approach involves solving transport equations of moments, where the weighting and abscissa values are calculated using the PD algorithm and where each iteration is performed again. Fox et al. [17] and Desjardins et al. [22] formulated the DQMOM approach based on the direct calculation of the weighting and the abscissa values [23]. The DQMOM contains an accurate description of the distribution function. The reconstruction is based on the method of maximum entropy. The coupling with the multi-fluid approach is comparatively easy due to the realizable direct calculation. The DQMOM approach [15] is also based on the direct transport of weights and abscissa values but also includes a sectional method [24]. The accurate description of vaporizing droplets and convective transport characterizes the sectional formulation of the moment method. The distribution function is determined within the respective section. Since the stability of moment methods depends on the number and the order of the moments in question, an increase in the moments leads to more extensive and sparse matrices, which in turn leads to stability problems in the calculation of the eigenvalues. The division into sections offers the advantage of increasing the number of moments [25]. The purpose of this article is to propose a model for droplet velocity. In the model, locally, only a single volume-averaged velocity has been considered for the droplets. The objective is to derive a droplet velocity model using the spray moments theory. The model employs an assumed gamma distribution through the derivation of droplet velocity. Model test cases are presented in Section 4. The comparisons are made with two sets of experimental data, including diesel and cryogenic liquid nitrogen spray before conclusions are drawn in the paper’s final quarter. The primary purpose of this work is to consider and model the droplet velocity. In the model, locally, only a single volume-averaged velocity has been considered for the droplets. The objective is to derive a droplet velocity model using the spray moments theory. The model employs an assumed gamma distribution through the derivation of droplet velocity. Model test cases are presented in Section 4 and compared with two experimental data sets, including diesel and cryogenic liquid nitrogen spray, before conclusions in the article’s final quarter.

2 Velocity moment (VM) model

When injecting a liquid spray into a gaseous flow (turbulent or otherwise), an exchange of motion occurs between the liquid and gaseous phases, generating a high shear in the flow. This phenomenon leads to the formation of vorticity around the liquid jet, which gradually breaks up into ligaments and droplets. Because it defines the transition from a dense liquid core to ligaments, these two states cannot be compared to spherical droplets as the formalism requires. In addition, to correctly solve the induction air generated by the spray, before the limitation of the calculation cost and the limits of the approach, the modeling of the formation of the diluted spray proves to be essential. The spray atomization processes can generate droplets of different sizes and shapes. Therefore, every single droplet has its velocity. That leads to an accurate prediction for subsequent operations such as a breakup, collision, and evaporation. The proposed model differs from the mono-dispersed model based on assigning a unified velocity for all droplets. In addition, the model differs from the Lagrangian approach, which is based on the graduate gradient for the droplet velocity by taking the effect of the gas phase induction. The VM model is settled by constructing the power relation between the droplet velocity and its diameter.

Consequently, when the droplet diameter is large and proportional to its mass, its velocity will be high compared with less droplet diameter. The prior analysis for the model is the selection of coefficients. The first coefficient is determined from the derivation and construction of the QMOMs. The user sets the second coefficient, but it can be modified during the calculation if the maximum velocity of the profile is greater than the Bernoulli velocity and to be chosen arbitrarily between { 0 , 1 } [26].

The parameter α is obtained by the constructed functional form to express the drop size distribution. D n denotes the normalization diameter. While the parameter β is obtained by comparison with experimental data – the value of the coefficient α , defining the thickness of the moment’s ratio. The model illustrates the evolution of the velocity profile at different times that correspond to varying values of the moment rate. For example, the droplet velocity is not zero at the edge of the surface, and the maximum speed varies (not equal to Bernoulli’s velocity). Therefore, the shape of the velocity profile evolves to ensure the conservation of the momentum.

For various moments, the equations of transport can be determined by multiplying the factor D k by the proposed distribution of drop size, a group of droplets with similar properties as would be solved in a multi-size Eulerian treatment, and integrating the product over the entire range of radius. Thus,

The aforementioned transport equation is the general form for the order of the differential operators, and the integration has been changed based on the physical meaning for the moment. The constructed transport equations for the moments of the number size distribution function n ( D ) with an extra transporting to their moment transport velocities.

and

By using these equations, the moment’s transport equation becomes:

Substituting the assumed relation that expresses u(D) from equation (1) in equation (4), the moment transport velocity takes the form of:

Substituting equation (6) in equation (5), the general moment transport equation becomes:

The source term, S k , interacts from this drop group with the phase of gas. Equation (7) shows that the existence of the moments M ( k + β ) needs closure to be solved, which is produced from the dependence of droplet velocity on its diameter, as explained in equation (1). A simple way to overcome the closure problem is explained after presenting the fourth-moment transport equation in terms of the volume of the droplet.

3 Numerical methods

Discretization of equations transforming these differential equations into a set of algebraic equations using derivative approximations. The methods used during numerical resolutions are finite differences, finite elements, and finite volumes [29]. In this study, the finite volume method is used. It entails solving partial differential equations on volumes surrounding each grid point. This method, even if it only provides first-order precision, has qualities that make it one of the most suitable for studying turbulent flows. Indeed, this one is very robust and makes it possible to treat equations comprising complex and nonlinear source terms.

Moreover, it has the advantage of satisfying mass conservation for each control volume. Finally, it can be used with a relatively coarse grid, which allows the implementation of codes at a reasonable cost. The finite volume method is a discretization method suitable for the numerical simulation of different conservation law equations. This method is widely used in phenomena studied by physicists and engineers, such as fluid mechanics, mass and heat transfer or chemical reactions. Some of the essential features of the finite volume method are similar to those of the finite element method: they can be used on arbitrary geometry, using a structured or unstructured grid, leading to robust designs. Another feature is the local conservative of digital streams, i.e., the stream is conserved from one cell to another [29]. This last characteristic makes the finite volume method very interesting when considering a problem for which flux is essential, such as in fluid mechanics, the simulation of semiconductor devices, mass and heat transfer, etc. The finite volume method is locally conservative because it is based on a “balance” approach; the local balance is written on each discretization cell, called the “control volume.” The divergence formula is then used to obtain an integral formulation of the fluxes on the control volume’s edge.

The in-house code solution is implemented in Gambit software and FORTRAN. The source codes are compiled with gfortran using the double-precision variables with a high optimization flag. The resulting algebraic system of linearized conservation and transport equations is solved with the in-house code SPRAYMOM. In order to solve the linear systems, such as the QMOMs system and the eigenvalue problems the product-difference algorithm methods is used.

An overview of the methods of numerical linear algebra has been used. For example, the CFD–grid represents a subroutine of linear solvers for large-scale scientific problems of parallel and single multi-grid methods with structured and unstructured computing grids and meshes, and supports various common programming languages. SPRAYMOM, in turn, contains a library of routines for calculating sparse matrices. These include iterative solvers, drills for preconditioning the matrices, and exercises for basic matrix operations. SPRAYMOM is written in the programming language FORTRAN. The LAPACK-95 software is also written in the programming language and includes routines for solving single, parallel, and linear systems of equations. The GFORTRAN module was created from f95-lapack for single and double methods. The solution methods of eigenvalue problems which are included the routines for the matrix factorization such as the LU or the QR method are solved by using the Cholesky decomposition for dense and diagonal-dominant matrices. The SPRAYMOM uses the “solve-data” routine based on the implicit QL and QR methods to solve eigenvalue problems to calculate the eigenvalues and eigenvectors of the asymmetric, tridiagonal matrix, where QL is algorithm for symmetric tridiagonal matrices, and QR is eigenvalue algorithm.

4 Results and discussion

For two-dimensional calculations, use the program Gambit–axial symmetric grids generated as shown in Figure 1. To transfer these grids correctly into the in-house CFD code, the axis of rotation must correspond to the X -axis. If the vector is selected correctly for gravity, this does not matter for further calculations. Exceptionally high is the number of grids in the chamber to choose in which spray droplets are expected to be spread out. Transition areas between fine and coarse meshing and wall areas can be further refined (adapted) later with different levels of compression ratio for the X -axis, as otherwise, the related discontinuities can occur at these locations.

Figure 1

The computational domain.

The number of cells can reach a multiple of the original scope. As can be seen, the grid is block structured in some places and unstructured in others. Spray droplets are described by their diameter, physical and chemical composition, velocity, and temperature. The droplet diameter and droplet velocity are the most essential calculated parameters and are the focus here. In the beginning, droplets were separated into groups whose diameters were between D-DD/2 and D+DD/2. The probability density functions (PDFs) of the droplet sizes of the spray droplets can be used here to classify the spray and count the number ( n ) of droplets restrained in each of these groups; then, a histogram is made.

In Figure 2, a comparison of drop size distribution results has been made between the experimental data of [30] and the current work. The red bar represents the practical work, and the blue represents the current simulation. The histogram of droplet sizes counts the number of droplets in each group. The QMOMs provides the capability to construct local droplet sizes and velocity distributions. It can be seen that a good matching between the gamma distribution (droplet size distribution) from the experimental work with the predictions. The main reason is that the gamma distribution is utilized using the Sauter mean radius ( r 32 ). It is mainly considered the essential parameter in the size distribution function, so it would be more appropriate to use the r 32 as the parameter implemented by the gamma distribution function. Figure 3 shows spray tip penetration as a function of start injection time for prediction cases and experimental data using the five grid resolutions defined in Table 1. It is clear from this comparison that the coarse grid (Grid A) importantly underpredicts spray penetration in the later stages of liquid droplet injection and does not acquire the precise behavior of spray penetration as a function of start injection time.

Figure 2

Comparison of drop size distribution results between current work and experimental data of [30] at 0.2 MPa.

Figure 3

Grid analysis test.

As the grids are refined (Grid B), the gradient of the spray tip penetration increases as injection time advanced, analogous to the experimental data behavior. In the other cases (Grid C and Grid D), the grid convergences are achieved with only little changes for further refinement but underpredicts spray penetration. Finally, the spray tip penetration slope is virtually over predicted for more improvement, as in the case Grid E, which is expected because there is severe over-resolution four times than the excellent agreement with the experimental data [31] (number of volumes are 7,799 as compared with 1925) (Figure 4).

Figure 4

Comparison of predicted spray penetration with the experimental data [31].

Spray tip penetration is defined as the distance between the outlet of the injector hole and the end of the spray. Its value is time dependent. Spray tip penetration has been studied for a long time by many researchers. Their experience involved injecting a liquid (from oil) under pressure in a re-watered chamber. The two variables were injection pressure and ambient gas pressure. It was noticed that increasing the injection pressure increased the spray penetration while, conversely, increasing the pressure of the ambient gas decreased it (Table 2).

The contour plots of fuel liquid mass shown in Figure 5 are presented to verify that the dominant liquid group occupies the core of the spray. Notably, most liquid mass also appears at the spray tip due to the droplets moving into the stagnant air being gradually slowed by drag force and caught up by the droplets behind, which experience a smaller drag force. The result of the high droplets velocity at the interface region between stagnant gas and spray droplets can change the droplets trajectory into the transverse direction and then flow backward because of the shear force effect. The newcomers then flow through the droplets at the interface and form the new tip, only to be slowed in turn and overtaken by droplets from behind once more.

Figure 5

Structure of the liquid volume fraction of the spray. (a) Uniform velocity, (b) relative velocity ad (c) velocity moment.

Figure 6 compares the experimental data of [30] and the current work. The solid cone injector atomizer produces the spray cone angle as a function of its injection pressure (operating conditions). The two curves explain that a similar trend is observed for spray cone angle variations in a constant frame of the atomizer. However, little difference in spray cone angle is seen between the predicted and measured experimental data at different ambient pressure cases (0.1–3.0 MPa).

Figure 6

Comparison between calculated and measured [30] initial spray cone angles at different injection pressure.

This comparison aims to observe the spray droplets’ dispersion, which gives a good indicator of the trajectory of spray velocity. Due to the possible atomization processes, primary or secondary, the size of the droplets formed by any injector can be well dispersed. Then, mathematical probability tools are used to characterize spray droplets, and the comparison of sprays is done based on objective quantitative criteria. For all types of fountains, it is possible to establish PDFs. Numerically based on volume, to define, for each droplet diameter, the probability of encountering such a droplet in the entire spray. It can be distinguished by using the numerical distribution. It PDF lists the number of droplets according to their size, and the volume distribution, denoted as n v ( V , x , t ) , represents the volume of liquid occupied by each diameter class. The volume distribution is always shifted to the right concerning the numerical distribution because the small droplets, even if numerous, represent a low volume of liquid. Figure 7 shows the unsteady behavior of the PDF.

Figure 7

Comparison between calculated probability density functions at different times, (a) velocity moment, (b) Relative velocity, and (c) Uniform velocity.

It can be seen that similar profiles are observed for cases of VM and UV. While much less homogeneous values can be seen in the relative velocity (RV) case, this is due to the gas velocity development. The effect of the droplet size-to-droplet velocity on the droplet size distribution is shown in Figure 8. It is observed that with the increase in the droplet size, the droplet velocity increases overall in the spray injection stages. Figure 7(a) shows the case of VM, where the distribution curves are shifted towards the smaller droplet sizes. This is due to the fact that the higher momentum of the dispersion and breakup processes results in better overall atomization and hence much smaller mass-mean diameter for the whole spray. The case of UV is shown in Figure 7(b).

Figure 8

Comparison between droplet velocity at different times, (a) velocity moment, (b) Relative velocity, and (c) Uniform velocity.

The constant velocity is assumed during the calculations, and the velocity profile lines increase during the elapsed time. The case of RV is shown in Figure 7(c). The effect of local gas phase velocity is considered, while the liquid is assumed to be uniformly distributed during the calculations. It should be noticed that the effect of drag interfacial force appears on the droplet velocity profile. At t = 0.15 ms, the values of the droplet velocity profile decreased to the lower range due to the drag force. Figure 9 shows th PDF with droplet velocity at different times. The figure shows the progressive dispersion of the spray droplets predicted by the model. In the VM case of VM at time = 0.01 ms, the lower values of PDF can be seen. That is because the mass fraction on the axis is close to 1, the density ratio is large, and the volume concentration of the corresponding liquid droplets is high. As time advanced, the maximum mass fraction decreased, and the profiles spread radially and approached the highest values for the PDF. In the UV case, the vertical lines refer to spray droplets dispersion. They start from lower PDF with lower droplet velocity and increase gradually. In some parts, the PDF remains at lower values due to not being fully atomized for the bulk of the liquid. In the RV case, the vertical lines refer to the spray droplet’s dispersion with the effect of the gas phase velocity. The difference is clear where the lower values are showing at time = 0.2 ms, where the difference between the liquid and gas velocities is maximum. For the rest, it can be seen that the maximum value of the profiles changes little throughout the dispersion of the spray velocity.

Figure 9

Comparison between calculated probability density function at different times and velocity profiles, (a) velocity moment, (b) Relative velocity, and (c) Uniform velocity.

Figure 10 displays the PDF with droplet diameter at different times for the three cases. At time = 0.01 ms, the RV case offers the maximum profile because, near the injector, the gas velocity is negative due to the induction of the gas phase velocity. However, the gas velocity has a minus sign and would be added to the droplet velocity. At time = 0.1 ms, all cases have approximately the same trend (Figure 11). At time = 0.2 ms, both VM and UV cases record the highest values. The reason is that the VM is assumed that modeling is therefore expressed as a function of droplet velocity and also depends on the speed of the droplet concerning the fluid. Figure 12 shows the unsteady droplet velocity profile for different droplet sizes. From the figure of the VM case at time = 0.01 ms, it can be seen that the biggest droplets have higher velocities than the smaller ones. At time = 0.2 ms, at the end of injection time and more atomization process, the largest droplets are accomplished with the maximum speed at a shorter time, and afterward, they break up into smaller droplets.

Figure 10

Calculated probability density function at different velocity profiles at other times, (a) Time = 0.01 s, (b) Time = 0.1 s, and (c) Time = 0.20 s.

Figure 11

Time evolution of velocity profiles with droplet diameter, (a) Time = 0.01 s, (b) Time = 0.1 s, and (c) Time = 0.20 s.

Figure 12

Time evolution of probability density function variation with droplet velocity profile different times, (a) Time = 0.01 s, (b) Time = 0.1 s, and (c) Time = 0.20 s.

Figure 13 shows the temporal variation of the PDF with droplet velocity at (a) 0.01 ms and (b) 0.2 ms. At time = 0.01 ms, both UV and RV cases assume a linear relation between PDF and droplet velocity. There is an increase in both the droplet velocity and the size during the injection period. At the same time, in the VM case, the PDF and droplet velocities vary nonlinearly. This result indicates that the droplet velocity decreases by drag force and atomization moves downstream. At time = 0.2 ms, in both UV and VM cases, the variation of the PDF and droplet velocity increases. In contrast, the RV case decreases due to the enlargement of the difference between the velocity of the two phases. Considering the momentum and turbulence modeling, the ( k − ε ) model is implemented to characterize the gas turbulent flow field. Figure 13 shows the axial gas phase velocity contours at time = 0.2 ms. UV and RV cases display the same structure. At the same time, developing the axial gas phase velocity structure for the VM case is better than other cases. This leads to a good induction for the gas phase from the stagnant condition.

Figure 13

Structure of the gas velocity at different velocity profiles.

5 Conclusion

The presented approach is used to develop the droplet velocity profile using QMOM. This model proposes the following:

1. A model for spray droplet velocity that can handle the entire size range generated during the atomization processes in the atmosphere.

2. The droplet velocity profile representation is defined using moments. The starting point uses the power–law relation, which is implemented in many spray droplet applications. The constructed functional form representing the drop size can be used to specify the first coefficient. The second coefficient’s value must be compared to experimental work as a continuous distribution function of diameter over the entire drop size range.

1. The gamma distribution function used to model the drop size distribution in a QMOMs indicates that the numerical results from the model agree with the experimental data for liquid nitrogen spray. This is a good distribution function to characterize cryogenic jets’ spray drop size distributions.

2. The application of distribution coefficients avoids inaccuracies at the zones of discontinuity for droplets of various sizes, and enables continuous integration across the droplet size spectrum.

3. The proposed formulation is used to predict spray tip penetration compared with experimental data and previous formulations for spray diesel.