Population Distribution in the Wake of a Sphere

The fluid physics of the heat and mass transfer from an object in its wake has much importance for natural phenomena as well as for many engineering applications. Here, we report numerical results on the population density of the spatial distribution of fluid velocity, pressure, scalar concentration and scalar fluxes of a wake flow past a sphere in the steady wake regime (Reynolds number 25 to 285). We find the population density to be well described by a Lorentzian distribution. We observe this apparently universal form both in the symmetric wake regime and in the more complex three dimensional wake structure of the steady oblique regime with Reynolds number larger than 225. The population density distribution identifies the increase in dimensionless kinetic energy and scalar fluxes with the increase in Reynolds number, whereas the dimensionless scalar population density shows negligible variation with the Reynolds number.


Introduction
The interactions between spherical bodies, such as, particles, bubbles and drops, and the ambient through which they move is a vast area of research, which has attracted attention over centuries in various scientific disciplines [1,2]. The flow past a sphere presents different regimes at different Reynolds number Re. The steady axisymmetric structure of a wake at low Reynolds number, up to Re ∼ 210 [3,4,5], followed by a steady oblique wake structure up to Re ∼ 280 [6,7], and an unsteady structure of the wake at higher Re [8, 9, etc.] had been studied both experimentally and numerically. The drag coefficient C D of a sphere, which varies with the roughness of the sphere surface and Re, was studied in detail both in experiments [10,11,12,13] and numerical simulations [14, 15, etc.]. At present, how the drag, lift and pressure coefficients vary both locally as well as globally with respect to the sphere is well known [5,15,16]. The two dimensional structures of the streamlines, vorticity and pressure contours along the orthogonal central planes through the sphere are also well known over various studies [5,6,17].
For many engineering applications and natural processes, the interaction between a sphere and the ambient also involves transport of various scalar species, either passively advected by the ambient flow or interacting actively with the flow through various physical processes, for example, through evaporation, buoyancy. The rate of scalar transport, in particular the convective heat transfer from spherical objects at various Re, has been investigated both numerically [17,18] and experimentally [19,20,21,22] to determine the heat transfer coefficient. Similar to the drag coefficient, attention was given to the dependence of the local Nusselt number (a ratio of the convective and the diffusive (conductive) heat transfer) on the sphere surface and its global average for different Re [17,22]. The profiles of the dimensionless temperature contours along the central orthogonal plane for various Re have also been described in the literature [17,23]. A coupled system involving an interplay between different scalars can also be present, for example, in case of the phase change during droplets evaporation or freezing resulting in heat and mass exchange with the ambient air. Such interaction has also been studied both experimentally [24,25] and numerically [26,23]. All these studies are mainly concerned with the average scalar flux at the surface of the sphere, which determines the mass and temperature change rate of the sphere. An overall description on the spatial structure of the wake, including the scalar concentrations and the convective fluxes for various Re was not fully explored.
Descriptive statistics on the spatial structure of the wake is of primary importance if the extent of the wake with certain properties needs to be quantified. Supersaturation in the wake of a precipitating cloud water droplets [27,28] for example requires a detailed analysis of the wake population. In this paper, we present a comprehensive numerical study on the details of the momentum and scalar transport in the wake of a sphere using a population density distribution for the steady axisymmetric and oblique wake regimes. A brief introduction of the numerical methods and computational details are described in Section 2. Results are presented and discussed in Section 3. Finally conclusions are given in Section 4.

Physical Model, Numerical Method and Boundary Conditions
We consider the flow which develops past a sphere, placed in incompressible viscous fluid with velocity u ∞ = (u ∞ , 0, 0), pressure p ∞ , and a constant density ρ. Together with the balances of mass and momentum, we consider also the transport of a passive scalar, that is any contaminant present in low concentration so that it does not influence the flow. Such dynamics is described in an Eulerian framework by an advectiondiffusion (AD) equation. If d s is the diameter of the sphere, θ the passive scalar concentration, θ s and θ ∞ are the scalar concentration on the surface of the sphere and in the external flow respectively, the problem can be suitably made dimensionless by using d p , u ∞ and θ s − θ ∞ as scales, and therefore by defining the dimensionless position, time, velocity, pressure, and scalar concentration as, Therefore, the dimensionless incompressible Navier-Stokes (NS) equations and the one-way coupled AD equation for the scalar are, where Re = u ∞ d s /ν is the Reynolds number (ν is the kinematic viscosity) and Sc = ν/κ θ is the Schmidt number, ratio between the kinematic viscosity and the scalar diffusivity κ θ . These equations are complemented by uniform flow boundary conditions far from the sphere (u * → (1, 0, 0), θ * → 0) and no slip boundary conditions on the surface of the sphere with a constant scalar concentration (u * = 0, θ * = 1). For sake of clarity, the * will be omitted in the following. These governing equations are numerically solved with the lattice Boltzmann method (LBM) [29,30]. A code is developed based on the open-source library, Palabos [31]. In LBM, the particle distribution function f (x, t) is governed by Here i is the index of the discrete velocity c, which defines the structure of lattice; x and t are the location of a lattice node and the time respectively. The collision operator Ω i (x, t) models the redistribution of the particle populations at each lattice node. In this study, we consider the Bhatnagar-Gross-Krook (BGK) collision operator [32], with which the population f i (x, t) relaxes towards its equilibrium state f eq i (x, t) according to the relaxation time scale τ . f eq i (x, t) and τ are defined as, Here w i is the weight; c s is the speed of sound. The macroscopic quantities, such as the density ρ and velocity u are moments of f i (x, t), according to respectively. For solving the fluid velocity field, the D3Q19 lattice is chosen as the non-linear momentum advection corrections are not very significant in the steady axisymmetric or oblique wake flows [33].
The one-way coupling between the fluid momentum ρu and the scalar concentration θ is solved by another LBM equation similar to equation 4, but with a distribution function g i (x, t) for the scalar. To recover the AD equation, the equilibrium distribution function g eq i (x, t) [34] and the relaxation time scale τ g are given as, ).
The scalar concentration θ is calculated according to θ = i g i (x, t) = i g eq i (x, t). Since only the zeroth and the first order moments of g i (x, t) are used to recover the AD equation from the LBM equation, a D3Q7 lattice is used for the scalar field [30].
The sphere is set in the origin of the reference frame, and the dimensionless domain is [−5, 20] × [−3.5, 3.5] × [−3.5, 3.5] in size (5 diameters upstream, 20 diameters downstream and 7 diameters in the transversal directions). The domain is discretized with a uniform Cartesian mesh with a grid size equal to 1/32 of the sphere diameter. Dirichlet and Neumann boundary conditions are considered for the inlet and outlet boundaries, respectively. For the lateral boundaries in transversal directions, periodic boundary conditions are applied. A second order extrapolation scheme, proposed   [35], is adopted for the curved boundary of the sphere.
The numerical setup is validated by comparing the drag coefficient, the length of the recirculating zone and the angle of separation with existing researches for the fluid velocity field. Tests have shown the mesh and domain independence for the flow around the sphere in the range of parameters considered. For example, in Figure  1(a), the drag coefficient C D obtained from our simulation is compared with empirical equations (equations 5 and 6) of Clift et al. (1978) [1] and with the numerical results of Johnson and Patel (1999) [6]. The drag coefficient deviates from the empirical equations maximum at Re = 25, with relative error 3.5%, which is further reduced with higher Re, e.g. less than 1% at Re = 200. Figure 1(b) presents the results of wake length L W along with numerical results of Johnson and Patel (1999) [6], Tomboulides and Orszag (2000) [5], and experimental data of Taneda (1956) [3], which reported transition to

Results on Spatial Structure of Steady Wake
Our work focuses on the wake behind a wet sphere in the steady axisymmetric regime (0 ≤ Re ≤ 220) and the steady oblique regime (225 ≤ Re ≤ 285). The difference in the overall features of these regimes can be appreciated from Figure 2, which visualizes the streamwise velocity u together with the contours of two advected scalar fields θ 1 and θ 2 of different scalar diffusivities in two perpendicular planes (z, x) and (y, x) passing through the center of sphere in parallel to the incoming flow. The Schmidt numbers for the scalars are 0.71 and 0.61, respectively, which correspond to the diffusivities of temperature and water vapor in air. The increase in Re features the thinning of the boundary layer, as well as a shrinking in the lateral extent of the wake and a stretching in the streamwise direction as in Figure 2 up to Re = 220. In the oblique regime, a tilt from the centerline (y = z = 0) along the (y, x) plane is observed, which is symmetric along (z, x) plane, see also [6,23]. This tilt in the oblique regime increases with Re until the wake becomes unstable and starts shedding vortices at Re ≥ 290. The apparent decrease in the streamwise length of the wake in the top panel of Figure 2 from Re = 225 to 275 is attributed to the tilting of the wake. The transport of any scalar θ is described by the same equation (3). The only difference lays in their Schmidt numbers, which governs their relative diffusivity. The different diffusivities govern the profiles of the scalars at the intermediate values of the dimensionless concentration, which shows differ in the external part away from the sphere boundary and in the far wake (for θ ∼ 0.2 to 0.4), as shown in Figure 2. Due to the Schmidt numbers, the gradient of θ 1 is less steep than the gradient of θ 2 . This feature, however, becomes less distinctive at higher concentrations of θ 1 and θ 2 near the sphere surface.
In order to provide a synthetic description of the flow field, we use a population density approach. For any variable, such as the longitudinal velocity component u, its population density distribution N (u) at a u 0 magnitude is defined as N (u 0 ) = dV u (u 0 )/du, where V u (u 0 ) is the volume of the region in which u is lower that u 0 . The distribution of u is shown in Figure 3 for three different Reynolds numbers (Re =75, 175 and 275). Figure 3  visible in the downstream zone, which has mostly negative p. Tilting is also observed in Figure 3(b) for Re = 275. Figure 3(c) and (d ) present the population density distribution N (u) of the longitudinal velocity component u in these two zones, computed along the entire orthogonal (z, x) and (y, x) planes of the computational domain, respectively. The distribution has been determined by dividing the range of u in 1000 bins, a resolution which allows for a smooth sample distribution while preserving its trend. In the upstream zone (bottom sets of curves in Figure 3(c,d )), N (u) shows a sharp decrease in population density as u decreases from the external ambient value of 1 towards the no-slip zero boundary condition at the sphere surface following a Lorentzian function in equation 7. Some sample population with u ≥ 1 is also observed which resembles the region of highest velocity magnitudes near the p ∼ 0 contour line. In order to create a visible scale separation, the N (u) of the downstream zone is shifted for the middle set of curves in Figure 3(c,d ). Negative values of velocity identify the recirculation zone behind the sphere. A large extent of the simulated wake can be well fitted by a Lorentzian distribution. The crescent like trend right after the ambient u = 1 is a result of the finite size of the simulation domain. Similar to the N (u) of the upstream zone, some sample population with u ≥ 1 is also observed in this downstream distribution, which are also coming from the p ∼ 0 region. N (u) of the entire plane is shifted for the top sets of curves in Figure 3(c,d ) with an amplification of its original magnitudes. As plotted in the insets, the two highest peaks at u ∼ 1 of the entire plane are the individual contributions from both the upstream and the downstream populations.
The Lorentzian or Lorentz-Cauchy distribution y(u; A, u c , b, y 0 ) is a single peak bell-shaped curve, defined as where y(u; A, u c , b, y 0 ) is the population density of samples of variable u, A is its integral over all possible values of u, u c is the position of its maximum where y takes the value 2A/(πb), with b being the width between its half maximums. Parameter y 0 is just an offset value. In the distribution of u, Figure 3(c,d ), a Lorentzian trend is observed in the intermediate range, which corresponds to the boundary layer and to the region external to the wake. An increase in N (u) is observed with increasing Reynolds numbers, indicating an increase in the dimensionless kinetic energy in this region. The out of plane tilting induced by the oblique wake at Re = 275 produces small spikes on top of an overall Lorentzian trend of the sample population along the (y, x) plane, as seen in Figure 3(d ). However, the oblique wake regime retains a symmetric structure along the (z, x) plane in our simulations but the out of plane tilting impacts the sample population. Therefore, N (u) in Figure 3(c) for Re = 275 only indicates a lower yet a smooth Lorentzian trend.
The existence of such a trend in the distribution of a variable indicates the existence of a matching region where the variable shows an algebraic variation from the values in the wake to the values in the external ambient. If the flow is axisymmetric and the flow structures are elongated in the streamwise direction, this variation is in the radial direction proportional to (y 2 + z 2 ) −1 (inverse of the square of the lateral distance from the axis). This algebraic matching region is not only present in the velocity field, but also in the associated pressure field and in the passively transported scalars. Figure 4 presents the spatial distribution of pressure p and transversal component of velocity v along the orthogonal (y, x) plane for various Re. In the axisymmetric regime, in Figure 4(a), the modulus of v is symmetric across the y = 0 plane but not the v. Similarly the modulus of w is also symmetric across the z = 0 plane in the axisymmetric regime, but not w. Complexity arises in the oblique regime, as neither p nor the modulus of v remains symmetric in the Figure 4(b). This is also seen in the population density distribution of v in Figure 4  v show dominance similar to Figure 4(b). In contrast to v, the positive and negative magnitudes of p are rather concentrated near the sphere respectively in the upstream and the downstream zones as in Figure 4(a,b). Similar to Figure 3(c,d ), the N (p) of the upstream zone (p ≥ 0 population) does not show significant variability with Re and exhibits a Lorentzian distribution. The N (p) of the downstream zone however shows local peaks at around p = −0.1, which marks the discontinuity in the sample population in Figure 4(a,b). A three dimensional spatial structure of the velocity components v and w for the oblique Re = 275 and the axisymmetric Re = 175 cases are shown in Figure 5, where the complexity in the oblique wake flow structure can be appreciated. Figure 6 presents the population density distribution of the scalar fields N (θ 1 ) and N (θ 2 ) across two central orthogonal planes (z, x) and (y, x) (similar to the previous Figures). Since the boundary conditions for the dimensionless scalars have a zero value in the ambient and a unit value on the sphere surface, their population density distribution shows the highest population around zero in Figure 6, followed by a domain induced crescent zone, and then a Lorentzian distribution in the intermediate values gradually approaching the surface unit value. The Lorentzian trend is again visible in the scalar population density, due to the similitude of the advection-diffusion equation of the scalars to the dynamics of momentum in regions with small pressure gradients. In the upstream region, the behaviour of velocity and scalars is very different due to the strong pressure gradient, while in the downstream region the difference is much milder. A closer look to the density distribution in the insets shows that the steady axisymmetric cases do not show a well distinguishable difference in the number density at different scalar magnitudes with the increase in Re, but only the threshold magnitude for the start of the Lorentzian trend increases. The shift in the threshold of Lorentzian trend is attributed due to the finite and a similar size of the simulation domain for all the Re cases and due to the shrink in the lateral extent of the wake but a stretch in the streamwise direction with increasing Re. The decrease in the sample population for the oblique cases in the left panel of Figure 6 for the orthogonal (z, x) plane is however due to the out of plane tilt of the wake which reduces the sample population. Whereas in the right panel for the orthogonal (y, x) plane, we see a step-wise perturbation on top of an overall Lorentzian trend in the oblique wake regime as a result of its tilt in this plane. Figure 7 presents the spatial distribution of the convective scalar fluxQ in the streamwise direction x, which is a product between θ and u. Spatial distribution ofQ along the orthogonal (y, x) plane in Figure 7(a,b) is someway different from the other flow quantities, since it shows highest positiveQ in the boundary layers and a negativeQ in the recirculating zone due to negative u. The non-symmetric spatial structure of the oblique (Re = 275) scalar flux is visible in Figure 7(b). The population density distribution N (Q) along the orthogonal (z, x) and (y, x) planes shows a different structure as expected. A Lorentzian trend is observed for a few limited sample populations, for example, for the samples between the white and pink contour lines in Figure 7

Discussions and Concluding Remarks
We present a detailed numerical analysis on the spatial structure of the wake flow using population density distribution for various Reynolds number in the steady wake regime. Incompressible Navier-Stokes equation is solved for the flow velocity and the one-way coupled advection-diffusion equations are solved for the scalars using the lattice Boltzmann method. The spatial evolution of various flow quantities, such as, longitudinal velocity component u, pressure p, passive scalar θ, convective scalar fluxQ in the wake of steady axisymmetric regime (Re ≤ 220) and oblique regime (225 ≤ Re ≤ 285) using a population distribution function N , shows a Lorentzian distribution which is proportional to the inverse of the square of the flow quantity (for example, N (p) ∝ p −2 ). This Lorentzian trend exhibits an algebraic decay in the number density of populations with different magnitudes of fluid quantities from the external ambient to the boundary layer in the wake and dominates the spatial distribution of the flow quantities outside the recirculating region. The horizontal components of fluid velocity, v and w, whereas show different spatial distributions not attributable to a Lorentzian one. Transition to the oblique wake regime at Re ≥ 225 in our simulations shows a complex three dimensional spatial evolution of the flow quantities, which also shows a Lorentzian trend. The population density distribution for the longitudinal velocity component u, shows an increase in its number density with increasing Re. Whereas the number density of the scalar populations remains the same for various axisymmetric Re. This feature however changes in case of the convective scalar flux, where an increase in its number density is observed again with an increase in Re.
Quantification of scalar transport in the wake of spherical objects is important for understanding various physical phenomena. For example, Bhowmick et al. (2020) [27], and also in Chouippe et al. (2019) [23] and Krayer et al. (2020) [28], scalar transport in the wake is used to understand the spatial distribution of supersaturation in the wake of precipitating cloud hydrometeors. By scaling the passive scalars as the temperature and the water vapor density fields around the droplets, we used three dimensional population distribution of supersaturation also in Bhowmick et al. (2020) [27] to quantify the supersaturated volume produced in the wake of precipitating cloud hydrometeors in presence of a sufficient temperature gradient in a slightly subsaturated cloudy ambient, which can activate cloud aerosols and thus contribute to the cloud life cycle.