Sub-Kolmogorov droplet dynamics in isotropic turbulence using a multiscale lattice Boltzmann scheme

The deformation and dynamics of a single droplet in isotropic turbulence is studied using a Lattice Boltzmann diffuse interface model involving exact boundary flow conditions to allow for the creation of an external turbulent flow. We focus on a small, sub-Kolmogorov droplet, whose scale is much smaller than the Kolmogorov length scale of the turbulent flow. The external flow field is obtained via pseudo-spectral simulation data describing the trajectory of a passive tracer in isotropic turbulence. In this way we combine the microscopic scale of the droplet and the macroscopic scale of the turbulent flow. The results obtained from this fully resolved model are compared to previous studies on sub-Kolmogorov droplet dynamics in isotropic turbulent flows, where an analytical model, which assumes the droplet shape to be an ellipsoid at all times, is used to describe the droplet deformation. Our findings confirm, that the hybrid pseudo-spectral Lattice Boltzmann algorithm is able to study droplet deformation and breakup in a regime well beyond the ellipsoidal approximation.


Introduction
The dynamics and breakup behaviour of immiscible droplets in laminar flows have been studied extensively in the literature, from applied studies involving concrete setups [4,5,6,7] to theoretical studies addressing their physical complexity [8,9,10,11,12]. A lot of attention has been dedicated to droplet deformation and breakup in stationary laminar flows [8,13,14,15,16,17,18,19,20,11,10,21,22], but less is known about the dynamics of droplets in homogeneous isotropic turbulent flows, which pose the challenge of solving a multi-physics problem, since the flow properties of the turbulent scales have to be accurately transferred to the scale of the droplet. We can differentiate between two interesting cases for the deformation of droplets in homogeneous and isotropic turbulent flows [23]: On the one hand, there are studies on droplet dynamics on scales larger than the Kolmogorov scale [24,25,26,27,28]. On the other hand, we investigate the dynamics and breakup statistics of sub-Kolmogorov droplets, i.e. droplets whose size is much smaller than the Kolmogorov scale [29,30,31,2,32,33]. The latter case is of particular interest to the authors, since the viscous stresses dominate over the inertial stresses, if the droplet size is smaller than the Kolmogorov scale of the outer flow [29]. In this regime the droplet is deformed similarly to the laminar case and the turbulence is only noticeable via the time-dependency of the solvent flow. Both [2] and [33] make use of the Maffettone-Minale model (MM-model) [3] to model the droplet in homogeneous isotropic turbulence, which only accounts for ellipsoidal deformations and thus can model breakup only via a cut off criteria. Even though there exist several computational studies on sub-Kolmogorov droplet dynamics and breakup in turbulent flows [29,30,31,2,32,33] there exists no experimental evidence of sub-Kolmogorov droplet break up [23], as such experiments are difficult to set up. However, there are experimental studies on breakup in flows through dilute random beds of fixed fibers [34] and breakup of a small single droplet in a chaotic advective flow [35]. These results [34,35] show that small droplets, whose size is comparable to the size of a sub-Kolmogorov droplet in a realistic turbulent flow, can break up in highly unsteady and time-dependent laminar flows, indicating that breakup of sub-Kolmogorov droplets is indeed possible and not merely a computational artefact. Moreover, our studies employ a novel multiscale algorithm, as we not only use the MM-model to measure the deformation of the droplet but compare it with fully resolved LBM (Lattice Boltzmann Method) simulations using a diffuse interface method, the Shan-Chen multicomponent model (SCM) [36,37].
But how can we couple LBM droplet dynamics to fully developed turbulence?
One way in doing so is to use a LBM hybrid boundary scheme, which is able to simulate an exact turbulent flow at the domain boundary. Exact flow boundary conditions for LBM are a very useful tool for multiscale physics simulations, as it is enabling scale separation in its very own nature. Exact flow boundaries for LBM were first developed by [49] for 2D lattices and then extended to 3D flows by [50,51,52] for perpendicular inlet and outlet flows. In [53] it is shown that exact boundary flow conditions can be implemented in 3D LBM simulations for the D3Q19 lattice. Even though the scheme presented in [53] seems rather suitable in modelling a sub-Kolmogorov droplet in homogeneous isotropic turbulence, we use a newly developed ghost boundary flow method. The main reason for this is that the method in [53] cannot determine both the density ρ b (x, t) and the velocity field v(x, t) at the lattice domain boundary. Since we require a constant density ratio of ρ d /ρ s , where the subscripts d and s denote the droplet and solvent phase respectively, for our hybrid simulations, we use the ghost node boundary scheme, which does not set constraints to either the density nor the velocity field at the boundary. Our numerical method can be seen as a novel hybrid approach to the sub-Kolmogorov droplet dynamics: The simulation domain can be of the same order as the droplet's size, see section 4, but should be ideally much larger than the droplet radius R, see section 5. The surrounding turbulent flow is used as an input via an open boundary condi-tion method [1], similar to the setups in [50,53]. This enables us to study the deformation and breakup behaviour of a sub-Kolmogorov droplet in fully developed isotropic turbulence, which we can compare to the results obtained via the widely used MM-model. We find that the fully resolved PS LBM hybrid simulation scheme is crucial to model sub-Kolmogorov droplet breakup in homogeneous isotropic turbulence, since it allows for non-ellipsoidal deformations close to breakup in contrast to the MM-model.

The MM-model and the Lattice Boltzmann algorithm
The phenomenological MM-model describes the coupling of droplet dynamics to a solvent flow. The droplet is always assumed to be ellipsoidal, so that we can describe it via a second rank tensor M ij , which is also referred to as the morphology tensor. Droplet deformation is characterised via the components of M ij (e.g. for an undeformed droplet M ij = δ ij ). The time evolution of M ij due to an external flow field is given by the MM equation [3]: where S ij and Ω ij are the symmetric and anti-symmetric parts of the velocity gradient G ij ≡ ∂ j v i respectively: are the third and second tensor invariants of M ij [3]. The degree of potential droplet deformability is given by the capillary number, which is defined as the ratio of viscous and interfacial forces of the droplet where G(t) ≡ G ij (t) ∞ is the time-dependent shear rate of the solvent flow, R the radius of the initially undeformed droplet, σ the surface tension between the droplet and the solvent and η s the dynamic viscosity of the solvent flow.
Unlike previous analytical approaches to model droplet deformation [9,54]  The droplet deformation in the MM-model is given by where L and W are the major and minor ellipsoidal axes respectively. Since we want to compare the LBM simulations with the MM-model, we need to make sure that we remain in the linear flow regime and check that our deformed LBM droplet is actually ellipsoidal at all times. The MM-model includes two parameters f 1 and f 2 , which are dependent on the viscous ratio χ = µ d /µ s , with µ d and µ s denoting the dynamic viscosities of the dispersed and solvent phase respectively, and the capillary number Ca [3]: In addition to the phenomenological MM-model we can use fully resolved Lattice Boltzmann simulations [55,56] to study sub-Kolmogorov droplet dynamics and breakup in homogeneous and isotropic turbulent flows. The Lattice Boltzmann Method (LBM) has been extensively used in the field of microfluidics, including extensions to accommodate non-ideal effects [57], coupling with polymer micro mechanics [58] and thermal fluctuations [47,48]. LBM has also been widely used for the modelling of droplet breakup behaviour [31,38,39,40,41,42,43,44,45,46]. In order to model multicomponent systems with the Lattice Boltzmann Model (LBM) we need to account for interfacial forces between different fluid components. This can be achieved with the Shan-Chen Multi-Component model (SCMC) [36,37], a diffuse interface model in the framework of the LBM. The hydrodynamical quantities, mass and momentum densities, can then be described as: The interaction at the fluid-fluid interface [59,60] is given by: where ρ σ (x, t) is the density field of the fluid component denoted by σ. G σ,σ is a coupling constant for the two phases σ and σ at position x and w i are the lattice isotropy weights. We use the same exact flow boundary conditions as outlined in [1]. In order to use arbitrary boundary values of the density ρ(x, t) and velocity u(x, t) fields of the solvent fluid we use ghost populations (or halos), which store the equilibrium distribution functions g eq i of the boundary density and velocity fields. The equilibrium distribution functions g eq i at the domain boundary are given by: with w i being the lattice weights for the set of lattice vectors c i , and ρ b (x, t) the density field at the domain boundary. Thus the ghost distributions update the boundary nodes during the LBM streaming step and effectively simulate an exact flow boundary given by the chosen density ρ b (x, t) and velocity u(x, t) fields of the outer fluid [1]. The streaming and collision steps are given by the Lattice Boltzmann equation: where Ω({g i (x, t)}) is the collision operator depending on the whole (local) set of lattice populations and ∆t is the simulation time step. For MRT (multi-relaxation time scale) the collision operator is linear and contains several relaxation times linked to its relaxation modes (depending on the lattice stencil) [61]. One relaxation time τ is directly linked to the kinematic viscosity ν in the which is one of the primary links between the LBM scheme and hydrodynamics [55,56]. Since the boundary scheme is not strictly mass conserving, we accept small mass fluctuations of both droplet and solvent, but keep the density ratio ρ d /ρ s constant, which is achieved by reinjecting mass into the droplet [62].

Simulation setup
We have introduced the ghost boundary flow method in section 3 as an exact flow boundary method: the boundary scheme enables us to enforce a density field ρ b (x, t) and u(x, t) at the boundaries at any LBM simulation time step Reynolds number of Re λ = 420. We consider droplets below the Kolmogorov scale, i.e. R η K , as we can limit ourselves to a linear velocity profile in this case [29]. The Kolmogorov scale is defined as where ν s is the kinematic viscosity and ε the energy dissipation of the solvent flow. The typical velocity increments of separation l of size |l| = l are which can be approximated by for l ∼ R 2 . Since the droplet's centre of mass remains fixed, i.e. we are in the frame of reference of the droplet, we obtain the following velocity relation in the case of a single phase flow: which is a linear velocity profile in the distance x from the droplet's centre of mass. With v (0) i (x, t) we denote the turbulent velocity profile imposed by the boundary, whereas v i (x, t) is the LBM obtained velocity profile given by equation (6). It is interesting to compare the numerical algorithm for an external homogeneous and isotropic turbulent flow in fully resolved LBM simulations to the study in Mukherjee et al. [28]. Instead of an exact boundary flow scheme Mukherjee et al. use a turbulent forcing scheme [28,62] to model the dynamics of droplet emulsions in homogeneous isotropic turbulence. In particular this method has the advantage that the boundary conditions can be chosen as being periodic, which does not require the LBM simulations to converge to a given solution as is the case for the exact boundary method we use. However, the droplets of the turbulent emulsions are larger than the Kolmogorov scale and thus the turbulent forcing scheme in [28] is not useful for our study of sub-Kolmogorov droplet dynamics. Following the setup of the exact boundary flow algorithm, we test whether our LBM scheme can relax to the imposed velocity field given in equation (14).

Validation results
Now that we have established the approximations we make to the velocity profile on the sub-Kolmogorov scale, see equation (14), we would like to create a pseudo-spectral LBM hybrid method, which takes the turbulent velocity gradients G ij (t) as an input. In order for the LBM model to yield accurate predictions, we need to simulate a fully developed turbulent flow at the domain boundaries, which then evolves into the centre of the domain. This is achieved by using both the ghost node boundary flow scheme described in section 2, which stores suitable equilibrium distribution functions on ghost nodes, and the pseudo-spectral values for G ij (t). The simulation setup is shown in figure 1 with the length of the simulation box L s and the diameter of the undeformed droplet turbulent flow field, see equation (14), to a given accuracy. Thus, recovering the velocity field in equation (14) effectively means recovering the turbulent velocity gradient G ij (t) of the pseudo-spectral DNS simulations, as the velocity profile is linear. We choose an LBM sampling time ∆t measured in lbu, which allows the LBM algorithm to relax to the velocity value required by equation (14).
In figure 3, we see that our LBM scheme reproduces indeed the turbulent signal G ij (t) we provided it with, since the original rescaled DNS signal for the G xx (t) component converges, with the value for G xx (t) obtained with LBM.
The rescaling for G ij (t) is with the one imposed at the boundary, see equation (14). The corresponding L 2 -error is provided in figure 4 and is defined via: where v(x, y, z, t) and v (0) (x, y, z, t)) are the velocity fields of the LBM simulation and the imposed turbulent linear velocity profile respectively. u (0) rms (t) is the rms value of v (0) (x, y, z, t). We would like to choose a minimal ∆t, as to optimise run time, we choose ∆t = 200 (lbu), since this value yields a time averaged global error of E (0) t < 0.1. This is a reasonable error threshold, since LBM simulations for droplet dynamics in oscillatory shear flows yield a similar maximal error [1]. However, we should remark, that the global error

Simulation results
Since we have established the validity of the hybrid PS-LBM algorithm in section 4, we can investigate droplet dynamics and breakup for a single droplet on the sub-Kolmogorov scale. Analogously to equation (4) we can define a measure for the deformation D LBM for our LBM simulations:   Therefore, it appears that the MM-model is not able to predict the correct critical capillary number Cacr for sub-Kolmogorov droplet breakup.
where Ω(t) is the time-dependent surface area of the deformed droplet and and 8. Figure 7 shows the elongated droplet for a maximal capillary number Ca max ≈ Ca cr , just before and after breakup. Interestingly, the droplet is now elongated and thus resembles droplets close to Ca cr in a pinch-off breakup process found in studies on both confined droplets in laminar flows [15,16,21,22] and computational studies on breakup of sub-Kolmogorov droplets [29,30,31], and breakup of droplets in porous media flows [34]. Figure 8 shows  [30] and [31], for Re = 1 and χ = 1.

Conclusion
The DNS-LBM hybrid scheme for sub-Kolmogorov droplets in fully developed homogeneous and isotropic turbulence has been used to model droplet deformation and breakup. The LBM results were also compared to MM-model predictions. We have found that LBM predicts a significantly lower critical capillary number Ca cr than the MM-model, see