# Dynamics of a droplet in shear flow by smoothed particle hydrodynamics

^{1}State Key Laboratory of Fluid Power and Mechatronic Systems, Department of Engineering Mechanics, Zhejiang University, Hanghzhou, China^{2}Department of Physics, Hangzhou Dianzi University, Hangzhou, China^{3}Hangzhou Shiguangji Intelligient Electronics Technology Co., Ltd., Hangzhou, China

The behavior of a droplet under shear flow in a confined channel is studied numerically using a multi-phase smoothed particle hydrodynamics (SPH) method. With an extensive range of Reynolds number, capillary number, wall confinement, and density/viscosity ratio between the droplet and the matrix fluid, we are able to investigate systematically the droplet dynamics such as deformation and breakup. We conduct the majority of the simulations in two dimensions due to economical computations, while perform a few representative simulations in three dimensions to corroborate the former. Comparison between current results and those in literature indicates that the SPH method adopted has an excellent accuracy and is capable of simulating scenarios with large density or/and viscosity ratios. We generate slices of phase diagram in five dimensions, scopes of which are unprecedented. Based on the phase diagram, critical capillary numbers can be identified on the boundary of different states. As a realistic application, we perform simulations with actual parameters of water droplet in air flow to predict the critical conditions of breakup, which is crucial in the context of atomization.

## 1 Introduction

The deformation and breakup of droplets in shear flow are ubiquitous in engineering applications. On microfluidic chips, droplets are utilized for microbial cultivation and material transport [1,2], and a thorough understanding of their dynamics in confined flows may improve the efficiency of production and transportation. In other environmental and industrial applications such as protection against harmful aerosols, ink-jet printing and atomization in nozzles [3–7], liquid droplets are typically in gas flows. Accordingly, a decent knowledge on their dynamics with a high density/viscosity ratio against the matrix fluid is significant. To this end, a comprehensive investigation on the dynamics of a droplet in shear flow, which involves a wide range of Reynolds number, capillary number, confinements of the wall, viscosity/density ratio between the two phases, is called for.

Since pioneering works by Taylor on droplet deformation in shear and extensional flows [8,9], enormous theoretical and experimental studies have been conducted. A series of works by the group of Mason [10–12] further studied the deformation and burst of droplets, and even depicted the streamlines inside and around the droplets. Chaffey and Brenner [13] extended a previous analytical approximation to a second order form, which is crucial for the non-elliptic deformation of a highly viscous droplet under large shear rate. Barthes-Biesel and Acrivos [14] expressed the solution of creeping-flow equations in powers of deformation parameters and applied a linear stability theory to determine the critical values for the droplet breakup. Hinch and Acrivos [15] investigated theoretically the stability of a long slender droplet, which is largely deformed in shear flow. However, early analytical works rarely considered effects of finite Reynolds number or wall confinements. In addition, numerous experimental studies have been conducted on the droplet deformation and breakup [16–19], where not only the effects of viscosity ratio between the droplet and the matrix fluid [20,21], but also wall confinements [22,23] have been taken into account.

With advance in computational science, numerical simulation has become a popular approach to study droplet dynamics in the past decades. Boundary integral method was among the first to be applied to study deformation of droplets in stationary and transient states [24], non-Newtonian droplets [25], and migration of a droplet in shear flow [26]. Moreover, Li et al. [27] employed a volume-of-fluid (VOF) method and Galerkin projection technique to simulate the process of droplet breakup. In the work of Amani et al. [28], a conservative level-set (CLS) method built on a conservative finite-volume approximation is applied to study the effect of viscosity ratio and wall confinement on the critical capillary number. In addition, lattice Boltzmann method (LBM) has been widely employed to study deformation, breakup and coalescence of droplets [29–33]; to model viscoelastic droplet [34] and surfactant-laden droplet [35]. We note that an interface tracing technique such as VOF, CLS, a phase-field formulation, or immersed boundary method is often necessary by a flow solver based on Eulerian meshes.

As a Lagrangian method, smoothed particle hydrodynamics (SPH) method has some advantages in simulating multiphase flows. Since different phases are identified by different types of particles, the interface automatically emerges without an auxillary tracing technique, even for a very large deformation. Moreover, inertia and wall effects can be taken into account straightforward, in contrast to theoretical analysis or the boundary integral method. Since its inception in astrophysics, SPH method has been largely developed and widely applied in various flow problems [36,37]. Morris [38] considered the surface tension based on a continuous surface force model and simulated an oscillating two-dimensional rod in SPH. Hu et al. [39] proposed a multi-phase model that handles both macroscopic and mesoscopic flows in SPH, where a droplet in shear flow was selected as a benchmark to validate the method. Other improvements and modifications have also been proposed for SPH in the context of multiphase problems [40–43]. Furthermore, a droplet or matrix flow with special properties can also be considered. For example, Moinfar et al. [44] studied the drop deformation under simple shear flow of Giesekus fluids and Vahabi [45] investigated the effect of thixotropy on deformation of a droplet under shear flow. Saghatchi et al. [46] studied the dynamics of a 2D double emulsion in shear flow with electric field based on an incompressible SPH method. There are also studies on colliding and coalescence process of droplets by SPH [47,48]. Simulation of bubbles in liquid is similar, but can encounter special challenges [49], due to the reverse density/viscosity ratio as that of droplet in gas.

Previously, simulations of multiphase flows by SPH method often investigated specific circumstances. Therefore, the objective of this paper is two fold: firstly, to simulate an extensive range of parameters to examine the SPH method for multiphase flows; secondly, to fill gaps of unexplored range of parameters and systematically investigate their influence on the droplet dynamics. The rest of the paper is arranged as follows: in Section 2, we introduce the multiphase SPH method and a specific surface tension model. We present validations and extensive numerical results in Section 3. We summarize this work after discussions in Section 4.

## 2 Methods

### 2.1 Governing equations and surface tension model

We consider isothermal Navier-Stokes equations with a surface tension for multiphase flow in Lagrangian frame

where *ρ*, **v** and *p* are density, velocity and pressure respectively. **F**_{b} is the body force, which is not considered in this study. **F**_{v}, **F**_{s} denote viscous force and surface tension at the interface between two phases, respectively.

Following previous studies of quasi-incompressible flow modeling [38], an artificial equation of state relating pressure to density can be written as

where *c*_{s} is an artificial sound speed and *ρ*_{ref} is a reference density. Theoretically, subtracting the reference density has no influence on the gradient of pressure, but it can reduce the numerical error of SPH discretizations for the gradient operator.

For a Newtonian flow, the viscous force **F**_{v} simplifies to

where *μ* is the dynamic viscosity. We assume surface tension to be uniform along the interface and do not consider Marangoni force. Therefore, the surface tension acts on the normal direction of the interface. Moreover, its magnitude depends on the local curvature as

where *σ*, *κ*, *δ*_{s} is a surface delta function and its discrete form shall be described later.

To describe the surface tension at the interface between two fluids, a continuous surface tension model is adopted. As a matter of fact, surface tension my be written as the divergence of a tensor **T** [50,51]

where

To represent a multiphase flow, we define a color function *c* and set a unique value for each phase, that is, *c*^{I} = 0 and *c*^{II} = 1 for the two phases, respectively. Apparently, the color function has a jump from 0 to 1 at the interface between phase I and II. Therefore, the unit normal vector can be represented by the normalized gradient of the color function as

and the surface delta function is replaced by the scaled gradient as

### 2.2 SPH method

In SPH, fluid is represented by moving particles carrying flow properties such as density, velocity and pressure. We largely follow the work of Hu and Adams [39] and provide a brief derivation here. Density of a particle is calculated by interpolating the mass of neighboring particles as

where mass *m*_{i} is constant for every particle. *W*_{ij} denotes a weight function for interpolation

where **r**_{ij} = **r**_{i} − **r**_{j} is a relative position vector from particle *j* to *i* and *h* is the smoothing length. We further define

to be an equivalent volume of particle *i* so that *V*_{i} = *m*_{i}/*ρ*_{i}.

The pressure gradient can be computed as

where *p*_{i} and *p*_{j} are obtained by Eq. 2. The viscous force can be calculated as

where **v**_{ij} = **v**_{i} − **v**_{j} is the relative velocity of particle *i* and *j* and

As suggested by Morris [38] and Hu et al. [39], a part of pressure contribution

to replace Eq. 6, where *d* is the spatial dimension. Combining Eqs 7, 8 and Eq. 14, we obtain

The gradient of color function between phase I and phase II can be calculated in SPH as

where *c*_{i} (or *c*_{j}) is initially assigned to be *c*^{I} or *c*^{II} according to which phase particle *i* (or *j*) consititutes. Substitute Eq. 16 into Eq. 15 to obtain stress tensor

Finally, the surface force term is calculated by the stress tensor using the SPH expression for divergence

It is simple to see that the discrete version of *δ*_{s} in SPH is

which has a finite support to remove the singularity and distributes the surface tension onto a thin layer of two fluids across the interface.

### 2.3 Computational settings

The quintic kernel is adopted as weight function

where *R* = *r*/*h* and *h* is the smoothing length. *ϕ* is a normalization coefficient which equals 1/120, 7/(478*π*) and 1/(120*π*) in one, two and three dimensions, respectively. We set *h* = 1.2Δ*x* with Δ*x* as the initial spacing distance between particles. This means that the support domain of the kernel function is truncated at 3.6Δ*x*, namely, the cutoff *r*_{c} = 3.6Δ*x*. According to our tests, a smoothing length of 1.2Δ*x* is almost optimal for an excellent accuracy while avoiding the pairing instability. A detailed discussion on this issue is referred to Price [52].

Since we adopt a weakly compressible formulation, the sound speed *c*_{s} should be large enough to restrict the density fluctuations. Based on a scale analysis, Morris et al. [38,53] suggested that

where Δ is the density variation and *U*, *L*, *F*, *κ* and *σ* are typical velocity, length, body force, curvature and surface tension coefficient, respectively. Accordingly, for multiphase flows the sound speed may be different for each phase. In all simulations, we set identical Δ ≤ 0.5% for each phase and calculate *c*_{s} accordingly.

At every time step, the minimal relative density is recorded among all particles, that is,

where particle *i* belongs to phase I and particle *j* belongs to phase II;

The explicit velocity-Verlet method is adopted for time integration and a time step is chosen appropriately for stability [38].

## 3 Numerical results

We consider a shear flow generated by two parallel walls with opposite velocity of magnitude *U*. Periodic boundaries apply in the *x* direction. The computational domain is with length *L* and height *H*. A circular droplet with radius *R*_{0} is initially located at the center of the computational domain, as shown in Figure 1A. The blue dashed lines represent periodic boundaries. The continuous phase has viscosity *μ*_{c} while the dispersed phase has viscosity *μ*_{d} = *λμ*_{c}.

**FIGURE 1**. Schematic representation of a droplet in shear flow. **(A)** Initial geometry and parameters. **(B)** Treatment of no-slip boundary condition.

No-slip boundary condition is applied at the wall-fluid interfaces using the method proposed by Morris [53]. As shown in Figure 1B, yellow particles represent fluid and blue particles represent the boundary wall. The wall particles and boundary move with velocity *v*_{b}, which depends on the shear rate. But when a wall particle *j* is in the support domain of a fluid particle *i*, a given velocity *v*_{j} = *d*_{j}/*d*_{i} (*v*_{b} − *v*_{i}) + *v*_{b} is used to calculate the viscosity force between particle *i* and particle *j* in Eq. 13.

Five dimensionless parameters that determine the deformation of the droplet are Reynolds number *R*_{0}/*H*, viscosity ratio *λ* = *μ*_{d}/*μ*_{c} and density ratio *α* = *ρ*_{d}/*ρ*_{c}, where *σ* is the surface tension coefficient, *ρ*_{d} and *μ*_{d} are density and viscosity of the dispersed fluid phase inside the droplet while and *ρ*_{c} and *μ*_{c} are for the continuous fluid phase, respectively.

In Section 3.1, we study the deformation for an intact droplet while considering the effects due to the five dimensionless numbers. In Section 3.2, we examine the breakup of the droplet. In Section 3.3, we summarize the droplet dynamics for both intact shape and breakup in phase diagrams. In Section 3.4, we demonstrate the deformation and breakup with physical parameters of a water droplet in air flow as an industrial application.

### 3.1 Droplet deformation

When the shear is mild, the droplet remains intact and deforms to arrive at a stable shape eventually. The degree of droplet deformation can be quantified by the Taylor deformation parameter *D* = (*A* − *B*)/(*A* + *B*), where *A* is the greatest length and *B* is the breadth of the droplet as shown in Figure 2. To validate our method, we first compare our results of transient deformations with that of Sheth and Pozrikidis using immersed boundary method within the finite difference method [54]. We follow their work to set *L* = *H* = 4*R*_{0} = 1, *ρ*_{d} = *ρ*_{c} = 1, *μ*_{d} = *μ*_{c} = 0.5 and adjust shear rate and surface tension. The two walls slide with velocities *x* = 2*R*_{0}/25 and *R*_{0}/25, corresponding to the droplet containing *N* = 484 and 1976 particles, respectively.

We present particle distributions and *D* as functions of time in Figure 3 for a typical simulation with *Re* = 0.125 and *Ca* = 0.45. We note that the deformation of the droplet may oscillate in time and its maximum elongation does not necessarily takes place at the steady state of a very long time.

**FIGURE 3**. Particle distribution for *H* =4*R*_{0}, *α* = *λ* =1, Re =0.125, *Ca* =0.45, *L* =4*R*_{0} at **(A)** initial configuration. **(B)** maximum elongation. **(C)** steady state. **(D)** Taylor deformation parameter as a function of time.

We further focus on the transient deformations in short time in Figure 4 so that we can compare our results with those of Sheth and Pozrikidis [54]. It can be readily seen that our results with low resolution Δ*x* = 0.02 or *N* = 484 already reproduce the reference very well for different Reynolds numbers and/or capillary numbers. As the reference is within a rather short time period, some interesting phenomenon such as oscillation of the Taylor deformation parameter *D* is not captured, as indicated for *Re* = 12.5 and *Ca* = 0.025 on Figure 4D.

**FIGURE 4**. Transient deformations of a droplet compared with results of Sheth and Pozrikidis [54] with *α* = *λ* =1, *H* =4*R*_{0}, *L* =4*R*_{0} and **(A)** Re =0.25, *Ca* =0.1,0.45. **(B)** Re =2.5, *Ca* =0.1,0.4. **(C)** Re =12.5, *Ca* =0.1,0.2. **(D)** Re =25, *Ca* =0.025,0.2.

To validate our method for vanishing Reynolds numbers, we calculate the stationary deformation and orientation of the droplet with respect to *Ca*. We follow Zhou and Pozrikidis [55] to set *L* = *H* = 8*R*_{0} = 2, *ρ*_{d} = *ρ*_{c} = 1, *μ*_{d} = *μ*_{c} = 0.5 and adjust shear rate and surface tension accordingly. The deformation parameter *D* and orientation *θ* (defined on Figure 2) as functions of *Ca* (up to *Ca* = 1) for *Re* = 0.01 are shown in Figure 5. Results for *Re* = 0.1 and 1 are also given for comparison, where droplet breakup already takes places at *Ca* ≳ 0.4 for *Re* = 1. The difference between the results of *Re* = 0.1 and *Re* = 0.01 is insignificant and they both resemble the results of boundary integral method for Stokes flow [55]. We can readily conclude that *Re* = 0.1 is small enough to approximate the Stokes flow. We further present the contours and streamlines for a typical evolution of droplet deformation at vanishing Reynolds number in Figure 6.

**FIGURE 5**. Effects of Re and *Ca* on 2D droplet deformation when *α* = *λ* =1, *H* =8*R*_{0}, *L* =8*R*_{0}. **(A)** Taylor deformation parameters. **(B)** Droplet orientation.

**FIGURE 6**. A typical evolution of deformation for an initially circular 2D droplet in shear flow: *α* = *λ* =1, Re =0.1, *Ca* =0.4, *H* =16*R*_{0}, *L* =16*R*_{0}. **(A)** Droplet deformation over time. **(B)** Streamlines: the color represents the magnitude of velocity and a red line indicates the droplet interface.

We commence to investigate the effects of confinement and set *L* = 16*R*_{0} to minimize the periodic artifacts. We first restrict out attention to *Re* = 0.1, *α* = 1 and *λ* = 1. Four ratios of confinement are considered: *H* = 2.4*R*_{0}, 4*R*_{0}, 8*R*_{0} and 16*R*_{0}. The deformation parameter as a function of *Ca* is shown in Figures 7A, B. As we can observe, a smaller distance of the two walls enhances the elongation of droplet and makes its long axis align more horizontally. As we relax the confinement, the relation between *D* and *Ca* becomes linear and the difference between *H* = 16*R*_{0} and *H* = 8*R*_{0} is already negligible. To further investigate the mechanism to the effects of wall confinement on the droplet deformation, we choose *Ca* = 0.4 plot the pressure around the droplet under four confinement ratio as shown in Figure 7C. We have subtracted a background pressure *p*_{0} = 1000*c*^{2}, where 1,000 is the initial density and *c* = 0.707 is used under *Ca* = 0.4. For *H* = 2.4*R*_{0} and *H* = 4*R*_{0} we plotted the entire flow field region. For *H* = 8*R*_{0} and *H* = 16*R*_{0}, we plot the middle region of 16*R*_{0} × 4*R*0. We can clearly see the presence of low-pressure regions at both ends of the droplet. Under different wall confinements, the pressure in these low-pressure regions does not differ significantly. On the sides of the droplet in the middle, there are high-pressure regions, and the pressure in these high-pressure regions increases with an increase in the wall confinements. These high-pressure regions on the sides cause the droplet to elongate and become more horizontal. We can observe that there is only a slight difference in the pressure field between *H* = 8*R*_{0} and *H* = 16*R*_{0}, so does the corresponding deformation of the droplet. A further exploration of why increasing the confinement causes an increase in pressure on sides of the droplet can be explained as follows: The proportion of the channel blocked by the droplet increases with the confinement. The outer flow is forced to bypass the droplet, resulting in a squeezing effect between the droplet and the wall. When the ratio *R*_{0}/*H* decreases, there is more space between the droplet and the wall to buffer this squeezing effect, dispersing the pressure over a larger area.

**FIGURE 7**. Effects of confinement ratio and *Ca* on 2D droplet deformation when *α* = *λ* =1, Re =0.1, *L* =16*R*_{0}. **(A)** Taylor deformation parameters. **(B)** Droplet orientation. **(C)** Pressure distribution around droplet under different confinements when *Ca* =0.4.

Furthermore, we simulate cases where the droplet and the matrix flow are two fluids with different physical properties. We first consider two fluids of the same density but with different viscosities. We choose a computational domain of 16*R*_{0} × 16*R*_{0} and set *Re* = 0.1, *α* = 1 and *λ* ranges from 0.1 to 10. Initial space Δ*x* among nearest particles is 2*R*_{0}/25 so a droplet contains 484 particles. The deformation parameter as a function of *Ca* is shown in Figure 8. As we can observe, the deformation increases as *λ* increases from 0.1 to 10. In this range of *λ*, a droplet with lower viscosity has a smooth inside circulation and fast reaction which can reduce the elongation [16,20].

**FIGURE 8**. Effects of viscosity ratio *λ* and *Ca* on 2D droplet deformation when *α* =1, Re =0.1, *H* =16*R*_{0}, *L* =16*R*_{0}. **(A)** Taylor deformation parameters. **(B)** Droplet orientation.

The other case is that fluids inside and outside the droplet have the same viscosity but different densities. The sound speed is chosen according to the ratio of initial density to balance the pressure

where *Ca* when *α* = 0.1. In this small Reynolds number regime (*Re* = 0.1), the density ratio has negligible influence and only the capillary number determines the droplet deformation.

**FIGURE 9**. Effects of density ratio *α* and *Ca* on 2D droplet deformation when *λ* =1, Re =0.1, *H* =16*R*_{0}, *L* =16*R*_{0}. **(A)** Taylor deformation parameters. **(B)** Droplet orientation.

In 3D simulations, the width of simulation box *W* is an additional computational parameter compared to the 2D simulations. To compare with analytical predictions or experiment data, the length and width of simulation box are numerical and should be large enough. One set of parameters of *Re* = 0.1, *Ca* = 0.2, *α* = *λ* = 1, *H* = 4*R*_{0} are selected and different length *L* and width *W* of simulation box are examined. According to our simulations, the deformation basically decreases with the increase of *L* and/or width *W*. We compare the Taylor deformation parameter *D* in steady states of our simulations with the analytical prediction of Shapira and Haber [56]. The differences between our results and analytical prediction under different *L* and *W* are plotted in Figure 10. It can be seen that when *L* is larger than 24*R*_{0} and *W* is larger than 8*R*_{0}, the results has little change with the increase of *L* and/or *W*. Figure 11 shows the steady deformation of 3D droplets in shear flow when *L* = 24*R*_{0}, *W* = 8*R*_{0}, *Re* = 0.1 and *α* = *λ* = 1 with different *Ca* and confinement in *H* direction, compared with theoretical predictions of Shapira and Haber [56] and experiment data of Sibillo et al. [57]. Our results agree well with both anlaytical and experiment references at *Ca* = 0.1 and 0.2, whereas are closer to the experimental data at *Ca* = 0.3. The deformation increases with the confinement ratio *R*_{0}/*H*, which has the same trend as for 2D cases.

**FIGURE 10**. Effects of box length *L* and width *W* on the Taylor deformation parameter *D* of 3D droplet deformation when *α* = *λ* =1, Re =0.1, *Ca* =0.2, *H* =4*R*_{0}, compared to the analytical prediction of [56]. **(A)** Each line represents the same *W/R*_{0}. **(B)** Each line represents the same *L/R*_{0}.

**FIGURE 11**. Effects of confinement ratio and *Ca* on 3D droplet deformation in shear flow when *α* =1 and Re =0.1, *L* =24*R*_{0}, *W* =8*R*_{0}, compared with predictions of Shapira and Haber [56] and experiment data of Sibillo et al. [57].

### 3.2 Droplet breakup

When the shear is strong, the droplet is over-stretched to break up. We find two patterns of breakup process under different viscosity ratios in simulations. As shown in Figure 12, when *α* = 1, *λ* = 0.2, *Re* = 0.1, *Ca* = 10, and *L* = *H* = 16*R*_{0}, a droplet is rotated and then stripped of its main body near the surface and gradually breaks apart. We call this breakup type A. This type is also found in the experiment study of Grace and they call it “tip streaming breakup” [20]. The conditions for type A breakup happening is exhibited in the next section. Figure 13 shows another set of typical snapshots of the droplet shape and flow fields in shear flow breaks when *α* = *λ* = 1, *Re* = 0.1, *Ca* = 0.9 and *L* = *H* = 16*R*_{0}. In this simulation, a droplet is stretched and its waist becomes slender and slender and finally breaks up. We call this breakup type B.

**FIGURE 12**. Breakup type **(A)** evolution of an initially circular 2D droplet breakup in shear flow with *α* =1, *λ* =0.2, Re =0.1, *Ca* =10, *H* =16*R*_{0}, *L* =16*R*_{0}. **(A)** Droplet shape. **(B)** Streamlines: The color represents the magnitude of velocity and red outlines in the background represent the droplet interface.

**FIGURE 13**. Breakup type **(B)** evolution of an initially circular 2D droplet breakup in shear flow with *α* = *λ* =1, Re =0.1, *Ca* =0.9, *H* =16*R*_{0}, *L* =16*R*_{0}. **(A)** Droplet shape. **(B)** Streamlines: The color represents the magnitude of velocity and red outlines in the background represent the droplet interface.

To encompass the breakup of a 3D droplet with a large elongation, we employ a rather long computational domain with *L* = 32*R*_{0}. Figure 14 shows the dynamics of the breakup with Re = 0.1, *H* = 2.857*R*_{0}, *Ca* = 0.46 and *α* = *λ* = 1. Left side are SPH particle distributions and right side are corresponding contour interfaces processed by SPH kernel interpolation into mesh cells. The color represents the magnitude of velocity. We adopt the same *Ca* and *R*_{0}/*H* as the experiment in creeping flow by Sibillo et al. [57]. The shape of the droplet in the breakup process of our simulation is very close to their experimental observation. Only a slight difference appears in the final stage: in the experiment, the droplet is divided into three main parts, while in our simulation the middle part continues to split into two smaller droplets. In contrast to the 2D case, a 3D droplet has a more slender shape before breaking up.

**FIGURE 14**. Evolution of an initially spherical 3D droplet breakup in shear flow when *α* = *λ* =1, Re =0.1, *Ca* =0.46, *H* =2.857*R*_{0}, *L* =16*R*_{0}: particle distribution (left) and interface (right).

### 3.3 Phase diagram

To clearly visualize the states of a droplet in different conditions, we consider a range of Reynolds numbers, capillary numbers, and confinements/density/viscosity ratios and summarize our simulation results into phase diagrams. Thereafter, we may estimate the critical capillary number *Ca*_{c} that segments the intact and breakup states and further investigate how it is influenced by other dimensionless parameters.

For *λ* = *α* = 1, we perform a group of 2D simulations with different Reynolds number Re = 0.01, 0.1, 1, 10 and confinement *H* = 2.4*R*_{0}, 4*R*_{0}, 8*R*_{0}, 16*R*_{0} with *Ca* ∈ [0.1, 1.1] and *L* = 16*R*_{0}. For a general overview, the states of the droplet are summarized in Figure 15. To get a clear view, we slice the phase diagram by two perspectives. Firstly, we divide results into groups of the same confinement to reveal the influence of Re on *Ca*_{c} as shown in Figure 16. Overall it is apparent that a higher Re reduces *Ca*_{c}. Three scenarios are special: under confinement *H* = 2.4*R*_{0}, 4*R*_{0} and 8*R*_{0}, we can not differentiate *Ca*_{c} between Re = 0.01 and 0.1.

**FIGURE 15**. Phase diagram of 2D droplets states under different confinement, Re and *Ca* when *α* = *λ* =1, *L* =16*R*_{0}.

**FIGURE 16**. Phase diagram of 2D droplets states under different confinement, Re and *Ca* when *α* = *λ* =1, *L* =16*R*_{0} and **(A)** *H* =2.4*R*_{0}. **(B)** *H* =4*R*_{0}. **(C)** *H* =8*R*_{0}. **(D)** *H* =16*R*_{0}.

From another perspective of *Ca versus* confinement ratio for each Re on Figure 17, we are not able to find a universal pattern. Under Re = 0.01, *Ca*_{c} decreases with *R*_{0}/*H* while under Re = 10, *Ca*_{c} increases with *R*_{0}/*H*. Whereas, under Re = 0.1 and 1, *Ca*_{c} has no monotonic relation with *R*_{0}/*H*.

**FIGURE 17**. Phase diagram of 2D droplets states under different confinement, Re and *Ca* when *α* = *λ* =1, *L* =16*R*_{0} and **(A)** Re =0.01. **(B)** Re =0.1. **(C)** Re =1. **(D)** Re =10.

Furthermore, we investigate effects of viscosity ratio *λ* = *μ*_{d}/*μ*_{c} ∈ [0.1, 10] on the droplet dynamics for Re = 0.1 and three confinement ratios *H* = 4*R*_{0}, 8*R*_{0} and 16*R*_{0}. The results are shown in Figure 18. For breakup type A, the droplet rotates and is stripped off as described in Section 3.2; Breakup type B represents that a droplet is stretched and breaks up in the middle. Under Re = 0.1, type A is observed only if the droplet has a much smaller viscosity compare to the matrix fluid. Overall, *Ca*_{c} decreases with the increase of *λ*. However, we notice a flatten trend or even a reverse trend with small difference for *Ca*_{c} from *λ* = 5 to *λ* = 10, as shown on the insets of Figure 18. According to the study of Karam et al. and Grace [16,20], a maximum transfer of energy takes place across an interface, which demands this trend.

**FIGURE 18**. Phase diagram of 2D droplets states under different confinement, viscosity ratios *λ* and *Ca* when *α* =1, Re =0.1, *L* =16*R*_{0}. **(A)** *H* =4*R*_{0}. **(B)** *H* =8*R*_{0}. **(C)** *H* =16*R*_{0}. **(D)** two patterns of breakup.

Due to highly computational cost in 3D, we only consider a moderate confinement *H*/*R*_{0} = 4 and perform a group of simulations to draw a phase diagram in the plane of *Ca* and Re, as shown in Figure 19. The size of the simulation box is *L* = 32*R*_{0}, *W* = 8*R*_{0}, *H* = 4*R*_{0}. As in 2D case, the critical *Ca*_{c} decreases with increasing Re in 3D, as shown in in Figure 16. However, the critical capillary number *Ca*_{c} in 3D case is significantly smaller than that of 2D case.

**FIGURE 19**. Phase diagram of 3D droplets states under different Re and *Ca* when *α* = *λ* =1, *H* =4*R*_{0}, *L* =32*R*_{0}, *W* =8*R*_{0}.

### 3.4 Water droplet in air flow

As one specific application, we employ our method to predict the breakup of a water droplet in shear flow of air. The critical capillary number or the shear rate determined is helpful to design an effective atomization device. Actual physical properties of water and air around 20°C are adopted: *ρ*_{d} = 998.2 *kg* ⋅ *m*^{−3}, *μ*_{d} = 1.0087 × 10^{−3} *Pa* ⋅ *s* and *ρ*_{c} = 1.205 *kg* ⋅ *m*^{−3}, *μ*_{c} = 1.81 × 10^{−5} *Pa* ⋅ *s* are set for water (dispersed) phase and air (continuous) phase, respectively; surface tension coefficient *σ* = 72.75 × 10^{−3} *N* ⋅ *m*^{−1} is set for the water-air interface.

We perform a relative large range of Reynolds numbers and depict a phase diagram on the plane of Re and *Ca* in logarithmic-logarithmic scales on Figure 20. This allows us to connect the results with the same droplet size and observe its behavior while changing Re and *Ca*. Points on each dotted line represent the droplet of the same radius, as marked in the figure. For example, we have a line of dynamics for the droplet with *R*_{0} = 10*μm* under shear rates of 1 × 10^{6}*s*^{−1}, 2 × 10^{6}*s*^{−1}, 5 × 10^{6}*s*^{−1}, 1 × 10^{7}*s*^{−1}, 2 × 10^{7}*s*^{−1}; another line of dynamics for the droplet with *R*_{0} = 100*μm* under shear rates of 5 × 10^{4}*s*^{−1}, 1 × 10^{5}*s*^{−1}, 2 × 10^{5}*s*^{−1}, 5 × 10^{5}*s*^{−1}, 1 × 10^{6}*s*^{−1}. Furthermore, we observe that if the Re is on the order of 100, the critical *Ca* for breakup is very sensitive to Re. We also perform a group of 3D simulations for a droplet with *R*_{0} = 50*μm* under shear rates of 1 × 10^{5}*s*^{−1}, 2 × 10^{5}*s*^{−1}, 5 × 10^{5}*s*^{−1}, 1 × 10^{6}*s*^{−1}, 2 × 10^{6}*s*^{−1}. The 3D results for the critical point of breakup is close to that of the 2D results.

**FIGURE 20**. Diagram for states of water droplets in air shear flow under different Re and *Ca* when *H* =16*R*_{0}, *L* =16*R*_{0}.

## 4 Conclusion and disucssion

In this study, we employed a multi-phase SPH method to simulate droplet deformation and breakup subjected to a simple shear flow in an extensive range of physical parameters. We performed both 2D and 3D simulations and validated them by benchmarks: transient deformations and steady shapes of droplets are compared with previous simulations, analytical derivations and experimental data. These results indicate that the method is reliable to simulate droplet dynamics in general. We wish to emphasize the convenience of SPH method in simulating multi-phase problems, as we can leverage on its Lagrangian nature and differentiate different phases by particle species. In addition, the algorithm and data structure for 2D and 3D simulations have tiny difference and therefore, it is a simple task to extend the code from 2D to 3D. Economical 2D simulations allow us to investigate a wide range of physical parameters in five dimensions, which serve as a guide to 3D realistic situations. From the results, we come to the following conclusions.

(1) A larger Reynolds number Re or capillary number *Ca* leads to a more considerable deformation of the droplet. The transient and steady-state deformations of the droplet in our study are in good agreement with the previous studies but beyond their time limits [54,55].

(2) Under low Reynolds number (Re = 0.1), a stronger confinement due to the walls enhances the steady-state deformation in both 2D and 3D simulations. When the walls are separated further apart, the Taylor deformation parameter is almost linear with respect to *Ca*. The influence of confinement on the deformation of a droplet has been studied by Shapira and Haber by a first-order analytical solution based on Lorentz’s reflection method. They proved that the walls do not influence the shape of deformed droplet but increases the deformation magnitude with a term of order

(3) The effects of wall confinement on the critical capillary number *Ca*_{c} are not universal under different Re. When Re = 0.1, a closer gap of walls reduces *Ca*_{c}. This is because a closer gap of walls increases the deformation as described above. But when Re is larger, the relation between *Ca*_{c} and the confinement ratio is unclear. From our observation, this non-monotonic relation results from an interplay of influences by the shear strength and the stability of the whole flow field. On the one hand, the shear stress transferred to the droplet from the wall is more pronounced in stronger confinement [56], thus closer walls reduce the *Ca*_{c}. On the other hand, the narrower channel reduces the instability of the flow and restricts droplet movements, thus increases the *Ca*_{c}.

(4) Under Re = 0.1 and the range of viscosity ratio *λ* causes a larger deformation. The effect of *λ* on *Ca*_{c} is not monotonic when *λ* > 1 and there is a minimum value of *Ca*_{c} between *λ* = 1 and *λ* = 10. The existence of a minimal *Ca*_{c} among different *λ* has also been found by previous experiment studies [16,20], when *λ* is about 1. The discrepancy between our results and the previous ones are attributed to the difference between 2D and 3D cases. At the same Re, the influence of density ratio on droplet deformation is much smaller compared with that of the viscosity ratio.

(5) As an application, a phase diagram obtained by actual physical parameters of water and air is depicted to predict the magnitude of shear rate for breaking a droplet of certain size, which is helpful in the designing atomization nozzles.

## Data availability statement

The original contributions presented in the study are included in the article, further inquiries can be directed to the corresponding author.

## Author contributions

KW: Formal Analysis, Investigation, Methodology, Validation, Writing–original draft, Writing–review and editing. HL: Formal Analysis, Supervision, Validation, Writing–review and editing. CZ: Conceptualization, Funding acquisition, Writing–review and editing. XB: Conceptualization, Formal Analysis, Funding acquisition, Supervision, Writing–review and editing.

## Funding

The authors declare financial support was received for the research, authorship, and/or publication of this article. The authors declare that this study received funding from the national natural science foundation of China under grant number: 12172330. This study also received funding from Hangzhou Shiguangji Intelligent Electronics Technology Co., Ltd, Hangzhou, China. The funder was not involved in the study design, collection, analysis, interpretation of data, the writing of this article, or the decision to submit it for publication.

## Conflict of interest

Author CZ was employed by Hangzhou Shiguangji Intelligient Electronics Technology Co., Ltd.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

## References

1. Anna SL. Droplets and bubbles in microfluidic devices. *Annu Rev Fluid Mech* (2016) 48:285–309. doi:10.1146/annurev-fluid-122414-034425

2. Dressler OJ, Casadevall i Solvas X, DeMello AJ. Chemical and biological dynamics using droplet-based microfluidics. *Annu Rev Anal Chem* (2017) 10:1–24. doi:10.1146/annurev-anchem-061516-045219

3. Liu J, Hao M, Chen S, Yang Y, Li J, Mei Q, et al. Numerical evaluation of face masks for prevention of covid-19 airborne transmission. *Environ Sci Pollut Res* (2022) 29:44939–53. doi:10.1007/s11356-022-18587-3

4. Lohse D. Fundamental fluid dynamics challenges in inkjet printing. *Annu Rev Fluid Mech* (2022) 54:349–82. doi:10.1146/annurev-fluid-022321-114001

5. Aydin O, Unal R. Experimental and numerical modeling of the gas atomization nozzle for gas flow behavior. *Comput Fluids* (2011) 42:37–43. doi:10.1016/j.compfluid.2010.10.013

6. Si C, Zhang X, Wang J, Li Y. Design and evaluation of a laval-type supersonic atomizer for low-pressure gas atomization of molten metals. *Int J Minerals, Metall Mater* (2014) 21:627–35. doi:10.1007/s12613-014-0951-4

7. Xu Z, Wang T, Che Z. Droplet deformation and breakup in shear flow of air. *Phys Fluids* (2020) 32:052109. doi:10.1063/5.0006236

8. Taylor GI. The viscosity of a fluid containing small drops of another fluid. *Proc R Soc Lond Ser A, Containing Pap a Math Phys Character* (1932) 138:41–8.

9. Taylor GI. The formation of emulsions in definable fields of flow. *Proc R Soc Lond Ser A, containing Pap a Math Phys character* (1934) 146:501–23.

10. Bartok W, Mason S. Particle motions in sheared suspensions: viii. singlets and doublets of fluid spheres. *J Colloid Sci* (1959) 14:13–26. doi:10.1016/0095-8522(59)90065-0

11. Rumscheidt F-D, Mason S. Particle motions in sheared suspensions xii. deformation and burst of fluid drops in shear and hyperbolic flow. *J Colloid Sci* (1961) 16:238–61. doi:10.1016/0095-8522(61)90003-4

12. Torza S, Henry C, Cox R, Mason S. Particle motions in sheared suspensions. xxvi. streamlines in and around liquid drops. *J Colloid Interf Sci* (1971) 35:529–43. doi:10.1016/0021-9797(71)90211-6

13. Chaffey CE, Brenner H. A second-order theory for shear deformation of drops. *J Colloid Interf Sci* (1967) 24:258–69. doi:10.1016/0021-9797(67)90229-9

14. Barthes-Biesel D, Acrivos A. Deformation and burst of a liquid droplet freely suspended in a linear shear field. *J Fluid Mech* (1973) 61:1–22. doi:10.1017/s0022112073000534

15. Hinch E, Acrivos A. Long slender drops in a simple shear flow. *J Fluid Mech* (1980) 98:305–28. doi:10.1017/s0022112080000171

16. Karam H, Bellinger J. Deformation and breakup of liquid droplets in a simple shear field. *Ind Eng Chem Fundamentals* (1968) 7:576–81. doi:10.1021/i160028a009

17. Flumerfelt RW. Drop breakup in simple shear fields of viscoleastic fluids. *Ind Eng Chem Fundamentals* (1972) 11:312–8. doi:10.1021/i160043a005

18. Stone HA, Bentley B, Leal L. An experimental study of transient effects in the breakup of viscous drops. *J Fluid Mech* (1986) 173:131–58. doi:10.1017/s0022112086001118

19. Guido S, Villone M. Three-dimensional shape of a drop under simple shear flow. *J rheology* (1998) 42:395–415. doi:10.1122/1.550942

20. Grace HP. Dispersion phenomena in high viscosity immiscible fluid systems and application of static mixers as dispersion devices in such systems. *Chem Eng Commun* (1982) 14:225–77. doi:10.1080/00986448208911047

21. Stone HA, Leal LG. The influence of initial deformation on drop breakup in subcritical time-dependent flows at low Reynolds numbers. *J Fluid Mech* (1989) 206:223–63. doi:10.1017/s0022112089002296

22. Vananroye A, Van Puyvelde P, Moldenaers P. Effect of confinement on droplet breakup in sheared emulsions. *Langmuir* (2006) 22:3972–4. doi:10.1021/la060442+

23. Vananroye A, Van Puyvelde P, Moldenaers P. Effect of confinement on the steady-state behavior of single droplets during shear flow. *J rheology* (2007) 51:139–53. doi:10.1122/1.2399089

24. Kennedy MR, Pozrikidis C, Skalak R. Motion and deformation of liquid drops, and the rheology of dilute emulsions in simple shear flow. *Comput Fluids* (1994) 23:251–78. doi:10.1016/0045-7930(94)90040-x

25. Toose E, Geurts B, Kuerten J. A boundary integral method for two-dimensional (non)-Newtonian drops in slow viscous flow. *J non-newtonian Fluid Mech* (1995) 60:129–54. doi:10.1016/0377-0257(95)01386-3

26. Uijttewaal W, Nijhof E. The motion of a droplet subjected to linear shear flow including the presence of a plane wall. *J Fluid Mech* (1995) 302:45–63. doi:10.1017/s0022112095004009

27. Li J, Renardy YY, Renardy M. Numerical simulation of breakup of a viscous drop in simple shear flow through a volume-of-fluid method. *Phys Fluids* (2000) 12:269–82. doi:10.1063/1.870305

28. Amani A, Balcázar N, Castro J, Oliva A. Numerical study of droplet deformation in shear flow using a conservative level-set method. *Chem Eng Sci* (2019) 207:153–71. doi:10.1016/j.ces.2019.06.014

29. Xi H, Duncan C. Lattice Boltzmann simulations of three-dimensional single droplet deformation and breakup under simple shear flow. *Phys Rev E* (1999) 59:3022–6. doi:10.1103/physreve.59.3022

30. van der Sman RG, Van der Graaf S. Emulsion droplet deformation and breakup with lattice Boltzmann model. *Comp Phys Commun* (2008) 178:492–504. doi:10.1016/j.cpc.2007.11.009

31. Farokhirad S, Lee T, Morris JF. Effects of inertia and viscosity on single droplet deformation in confined shear flow. *Commun Comput Phys* (2013) 13:706–24. doi:10.4208/cicp.431011.260112s

32. Komrakova A, Shardt O, Eskin D, Derksen J. Lattice Boltzmann simulations of drop deformation and breakup in shear flow. *Int J Multiphase Flow* (2014) 59:24–43. doi:10.1016/j.ijmultiphaseflow.2013.10.009

33. Huang B, Liang H, Xu J. Lattice Boltzmann simulation of binary three-dimensional droplet coalescence in a confined shear flow. *Phys Fluids* (2022) 34:032101. doi:10.1063/5.0082263

34. Wang D, Tan DS, Khoo BC, Ouyang Z, Phan-Thien N. A lattice Boltzmann modeling of viscoelastic drops’ deformation and breakup in simple shear flows. *Phys Fluids* (2020) 32:123101. doi:10.1063/5.0031352

35. Zong Y, Zhang C, Liang H, Wang L, Xu J. Modeling surfactant-laden droplet dynamics by lattice Boltzmann method. *Phys Fluids* (2020) 32:122105. doi:10.1063/5.0028554

36. Monaghan J. Smoothed particle hydrodynamics and its diverse applications. *Annu Rev Fluid Mech* (2012) 44:323–46. doi:10.1146/annurev-fluid-120710-101220

37. Ye T, Pan D, Huang C, Liu M. Smoothed particle hydrodynamics (sph) for complex fluid flows: recent developments in methodology and applications. *Phys Fluids* (2019) 31:011301. doi:10.1063/1.5068697

38. Morris JP. Simulating surface tension with smoothed particle hydrodynamics. *Int J Numer Methods Fluids* (2000) 33:333–53. doi:10.1002/1097-0363(20000615)33:3<333::aid-fld11>3.0.co;2-7

39. Hu XY, Adams NA. A multi-phase sph method for macroscopic and mesoscopic flows. *J Comput Phys* (2006) 213:844–61. doi:10.1016/j.jcp.2005.09.001

40. Wang Z-B, Chen R, Wang H, Liao Q, Zhu X, Li S-Z. An overview of smoothed particle hydrodynamics for simulating multiphase flow. *Appl Math Model* (2016) 40:9625–55. doi:10.1016/j.apm.2016.06.030

41. Zhang M. Simulation of surface tension in 2d and 3d with smoothed particle hydrodynamics method. *J Comput Phys* (2010) 229:7238–59. doi:10.1016/j.jcp.2010.06.010

42. Tartakovsky AM, Panchenko A. Pairwise force smoothed particle hydrodynamics model for multiphase flow: surface tension and contact line dynamics. *J Comput Phys* (2016) 305:1119–46. doi:10.1016/j.jcp.2015.08.037

43. Yang Q, Yao J, Huang Z, Asif M. A comprehensive sph model for three-dimensional multiphase interface simulation. *Comput Fluids* (2019) 187:98–106. doi:10.1016/j.compfluid.2019.04.001

44. Moinfar Z, Vahabi S, Vahabi M. Numerical simulation of drop deformation under simple shear flow of giesekus fluids by sph. *Int J Numer Methods Heat Fluid Flow* (2022) 33:263–81. doi:10.1108/hff-01-2022-0067

45. Vahabi M. The effect of thixotropy on deformation of a single droplet under simple shear flow. *Comput Math Appl* (2022) 117:206–15. doi:10.1016/j.camwa.2022.04.023

46. Saghatchi R, Ozbulut M, Yildiz M. Dynamics of double emulsion interfaces under the combined effects of electric field and shear flow. *Comput Mech* (2021) 68:775–93. doi:10.1007/s00466-021-02045-x

47. Hirschler M, Oger G, Nieken U, Le Touzé D. Modeling of droplet collisions using smoothed particle hydrodynamics. *Int J Multiphase Flow* (2017) 95:175–87. doi:10.1016/j.ijmultiphaseflow.2017.06.002

48. Xu X, Tang T, Yu P. A modified sph method to model the coalescence of colliding non-Newtonian liquid droplets. *Int J Numer Methods Fluids* (2020) 92:372–90. doi:10.1002/fld.4787

49. Zhang A, Sun P, Ming F. An sph modeling of bubble rising and coalescing in three dimensions. *Comp Methods Appl Mech Eng* (2015) 294:189–209. doi:10.1016/j.cma.2015.05.014

50. Brackbill JU, Kothe DB, Zemach C. A continuum method for modeling surface tension. *J Comput Phys* (1992) 100:335–54. doi:10.1016/0021-9991(92)90240-y

51. Lafaurie B, Nardone C, Scardovelli R, Zaleski S, Zanetti G. Modelling merging and fragmentation in multiphase flows with surfer. *J Comput Phys* (1994) 113:134–47. doi:10.1006/jcph.1994.1123

52. Price DJ. Smoothed particle hydrodynamics and magnetohydrodynamics. *J Comput Phys* (2012) 231:759–94. doi:10.1016/j.jcp.2010.12.011

53. Morris JP, Fox PJ, Zhu Y. Modeling low Reynolds number incompressible flows using sph. *J Comput Phys* (1997) 136:214–26. doi:10.1006/jcph.1997.5776

54. Sheth KS, Pozrikidis C. Effects of inertia on the deformation of liquid drops in simple shear flow. *Comput Fluids* (1995) 24:101–19. doi:10.1016/0045-7930(94)00025-t

55. Zhou H, Pozrikidis C. The flow of suspensions in channels: single files of drops. *Phys Fluids A* (1993) 5:311–24. doi:10.1063/1.858893

56. Shapira M, Haber S. Low Reynolds number motion of a droplet in shear flow including wall effects. *Int J multiphase flow* (1990) 16:305–21. doi:10.1016/0301-9322(90)90061-m

Keywords: droplet, multiphase flow, surface tension, shear, SPH

Citation: Wang K, Liang H, Zhao C and Bian X (2023) Dynamics of a droplet in shear flow by smoothed particle hydrodynamics. *Front. Phys.* 11:1286217. doi: 10.3389/fphy.2023.1286217

Received: 31 August 2023; Accepted: 29 September 2023;

Published: 12 October 2023.

Edited by:

He Li, University of Georgia, United StatesReviewed by:

Kaixuan Zhang, Nankai University, ChinaJiayi Zhao, University of Shanghai for Science and Technology, China

Copyright © 2023 Wang, Liang, Zhao and Bian. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Xin Bian, bianx@zju.edu.cn