Abstract
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 framewhere ρ, v and p are density, velocity and pressure respectively. Fb is the body force, which is not considered in this study. Fv, Fs 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 aswhere cs 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 Fv simplifies towhere μ 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 aswhere σ, κ, are surface tension coefficient, curvature and unit normal vector to the concave side, respectively; δ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, cI = 0 and cII = 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 asand 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 aswhere mass mi is constant for every particle. Wij denotes a weight function for interpolationwhere rij = ri − rj is a relative position vector from particle j to i and h is the smoothing length. We further defineto be an equivalent volume of particle i so that Vi = mi/ρi.
The pressure gradient can be computed aswhere pi and pj are obtained by Eq. 2. The viscous force can be calculated aswhere vij = vi − vj is the relative velocity of particle i and j and is the distance between them.
As suggested by Morris [38] and Hu et al. [39], a part of pressure contribution is removed to avoid attractive force and improve the stability of the interactions between SPH particles. Therefore, we employto replace Eq. 6, where d is the spatial dimension. Combining Eqs 7, 8 and Eq. 14, we obtainThe gradient of color function between phase I and phase II can be calculated in SPH aswhere ci (or cj) is initially assigned to be cI or cII 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 iswhich 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 functionwhere 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 rc = 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 cs should be large enough to restrict the density fluctuations. Based on a scale analysis, Morris et al. [38,53] suggested that should be comparable to the largest ofwhere Δ 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 cs 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; , are initial densities for the two phases, respectively. Thereafter, , are subtracted as reference density for each phase in Eq. 2 to compute the particle pressure. This operation is performed to reduce numerical errors in calculating the pressure gradient while still keeping repulsive forces between particles.
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 R0 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
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 vb, 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 vj = dj/di (vb − vi) + vb 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 , Capillary number , confinement ratio R0/H, viscosity ratio λ = μd/μc and density ratio α = ρd/ρc, where is the shear rate, σ 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 = 4R0 = 1, ρd = ρc = 1, μd = μc = 0.5 and adjust shear rate and surface tension. The two walls slide with velocities to generate a clockwise rotation of the droplet. Two resolutions are considered for particles initially placed on squared lattice: Δx = 2R0/25 and R0/25, corresponding to the droplet containing N = 484 and 1976 particles, respectively.
FIGURE 2
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
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
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 = 8R0 = 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
FIGURE 6
We commence to investigate the effects of confinement and set L = 16R0 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.4R0, 4R0, 8R0 and 16R0. 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 = 16R0 and H = 8R0 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 p0 = 1000c2, where 1,000 is the initial density and c = 0.707 is used under Ca = 0.4. For H = 2.4R0 and H = 4R0 we plotted the entire flow field region. For H = 8R0 and H = 16R0, we plot the middle region of 16R0 × 4R0. 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 = 8R0 and H = 16R0, 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 R0/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
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 16R0 × 16R0 and set Re = 0.1, α = 1 and λ ranges from 0.1 to 10. Initial space Δx among nearest particles is 2R0/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
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 pressurewhere , and , are sound speeds and reference densities used for fluids outside and inside the droplet. As shown in Figure 9, the difference between deformations of droplet under density ratio 0.1–10 is very small except obvious lower inclination at small 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
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 = 4R0 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 24R0 and W is larger than 8R0, 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 = 24R0, W = 8R0, 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 R0/H, which has the same trend as for 2D cases.
FIGURE 10
FIGURE 11
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 = 16R0, 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 = 16R0. 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
FIGURE 13
To encompass the breakup of a 3D droplet with a large elongation, we employ a rather long computational domain with L = 32R0. Figure 14 shows the dynamics of the breakup with Re = 0.1, H = 2.857R0, 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 R0/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
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 Cac 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.4R0, 4R0, 8R0, 16R0 with Ca ∈ [0.1, 1.1] and L = 16R0. 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 Cac as shown in Figure 16. Overall it is apparent that a higher Re reduces Cac. Three scenarios are special: under confinement H = 2.4R0, 4R0 and 8R0, we can not differentiate Cac between Re = 0.01 and 0.1.
FIGURE 15
FIGURE 16
From another perspective of Caversus confinement ratio for each Re on Figure 17, we are not able to find a universal pattern. Under Re = 0.01, Cac decreases with R0/H while under Re = 10, Cac increases with R0/H. Whereas, under Re = 0.1 and 1, Cac has no monotonic relation with R0/H.
FIGURE 17
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 = 4R0, 8R0 and 16R0. 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, Cac decreases with the increase of λ. However, we notice a flatten trend or even a reverse trend with small difference for Cac 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
Due to highly computational cost in 3D, we only consider a moderate confinement H/R0 = 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 = 32R0, W = 8R0, H = 4R0. As in 2D case, the critical Cac decreases with increasing Re in 3D, as shown in in Figure 16. However, the critical capillary number Cac in 3D case is significantly smaller than that of 2D case.
FIGURE 19
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−3Pa ⋅ s and ρc = 1.205 kg ⋅ m−3, μc = 1.81 × 10−5Pa ⋅ s are set for water (dispersed) phase and air (continuous) phase, respectively; surface tension coefficient σ = 72.75 × 10−3N ⋅ 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 R0 = 10μm under shear rates of 1 × 106s−1, 2 × 106s−1, 5 × 106s−1, 1 × 107s−1, 2 × 107s−1; another line of dynamics for the droplet with R0 = 100μm under shear rates of 5 × 104s−1, 1 × 105s−1, 2 × 105s−1, 5 × 105s−1, 1 × 106s−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 R0 = 50μm under shear rates of 1 × 105s−1, 2 × 105s−1, 5 × 105s−1, 1 × 106s−1, 2 × 106s−1. The 3D results for the critical point of breakup is close to that of the 2D results.
FIGURE 20
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 [56]. The experiment data of Sibillo et al. illustrate satisfactory agreement with the predictions of Shapira and Haber except for the droplet being within a small gap, where the reflection analysis is expected to fail [57]. Our 3D simulation results resemble the whole set of experiment data even when the droplet is within the small gap, which suggests the method as an applicative tool for more realistic situations in microfluidics.
(3) The effects of wall confinement on the critical capillary number Cac are not universal under different Re. When Re = 0.1, a closer gap of walls reduces Cac. This is because a closer gap of walls increases the deformation as described above. But when Re is larger, the relation between Cac 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 Cac. On the other hand, the narrower channel reduces the instability of the flow and restricts droplet movements, thus increases the Cac.
(4) Under Re = 0.1 and the range of viscosity ratio , a higher λ causes a larger deformation. The effect of λ on Cac is not monotonic when λ > 1 and there is a minimum value of Cac between λ = 1 and λ = 10. The existence of a minimal Cac 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.
Statements
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.
AnnaSL. Droplets and bubbles in microfluidic devices. Annu Rev Fluid Mech (2016) 48:285–309. 10.1146/annurev-fluid-122414-034425
2.
DresslerOJCasadevall i SolvasXDeMelloAJ. Chemical and biological dynamics using droplet-based microfluidics. Annu Rev Anal Chem (2017) 10:1–24. 10.1146/annurev-anchem-061516-045219
3.
LiuJHaoMChenSYangYLiJMeiQet alNumerical evaluation of face masks for prevention of covid-19 airborne transmission. Environ Sci Pollut Res (2022) 29:44939–53. 10.1007/s11356-022-18587-3
4.
LohseD. Fundamental fluid dynamics challenges in inkjet printing. Annu Rev Fluid Mech (2022) 54:349–82. 10.1146/annurev-fluid-022321-114001
5.
AydinOUnalR. Experimental and numerical modeling of the gas atomization nozzle for gas flow behavior. Comput Fluids (2011) 42:37–43. 10.1016/j.compfluid.2010.10.013
6.
SiCZhangXWangJLiY. 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. 10.1007/s12613-014-0951-4
7.
XuZWangTCheZ. Droplet deformation and breakup in shear flow of air. Phys Fluids (2020) 32:052109. 10.1063/5.0006236
8.
TaylorGI. 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.
TaylorGI. 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.
BartokWMasonS. Particle motions in sheared suspensions: viii. singlets and doublets of fluid spheres. J Colloid Sci (1959) 14:13–26. 10.1016/0095-8522(59)90065-0
11.
RumscheidtF-DMasonS. Particle motions in sheared suspensions xii. deformation and burst of fluid drops in shear and hyperbolic flow. J Colloid Sci (1961) 16:238–61. 10.1016/0095-8522(61)90003-4
12.
TorzaSHenryCCoxRMasonS. Particle motions in sheared suspensions. xxvi. streamlines in and around liquid drops. J Colloid Interf Sci (1971) 35:529–43. 10.1016/0021-9797(71)90211-6
13.
ChaffeyCEBrennerH. A second-order theory for shear deformation of drops. J Colloid Interf Sci (1967) 24:258–69. 10.1016/0021-9797(67)90229-9
14.
Barthes-BieselDAcrivosA. Deformation and burst of a liquid droplet freely suspended in a linear shear field. J Fluid Mech (1973) 61:1–22. 10.1017/s0022112073000534
15.
HinchEAcrivosA. Long slender drops in a simple shear flow. J Fluid Mech (1980) 98:305–28. 10.1017/s0022112080000171
16.
KaramHBellingerJ. Deformation and breakup of liquid droplets in a simple shear field. Ind Eng Chem Fundamentals (1968) 7:576–81. 10.1021/i160028a009
17.
FlumerfeltRW. Drop breakup in simple shear fields of viscoleastic fluids. Ind Eng Chem Fundamentals (1972) 11:312–8. 10.1021/i160043a005
18.
StoneHABentleyBLealL. An experimental study of transient effects in the breakup of viscous drops. J Fluid Mech (1986) 173:131–58. 10.1017/s0022112086001118
19.
GuidoSVilloneM. Three-dimensional shape of a drop under simple shear flow. J rheology (1998) 42:395–415. 10.1122/1.550942
20.
GraceHP. 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. 10.1080/00986448208911047
21.
StoneHALealLG. The influence of initial deformation on drop breakup in subcritical time-dependent flows at low Reynolds numbers. J Fluid Mech (1989) 206:223–63. 10.1017/s0022112089002296
22.
VananroyeAVan PuyveldePMoldenaersP. Effect of confinement on droplet breakup in sheared emulsions. Langmuir (2006) 22:3972–4. 10.1021/la060442+
23.
VananroyeAVan PuyveldePMoldenaersP. Effect of confinement on the steady-state behavior of single droplets during shear flow. J rheology (2007) 51:139–53. 10.1122/1.2399089
24.
KennedyMRPozrikidisCSkalakR. Motion and deformation of liquid drops, and the rheology of dilute emulsions in simple shear flow. Comput Fluids (1994) 23:251–78. 10.1016/0045-7930(94)90040-x
25.
TooseEGeurtsBKuertenJ. A boundary integral method for two-dimensional (non)-Newtonian drops in slow viscous flow. J non-newtonian Fluid Mech (1995) 60:129–54. 10.1016/0377-0257(95)01386-3
26.
UijttewaalWNijhofE. The motion of a droplet subjected to linear shear flow including the presence of a plane wall. J Fluid Mech (1995) 302:45–63. 10.1017/s0022112095004009
27.
LiJRenardyYYRenardyM. Numerical simulation of breakup of a viscous drop in simple shear flow through a volume-of-fluid method. Phys Fluids (2000) 12:269–82. 10.1063/1.870305
28.
AmaniABalcázarNCastroJOlivaA. Numerical study of droplet deformation in shear flow using a conservative level-set method. Chem Eng Sci (2019) 207:153–71. 10.1016/j.ces.2019.06.014
29.
XiHDuncanC. Lattice Boltzmann simulations of three-dimensional single droplet deformation and breakup under simple shear flow. Phys Rev E (1999) 59:3022–6. 10.1103/physreve.59.3022
30.
van der SmanRGVan der GraafS. Emulsion droplet deformation and breakup with lattice Boltzmann model. Comp Phys Commun (2008) 178:492–504. 10.1016/j.cpc.2007.11.009
31.
FarokhiradSLeeTMorrisJF. Effects of inertia and viscosity on single droplet deformation in confined shear flow. Commun Comput Phys (2013) 13:706–24. 10.4208/cicp.431011.260112s
32.
KomrakovaAShardtOEskinDDerksenJ. Lattice Boltzmann simulations of drop deformation and breakup in shear flow. Int J Multiphase Flow (2014) 59:24–43. 10.1016/j.ijmultiphaseflow.2013.10.009
33.
HuangBLiangHXuJ. Lattice Boltzmann simulation of binary three-dimensional droplet coalescence in a confined shear flow. Phys Fluids (2022) 34:032101. 10.1063/5.0082263
34.
WangDTanDSKhooBCOuyangZPhan-ThienN. A lattice Boltzmann modeling of viscoelastic drops’ deformation and breakup in simple shear flows. Phys Fluids (2020) 32:123101. 10.1063/5.0031352
35.
ZongYZhangCLiangHWangLXuJ. Modeling surfactant-laden droplet dynamics by lattice Boltzmann method. Phys Fluids (2020) 32:122105. 10.1063/5.0028554
36.
MonaghanJ. Smoothed particle hydrodynamics and its diverse applications. Annu Rev Fluid Mech (2012) 44:323–46. 10.1146/annurev-fluid-120710-101220
37.
YeTPanDHuangCLiuM. Smoothed particle hydrodynamics (sph) for complex fluid flows: recent developments in methodology and applications. Phys Fluids (2019) 31:011301. 10.1063/1.5068697
38.
MorrisJP. Simulating surface tension with smoothed particle hydrodynamics. Int J Numer Methods Fluids (2000) 33:333–53. 10.1002/1097-0363(20000615)33:3<333::aid-fld11>3.0.co;2-7
39.
HuXYAdamsNA. A multi-phase sph method for macroscopic and mesoscopic flows. J Comput Phys (2006) 213:844–61. 10.1016/j.jcp.2005.09.001
40.
WangZ-BChenRWangHLiaoQZhuXLiS-Z. An overview of smoothed particle hydrodynamics for simulating multiphase flow. Appl Math Model (2016) 40:9625–55. 10.1016/j.apm.2016.06.030
41.
ZhangM. Simulation of surface tension in 2d and 3d with smoothed particle hydrodynamics method. J Comput Phys (2010) 229:7238–59. 10.1016/j.jcp.2010.06.010
42.
TartakovskyAMPanchenkoA. Pairwise force smoothed particle hydrodynamics model for multiphase flow: surface tension and contact line dynamics. J Comput Phys (2016) 305:1119–46. 10.1016/j.jcp.2015.08.037
43.
YangQYaoJHuangZAsifM. A comprehensive sph model for three-dimensional multiphase interface simulation. Comput Fluids (2019) 187:98–106. 10.1016/j.compfluid.2019.04.001
44.
MoinfarZVahabiSVahabiM. 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. 10.1108/hff-01-2022-0067
45.
VahabiM. The effect of thixotropy on deformation of a single droplet under simple shear flow. Comput Math Appl (2022) 117:206–15. 10.1016/j.camwa.2022.04.023
46.
SaghatchiROzbulutMYildizM. Dynamics of double emulsion interfaces under the combined effects of electric field and shear flow. Comput Mech (2021) 68:775–93. 10.1007/s00466-021-02045-x
47.
HirschlerMOgerGNiekenULe TouzéD. Modeling of droplet collisions using smoothed particle hydrodynamics. Int J Multiphase Flow (2017) 95:175–87. 10.1016/j.ijmultiphaseflow.2017.06.002
48.
XuXTangTYuP. A modified sph method to model the coalescence of colliding non-Newtonian liquid droplets. Int J Numer Methods Fluids (2020) 92:372–90. 10.1002/fld.4787
49.
ZhangASunPMingF. An sph modeling of bubble rising and coalescing in three dimensions. Comp Methods Appl Mech Eng (2015) 294:189–209. 10.1016/j.cma.2015.05.014
50.
BrackbillJUKotheDBZemachC. A continuum method for modeling surface tension. J Comput Phys (1992) 100:335–54. 10.1016/0021-9991(92)90240-y
51.
LafaurieBNardoneCScardovelliRZaleskiSZanettiG. Modelling merging and fragmentation in multiphase flows with surfer. J Comput Phys (1994) 113:134–47. 10.1006/jcph.1994.1123
52.
PriceDJ. Smoothed particle hydrodynamics and magnetohydrodynamics. J Comput Phys (2012) 231:759–94. 10.1016/j.jcp.2010.12.011
53.
MorrisJPFoxPJZhuY. Modeling low Reynolds number incompressible flows using sph. J Comput Phys (1997) 136:214–26. 10.1006/jcph.1997.5776
54.
ShethKSPozrikidisC. Effects of inertia on the deformation of liquid drops in simple shear flow. Comput Fluids (1995) 24:101–19. 10.1016/0045-7930(94)00025-t
55.
ZhouHPozrikidisC. The flow of suspensions in channels: single files of drops. Phys Fluids A (1993) 5:311–24. 10.1063/1.858893
56.
ShapiraMHaberS. Low Reynolds number motion of a droplet in shear flow including wall effects. Int J multiphase flow (1990) 16:305–21. 10.1016/0301-9322(90)90061-m
57.
SibilloVPasquarielloGSimeoneMCristiniVGuidoS. Drop deformation in microconfined shear flow. Phys Rev Lett (2006) 97:054502. 10.1103/physrevlett.97.054502
Summary
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
Volume
11 - 2023
Edited by
He Li, University of Georgia, United States
Reviewed by
Kaixuan Zhang, Nankai University, China
Jiayi Zhao, University of Shanghai for Science and Technology, China
Updates
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
Disclaimer
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.