Direct Numerical Simulation of Water Droplets in Turbulent Flow

Details on the fall speeds of raindrops are essential in both applications and natural events, such as rain-rate retrieval and soil erosion. Here, we examine the influence of turbulence on the terminal velocity of two water drops of different sizes. For the first time, computations of droplets in turbulent surroundings are conducted with a direct numerical simulation code based on a volume of fluid method. Both the drop surface deformation and internal circulation are captured. The turbulence intensity at the inflow area, as well as the turbulence length scale are varied. In turbulent flow, we find a decline in the terminal velocities for both drops. Based on the study of the wake flow characteristics and drop surface deformation, the decrease in the terminal velocity is found to be directly linked to a shortening of the wake recirculation region resulting from an earlier and more drastic increase in the turbulence kinetic energy in the shear layer. The turbulent surroundings trigger substantial rises in the drop axis ratio amplitude and a slight increase in the drop oscillation frequency, but barely influence the time-averaged drop axis length.


Introduction
Details on the fall speeds of raindrops are essential in numerous applications. These include using areal detectors and Doppler radar to estimate drop size distributions [1,2], employing radars to retrieve rainfall rates [3][4][5], modelling of coalescence and collisional breakup processes [6], soil erosion due to the impact of natural rainfall [7][8][9] and telecommunications [10]. In these applications and many others, raindrops of various sizes have been assumed to fall at their terminal velocities v t . The measurements by Gunn and Kinzer of v t of water droplets are regarded as the standard and are still generally used today [11].
However, several recent studies employing a variety of measuring techniques provide proof of the substantial fall speed discrepancies of raindrops under natural rains conditions and challenged the presumption of raindrops falling at their terminal speed. Löfler-Mang and Joss [12] discovered a substantial scattering of raindrop speeds around the terminal velocity v t curve. They associated the scattering with the turbulent air close to the ground and instrumental limitations. Testik et al. [13] found a reduction of around 10% in the fall speed of 1.9 mm raindrops from v t throughout calm wind conditions. They proposed that the increased drag force on the raindrops was triggered by large amplitude oscillations. Using 2D cloud and precipitation probes, Montero-Martinez et al. [14] identified the presence of super-terminal raindrops with sizes less than 0.7 mm during strong rainfall events. An explanation for the presence of super-terminal drops is that large raindrops have broken up into smaller ones and retain the speed of the original larger drop. The existence of super-terminal raindrops was also detected and confirmed by the following work of Larsen et al. [15] and Montero-Martinez and Garcia-Garcia [16]. Niu et al. [17] discovered that the fall velocities of raindrops covered a wide range for all sizes of raindrops in both stratiform and convective rains. Small raindrops (<1.4 mm) exhibit extensive positive mean deviations from their terminal velocities and larger drops with weaker negative deviations from v t . They associated the results with the influence of lower air density and a typical result of combined effects (e.g., turbulence/organised air motions, coalescence/breakup, instrumental errors). Using a C-band polarimetric radar and 2D Video Disdrometer (2DVD), Thurai et al. [18] identified a robust negative skewness and a wide broadening of the fall velocity distribution for 3 mm raindrops in a rain event dominated by an embedded convective line. They suggested that asymmetric horizontal-mode oscillations may result in an increase of the drag and a decrease of drop fall speeds. Montero-Martinez and Garcia-Garcia [16] detected sub-terminal raindrops smaller than 2 mm in size. During windy periods, there exists a larger portion of sub-terminal drops than in calm wind periods. Drop oscillations, the turbulence effect and transverse drift enhancement and deviations in the vertical trajectory were given as potential influential factors for the presence of sub-terminal raindrops. Subsequently, Bringi et al. [19] discovered that during heavy rainfall with strong wind gusts and high-intensity turbulence, the fall velocity distributions for drops of 1.3 mm, 2.0 mm and 3.0 mm were significantly broadened with negative skewness. In the case of 1.3 mm and 2.0 mm drops, the decrease in the mean fall speeds are stronger, with a maximum decrease around 25-30% below v t . With higher turbulence strength, the decrease in mean fall speeds was found to be more significant.
In the present study, we examine the influence of turbulence on falling water drops through Direct Numerical Simulation (DNS). Stout et al. [20] simulated the motion of particles within isotropic turbulence. To predict the fluid velocity, a stochastic model based on the Markov chain method [21] was used. They showed that nonlinear drag effects may result in a slowing of the settling velocities for large, heavy particles. The reduction in the mean settling velocity of particles around 2 mm can reach more than 35%. The reduction in the mean settling velocity increases with increasing particle Reynolds numbers. The Markov chain method is efficient in simulations of turbulence corresponding to prescribed turbulence statistics, but cannot provide accurate quantitative details of the turbulence, e.g., coherent structures. DNS does not suffer from this issue and can adequately resolve the features of the flow field because the grid resolution of a DNS goes down to the Kolmogorov length scale. In addition, during the falling process, the shape altering oscillations for drops larger than 1 mm in size need to be considered. We therefore conducted numerical simulations in the present study with our DNS code FS3D, which applies the volume of fluid (VOF) method with PLIC to account for the surface detection of a droplet.
The primary goal of this work is to understand how turbulence influences v t of a 3 mm and a 2 mm water drop. We evaluated the influence of turbulence on the drop oscillation behaviour and the flow field over the drop (e.g., size of the wake recirculation region, base suction pressure). Within all our simulation cases, v t of both water drops decreases in turbulent surroundings. We found this relates directly to the shortening of the wake recirculation region.
The following provides a brief discussion of the physics regarding the dynamic behaviour of a falling particle and a raindrop.

Flow over a Spherical Particle
A spherical particle falling in still air undergoes acceleration until the steady-state is reached when the aerodynamic drag force is in balance with the gravitational and buoyancy forces acting on the particle. The particle then falls at its terminal velocity. Flow over particles of different sizes and those travelling at different speeds may differ significantly, thereby, resulting in large variations of the drag coefficient C d , with C d = 4Dg∆ρ/(3ρ air v 2 rel ) [22]. Here, D donates the particle diameter, g the gravitational acceleration, ∆ρ the density difference between the particle and the air, v rel the relative velocity between the air and the particle and ρ air the density of the air. The flow field can be characterised by the Reynolds number (Re = ρvD/µ), which directly relates to C d . Here, ρ, v and µ are the density of the fluid, the velocity of the particle and the dynamic viscosity of the fluid, respectively.
Clift et al. [22] recommended a standard drag curve for a sphere, where C d is plotted as a function of Re. Previous studies provided a detailed discussion including the rich physics phenomenon of the flow field over a sphere at different Reynolds numbers. These include the work of Clift et al. [22], Johnson and Patel [23], Sakamoto and Haniu [24], Leweke et al. [25], Mittal [26], Kim and Durbin [27], Yun et al. [28] and Van Dyke [29].

Hydrodynamics of Raindrops
Falling raindrops can vary in diameter between 100 µm and several millimetres; a typical raindrop lies between 1 mm and 3 mm [30]. Smaller raindrops (<1 mm) are usually spherical as they fall through the air [31], mainly due to high surface tension forces. Larger raindrops undergo relatively low amplitude oscillations and take the shape of an oblate spheroid with a flattened base. The base becomes more flattened as the drop diameter increases. In the work of Szakall et al. [30], time averaged axis ratios for a variety of differently-sized drops, obtained from a number of experiments, were summarised.
Drop oscillation and surface deformation are primarily triggered by the aerodynamic forces acting on the surface of a raindrop [32]. Small raindrops oscillate with high frequencies and negligibly small oscillation amplitudes. The oscillation amplitude increases with increasing drop sizes, and the oscillation becomes noticeable for drops lager than 1 mm in size (Re ≥ 300) [32]. As the drop size increases, the oscillation frequency decreases significantly, with 1 mm drops having a frequency of 300 Hz and 5 mm drops having a frequency of around 30 Hz [30,32]. Continuous changes of the pressure and drag induced by the vortex shedding starting at an onset Reynolds number of about 300 was proposed as one intrinsic mechanism for the driving of continuous oscillation of drops ≥1 mm. According to Tokay and Beard [33], elements such as wind shear, turbulence and drop collisions were not sufficient to maintain the oscillations of raindrops against viscous dissipation.
Another essential characteristic of a drop is the internal circulation, which is initiated by the external flow fields. The internal circulation depends on the size of the drop and interacts with both the drop oscillation and the surface deformation. Pruppacher and Beard [34] observed a well-developed internal circulation and an under-developed secondary recirculation within a water drop that falls at its terminal velocity (Re = 135). The secondary recirculation circulates in the opposite direction and is induced by the reverse flow in the wake region [34]. For drops larger than 1 mm in size, the organised internal circulation tends to be disrupted. It forms and breaks up periodically as a result of the drop oscillation and surface deformation induced by vortex shedding [35]. Szakall et al. [36] stated that the opposite recirculation within large raindrops becomes more energetic with increasing Reynolds numbers and may ultimately interact with the primary circulation. Hence, the flow within the drops becomes chaotic; the drops are internally mixed. LeClair et al. [35] proposed that the internal circulation barely influences the shape of a drop smaller than 5 mm in size and that the influence of the internal circulation on the drag of a water drop ≤1 mm falling in air is negligible. Their numerical analysis also displayed a backward shift of the flow separation on the drop surface with 10 ≤ Re ≤ 300 as a result of the effect of the internal circulation, which is driven by the tangential stress from the airflow.

Drag Coefficients and Terminal Velocities of Water Drops
As we have discussed, small water drops of less than 1 mm in size retain basically a spherical shape, and the drag force on the drop is barely influenced by the regular internal circulation until the Reynolds number reaches around 300. These coincide with the minimum difference in the drag coefficients C d of a rigid sphere (Re ≤ 300) and a water drop smaller than 1 mm in size at its terminal velocity [31], as shown in Figure 1. As Re reaches values larger than 300, water drops begin to oscillate, deformation of the drop is more obvious, and the internal circulation becomes chaotic at higher Re. Accordingly, the divergence in C d intensifies progressively with increasing Reynolds number (larger drops). Terminal velocities, as well as the corresponding Reynolds numbers, for water drops of various sizes are displayed in Figure 2, according to the measurements of Gunn and Kinzer [11].  [11]) in comparison to the drag curve of rigid spheres suggested by Clift [22].

Numerical Methods
We conducted direct numerical simulations of a water drop falling in turbulent flow using the DNS multiphase code FS3D [37]. FS3D solves the incompressible Naiver-Stokes equations. For the description of multiphase flows and the interfaces between immiscible phases, the VOF method [38] using a Piecewise Linear Interface Calculation (PLIC) is applied [39].
The flow field inside a drop and its surrounding air are calculated using the conservation equations for the mass and momentum of incompressible flow, Here, ρ donates the density, t the time, u the velocity vector, p the pressure, µ the dynamic viscosity and g the gravity vector. The surface tension force f st is modelled with the balanced-force Continuum Surface Force (CSF) method combined with Height Function (HF) curvature technique proposed by Popinet [40]. The pressure field is computed by solving a Poisson equation implicitly, using a multi-grid algorithm [41].
The VOF variable f is used to identify the location of the gas and liquid phases and is defined as: for interfacial cells, 0 in the continuous phase. ( The transport equation of the variable f reads: To suppress numerical diffusion of the liquid phase, the PLIC method proposed by Rider and Kothe [39] is used to reconstruct the interface after each time step. Within each interface cell, the reconstruction is done by the approximation of a normal vector n of the interface where n = −∇ f / ∇ f . The position of the reconstructed plane corresponds to the VOF variable f in each interface cell and is solved through iterations [41].
Following the one-field formulation, the material properties of the liquid and gas phase are calculated as follows, where the subscripts g and l indicate the gas phase and the liquid phase, respectively. The finite volume method based on a MAC staggered grid is employed for the spatial discretisation [42]. Hence, all scalars are stored at cell centres and velocity vectors at cell faces. A first-order explicit Euler scheme is used for the time integration. The maximum time step sizes are restricted by the Courant-Friedrichs-Lewy (CFL) number.
The random spot method of Kornev and Hassel [43] is used to generate anisotropic turbulence, with a freely decaying energy spectrum, within the flow field. The idea is to place a collection of random rotating vortons with a specific inner velocity distribution on the front edge of the inflow area, which then move across the inflow plane at a uniform velocity. An application of the turbulence generator applying this method in FS3D can be found in the work of Zhu et al. [44].

Computational Setup
The falling of a water drop in the air flow was computed with an equidistant Cartesian mesh in a 3D rectangular domain. The dimension of the domain is 19.2D × 12D × 12D, with D representing the equivalent drop diameter. The resolution of the drop diameter is about 26.7 cells; the whole domain is resolved by 52.4 × 10 6 computational cells. Figure 3 shows the computational domain in the form of a 2D slice across the centre of a falling water drop. To minimise the computational effort, a grid that moves in accordance with the falling drop was used. This means that the water drop is positioned at the centre of the yz-plane for the whole falling process with a constant distance of 6D from the inflow and a wake length of 13.2D. The domain width in both z-and y-directions is 12D. This setup was chosen out of four separate series of convergence tests as to guarantee that for all simulation cases, both the grid resolution and the domain size (e.g., wake length, domain width, distance of the drop centre from the inlet plane) have a negligible influence on the fall velocity of a water drop in turbulent air flow (see Appendix A).
As seen in Figure 3, the gravitational acceleration g points in the direction of the negative x-axis. We applied the no slip wall condition at the inflow and continuous boundary conditions (Neumann) on the rest of the domain boundaries. In the case of still air, the air flow at 293.15 K enters the computational domain with a uniform velocity. In the case of a water drop falling across turbulent air flow, the random spot method by Kornev and Hassel [43] was applied for the generation of a turbulent inflow. The grid used is fine enough to resolve the smallest turbulent scales. For a 3 mm water drop falling at its terminal velocity, the ratio between the cell size L c and the Kolmogorov length scale λ k = L t /Re t 3/4 is around 2.3 with Re t = ρ air u 2 L t /µ air , the turbulence length scale L t = 0.75 mm and a turbulence intensity Tu = 10%. In Table 1, the relevant physical properties of the water and air are summarised, with σ, µ and ρ, denoting the surface tension, the dynamic viscosity and the density, respectively. Table 1. Physical properties of the air and water drop (20 • C) within the simulations.

Results
In this section, we present the simulation results regarding how turbulence influences a 3 mm and 2 mm falling water droplet. We will also show that the results for the still air cases match quite well with previous studies, in terms of the terminal velocities, drop oscillation behaviour and wake flow characteristics. In cases of turbulent inflow, we set the turbulence intensity to either 5% or 10% and the turbulence length scale L t to 1/4D or 1/2D. Tables 2 and 3 summarise the v t measured by Gunn and Kinzer [11] and the current simulation results from FS3D. In Table 3, with regard to the name of each test case, D indicates the equivalent drop diameter, L the ratio of the turbulence length scale and the drop diameter and T the turbulence intensity. Hence, D2L050T10denotes, for instance, the case of a 2 mm water drop falling in a turbulent surroundings. At the inflow area, the turbulence length scale is set to 1 mm, and the turbulence intensity is 10%.  Results from the simulations were all averaged over a time period of 0.5 s after the fall velocity of the drop stabilised. The trend for v t of both the 3 mm and 2 mm water drops matches well with the measured values, as in Tables 2 and 3. Figure 4 illustrates the v t of a 3 mm water drop according to Gunn and Kinzer [11], along with the evolution of the fall velocity of the simulated 3 mm water drop that falls in still air.  In most cases (except D3T0), the water drop was initially set to fall at a velocity close to the expected v t . This approach saves the computational time required for a falling water drop to finally stabilise at its v t .

Terminal Velocity
In turbulent air flow, v t of both water drops is reduced. The fluctuations of the drop fall velocities are comparable in all conducted cases. The maximum reduction in v t reaches 7.1% and occurs in the case D3L050T10. The results in Table 3 suggest that in the range under investigation, a larger turbulence length scale and higher turbulent intensities at the inflow area tend to cause a stronger reduction of v t , as well as a more intense rise in C d . In addition, in terms of v t , turbulence has more influence on the 3 mm water drop with a higher Re than on the 2 mm water drop.  Figure 6 shows an example of the temporal variation of the surface area ratio A dr /A sp , the normalised major (horizontal or lateral) axis lengths 2a/D, 2a /D and minor (vertical) axis lengths 2b/D of a 3 mm water drop falling in stagnant air over several oscillation periods. A dr and A sp denote the surface area of a drop and that of a sphere of equivalent diameter, respectively. 2a/D was measured in the y-direction (orange curve) and 2a /D in the z-direction (green curve). The variation in A dr /A sp and 2b/D fits well. As the minor axis reaches its minimum, the drop surface area reaches its maximum and vice versa. In Figure 7, we demonstrate the surface deformation across one oscillation period of a 3 mm water drop falling in still air. During the falling process, the water drop is flattened and takes an oblate shape. The deformation of the water drop at the instants A, B, and C displays a good match with the photos for the falling rain drops of 2.7 mm and 3.45 mm taken by Pruppacher and Beard [34]. In addition, we compared our results regarding the oscillation frequency, the axis ratio amplitude and the mean axis ratio with the measurements and theoretical models from previous studies (Figures 2-4 from the review of Szakall et al. [30]) and found a quite good agreement for the simulated 3 mm and 2 mm water drops that fall in still air (see Appendix B).  Table 4 summarises the time-averaged drop axis ratio, horizontal axis length 2ā, oscillation frequency and the temporal axis ratio amplitude for all conducted cases. Turbulence exerts only a minimum effect on both the drop axis ratio and horizontal axis length. The largest increase in both parameters is less than 1%, which occurs in the cases D3L050T10 and D3L050T05. Thus, we suggest that the mean effect of 2ā on the pressure drag is invariant to moderate turbulence. In turbulent surroundings, the increase in the axis ratio amplitude is obvious particularly in cases such as D3L050T10 and D2L050T10, both with an increase of over 50%. This means a higher temporal variation of the oscillation amplitude for a water drop that travels across a flow with a larger turbulence length scale (comparable to the size of the drop) and higher turbulence intensity. However, this large variation in the oscillation amplitude was not sufficient for an abrupt fall velocity variation (see Figure 5). In turbulent air flow, the oscillation frequencies either increase slightly or remain at the same values. The increase is more obvious for the cases with a flow field of a larger turbulence length scale and higher turbulence intensity.  Figure 8 compares the instantaneous velocity field of a 3 mm water drop falling in stagnant air (D3T0) with that in turbulent air flow (D3L050T10). The flow separates from the drop surface near a point at which the velocity in the stream-wise direction U x turns from a positive to a negative value. This is illustrated by the colour change from green to blue near the drop surface. The boundary layer thickness is indicated by the change in colour from red to yellow near the drop surface, which occurs at about 0.95 v t . Downstream of the drop base, the blue colour, where U x adopts negative values, characterises the wake recirculation. As shown in Figure 8, the turbulent surroundings do not seem to have an obvious influence on both the flow separation point and the boundary layer thickness, but it has a strong influence on the wake recirculation area. The water drop falling in stagnant air shows a recirculation area that is about twice the length of that of the water drop in turbulent air flow. A comparison of the instantaneous flow field including the velocity field along with the streamlines, the vorticity field and the pressure field between both cases is depicted in Figures 9 and 10. The vortices shown in both figures visualise the wake recirculation region of the water drop, which ends where the streamlines merge. The colour change downstream the wake recirculation clearly illustrates the recovery of flow velocity in both cases, where the tortuosity of the streamlines is weakened.  Compared to the case in stagnant air, the streamlines in the case D3L050T10 are crooked in the inflow region and far field. In both cases, the vorticity magnitudes increase significantly on the front half of the drop surface, which then decrease continuously along the shear layers. At the rear end of the recirculation region, the flow demonstrates high pressure losses and strong vorticity. We explored the vortex shedding process of both cases and discovered that the turbulent surroundings triggers an earlier onset of the shear layer instability, as well as a more irregular, unsymmetrical and much smaller wake recirculation region. In the case D3L050T10, the vortex is formed near the drop surface immediately after the flow separates. In the case of stagnant air, the main vortices within the shear layer mostly start to form at around 1.0D downstream of the flow separation point. We observe increasing pressure losses at these vortices as they travel downstream. The detachment of a vortex within the shear layer can be recognised by the saddle point of the streamlines beside the vortex, where the pressure loss and the vorticity magnitude in the vortex core reach local maxima. With a drastic reduction in the size of the wake recirculation, the strong pressure loss draws nearer to the base of the water drop. The base suction pressure -C pb is thus higher in the case D3L050T10 than that in the still air case. Similar changes in the near wake region were found for the cases with a less intense turbulent inflow, as well as in the cases with the 2 mm water drop. Hence, we think that this might be one potential influential factor that relates to the reduction of the terminal velocity of a water drop.

Mean Flow Parameters
In this section, we present the quantitative details of the flow field characteristics. We will first show that our results for the still air case are in good agreement with previous studies. Table 5 presents a summary of the time-averaged flow parameters, which include the flow separation angleθ s , recirculation bubble sizeL/D (see Figure 11), base pressure coefficient −C pb = (p b − p ∞ )/(1/2ρu 2 ), vortex formation lengthL v /D, minimum pressure coefficient on the wake axis −C pmin and the Strouhal number St = f rD/U. Here, p b denotes the time averaged pressure at the droplet basis and f r the vortex shedding frequency. Table 5. Comparison of statistical flow parameters for the simulated water drops falling in stagnant air and turbulent surroundings, e.g., the flow separation angle, base suction pressure, length of the wake recirculation, minimum pressure coefficient, vortex formation length and Strouhal numbers. To attain converged statistics, we averaged the flow field of the 3 mm and 2 mm water drops with results from 400 time steps for a time period of 0.28 s and 0.  To the best of our knowledge, there exists no measurement data for the flow characteristics (e.g., flow separation, size of the recirculation bubble) of a water drop at Re ≈ 1000. Since both the 3 mm and 2 mm water drops still retain a sphere-like shape in spite of the surface deformation, we may expect the flow in the near wake of a sphere to be similar to that of a water drop at the same Re. Hence, we conducted the simulation of a sphere falling in air with Re = 1000 for the validation. The computation was done using FS3D with a recently developed method for arbitrarily-shaped rigid bodies [45], which was specifically designed to account for the flow velocity within each interface cell of a rigid body.
The simulated time-averaged flow separation occurs at 101.7 • from the stagnation point of the sphere, and the length of the recirculation area is found to be 1.71D. Both values match very well with the results from Tomboulides and Orszag [46]. According to their work, the flow separates at 102 • , and the length of the recirculation region is about 1.7D. In the case D2T0 (Re ≈ 1000), the flow separation occurs earlier on the drop surface withθ s = 90.6 • . We suspect the deformation of the water drop to be one of the main reasons resulting in an earlier onset of the flow separation. Note that, during the simulation, the water drop displays a minimum shift of the drop centre relative to its original position through the falling process. Since the calculation of the separation angle is based on the original position of the drop centre, this has brought some uncertainties. We estimated the shift of the water drop through overlapping the time-averaged pressure distribution of the drop onto the contour of a 2 mm sphere and found an uncertainty of +2.3 • in the separation angle, meaning that in the case D2T0, the flow separation may take place at an angle of up to 92.9 • .
To validate the velocity field in the near wake region, we compared the stream-wise velocity U x /U on the wake centre-line with previous simulation results (Re = 1000) from Tomboulides and Orszag [46], as well as with the experimental data (Re = 960) by Wu and Faeth [47]. Figure 12 shows a good agreement between the results from FS3D and those from the previous studies. The length of the wake recirculation bubbleL/D of the 2 mm water drop is 2.5D, which is about 47% longer than that of the sphere. This matches well with the displacement of U x /U in the downstream direction, as shown in Figure 12. In the case of a falling water drop, the drop surface deformation and drop oscillation will have some influence on the length and width of the wake recirculation region. Furthermore, we suspect an increase in the width of the shear layer, which wraps around the water drop, as a result of the earlier onset of flow separation. This should eventually lead to a largerL/D. The large-eddy simulation (LES) conducted by Yun et al. and the DNS by Rodriguez et al. [28,48] show comparable results forL/D andθ s in the case of a flow over a sphere at Re = 3700.L/D was predicted to be 2.  Through our observation, the Strouhal number St of the 2 mm water drop falling in still air was found to be significantly lower than that of the 3 mm water drop. Moreover, the vortex shedding process of these two cases has been noted to be slightly different. In the case of the 3 mm water drop, vortices are often observed travelling further downstream after detaching from the shear layer, while in the case of the 2 mm water drop, vortices are most likely to be drawn to the centre-line of the wake after their detachment, which then disappear shortly after. A comparison of St for the wake of simulated water drops and that of spheres at different Re is shown in Figure 13. St of the simulated sphere obtained through observation is 0.26 and lies between the low-mode St associated with the large-scale wake instability and the high-mode St for the small-scale shear layer instability measured in the previous studies. At the same Re, St appears to be higher in the rigid sphere case than in the case of a water drop. [49] studied the wake of a drop in a liquid system with Re ranging from 350 to 500. They also reported the St of the drop, which were lower than those of a rigid sphere. Figure 13. Comparisons of the Strouhal numbers for the simulated 3 mm and 2 mm water drops in stagnant air, the simulated sphere at Re = 1000 and the spheres from the previous studies. The figure is reproduced from [28] and modified with our simulation results, with the permission of AIP Publishing.

Magarvey and Bishop
Turbulence has a great influence on the vortex shedding process in the shear layer, as well as on the flow characteristics in the near wake of the water drop. In turbulent cases, the formation and shedding of vortices take place at an earlier point in the shear layer and at higher frequencies.
Cases D3L050T10 and D2L050T10 exhibit the earliest vortex formation near the drop surface; cases D3L025T05 and D2L025T05 display the lowest level of the displacement of vortex formation in the upward direction relative to the still air cases.
As in Table 5, the turbulent surroundings cause a backward movement of the flow separation pointsθ s , an increase in the critical wake flow parameters St, −C pb and −C pmin , as well as a shortening inL v /D andL/D. It is clear that in stronger turbulence conditions, where L t = D/2 or Tu = 10% at the inflow, there is a more intense delay of the flow separation and a stronger rise in St.
Changes in the other four wake parameters all display a similar tendency that the reduction of L v /D andL/D, as well as the corresponding increase of −C pmin and −C pb become more significant in stronger turbulence conditions, which is consistent with the rise in C d . Figure 14 displays the strong relation between these four parameters amongst all simulation cases. As illustrated in the figure, the reduction inL v /D andL/D corresponds to a proportional increase in −C pmin and −C pb , respectively. More interestingly, for each case, the ratio between −C pmin and −C pb falls between 1.29 and 1.33, with a maximum ratio discrepancy of around 3%. In addition, the figure indicates less impact from the turbulence on the 2 mm water drop than on the 3 mm drop, where the shortening of L/D ranges from 16.4% to 45.6% with the 2 mm drops; with the 3 mm drops, from 37.2% to 56.4%. The increase in −C pb lies between 9.6% and 26% in the 2 mm drop case and between 26% and 40% in the 3 mm drop case. Figure 15 illustrates a strong connection betweenL/D and C d . In the case of turbulent surroundings, C d increases linearly with the decrease ofL/D. This however does not occur in the still air case, where the C d of both 3 mm and 2 mm water drops is located higher than the elongation of the dotted lines. This indicates that other factors, which tend to reduce C d in the case of a turbulent surroundings, may exist, and especially so for a 3 mm water drop.  Table 5).  Table 5).
In Figure 16, C d of the simulated water drops is plotted as a function of St. These data points lie close to the straight line with a slope of 0.299 (3 mm drop) and 0.176 (2 mm drop).  Table 5).
In the following section, we examine the distribution of the time-averaged Reynolds stresses (e.g., the stream-wise Reynolds stress u u /u 2 ∞ , the transverse Reynolds stress v v /u 2 ∞ and the Reynolds shear stress u v /u 2 ∞ ) and turbulent kinetic energy TKE in the flow field, as these further reveal the influence of turbulence on the shear layer and wake instabilities, which are closely associated with the vortex shedding process, recirculation in the near wake, as well as the base suction pressure. Figure 17 displays u v /u 2 ∞ of flow past a sphere at Re = 1000 simulated with FS3D. The result is in good agreement with that computed by Poon et al. [50]. Figures 18 and 19 depict the distribution of u u /u 2 ∞ , v v /u 2 ∞ , u v /u 2 ∞ and TKE for the 3 mm and 2 mm water drops, respectively. Compared to the sphere case, the Reynolds stress near the drop surface is increased, and we think that this relates more to the dynamic drop oscillation. For all conducted cases, the peak value of u u /u 2 ∞ appears at the shear layer and lies stream-wise between the end of the vortex formation region and wake recirculation region. Regions of stronger v v /u 2 ∞ are seen near the centre-line of the wake. The maximum of v v /u 2 ∞ is located directly behind the end of the recirculation area. The field of v v /u 2 ∞ and that of w w /u 2 ∞ (not displayed) downstream the drop base are found to be quite similar. The contours of the Reynolds stresses, in accordance with the shortening of the wake recirculation area, are thus shifted upwards closer to the drop base. Cases D3L050T10 and D2L050T10 display the largest displacement upwards for the 3 mm and the 2 mm drop case, respectively. In addition, in the cases of higher turbulence intensities and a larger turbulence length scale, TKE reaches higher values at the shear layer at an earlier stage due to stronger velocity fluctuations between the shear layer and the free stream. The stronger and faster increase in the TKE level corresponds to a more obvious displacement of Reynolds stress components, as well as a more significant reduction inL/D, as seen in Figures 18 and 19.
In cases D3L050T10 and D2L050T10, TKE reaches the highest value of around 0.1 and 0.08 at 0.85-1.1D and 1.15-1.25D at the shear layer downstream of the drop base, respectively. This implies that there is still room for a further reduction ofL/D in the case of a stronger turbulent air flow, where C d continues to increase until the boundary layer and the shear layer become fully turbulent and drag crisis sets in.  Figure 20a depicts the time-averaged streamlines for the 3 mm and 2 mm drop falling in stagnant air. The symmetrical wake recirculation area in both cases indicates good statistical convergence. In contrast to the 3 mm water drop, the 2 mm water drop is less flattened; the recirculation area appears to be shorter and narrower. Inside the 3 mm water drop, the streamline pattern illustrates a pair of well-developed circulation. The flow inside the 2 mm water drop displays a pair of well-organised secondary recirculation opposing the primary circulation.

Internal Circulation
Instantaneous streamlines of the 2 mm drop show that as the drop oscillates and gets more flattened with increasing major axis length, the core of the secondary recirculation is shifted towards the drop surface near the flow separation. As the major axis reaches its maximum, the secondary recirculation is disrupted and tends to break down. Inside the 3 mm water drop, the opposing circulation rarely appears and disappears shortly after as a result of the mixing effect of the near-wake flow. Figure 20b,c displays an example of the streamlines for 3 mm and 2 mm water drops at different time steps, respectively.

Conclusions
The present work investigates how the terminal velocity of a 2 mm (Re ≈ 900) and 3 mm (Re ≈ 1600) water drop are influenced by the turbulent surroundings by means of direct numerical simulations. We examined the turbulence effects on both the dynamic drop oscillation and flow field characteristics. At the inflow area, the turbulence intensity Tu was set to vary between 0%, 5% and 10%; with the turbulence length scale L t between a quarter or a half of the drop diameter. In the cases of still air, the simulated terminal velocity, as well as the drop oscillation behaviour (e.g., the drop oscillation frequency, drop axis ratio and axis ratio amplitude) are in good agreement with the results from previous studies. The characteristics of the flow field were validated through the DNS of the flow over a rigid sphere at Re = 1000. The flow separation point, the size of the wake recirculation region, the Strouhal number, the stream-wise velocity along the centreline of the near wake, as well as the Reynolds shear stresses all match well with previous studies.
The DNS results show a decrease in the terminal velocities v t for both water drops falling in turbulent air flow. The drag coefficients C d of both water drops increase. A stronger turbulent air flow (higher Tu and a larger L t ) at the inflow region results in a greater decrease in v t .
The study of the drop oscillation behaviour and the characteristics of flow over the water drop reveals that the rise of C d is strongly associated with the decrease in the length of the wake recirculation regionL/D, as well as the rise in the base suction coefficient −C pb . In terms of the drop surface deformation, the time-averaged drop axis ratio and length of the drop major axis are barely influenced by the turbulence. In the cases of stronger turbulence surroundings, the increase of the drop oscillation frequency becomes obvious; the rise in the axis ratio amplitude becomes much more significant. However, the fluctuations of the drop fall velocities are still comparable in all simulation cases. A great impact from the turbulent air flow on the near wake flow region is revealed by the instantaneous streamlines, which clearly display the increased irregularity and shortening of the recirculation area in the near wake region. Accordingly, the pressure loss at the rear end of the wake recirculation is shifted upstream.
Regarding the influence of the turbulence on the statistical flow parameters, we found an increase in the flow separation angle, base suction pressure −C pb , minimum pressure coefficient −C pmin and Strouhal number St and a decrease of the wake recirculation lengthL/D and vortex formation lengthL v /D. The impact from a stronger turbulent air flow is more intense on these parameters. Strong correlations between the wake parameters −C pmin , −C pb ,L v /D andL/D are displayed.
For all simulation cases, the ratio between −C pmin and −C pb lies near a constant value (around 1.3). The reduction inL v /D andL/D corresponds to a proportional increase of −C pmin and −C pb .
Furthermore, in the cases of turbulent surroundings, C d increases linearly with a decrease of L/D. The turbulent air flow results in an earlier and more frequent vortex shedding. The increase in the Reynolds stresses, as well as in the turbulent kinetic energy TKE takes place earlier at the shear layer and is more significant in the cases of stronger turbulent ambient flow. In accordance with the shortening of the wake recirculation area, the contours for the shear stresses and TKE are shifted upwards closer to the drop base.
In addition to the above-mentioned observations, the time-averaged streamlines display a pair of well-organised circulation inside the 3 mm drop and an extra pair of secondary recirculation inside the 2 mm drop in the still air case.
While the current study involves only cases with a limited range of turbulence levels, it is shown that with stronger turbulent air flow, there is still room for a further increase in TKE at the shear layer, shortening of the recirculation area and consequently an increase of C d until the drag crisis sets in. Regarding the simulation results from FS3D, a couple of interesting questions remain unanswered. These include the limit ofL/D, L t and Tu for the presence of drag crisis, the impact of ambient turbulence on the drag and the wake flow of a smaller falling water droplet (with a lower Re). Whether the linear relation between C d andL/D still holds with stronger turbulent air flow also requires further investigations. In addition, it remains to be explored whether or at which turbulence level the influence of the internal circulation/recirculation on the near wake would become strong enough to affect the fall velocity of a drop.    Figure A1 summaries the results, the velocity ratios between the simulated fall velocities and v t measured in still air of all cases, from these four studies. Based on the convergence of these results, we chose a grid resolution of 26.7 cells per drop diameter for the further simulations. With this resolution, we resolved the flow field until the Kolmogorov length scale (see Section 3). The wake length and domain width were selected as 13.2D and 12D, respectively. We suggest that the deviations of the velocity ratio in the cases d1-d4 are more related to the spacing for the pressure field around the flow stagnation point. A large spacing should provide a more reasonable solution of the flow field. Thus, we selected the position of the drop at a stream-wise direction of 6D, with consideration of the decay in the turbulent air flow from the inlet. Figure A2 shows a comparison of the simulated axis ratios, axis ratio amplitudes and drop oscillation frequencies for the 3 mm and 2 mm water drops in still air with results from previous studies. Values from the simulations were averaged over at least 0.5 s after the fall velocity of the water drop stabilised. The axis ratio amplitude was calculated with the difference between the maximum and the minimum temporal axis ratio during the drop oscillation periods. As shown in the figure, the simulated axis ratios and oscillation frequencies of both drops fit quite well with the results from previous studies. The axis ratio amplitudes of both water drops lie near most of the data from previous measurements.

Appendix B. Validations for the Drop Deformation
(a) (b) (c) Figure A2. Comparison of the axis ratios, the axis ratio amplitudes and the oscillation frequencies of the simulated 3 mm and 2 mm water drops in stagnant air: (a) Axis ratio. (b) Axis ratio amplitude. (c) Oscillation frequency. The figure is adapted from [30].