{'fraction_non_alphanumeric': 0.036938569989929504, 'fraction_numerical': 0.015347432024169185, 'mean_word_length': 4.479143676892519, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Lattice Boltzmann simulations of droplet breakup in confined and time-dependent flows

We study droplet dynamics and breakup in generic time-dependent flows via a multicomponent Lattice Boltzmann algorithm, with emphasis on flow start up conditions. We first study droplet breakup in a confined oscillatory shear flow via two different protocols. In one set up, we start from an initially spherical droplet and turn on the flow abruptly ("shock method"); in the other protocol, we start from an initially spherical droplet as well, but we progressively increase the amplitude of the flow, by allowing the droplet to relax to the steady state for each increase in amplitude, before increasing the flow amplitude again ("relaxation method"). The two protocols are shown to produce substantially different breakup scenarios. The mismatch between these two protocols is also studied for variations in the flow topology, the degree of confinement and the inertia of the fluid. All results point to the fact that under extreme conditions of confinement the relaxation protocols can drive the droplets into metastable states, which break only for very intense flow amplitudes, but their stability is prone to external perturbations, such as an oscillatory driving force.


I. INTRODUCTION
Fluid dynamics phenomena, involving droplet dynamics, deformation, and breakup, are prominent in the field of microfluidics and even in general complex flows at larger scales. Beyond the practical importance in a variety of concrete applications [1][2][3][4], they are also relevant from the theoretical point of view, due to the complexity of the physics involved [5][6][7][8][9]. Droplet deformation is characterized via the capillary number, where η s is the dynamic viscosity of the solvent, R the radius of the initially undeformed spherical droplet, σ the surface tension and G the shear rate intensity [6,10]. The value of Ca at break up is denoted by Ca cr , the critical capillary number. A lot of attention has been dedicated to droplet deformation and breakup in stationary flows [5,11,12], and, in particular, the effect of the degree of confinement on the flow dynamics [13][14][15]. The degree of confinement is parametrized by the ratio α = 2R/L, where L denotes the shear wall separation. Confinement is frequently encountered in experimental setups of droplet dynamics in simple shear flows [13][14][15][16][17][18][19][20][21][22][23][24] and can be enhanced by changing α. There are some theoretical models which were developed to capture the experimental phenomenology of confined droplet dynamics, analytical models [25,26], which extended the theory by Taylor [6,27], and phenomenological models [28][29][30]. The validity of the analytcial models were verfied in Ref. [31] and the phenomenological models in Ref. [32]. Of particular interest are the results in Ref. [14], which show that, for non vanishing α breakup differs substantially from the unconfined shear case both qualitatively and quantitatively for all viscosity ratios χ = η d /η s , where η d,s is the dynamic viscosity of the droplet (d) or solvent (s) phase. Additionally, the dependency of the critical capillary number Ca cr on the droplet's inertia is a central area of interest [8,, with the degree of inertia being given by the Reynolds number, where ν s is the kinematic viscosity of the solvent. Furthermore, breakup is influenced by the start up conditions, as demonstrated in experimental and theoretical studies [54][55][56][57][58]. This phenomenon is rather subtle and can have different effects depending on the protocol in use. The dependency on the rate of increase of the shear rate G was confirmed by [54] via supporting calculations based on the model by Taylor [6]. A theoretical model developed by Hinch et al. [55] shows that stable droplet equilibria below the critical capillary number Ca cr breakup are only possible for a sufficiently low increase in G. Furthermore, Renardy [58] has shown that although these stable equilibria require a slow increase in the shear rate G they are unique and do not depend on the rate of change of G. We stress that even though the effect of start up conditions on break up has been investigated [54][55][56][57][58], the role of confinement with varying start up conditions on droplet dynamics and break up is not clear. Moreover, it is unclear how break up is affected, if the flows are time-dependent [59][60][61][62][63]. The aim of the present paper is to take a step further in this direction. With the use of numerical simulations we show that at capillary numbers close to breakup, confinement allows for the existence of a metastable flow configuration next to the solution of the Stokes equation found in Ref. [58]. This metastable state is prone to perturbations and collapses to the Stokes solution, if we have a time-dependent flow with a sufficiently large shear frequency. It should be stressed that this result is unique to the case of a confined droplet in an oscillatory shear, as this metastable configuration is not present neither for an unconfined droplet in an oscillatory shear flow nor in the case of an oscillatory elongational flow. Our studies can be seen as an extension to Refs. [34,58], where the influence of inertia on droplet breakup was studied, whereas we deal with time-dependent cases, where the temporal rate of change of the shear intensity is comparable to the droplet relaxation time, This work is a follow up study of [63], where stable time-dependent droplet dynamics was investigated via a multicomponent lattice Boltzmann scheme and a phenomenological model [28,29]. It was found that droplet deformation depended strongly on an external timescale, the oscillation frequency of an oscillatory shear flow, for a confined droplet. For relatively large oscillation periods close to the value of t d the droplet is hardly deformed by the solvent shear flow, which was described as the "transparency effect" in Ref. [63]. The findings in Ref. [63] have been validated by comparing the lattice Boltzmann results to the results obtained via a phenomenological droplet deformation model, the Maffettone-Minale model [28,29]. This paper is organized as follows: Sec. II gives a brief overview on the lattice Boltzmann algorithms and models in use. In Sec. III, we outline the general details of droplet break up with an emphasis on confined systems and simple shear flows. In Sec. IV, we investigate break up in a time-dependent (oscillatory) shear flow under strong confinement. A mismatch between two protocols, involving different start up conditions of the flow, leads us to investigate break up conditions under the influence of inertia (Sec. V) and the effect of confinement (Sec. VI) Moreover, we check whether the protocol mismatch depends on the flow topology (Sec. VII).

II. LATTICE BOLTZMANN ALGORITHMS AND METHODS
We use lattice Boltzmann simulations [64,65] to study droplet break up in confined and time-dependent shear and elongational flows. The lattice Boltzmann method (LBM) has been extensively used in the field of microfluidics, including extensions to accommodate nonideal effects [66], coupling with polymer micro-mechanics [67] and thermal fluctuations [68,69]. LBM has also been used widely for the modeling of droplet break up behavior [41,44,[70][71][72][73][74][75][76][77]. 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 multicomponent model (SCMC) [78,79], a diffuse interface model in the framework of the LBM. The hydrodynamical quantities, mass and momentum densities, can then be described as: where g σ i (x, t) denotes the populations in the LBM model for the fluid component σ and c i are the lattice velocities. For example, for a two component system with species A and B the index σ can take the values σ = A and σ = B. The interaction at the fluid-fluid interface [80,81] 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 open flow boundary conditions as outlined in Ref. [63]. 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 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 simulations domain boundary. Thus the ghost distributions update the boundary nodes during the LBM streaming step and effectively simulate an open flow boundary given by the chosen density ρ b (x, t) and velocity u(x, t) fields of the outer fluid [63]. 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 timescale) the collision operator is linear and contains several relaxation times linked to its relaxation modes (depending on the lattice stencil) [82]. One relaxation time τ is directly linked to the kinematic viscosity ν in the system which is one of the primary links between the LBM scheme and hydrodynamics [64,65]. The boundary scheme described here is not strictly mass conserving, so we correct the local population mass densities To cure mass conservation [83][84][85]. This is not the case in unconfined system, where we can accept small mass fluctuations of both droplet and solvent, but have to reinject mass into the droplet [86].

III. SIMULATION SETUP AND DEFINITIONS
In this section we define what we mean when we speak of droplet break up and characterize the simulation setups. We deal with both a confined droplet in a simple shear flow and an unconfined droplet in a uniaxial extensional (elongational) flow. The velocity gradient matrix for both shear and elongational flows is given by where ∇v = G and β is a parameter characterising the flow type. The shear flow setup is equivalent to the one used in Ref. [63] with β = 0 in Eq. (9) except that the flow is unconfined and elongational with an oscillatory velocity gradient amplitude G(t) given by Eq. (9) with β = 1. Droplet deformation can be characterized by the capillary number Ca. In the case of a shear flow including confinement the shear rate is given by with L z being the channel width responsible for the droplet confinement and u 0 being the maximum wall velocity amplitude. This definition may also be extended to time-dependent shear flows [63] In accordance with [14] we define the critical capillary number Ca cr as the value of Ca for which an initially spherical droplet breaks up, which is achieved by a sudden increase in the shear rate amplitude G. We refer to this break up protocol as the Shock Method. In addition we can gradually increase the shear rate G starting from a value for which the droplet is only marginally deformed [54,55,58]. A fixed increase ∆G (or ∆u 0 in the case of Eq. (11)) is equivalent to a fixed increment rate ∆Ca for the capillary number. This way the droplet and the solvent flow are given more time to relax to their respective equilibrium distributions at specific Ca. We call this protocol the Relaxation method. It should be stressed that breakup in the relaxation method has a small dependency on ∆Ca. If ∆Ca is very large, e.g. ∆Ca ∼ Ca cr , the value for Ca cr will be the same as the one obtained through the shock method. Thus, the ∆Ca has to be chosen sufficiently small enough for the relaxation method to work. Essentially, the relaxation method captures the deformation history of the droplet before breakup with an accuracy given by ∆Ca contrary to the shock method. The relaxation method is especially important for droplet dynamics in palatially evolving shear flows in the case of a smoothly varying local shear both spatially and temporarily. A variation of the relaxation method for time-dependent oscillatory flows, i.e. where the shear amplitude G(t) = G 0 cos(ωt), is to consider the flow and droplet configuration at a capillary number Ca close to Ca cr and then to increase the oscillatory shear frequency ω f = ω/(2π) until break up, starting from the stationary case of ω f = 0. As in Ref. [63] we use a dimensionless frequency ω f t d in our discussion, where t d is the droplet relaxation time defined in Eq. (3). In the presence of a flow with nonzero frequency ω f t d , we focus on Ca max , which denotes the maximum value of the time-dependent capillary number Ca(t) over one oscillatory cycle [63]. An instance of droplet break up in an oscillatory simple shear flow is depicted in Fig. 1. The droplet is oscillating between two maximally elongated states for Ca < Ca cr and breaks up during the flow build up for Ca > Ca cr in the case of the shock method. The droplet elongation is characterized by the droplet length L(t), which is defined as the longest axis of the elongated droplet, and L cr denotes the droplet length in the critical case Ca ≥ Ca cr . The time evolution of L(t) is also shown for the two cases Ca < Ca cr and Ca > Ca cr in Fig. 1, which shows that break up occurs at around t = 17000 lbu with lbu denoting lattice Boltzmann Units. In all simulations in this article the viscous ratio χ ≡ 1 and the density ratio ρ d /ρ s ≡ 1. If not explicitly stated otherwise, then the confinement ratio for simple shear flows α ≡ 2R/L z , where R is the radius of the spherical undeformed droplet and L z the channel width, is set to α = 0.75.

IV. DROPLET BREAK UP IN AN OSCILLATORY SHEAR FLOW
Similarly to Ref. [63] we consider a droplet in a confined oscillatory simple shear flow; see Fig. 2. The setup is shown in Fig. 1 with a confinement ratio α = 0.75 and a time-dependent shear rate G(t) = 2u 0 /L z cos(2πω f t), where ω f is the frequency of the outer oscillatory flow [59][60][61][62][63]. Our main focus is the dependency of Ca cr on the normalized shear frequency ω f t d of the oscillatory outer flow. Droplet dynamics in oscillating flows may feature a so called transparency effect [63], which states that the droplet is hardly deformed if ω f t d ∼ 0.1, i.e. the timescale of the oscillating shear flow 1/ω f is of the similar order as the droplet relaxation timescale t d . The droplet dynamics are hardly influenced by the shear frequency for ω f t d ∼ 10 −4 and the transparency effect is noticeable for ω f t d ∼ 10 −2 and higher frequencies, which leads to a sudden increase in the critical capillary number. To stay in tune with experimental results [13,14,21,22], we limit the range of the critical capillary number close to Ca cr ∼ 1.0. In Fig. 3 we can see that the droplet break up behavior is significantly different for our two LBM simulation protocols, the shock and relaxation method. The shock method implies that droplet break up is independent of the oscillatory shear frequency ω f t d , significant changes in Ca cr only occur close to the transparency effect region at high frequencies (ω f t d ∼ 10 −2 ). The relaxation method is of a different nature: first of all Ca cr in the low-frequency region (ω f t d ∼ 10 −4 ) is larger than the values obtained with shock method (see also Sec. VI). Moreover, for intermediate frequencies ω f t d ∼ 5 × 10 −3 we observe that break up occurs at a significantly smaller Ca cr than in the low-frequency range and is now of a comparable value to Ca cr obtained via the shock method. The mismatch between the two protocols in the low-frequency regime in Fig. 3 is in disagreement with previous studies of start up conditions of droplet break up in confined simple shear flows [14,58]. However, the shock method produces results in accordance with the literature [14], as the dashed line in Fig. 3 indicates. It should also be noted, that the destabilization of the "relaxation branch" is rather sudden and takes place at very small ω f t d . This suggests that the protocol mismatch is due to metastable solution (relaxation method) existing next to a stable solution (shock method) in the low-frequency range ω f t d ≤ 0.02. The protocol mismatch seems rather puzzling: according to Renardy [58] the solution should be unique. However, our setup differs in a few points from the one in Renardy [58]. First of all, the droplet is strongly confined (α = 0.75) in our setup (see Fig. 1), which could have a strong effect on the values Ca cr for varying start up conditions. Moreover, inertia might stabilize the droplet in the case of the relaxation method. Therefore, the protocol mismatch might disappear in the Stokes limit. In addition, one may also wonder what is the effect of flow topology, as an inherently different flow field might lead to a similar protocol mismatch. Given these considerations, in the following sections, we will investigate the cause of the mismatch by considering both inertial effects, as is the case in Ref. [58], (see Sec. V) and the importance of confinement in stationary shear flows (see Sec. VI). Regarding the importance of flow topology, we investigate time-dependent break up in an elongational flow in Sec. VII.

V. INERTIAL EFFECTS
In Ref. [58] it is shown that the solution of the Stokes equation in confined simple shear flows is unique and does not depend on neither the initial conditions of the droplet nor the solvent flow configuration. Thus, one might think that the protocol mismatch might be due to inertial effects and would disappear, if we were close to the time-dependent Stokes limit of Re ≡ 0. Interestingly, the LBM formalism allows us to directly set Re = 0, as we can eliminate the nonlinear terms in the equilibrium distribution functions in the LBM algorithm, Eq. (6), which leads us to a modified equation ( 12), accounting only for the linear terms in the velocity field u(x, t). It should be remarked that only the nonlinearites of the Navier-Stokes equation are removed in this way, since the inertia embedded in the time derivative of the velocity field u(x, t) does not disappear and may still play a role during the non steady break-up process. Inertial effects tend to stabilize the droplet [56,57] for low Re < 1, whereas Ca cr ∼ 1/Re for large Re > 10 [34]. This suggests, that the stabilization effect of low Re are responsible for the protocol mismatch, which consequently should disappear in the Stokes limit Re = 0. We investigate the dependency of Ca cr on Re, as shown in Fig. 4. For the case Re = 0 we use only the linear terms of the equilibrium distribution functions given by  There is a mismatch between the predictions of the two break up protocols. Whereas droplet break up is largely independent in the case of the shock method, except for the asymptotic behavior in the high-frequency region, the relaxation method in the low-frequency limit predicts a higher Cacr than the ones of the shock method. This mismatch is investigated in the article. The error bars are estimated via steps in the critical capillary number ∆Ca. Both curves are interpolated via bezier curves.
The simulations are carried out for a stationary shear flow, with the setup described in Fig. 1. We can see that the mismatch between the break up protocols, does not depend on inertia and is even present in the Stokes limit of Re = 0. We conclude that the mismatch between the two break up protocols is not influenced by any stabilization effects of inertia [56,57] for the given range of Reynolds numbers Re ∼ 0.0, . . . , 1.5.

VI. CONFINEMENT EFFECTS
We now focus on both confinement and start up conditions in the shear rate amplitude G for droplet break up in a stationary shear flow. The setup is once again the one in Fig. 1, a confined droplet in a stationary (ω f t d = 0) shear flow, but now we vary the confinement ratio α and, in the case of the relaxation method, the rate of change of the shear amplitude G, resulting in increments of the capillary number ∆Ca. Our results are summarized in Fig. 5. We can see, as was shown in Ref. [58], that the critical capillary number Ca cr is independent of the start up conditions for low confinement ratios (α ≤ 0.5), as both the shock method and the relaxation method yield the same results with respect to the simulation errors. However, if the droplet is strongly confined (α ≥ 0.6), then the two methods yield very different results, with the Ca cr predicted by the relaxation method being substantially larger than the one predicted by the shock method. It should be noted, that Ca cr is independent of ∆Ca, given that ∆Ca is small enough, which can be seen from Fig. 5, where the values of Ca cr overlap in respect to their error ranges for different ∆Ca and the same α. Fig. 6 shows the length of the elongated droplet as a function of the LB simulation time for the different shear start up methods: we can see that for the shock method droplet break up occurs soon after the maximal elongation, whereas for the relaxation the droplet experiences a sequence of maximal extensions and subsequent retractions after breaking up for a given Ca cr at its critical length L cr (t). We conclude that both a slow start up of the outer flow (relaxation method) and a strong confinement of the droplet (α ≥ 0.6) are necessary for the mismatch reported in Fig. 3 in the low-frequency limit. The eventual collapse of the relaxation method solution on to the one found by the shock method suggests, that the relaxation method branch in the low-frequency limit in Fig. 3

VII. FLOW TOPOLOGY
We now investigate the protocol mismatch in terms of the flow topology. Instead of an oscillatory shear flow, we consider break up in an elongational (or uniaxial extensional) flow; see Fig. 7. This flow is by its very nature unconfined, so we would expect to not see a mismatch, as is the case for α = 0 in the case of the confined shear flow; see Sec. VI. The results are shown in Fig. 8. Interestingly, a mismatch between the two droplet protocols is absent and the predictions agree well with each other in terms of their respective errors. This shows that strong confinement (α ≥ 0.75) is necessary for the existence of the protcol mismatch shown in Fig. 3. Moreover, Fig. 8, shows that droplet break up in an oscillatory elongational flow is frequency dependent, with an exponential dependence between the oscillation frequency ω f t d and the critical capillary number Ca cr . The low-frequency limit matches the stationary flow predictions of [39].

VIII. CONCLUSIONS AND OUTLOOK
We have shown that the interplay of varying start up conditions and strong confinement ratios can lead to qualitatively and quantitatively different droplet break up conditions in stationary shear flows, unlike the stable equilibria found for varying start up conditions [58] or the ones found for varying degrees of confinement [14]. Having investigated the effects of inertia, confinement and flow topology, we conclude that the protocol mismatch between the shock and the relaxation method are due to a high degree of confinement for a droplet in a shear flow (α = 0.75). However, the break up solution found via the relaxation method is only metastable, since it becomes unstable in the case of a time-dependent, oscillatory shear flow. The protocol mismatch is thus solely due to an extra metastable solution in a strongly confined shear flow and disappears in the presence of small perturbations (e.g., amplitude variations in an oscillatory shear flow) in accordance with the uniqueness of the Stokes solution [14,58]. We have also shown the dependency of the critical capillary number Ca cr on the normalized oscillation frequency ω f t d in both oscillatory shear and elongational flows. In the case of the elongational flow, Ca cr increases with increasing ω f t d , whereas no simple functional dependence can be found for the oscillatory shear flow, since Ca cr also depends on the flow start up and degree of confinement. We should stress again that the results presented in this work are only valid for χ = 1, since the viscosity ratio influences the breakup of a confined droplet [9,14,21,24]. On the one hand, for small viscosity ratios χ ≈ 0.3 the confined shear flow stabilizes the droplet and breakup is more difficult to occur than for χ = 1. On the other hand for large viscosity ratios χ ≈ 5.0 the confined droplet is destabilized and breakup is more likely to happen than for χ = 1. [14,24]. It would be interesting to see whether the metastable solution can be found in an experimental setup or whether it is too prone to perturbations to manifests itself.