Two-way Euler-Lagrange Simulation on Mid-ﬁeld Spray under Multi-physics Control

Spray formed by the process of high-speed gas shooting out of a round nozzle, along with liquid component teared, stretched and atomized into tiny droplets, exists extensively in natural environment and industrial activities. Many experimental and numerical works have been done for better understanding and estimating the complicated multiphase dynamics. The overarching aim of this report includes (i) simulating high-ﬁdelity turbulent sprays, (ii) validating mid-ﬁeld transport of polydisperse droplets and (iii) testing various artiﬁcial manipulation on the spray by external multiphysics control such as acoustic/pressure wave, electro-static ﬁeld and controlled introduction of swirl.


Introduction
The spatial development of a liquid spray can be conceptually divided into two parts: (i) the near field, where the injected liquid field still remains somewhat coherent but tends to distort, stretch, and breakdown into droplets due to the inter-phase shear forces and (ii) the middle field, where the liquid phase has completely broken into tiny spherical droplets and no more further fragmentation occurs. The Euler-Lagrange method is ideal in simulating this fully-dispersed mid-field. Specifically, each droplet becomes an individual Lagrangian point, whose instantaneous trajectory is traceable. To simplify the problem, we ignore the inner flow and surface tension of droplets and consider them as solid rigid particles of constant density but diverse diameter.
The inter-phase aerodynamic force and heat flux are coupled by point-particle models. In a two-way coupled simulation, the aerodynamic force and heat transfer exerted on droplets are fed back to local Eulerian gas field with a Gaussian spreading. The aerodynamic force consists of two main contributions: (i) fluid-droplet force, which is commonly the main force that results from the local macroscopic flow, and (ii) fluid-mediated droplet-droplet force, which acts as a perturbation force that takes microscopic wakes of neighboring droplets into account. The latter contribution can be significant to droplet dynamics even at the macroscopic level, especially under dense droplet-laden condition. However, spray experiments show that liquid volume fraction within the middle field stays quite low (usually < 10 −3 ). In this dilute limit, the droplet-droplet force, as well as collision and coalescence, are negligible.
There have been a number of computational and experimental studies focused on better understanding of the spray processes and the numerical studies have used fully-resolved, Euler-Lagrange and Euler-Euler approaches [1,2,3,4,5,6,7,8]. In this work, we follow our previous two-way Euler-Lagrange simulations on mid-field spray [9] and focus on testing the effect of external acoustic field. The standing wave of the acoustic pressure applies a second-order radiation force on droplets that could redistribute the spray dispersion. However, the acoustic effect varies with the droplet size, Reynolds number, acoustic energy, spatial location, and many other factors. To evaluate the acoustic effect in this complex system, we perform spray simulations within an axisymmetric acoustic field and analyze the statistics of droplet dispersion.
The organization of the rest of the paper is as follows. In Section "Simulation setup", we describe the simulation parameters and schematics. Then, Section "Governing equations" presents the governing equations of the turbulent flow and the droplets. Subsequently, we compare the difference in spray features without and with the acoustic control in Section "Results and discussion". Finally, we conclude in Section "Conclusion".

Simulation setup
This simulation is based on a corresponding experiment elaborated in [9]. We choose the length scale as the nozzle diameter d 0 * = 0.01 m and the velocity scale as the extruding velocity of the gas-phase U 0 * = 76.7 m/s, where the subscript "*" means a dimensional quantity (any variable without a subscript "*" is non-dimensional). The time scale is consequently T * = d 0 * /U 0 * = 1.30 × 10 −4 s. The kinematic viscosity of air at the laboratory temperature of 20 • C is ν f * = 1.52×10 −5 m 2 /s that yields a Reynolds number of Re = U 0 * d 0 * /ν f * ≈ 50, 000. The density ratio between water ρ p * = 998 kg/m 3 and air ρ f * = 1.20 kg/m 3 is ρ p = ρ p * /ρ f * ≈ 832.
The subscripts "f" and "p" denote the continuous phase (air flow) and the dispersed phase (water droplets), respectively. In the fully-dispersed mid-field, we track every individual droplet as a solid sphere (Lagrangian particle) of unique diameter within d p * ∈ [1.2, 250] µm. Droplet evaporation is ignored since the time span of each fast-moving droplet within the simulation domain is short. The droplet coalescence is unimportant since the droplet volume fraction is always lower than 10 −3 . Secondary break-up is also negligible in the mid-field. For notation of convenience, we will refer to the gas phase as fluid and water droplets as particles. Figure 1 shows the schematics of the simulation. The physical prototype is a water-air twophase round jet spray ejected from a co-axial air-blast nozzle (the bottom outlet of the nozzle is of non-dimensional diameter d 0 = 1 positioned at the origin of the Cartesian coordinate). Once ejected, the gas phase forms a turbulent round jet flow that was solved by LES (The two gray-scale contours on the back and left planes are the instantaneous velocity magnitude and λ ci respectively on the center slice). In the near-field 0 < z 10, the initially intact water phase is torn up by the shear stress of the high-speed air flow, and atomizes into numerous droplets of a wide range of size. Subsequently in the mid-field z 10, the water-phase is considered to be fully-dispersed that has completely broken into millions of spheroidal droplets. The blue dot cloud is an instantaneous distribution of the entire Lagrangian droplets. The two red dashed lines mark the radial boundaries of the droplet distribution (the averaged axial number flux on the boundary is about 1% of the value on the centerline), and the angle between them was computed as 10.54 • . By counting the droplets across the plane z = 50 over t ≈ 1000 nondimensional time, the PDF of droplet distribution is plotted as the bivariate histogram on the bottom plane with the peak value about 0.0214. It should be noted that this PDF is sharper than the Gaussian-like distribution. All of the above single-phase and two-phase simulations have been rigorously validated with the experimental measurements in [9]. Based on the above Euler-Lagrange LES of the spray, we further apply an axisymmetric acoustic wave that is homogeneous in the z-direction and effective within z > 0 (illustrated as the yellow semi-transparent volume). The cross-section of the acoustic force F ac at z = 20 field is shown as a semi-transparent contour. The radial variations of the acoustic force (red, vertical) and the acoustic pressure (blue, horizontal) are plotted on the right hand side. The blue solid and dashed lines are the fluctuating range of the pressure standing wave. Within one wave length λ = 60, the acoustic force field can be uniformly divided into four regions. Because the acoustic force points in the direction from the anti-node (solid, vertical) to the pressure-node (dashed, vertical), the moving direction of droplets (black sphere) within different regions are shown as the short arrows. Therefore, spray droplets within region I and III (respectively, II and IV) will radially expand (respectively, concentrate) in the radial direction.
The original simulation is performed in a rectangular box of size L x × L y × L z = 50 × 50 × 60 and split into N x × N y × N z = 20 × 20 × 64 spectral elements. Detailed settings of the domain, mesh, boundary conditions and the solver are elaborated in [9] (refer to the Case "TL2").The original spray simulations without acoustic control have been thoroughly validated by the corresponding experiment. In this work, we apply the acoustic wave directly on the original spray simulations. To simplify the acoustic field, we assume the spatial variation of the acoustic force to be expressed as It should be noted in the present configuration the centerline r = 0 is an anti-node and so is r = 30 and 60, with the pressure nodes in between. In the above equation the acoustic force F ac is solely in the radial direction. F ac,m = Y ff st EA p is the fluctuation magnitude and A p is the droplet cross-sectional area. The wave length is chosen as λ = 60 so that the spray droplets are mainly in the expansion region. The acoustic energy is chosen as E = 0.1. We assume the acoustic effect on deformation and break-up of droplets to be negligible. We also ignore other less important factors such as the influence of droplet clustering and wave reflection on the original acoustic field.

Governing equations
By defining a phase indicator function (I = 1 in the fluid volume, I = 0 in the particle volume) and the Gaussian distribution G(x) = 1 ( √ 2πσ) 3 exp |x| 2 2σ 2 as the mesoscopic filter, the volume fractions of the fluid φ f and particle φ p phases and the general filtered (volume-averaged) quantity ϕ are expressed as where ϕ can be any Eulerian field of the fluid phase such as velocity, pressure and temperature, while [ϕ] mean the corresponding fully-resolved field solved by direct numerical simulation (DNS). For the large-scale simulations at such a high Reynolds number in this work, we will solve the volume-averaged mesoscopic flow. Therefore, the governing equations of the fluid flow become where u and P are actually the Gaussian-filtered large-eddy flow velocity and pressure. ν t (x, t) is the turbulent kinematic viscosity calculated by the dynamics Smagorinsky sub-grid stress is the Eulerian feedback force from the N p droplets within the computational domain. Although the summation is formally written to be over all the droplets, the actual summation will include only the few droplets that happen to be near the point x where the Gaussian is larger than a set small threshold. Here x p,i and F hyd,i are the center position and the hydrodynamic force of the i th particle.
The dispersed phase is solved by time-integrating D Dt where m p,i = ρ p π 6 d 3 p,i is the particle mass. The total force on each droplet is the summation of the hydrodynamic force and the acoustic force. Since the density ratio ρ p 1, we assume the total aerodynamic force to be the quasi-steady force with the finite Reynolds correction Φ (Re p,i ) = 1+0.15 Re 0.687 p,i , where Re p,i = |u (x p,i )+u (x p,i )−v i |d p,i /ν f is the particle Reynolds number. In order to account for the effect of the sub-grid turbulent kinetic energy in the droplet motion, the sub-grid perturbation flow velocity u (x p,i ) at the center of each Lagrangian droplet is modeled by the Langevin equation [10,11] as where ε = ν t |S| 2 is the dissipation rate that depends on the filtered strain rate tensor S, ∆t is the simulation time step and C s is the dynamic Smagorinsky coefficient. Parameters C 0 = 2.1 and C Y = 0.039 come from previous works [12,13] and ξ ∼ N (0, 1) is a Gaussian-distributed random variable.
The acoustic force on the finite-sized spherical droplet in the viscous compressible medium has been derived [14,15] using the perturbation method. In an one-dimensional standing acoustic wave field, the net first-order acoustic force on the droplet is zero, while the secondorder radiation force is nonzero even after time averaging. This second order force is parallel to the pressure wave transmission direction and always points from anti nodes to the nearest pressure node. The force exerted on an individual droplet can also be expressed as where E is the acoustic energy density, A p = πd 2 p /4 is the cross-section area of the droplet, λ is the acoustic wave length, and h is the distance between the droplet and the nearest anti-node. The coefficient and it dependent on flow Reynolds number, wave number (k 0 = 2π/λ) and droplet diameter as shown in Figure 2. α n and β n are the real and imaginary parts of the n-th order perturbation of the scattering coefficient S n for the medium fluid outside the droplet [14,15]. This radiation force model for standing wave is valid if the acoustic wavelength is much larger than the droplet diameter and the flow diffusive thickness λ d p , δ.  In order to study how the external acoustic wave alters the original spray dispersion, we define the acoustic time scale τ ac = m p U 0 /F ac accounting for the particle responding time to the acoustic force. In comparison with the hydrodynamic force, Figure 3 plots the ratio of the acoustic time scale (multiplied by the acoustic energy) over the particle time scale τ p = ρ p d 2 p Re /18. In general, the acoustic force drives particles away from their original trajectories following the jet flow, while the hydrodynamic force correspondingly obstructs droplets from the external traction. The overall balance between these two counter effects, which is reflected by the ratio of the two time scales, dictates the consequent droplet dispersion. Within the particle diameter range in the simulations d p ∈ [1.2 × 10 −4 , 25 × 10 −2 ], the time scale ratio generally decreases exponentially with respect to the particle diameter, although there is an intermediate transitional regime that the time scale ratio inversely increases slightly. The Reynolds number and the acoustic wave length almost make no difference on the general shape of the curve, but only horizontally shift it to the left and right, respectively. This plot has been generated for unit acoustic energy. It should be noted that the acoustic force can be obtained by scaling with the actual acoustic energy. In addition, the large-size droplets can overlook the time scale ratio or the force balance, as long as their inertial effect is overwhelming. Figure 4 shows the former stationary state without the acoustic force and the acoustic controlled spray at t = 753 non-dimensional time after adding the acoustic force. On the one hand, the simulation without acoustic forcing illustrated in Figure 4(a) has been carefully validated by the former study [9], in which the droplet number flux, mean diameter, and axial velocity variations on the radial direction on five downstream planes have been proven to be highly consistent with the experiment with exactly the same input parameters. On the other hand, the acoustic model derived analytically [14,15] accurately evaluates the radiation force on droplets. Therefore, the current result shown in Figure 4(b) as the combination of the former two studies are of high reliability.  It is evident that the droplet cloud has been expanded in the radial direction. In the original simulation without the acoustic wave, the global droplet cloud is a thinner cone shape with more small (green) droplets inside while the larger (red) outside. After turning on the acoustic force, Both the large and small droplets are moved towards the spray edge so that the whole spray gets thicker. We notice that at the new spray edge, although the small droplets are still fewer than the large ones, much more middle size (yellow) droplets begin to show up.

Results and discussion
The second feature is that the droplet size distribution at the spray edge is varying in the downstream direction. There are more small droplets in the upstream edge, while the downstream edge is the opposite. This variation is mainly due to the inertial difference for different sizes. Given the same external force, smaller droplets are easier to be dragged away from their original trajectory, while the large droplets change their direction only slowly and follow a trajectory of much larger radius of curvature.
The last qualitative finding is that the original cone shape of the droplet cloud is not proportionally expanded, but now becomes more like a cylinder. This is partially because of the different responses and resulting different trajectories of the varying sized droplets. The other reason is that the acoustic field is homogeneous in the axial direction, and the acoustic force decays to zero at r = 15 at all z.
For an individual droplet, its dynamics are influenced by many factors. Generally speaking, the consequent trajectory it traces is the result of reconciliation among the flow force, acoustic effect, and its own inertia. But the spatially varying acoustic force and the temporarily varying turbulent eddies make the trajectory stochastic and undetermined. In comparison, the constant inertial effect is a much simpler factor purely dependent on the diameter. The droplet diameter determines both its reaction time to the acoustic and flow fields. For small droplets, the overall acoustic effect essentially depends on the balance between the acoustic force and the counteracting drag force, while for the large droplets, their motion depends on the balance between the acoustic force and droplet inertia. Figure 5 shows the normalized droplet number flux. The overall number flux in Figure 5a obviously decreases near the centerline but increases near the spray edge, and this change keeps amplifying in the downstream direction. Based on the value of normalized number flux near the centerline r = 0 with the acoustic force for the eight size levels discussed in [9] first decreases from level 1 to level 3, then levels off and even increase slightly from level 3 to level 6, but finally again continue to increase from level 6 to level 8. Especially for level 8, we find that the downstream curves start away from the centerline r = 0, which mean there are no longer such large size droplets at near the centerline since all of they are driven away by the acoustic force. Figure 5. Plots of normalized number fluxes of (a) overall, (b) to (i) are droplets of size level 1 to 8. "aco." means the case with the acoustic force, and "ori." means the original simulation without the acoustic force.
The mean diameters in Figure 6 shows that more small size droplets move to the outer annulus than the large size droplets so that the mean diameter increases in small radius but decreases in large radius, but this situation gets diminished in the downstream location. As discussed in the governing equation subsection, larger size droplets get expanded by the acoustic force much slower due to their stronger inertia, even though their acoustic-to-particle time scale ratio are smaller. Therefore, the expansion of larger size droplets will be more revealing in the more downstream position. Moreover, due to the radial variation of the acoustic force, we find the mean diameter has already got decreased near the centerline r = 0 due to the escape of large size droplets.

Conclusions
The mid-field spray under the acoustic control is a complex system that involves a number of factors and varies dramatically in different specific configurations. This work starts with a simple acoustic field and tries to investigate the inherent mechanisms that determines the droplet dispersion. We follow the original two-way Euler-Lagrange spray simulations without the additional control presented in [9] and compare the changes after adding the acoustic control, in which the acoustic wave length is chosen wide enough to ensure almost all droplets are within the expending region. Finally, we find that the droplet cloud is greatly expanded influenced by the acoustic control, and the new droplet could is no longer the cone shape but more like a cylinder.
We have analyzed the droplet dispersion by calculating the droplet number flux and mean diameter. Detailed statistics show that the acoustic effect is generally increasing with respect to the droplet diameter. However, the acoustic effect is not monotonic. It decreases somewhat for mid-size droplets, which can be seen for the λ = 60 curve in Figure 3. The acoustic influence is also varying in the axial direction. This effect is earlier to influence the small droplets than the large droplets.