Probabilistic derivation of droplet velocity using quadrature method of moments

: Accurately prediction of dispersion and polydispersity of droplet ﬂ ow is not a trivial task due to the complex behavior of the droplet size distribution (DSD) and the strong state of the instantaneous velocity of a droplet on the shape and size of the droplet. Describing the distribution of sizes and velocities of droplets initially formed in sprays is an essential piece of information needed in spray modeling, which is used to de ﬁ ne the initial state of the spray droplets in the downstream two-phase ﬂ ow ﬁ elds ’ predictive computations. The predictive model for the droplet size and velocity distributions in sprays is formulated as the droplet ’ s velocity magnitude has a power – law rela-tionship with the droplet in this study. The present model incorporates the deterministic and stochastic aspects of the spray formation process. The quadrature method of moments (QMOM) is applied to solve numerically the transport equations of the probability density function coupled with conserved source terms incompressible Navier-Stokes equations for the liquid phase. The sub-models are connected by dif-ferent source terms signifying the liquid-gas interaction. Equations of transport for spray moments are derived from DSD, and closure is attained using a gamma distribution. The integer spray moments concerning the volume are used to construct the continuous distribution of QMOM. In contrast, the velocity moments are used to determine the droplet velocity as a constant function of the droplet diameter. The model is ﬁ rst applied to simulate a diesel spray tip penetration under nonreactive conditions with di ﬀ erent droplet velocity pro ﬁ les to validate the approach with experimental data. An additional case of a liquid nitrogen spray is applied to show the gamma distribution ’ s ability to describe spray drop size distribution. The model generally predicts reasonable agreement with the experimental for both cases.

Abstract: Accurately prediction of dispersion and polydispersity of droplet flow is not a trivial task due to the complex behavior of the droplet size distribution (DSD) and the strong state of the instantaneous velocity of a droplet on the shape and size of the droplet. Describing the distribution of sizes and velocities of droplets initially formed in sprays is an essential piece of information needed in spray modeling, which is used to define the initial state of the spray droplets in the downstream two-phase flow fields' predictive computations. The predictive model for the droplet size and velocity distributions in sprays is formulated as the droplet's velocity magnitude has a power-law relationship with the droplet in this study. The present model incorporates the deterministic and stochastic aspects of the spray formation process. The quadrature method of moments (QMOM) is applied to solve numerically the transport equations of the probability density function coupled with conserved source terms incompressible Navier-Stokes equations for the liquid phase. The sub-models are connected by different source terms signifying the liquid-gas interaction. Equations of transport for spray moments are derived from DSD, and closure is attained using a gamma distribution. The integer spray moments concerning the volume are used to construct the continuous distribution of QMOM. In contrast, the velocity moments are used to determine the droplet velocity as a constant function of the droplet diameter. The model is first applied to simulate a diesel spray tip penetration under nonreactive conditions with different droplet velocity profiles to validate the approach with experimental data. An additional case of a liquid nitrogen spray is applied to show the gamma distribution's ability to describe spray drop size distribution. The model generally predicts reasonable agreement with the experimental for both cases.

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 twofluid 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][8][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][15][16]. Fox et al. [17] introduced the socalled 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][18][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 volumeaveraged 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.

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 datathe 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 ) 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.

Transport equations for moments based on droplet volume
The geometry of the flow is assumed to be dilute dispersed flows of spherical droplets. The droplet volume is expressed in terms of its diameter as = V D π 6 3 . Here, a modification will be made to the number size distribution function in terms of droplet volume, which is denoted by In addition, the number of droplets per unit volume in the interval 2 , the following relation becomes: From equation (1), the velocity of droplets, u V ( ), is illustrated in terms of droplet volume as follows: Transforming the equation of population balance in terms of the volume size distribution becomes: However, from the aforementioned equation, it can be possible to convert the transport equations for moments as follows: ) . For such cases, closure relations are needed to be computed + M k β . In brief, the momentum equation is written as follows:

Gamma distribution
The general gamma number size distribution is given by and where k Γ( ) is the gamma function, and r 32 is the Sauter mean radius of the number size distribution of the drops.
= ∕ r M M 32 3 2 defines this. For numerical calculations, the gamma function can be approximated by [27]: with an error of at most 1% for values of > M 0 is calculated from (6) by setting = i 0. For all > k 0, the gamma distribution is defined. On the other hand, several constraints should be applied in practice, arising from the submodels employed, particularly the drag, breakup, and collision models used.

Cumulative gamma function
To evaluate, values of the cumulative gamma function must be tabulated for general matters of > k 0 and xc in the range ≤ ≤ ∞ x 0 c . A series of partial integrations leads to the recursion relation: When ≤ − ≤ k m 1 2 , this recursion is stopped because the values of Γ c then vary over a much smaller range than for any other unit integer range [28].

Determination of the polydispersed phase velocity
The value of the parameters β and α needs to be specified from the: With the gamma distribution in equation (18), it is simple to demonstrate the identity: thus, Taking = k 2 in equation (19), and with the use of Eqs. (20) and (21), leads to The parameter β needs to be limited to a relatively small decimal range. However, the content of drop velocity predicted locally in the spray may significantly exceed the actual spread. To set an appropriate restricting value for β.

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. SPRAY-MOM 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.

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.
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.
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).
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 rewatered 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 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).
This comparison aims to observe the spray droplets' dispersion, which gives a good indicator of the trajectory of  Figure 4: Comparison of predicted spray penetration with the experimental data [31]. Gas density (kg/m 3 ) 13.8 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 nv 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.
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 sizeto-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 massmean diameter for the whole spray. The case of UV is shown in Figure 7 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 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 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.

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.
The advantages of the new representation of the droplet velocity profile over previous models are as follows: (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.