PHYSICAL REVIEW E 107, 035106 (2023)

Effect of droplet deformability on shear thinning in a cylindrical channel

Danilo P. F. Silva ,* Rodrigo C. V. Coelho , Margarida M. Telo da Gama , and Nuno A. M. Araújo
Centro de Física Teórica e Computacional, Faculdade de Ciências, Universidade de Lisboa, 1749-016 Lisboa, Portugal
and Departamento de Física, Faculdade de Ciências, Universidade de Lisboa, 1749-016 Lisboa, Portugal

(Received 10 August 2022; revised 20 December 2022; accepted 28 February 2023; published 21 March 2023)

Droplets suspended in fluids flowing through microchannels are often encountered in different contexts and
scales, from oil extraction down to microfluidics. They are usually flexible and deform as a product of the
interplay between flexibility, hydrodynamics, and interaction with confining walls. Deformability adds distinct
characteristics to the nature of the flow of these droplets. We simulate deformable droplets suspended in a
fluid at a high volume fraction flowing through a cylindrical wetting channel. We find a discontinuous shear
thinning transition, which depends on the droplet deformability. The capillary number is the main dimensionless
parameter that controls the transition. Previous results have focused on two-dimensional configurations. Here
we show that, in three dimensions, even the velocity profile is different. To perform this study, we improve and
extend to three dimensions a multicomponent lattice Boltzmann method which prevents the coalescence between
the droplets.

DOI: 10.1103/PhysRevE.107.035106

I. INTRODUCTION ratio between the droplet and the surrounding (continuous)

fluid [7]. When the viscosity ratio is equal to unity, droplets
Particles and droplets suspended in fluids flowing through
accumulate at the center of the channel, while for higher
channels are often encountered in industrial and natural pro-
values of the viscosity ratio the droplets may accumulate at
cesses. The flow of suspensions, microfluidic devices, and
the walls, leaving a droplet depleted zone in the middle. A
drug delivery are typical applications [1]. For instance, mi-
more detailed study on particle migration and concentration
crofluidics has matured and provided new techniques for the
was conducted using Stokesian dynamics in three dimensions
control of soft suspensions. For example, this has allowed re-
[8], which allows accessing much smaller Reynolds numbers
searchers to simulate microgeometries with enhanced control
(Re  1). However, it appears to be restricted in terms of
of droplets for a variety of applications such as the creation
applications since it is mainly limited to the study of spherical
of emulsions and enhancement of mixing reagents, among
particles in Stokes flow.
others [2].
The lattice Boltzmann method (LBM) provides a way to
Generally, these complex suspensions consist of soft de-
include surface tension. The LBM is commonly used to study
formable particles which exhibit different flow behavior when
deformable particles under different conditions [9–12]. For in-
compared to their hard counterparts [3,4]. Systems where
stance, it has been used to investigate how three-dimensional
particles deform due to hydrodynamic interactions have not
(3D) particle migration and distribution are affected by de-
been systematically studied, especially in three dimensions,
formability and inertia at different Reynolds numbers, in
and a deeper understanding of the key physical mechanisms
particular how distinct flow focusing emerges at increasing
underlying their single and collective behavior is needed.
Reynolds number for strongly deformable particles [13]. In
There has been a growing effort to model the role of
those limits, a nonmonotonic behavior of the apparent viscos-
deformability in the dynamics of soft suspensions in chan-
ity of the suspension was reported. The LBM is an efficient
nels. Several numerical methods have been developed. A
algorithm allowing for extreme simulations such as with bil-
well-known method is the boundary integral method (BIM)
lions of grid points [14].
and its variants, especially used in the study of multiphase
The flow of suspensions exhibits a wide variety of non-
problems [5]. A major challenge to the use of these methods
Newtonian behavior such as shear thinning, which has been
is how to accurately and effectively include surface tension.
observed in experimental studies of hard suspensions such as
However, the BIM is applicable to certain types of flows such
glass balls in oil [15] and cornstarch particles in water [16].
as Stokes or potential flows. A popular variant of the BIM
A striking phenomenon is the observation of shear thinning in
is the method of interfacial dynamics [6], which is suitable
noncolloidal suspensions reported in numerous experiments
for the simulation of incompressible soft particles. Using this
[17–20]. In these experiments, particles larger than 20 µm
method, it has been shown that ordered droplets driven by
were used and the effects of particle migration, sedimentation,
a pressure difference in a periodic channel accumulate at
and confinement were investigated where Brownian effects
different regions of the channel depending on the viscosity
are negligible [21].
Previous experimental studies have shown that shear thin-
ning in soft suspensions of droplets depends strongly on the
* size of the droplets [22–24]. These experiments typically
[email protected]

2470-0045/2023/107(3)/035106(9) 035106-1 ©2023 American Physical Society

SILVA, COELHO, DA GAMA, AND ARAÚJO PHYSICAL REVIEW E 107, 035106 (2023)

involve oil droplets in water studied over a large range of TABLE I. Velocity vectors √ and weights for the D3Q41 lattice.
volume fractions. In parallel with experiments, shear thinning The speed of sound cs2 is 1 − 2/5.
has also been reported in emulsions through numerical simu-
lations, for example, using interface tracking methods [25]. ci wi
The collective dynamics of a large number of deformable √
(0,0,0) 2(5045 −√1507 10)/2025
particles in microchannels was also investigated [26] via 3D (±1, 0, 0), (0, ±1, 0), (0, 0, ±1) 377/(5 10)√− (91/40)
simulations using an immersed boundary method. Addition- (±1, ±1, 0), (±1, 0, ±1), (0, ±1, ±1)
ally, the immersed boundary method has been used to address √− 17 10)/50
(55
(±1, ±1, ±1) (233 10 −√730)/1600
more complex situations such as the implementation of mem-
(±3, 0, 0), (0, ±3, 0), (0, 0, ±3) (295 − 92√ 10)/16200
brane viscosity in fluid-filled bodies like capsules [27,28].
(±3, ±3, ±3) (130 − 41 10)/129600
The trajectories and the plug-flow profile of the particles were
analyzed as a function of the deformability and channel size.
Other works [29,30] investigated the effect of elasticity and
confinement of capsules such as red blood cell suspensions deformability is studied for a dense suspension flow in a 3D
for different cell rigidity. They observed shear thinning be- cylindrical channel. We highlight the differences and similar-
havior and reported that at high enough capillary numbers the ities between the results in 2D and 3D flows. In Sec. IV we
effective viscosity of the suspension converges to the solvent make some final observations. The Appendix contains numer-
viscosity. These studies suggest that a key parameter when ical validation [35].
investigating the flow of deformable particles is the capillary
number, defined as the ratio of viscous forces to interfacial
tension forces. II. METHOD
Recently, a discontinuous shear thinning behavior was re- We extended to three dimensions a multicomponent model
ported in the flow of deformable droplets in two dimensions with frustrated coalescence. For flexible droplets, the motion
[31,32]. This behavior is associated with a nonequilibrium of the fluid is represented by a set of distribution functions
transition between a soft and a hard phase, which depends fk,i (x, t ) at position x and time t, where the subscripts k and
on the area fraction of droplets in suspension and the applied i denote the fluid component and discrete velocity directions,
pressure difference. Strong shear thinning behavior was re- respectively. The time evolution of fk,i (x, t ) is given by the
ported at a higher area fraction and it was further revealed that discrete Boltzmann equation
when the area fraction is higher than 0.5, the shear thinning
becomes discontinuous, i.e., there is a jump in the viscosity 1 eq
at a critical value of the forcing. In addition, the velocity (in fk,i (x + ci , t + 1) − fk,i (x, t ) = − [ fk,i (x, t ) − fk,i (x, t )]
τk
the flow direction) was measured and it was reported that the
velocity of the droplets remains close to that of the continuous + Fk,i , (1)
fluid throughout the channel. Simulations were performed
using a 2D hybrid LBM. A more recent study [33] inves- where τk is the relaxation time for each component and Fk,i
tigated and confirmed the robustness of this discontinuous is the forcing term. Unless otherwise stated, we express our
behavior. This study reported additional discontinuous jumps results in lattice units (l.u.), which means that the lattice spac-
in viscosity over a larger parameter range, showing how the ing x and the time step t are equal to one. The equilibrium
discontinuous shear thinning is preserved at lower values of distribution is given by
the confinement, defined as the ratio between channel height  
and droplet radius. To do so, they simulated a large 2D sys- eq ueq · ci (ueq · ci )2 (ueq )2
fk,i = ρk wi 1 + + − , (2)
tem with up to 500 droplets (approximately 0.85 of the area cs2 2cs4 2cs2
fraction).
The effect of the viscosity ratio as well as the ratio where cs is the speed of sound. The lattice used in the stream-
of droplet to channel radius was also investigated. While ing step is the D3Q41 lattice (see Table I for the vectors ci
discontinuous shear thinning behavior of droplets has been and weights wi ). As will be discussed, using higher-order
investigated over a wide parameter range in two dimensions, lattices to solve the Boltzmann equation reduces the spurious
these effects have not yet been identified and tested in 3D ge- velocities. Here ueq is the equilibrium velocity given by
ometries. In particular, it is elusive how deformability affects
 ρk uk
the discontinuous shear thinning and the overall fluid flow in k τk
3D systems. In this paper we simulate droplets in 3D flows u eq
=  ρk , (3)
and investigate the discontinuous shear thinning behavior as k τk

a function of the surface tension, which is related to droplet

deformability. We improve and extend to three dimensions a where ρk uk is the kth component of momentum. The Guo
multicomponent LBM with frustrated coalescence [34]. We forcing scheme [36] was adopted to implement the forces
analyze the shear thinning for different flow conditions and acting in a fluid as it yields a viscosity-independent surface
droplet deformability and discuss the main properties of the tension. Thus,
velocity profiles.    
1 ci − ueq (ci · ueq )ci
The paper is organized as follows. The lattice Boltzmann Fk,i = 1 − wi + · Fk , (4)
model is described in Sec. II. In Sec. III the effect of droplet 2τk cs2 cs4

035106-2
EFFECT OF DROPLET DEFORMABILITY ON SHEAR … PHYSICAL REVIEW E 107, 035106 (2023)

TABLE II. Velocity vectors and weights for the D3Q39 lattice.
The speed of sound cs2 is 2/3.

ci wi
(0,0,0) 1/12
(±1, 0, 0) 1/12
(±1, ±1, ±1) 1/27
(±2, 0, 0) 2/135
(±2, ±2, 0) 1/142 FIG. 1. Spurious velocities (arrows) around a circular interface
(±3, 0, 0) 1/1620 (red) of radius r = 20 at equilibrium for Gk,1 = −7.9 and Gk,2 = 4.9.
The largest spurious velocities are (a) 0.017 l.u. for the D3Q19 lattice
and (b) 0.0092 l.u. for the D3Q41 lattice. The velocity field is on the
where Fk is the sum of all the forces. The equilibrium velocity same scale in both panels. The ratio of the density inside and outside
is the same as in Eq. (2) [36]. Also, the droplet is 1.133 for D3Q19 and 1.062 for D3Q41.


40 
40
Fk
ρk = fk,i , ρk uk = fk,i ci + . (5) fluid. We use Gks = −0.35, resulting in a contact angle of 90◦
i=0 i=0
2
(neutral wetting).
The barycentric velocity u of the fluid mixture, i.e., the phys- The total force Fk in Eq. (4) is the sum of the external body
ical velocity, is given by forces Fkb (e.g., gravity, which we neglect here), the internal
 forces, and the solid boundary interaction force: Fk = Fkb +
ρk uk 
u= k , ρ= ρk , (6) Fkc + Fkr + Fks . For simplicity, we assume that GA,1 = GB,1 ,
ρ k GA,2 = GB,2 , and GA,B = GB,A . The parameters GA,1 , GB,2 ,
where ρ is the total density of the fluid mixture. etc., are the strength coefficients of the interaction forces.
There are three internal forces acting on the fluid: an intra- More specifically, a positive (or negative) strength coefficient
component short-range attractive force, a midrange repulsive represents a repulsive (or attractive) interaction. To have the
force, and a repulsive force between different components. desired repulsive and attractive forces Gk,1 < 0, Gk,2 > 0,
The midrange force prevents droplet coalescence and for that and Gk,k > 0. In the competing force Fkc , the attractive force
we use the D3Q39 lattice [37] (see Table II). For the attractive must overcome the repulsive force to form droplets, so we
short-range and repulsive forces between components, we use set |Gk,1 | > |Gk,2 |. Details on how the strength parameters
the D3Q41 lattice [37]. The intracomponent forces acting on Gk,1 , Gk,2 , and Gk,k and reference density ρ0 are related to
the same component are the surface tension can be found in [39].
Multiphase and multicomponent LBMs suffer from spuri-

40
ous velocities caused by an imbalance between the interaction
Fkc = − Gk,1 ψk (x) wi ψk (x + ci )ci forces. These velocities can increase if the viscosity ratio
i=0
deviates much from one or if the surface tension is high.

38 Thus, the simulations might become unstable for realistic
− Gk,2 ψk (x) p j ψk (x + c j )c j , (7) physical parameters. Previous works have demonstrated that
j=0 by increasing the isotropy of the lattice used to calculate the
interaction forces in the pseudopotential models it is possible
where ψk (ρk ) = ρ0 (1 − e−ρk /ρ0 ) is a pseudopotential with a
to reduce the spurious velocities [40–42]. In two dimensions
uniform reference density ρ0 . Additionally, Gk,1 and Gk,2
the lattice used for the midrange repulsive force was D2Q25
are the self-interaction forces within each component. The
[43,44]. In three dimensions the equivalent of the D2Q25
repulsive forces between the components are implemented as
lattice would be a D3Q125 lattice, which includes all the
usual [38]:
neighbors in the first and second belts. However, such a lat-
 
40 tice would be computationally expensive. Instead, we use the
Fkr = −ρk (x) Gkk wi ρk (x + ci )ci . (8) D3Q39 lattice for the midrange force (see Tables I and II).
k i=0 A nonphysical effect causes droplets to become glued in the
The adhesion force with the solid boundary is given by direction of the diagonal vectors of the lattice. We also observe
this effect in the usual 2D models [39,41,45]. Peng et al. [46]

solid
showed that using lattices with higher isotropy in the stream-
Fks = −Gks ρk (x) wi s(x + ci )ci , (9) ing process significantly reduces the spurious velocities. In
i=0 our simulations, we use an eighth-order isotropic lattice in the
where s(x) is the switch function, which takes the values 0 streaming step. By contrast, the D3Q19 lattice has only fourth-
and 1 for fluid and solid nodes, respectively, and Gks is the order isotropy. In Fig. 1 we compare the spurious velocities
interaction strength between fluid component k and the solid when these two lattices are used in the streaming step. We note
boundary. The subscript k represents the component A or B. that it reduces spurious velocities by nearly one-half. Notice
Additionally, to the force (9), other forces [Eqs. (8) and (7)] in Fig. 1 that spurious velocities are stronger in the diagonal
are applied to the solid nodes by using a virtual solid density, direction of the grid, which explains why the droplets stick to
which in our case is the initial density of the surrounding each other at these points but not along the grid axes.

035106-3
SILVA, COELHO, DA GAMA, AND ARAÚJO PHYSICAL REVIEW E 107, 035106 (2023)

will coalesce immediately due to the presence of diffusive in-

terfaces. Thus, we allow a small distance between the droplets
(three to five lattice points) to avoid this. So, while not fully
jammed, the droplets are in close contact with each other. Our
simulations are parallelized using OPENMP and we use eight
cores for our simulations, each with 2.3 GHz. Our simulations
deliver approximately 1.8 × 106 lattice updates per second.
The length and radius of the cylinder are 95 and 47.5 l.u.,
respectively. We impose no-slip boundary conditions in ve-
locity and a neutral wetting boundary condition which creates
a layer of wetting droplets along the wall. The fluid is driven
by a fixed external force. Figure 2 shows a snapshot of the
initial setup.
To study the effect of droplet deformability we first vary
the external force F in the z direction and measure the effec-
tive relative viscosity μr = μ/μ0 , where μ is the (apparent)
dynamic viscosity of the fluid with droplets and μ0 is the
dynamic viscosity of the fluid without droplets. The effective
viscosity in terms of the volume flow rate

Q0
μr = (10)
Q

is measured, where Q0 is the flow rate of the fluid without

droplets and Q is theflow rate with droplets. The flow rate
Q is defined as Q = A vz dA, where vz is the axial velocity
and A is the cross-sectional area. This measurement of μr is
obtained for different values of surface tension γ (by vary-
FIG. 2. Initial setup. (a) Here 105 droplets are arranged in a
ing ρ0 ). The surface tension may be used as a measure of
hexagonally non-close-packed system with a volume fraction of
the droplet deformability; however, since we compare forces
0.61. A force F is imposed in the z direction and varied. The walls
resulting from fluid motion with forces resulting from surface
have a neutral wetting condition [47]. The droplets (red) are inside
the cylindrical channel (gray). The system is periodic in the direction
tension, the capillary number Ca = μ0 v/γ is a more appro-
of the flow. The parameters for the interactions are Gk,1 = −7.9 and priate dimensionless parameter, where v is the characteristic
Gk,2 = 4.9 and the viscosity ratio of the droplet and the continuous velocity (taken as the maximum velocity of the fluid). We vary
fluid is set to one. (b) Top view. (c) Front view. it by changing the surface tension. The flow rate is measured
along a plane perpendicular to the flow (XY ). We measure
the flow velocity on each lattice node along this plane to
Regarding computational cost, in the present model the obtain the overall flow rate Q at a particular time step. Due
number of droplets does not increase the computational cost to the discretization of the droplets on the lattice, the value of
at a given system size. This is because the repulsion between Q fluctuates around an average value, as the intersection of the
droplets is based on the density of each lattice Boltzmann flowing droplets with the pane is not fixed. Then Q is averaged
node. In other models [27,28], an increase in the number over subsequent instants of time (over 15 frames with 5000
of immersed flexible fluid-filled bodies leads to an increase iterations between them) in the steady state. We study the
in computation, which is significant for a larger number of relationship between effective viscosity and surface tension.
immersed objects. However, these do not suffer from spurious We plot the dependence of μr on the surface tension in Fig. 3.
velocities as in the present model. They describe funda- In Fig. 3 we see the discontinuous behavior which happens at
mentally different mechanisms due to the different boundary larger values of F as we increase the surface tension. In ad-
conditions on the droplet and the capsule. While the cap- dition, we see that the viscosity curves collapse when plotted
sule has a fluid-solid boundary, the droplet has a fluid-fluid against the capillary number. The discontinuous drop of the
boundary. viscosity suggests shear thinning. For soft particles, deforma-
bility has been shown to promote shear thinning analogous to
cells flowing in microenvironments [48]. We report data col-
III. FLOW OF DROPLETS
lapse in both regimes by taking the maximal velocity in each
To understand the effect of droplet deformability on the respective regime (before or after shear thinning) as shown
overall flow of concentrated suspensions, we simulate a in Fig. 3(b). Shear thinning is also observed in non-Brownian
densely packed system of 105 droplets of radius r = 8 lat- hard-sphere suspensions [49] and the decrease in viscosity has
tice units in a 3D cylindrical channel with periodic boundary been studied under more controlled conditions for those par-
conditions in the z direction (see Fig. 2). The droplets are ticles [17,50]. The results are nearly independent of the size
not close packed, resulting in a total volume fraction of 0.61. of the channel [35]. For a channel with 170 nodes, the largest
Initially, the droplets cannot be close packed; otherwise they difference between the two curves of relative viscosity is less

035106-4
EFFECT OF DROPLET DEFORMABILITY ON SHEAR … PHYSICAL REVIEW E 107, 035106 (2023)

FIG. 3. Discontinuous shear thinning for different values of γ . (a) Discontinuous shear thinning occurs for progressively larger values of
the external force which drives the flow. The black arrow indicates the threshold viscosity μ∗r , i.e., the relative viscosity before shear thinning.
(b) Relative viscosity μr (normalized by the threshold viscosity μ∗r ) as a function of the capillary number Ca. The transition is at approximately
0.025. Simulations were carried out for different values of surface tension γ .

than 2%. This information can be found in the Supplemental in flow is however more complex. The flow in this channel
Material [35], which includes Refs. [51,52]. is described by small characteristic dimensions and veloci-
We define F∗ as the threshold force applied to the system ties. Consequently, the flow is characterized by low Reynolds
necessary to trigger discontinuous shear thinning (see Fig. 4). number Re and is laminar. Low-Re flows also indicate that
We observe that at higher values of γ the force required to the viscous forces are relevant. While the flow is laminar, the
trigger discontinuous shear thinning F∗ increases in a roughly regime is not Stokesian. Viscous forces usually arise due to
linear fashion until it saturates. This is shown in Fig. 4(a). In friction, for example, near the walls. The presence of both sus-
fact, as we increase γ the droplets become less flexible and pended and adhered droplets in a 3D channel leads to distinct
eventually reach a state similar to that of hard particles. We ex- profiles. In particular, the droplets adhered at the channel wall
pect in that limit shear thickening, as reported, for example, in offer resistance to the flow of other droplets up to a certain
Ref. [15] for hard suspensions (cornstarch in water and glass point.
spheres in oil) between parallel plates (rotating top plate). We plot the velocity distribution as a function of
Figure 4(b) shows that as the droplets become less deformable the nondimensional radius r/R. We measure the velocities
(increasing γ ), the threshold viscosity μ∗r increases. of the continuous phase and of the droplet phase separately
The simulations are performed for a fixed value of the before (high μr ) and after (low μr ) the discontinuous shear
ratio between the channel diameter and the droplet radius. In thinning. An average in space is taken between the radial
two dimensions [33], it is reported that the effective relative distances r and r + r, at a fixed r in steady-state flow.
viscosity μr increases with this ratio. We do not vary this The velocities of the droplet phase vd (r) and of the contin-
ratio, but we do not expect that the effect of confinement uous phase v f (r) are identified by means of the density field
is qualitatively different in three dimensions. The fully de- for the droplet component. A threshold of 0.6 is considered,
veloped velocity profile across a cylindrical channel with no which is half of the density of the droplets, 1.2. Regions
droplets can be computed analytically from the Navier-Stokes with density above this threshold are considered part of the
equations. For a cylindrical channel with a circular cross droplets, while those below are part of the continuous fluid.
section of radius R, the solution is a parabolic velocity distri- The diffusive interface has a thickness of 2.84 l.u. The r
bution [53]. The velocity distribution for droplets suspended value is 6 l.u. for droplets and 4 l.u. for continuous fluid. When
shear thinning has occurred, the v f (r) and vd (r) are similar
and are consistent with plug flow behavior. Before shear thin-
ning the flow is no longer well described as plug flow (as
reported in two dimensions [32]), i.e., the radial profile of
the continuous fluid velocity defined as v f (r) = v f (r, θ , z)
in cylindrical coordinates is higher than that of the droplets
phase vd (r) = vd (r, θ , z). The average is taken along z and
angle θ . After shear thinning v f (r) ≈ vd (r). We also observe
that for a low number of droplets, no shear thinning occurs
(plot not shown) as in two dimensions. The situation before
shear thinning, however, contrasts with that of a 2D channel
FIG. 4. (a) Threshold force F∗ required to trigger discontinuous where the system flows slowly and is in a nearly jammed
shear thinning as a function of γ . Here F∗ increases roughly linearly state. While previous studies in two dimensions [31,33] have
with the droplets’ surface tension. (b) Corresponding threshold vis- not reported differences between the fluid and droplet ve-
cosity μ∗r as a function of γ . locity before shear thinning, we observe differences in a 3D

035106-5
SILVA, COELHO, DA GAMA, AND ARAÚJO PHYSICAL REVIEW E 107, 035106 (2023)

FIG. 5. Velocity profile given by the dimensionless Reynolds

number Re = VzνD , where D is the diameter of the channel, ν is the
kinematic viscosity, and Vz  is the average velocity in the z direction.
The capillary numbers are Ca = 0.017 and 0.0314. The droplet and
fluid velocity are plotted separately. FIG. 7. Relative velocity fields and droplet shape cross
section (a) before and (b) after discontinuous shear thinning. The
capillary numbers are (a) Ca = 0.017 and (b) Ca = 0.0314. The
geometry. The velocity of the continuous fluid v f (r) is sig- relative velocity is obtained by subtracting the velocity at the center
nificantly higher than vd (r), as shown in Fig. 5. We point out of the droplet and it is on the same scale in both images. The cross
that for a 3D channel there is space between the droplets (see section is obtained along the ZY plane. Large black arrows indicate
Fig. 6) for the continuous fluid to flow in contrast to the 2D the direction of flow. (c) Relative position of the sample droplet in
case with a similar arrangement. While the velocity is con- question marked as A (front view, i.e., XY plane).
tinuous at the interface, the gaps between some droplets are
large enough [see Fig. 2(c)] so that before shear thinning the
average continuous fluid flows v f (r) with a higher velocity in the case of a droplet), in particular, the difference between
than the average velocity of the droplets phase vd (r). In the droplet and continuous fluid velocities. The flow of the
Fig. 5, before shear thinning v f (r) and vd (r) are 0.000 88 continuous fluid between 2D droplets is affected by geometri-
and 0.000 21, respectively. Thus, the continuous fluid flows cal restrictions caused by narrow passages and dead-end pores
4.2 times faster than the droplets for this particular force. between droplets [54,55], especially at high area fractions.
As we increase the forcing, the fluid eventually pushes By contrast, 3D geometries with droplets allow significant
the droplets, leading to the discontinuous shear thinning. The interconnected space between the droplets and therefore the
droplets and fluid then flow with matching velocities. We continuous fluid can flow through it. The importance of this
highlight that, although simulations in 2D geometries can effect and the role it plays in the hydrodynamics of droplets
capture the essential flow qualities such as discontinuous are not well understood.
shear thinning and plug flow behavior, 3D geometries give The deformation of the droplets also causes distinct in-
more realistic results (in two dimensions a circular cross ternal flows. Figures 7(a) and 7(b) show cross sections of
section would be representing a cylinder and not a sphere like a nondeformed droplet before and a deformed droplet after
shear thinning, respectively. The droplet in question is next to
the wetting layer of droplets. Droplet deformation is naturally
associated with shear thinning as the deformation promotes
the flow of droplets. In line with the results in two dimensions
[31], not only do the droplets deform but internal currents are
set up within the deformed droplets. Notice that the internal
flow is almost nonexistent before shear thinning when the
droplets are hardly deformed.

IV. CONCLUSION
We have simulated the flow of droplets in a 3D channel. For
progressively higher values of the external force, we observed
discontinuous shear thinning. We studied in this transition
the effect of the surface tension, which is proportional to the
droplet deformability. At higher surface tension the droplets
FIG. 6. Streamlines (blue) for a central region of the channel just are less deformable and thus larger values of force are required
before shear thinning, which indicates that the continuous fluid flows for discontinuous shear thinning to occur. We observed that
between droplets (red). The system has Ca = 0.025. this transition occurs at a given capillary number. We also

035106-6
EFFECT OF DROPLET DEFORMABILITY ON SHEAR … PHYSICAL REVIEW E 107, 035106 (2023)

FIG. 8. Pressure difference p as a function of the inverse of the

droplet radius r. We simulate a static droplet for a fixed ρ0 (reference
density) and let it reach the steady state. The symbols represent the
pressure difference inside and outside the droplet in our simulations.
The solid lines are linear fits corresponding to Eq. (A1), where the
slope is γ .

noticed that μr increases with the surface tension. This is to be

expected as the droplets become less deformable. Simulations
for low volume fractions do not exhibit any discontinuous
shear thinning in line with 2D studies. We analyzed the ve-
locity profiles before and after shear thinning and noticed that
after shear thinning the velocity profile is that of plug flow
similar to 2D flows. However, closer inspection of the velocity
profile before shear thinning revealed an essential feature that
distinguishes our results from those in two dimensions. We
noticed that the fluid and droplet velocities are different, with FIG. 9. Disjoining pressure as a function of the distance
the fluid velocity being higher and that of the droplets being h between two flat interfaces with fixed Gk,1 = −3.5 and Gk,2 = 2.5
near zero. Additionally, we were able to observe internal cur- achieved with competing forces Fc given by Eq. (7), which provides
rents in the droplets after shear thinning. A possible extension the repulsive force to avoid coalescence. A positive value indicates
of this study is to analyze how the viscosity ratio (between the repulsion, while a negative value indicates attraction. (a) Disjoining
droplets and the fluid) affects the overall flow. pressure for different values of the reference density ρ0 . (b) Dis-
To carry out this study we extended to three dimensions joining pressure for different values of the viscosity ratio M with
a previous multicomponent LBM that prevents coalescence ρ0 = 1.00.
between the droplets. We measured the disjoining pressure,
which is independent of the viscosity ratio. We were able to
reduce spurious velocities which cause unphysical effects resistance (the extra resistance due to the presence of an object
in the simulations such as sticking droplets and coales- in a channel) has been suggested as a parameter for character-
cence when a collision occurs in certain directions. This was izing such flows [56]. It would be interesting to calculate the
achieved by using higher-order lattices in the streaming step. dependence of the hydrodynamic resistance on the capillary
We note some prospects for experimental studies of de- number in experiments with a setup similar to the one pre-
formable droplets in microchannels. Recently, hydrodynamic sented in this paper.

TABLE III. Surface tension γ for different reference densities ρ0

in lattice units. The values of γ are obtained from linear fits in Fig. 8.

ρ0 γ
1.000 0.026
1.025 0.023
1.050 0.021
1.075 0.019
1.100 0.017
1.150 0.012 FIG. 10. Snapshots of a collision between two droplets illustrat-
ing the frustrated coalescence.

035106-7
SILVA, COELHO, DA GAMA, AND ARAÚJO PHYSICAL REVIEW E 107, 035106 (2023)

ACKNOWLEDGMENTS LBM simulations with different values for the droplet radius
were performed. As seen in Fig. 8, the pressure difference
We acknowledge financial support from the Portuguese
p increases linearly with the inverse of the droplet radius
Foundation for Science and Technology under Contracts No.
r (curvature). To vary the surface tension, we fix the param-
PTDC/FIS-MAC/28146/2017 (LISBOA-01-0145-FEDER-
eters Gk,1 and Gk,2 and change the uniform reference density
028146), No. UIDB/00618/2020, No. UIDP/00618/2020,
ρ0 (effectively varying the interaction forces). The solid lines
and No. 2020.08525.BD.
represent results calculated from Laplace’s law and the sym-
There are no conflicts of interest to declare.
bols are obtained from the simulations. The slope of the solid
lines is γ . The simulation results are consistent with Laplace’s
APPENDIX: NUMERICAL TESTS law. The values of γ are shown in Table III. The ratio of the
In this Appendix we analyze the properties of the 3D total densities inside and outside the droplets for the values
multicomponent method. We measure the surface tension, of ρ0 = 1.000 and 1.150 are 1.058 and 1.062, respectively.
which is related to the deformability of the droplets, and When the parameter GAB exceeds a certain threshold it gives
we analyze the disjoining pressure. To measure the surface rise to stable interfaces between fluids A and B with positive
tension, we perform the Laplace test in two dimensions. At surface tension. However, when the droplets approach each
equilibrium, the curved interface of a static droplet increases other, a thin film is formed, leading to coalescence. The phase
the pressure inside the droplet. The radius, pressure differ- separation promotes negative disjoining pressure and thus the
ence, and surface tension must satisfy the Young-Laplace competing interactions prevent the coalescence of neighbor-
equation [38] ing droplets and give rise to positive disjoining pressure as
seen in Fig. 9. This mechanism has been used in other stud-
γ
p = pin − pout = , (A1) ies to simulate noncoalescing droplets [34,57,58]. We further
r show this mechanism in action for 3D droplets, in Fig. 10. The
where r is the radius of the droplet, pin and pout are the disjoining pressure is also independent of the viscosity ratio
pressures inside and outside the droplet, respectively, and M defined as the viscosity of the droplet over the viscosity of
γ is the surface tension. To test Laplace’s law, a series of the continuous fluid [59].

[1] Q. Xu, M. Hashimoto, T. T. Dang, T. Hoare, D. S. Kohane, [18] S. Dai, E. Bertevas, F. Qi, and R. I. Tanner, J. Rheol. 57, 493
G. M. Whitesides, R. Langer, and D. G. Anderson, Small 5, (2013).
1575 (2009). [19] Y. Lin, G. W. H. Tan, N. Phan-Thien, and B. C. Khoo, J. Non-
[2] S. Sohrabi, N. Kassir, and M. K. Moraveji, RSC Adv. 10, 27560 Newtonian Fluid Mech. 212, 13 (2014).
(2020). [20] A. Singh and P. R. Nott, J. Fluid Mech. 490, 293 (2003).
[3] N. Koumakis, A. Pamvouxoglou, A. S. Poulos, and G. [21] N. A. M. Araújo, L. M. C. Janssen, T. Barois, G. Boffetta, I.
Petekidis, Soft Matter 8, 4271 (2012). Cohen, A. Corbetta, O. Dauchot, M. Dijkstra, W. M. Durham,
[4] K. van der Vaart, Y. Rahmani, R. Zargar, D. Bonn, and P. Schall, A. Dussutour, S. Garnier, H. Gelderblom, R. Golestanian, L.
J. Rheol. 57, 1195 (2013). Isa, G. H. Koenderink, H. Löwen, R. Metzler, M. Polin, C. P.
[5] T. Hou, J. Lowengrub, and M. Shelley, J. Comput. Phys. 169, Royall, A. Łarić et al., Soft Matter 19, 1695 (2023).
302 (2001). [22] Y. Otsubo and R. K. Prud’homme, Rheol. Acta 33, 303
[6] C. Pozrikidis, J. Comput. Phys. 169, 250 (2001). (1994).
[7] H. Zhou and C. Pozrikidis, Phys. Fluids 6, 80 (1994). [23] R. Pal, AIChE J. 42, 3181 (1996).
[8] P. R. Nott and J. F. Brady, J. Fluid Mech. 275, 157 (1994). [24] R. Pal, J. Colloid Interface Sci. 225, 359 (2000).
[9] B. Dünweg and A. J. C. Ladd, Adv. Polym. Sci. 221, 89 (2009). [25] M. Loewenberg and E. J. Hinch, J. Fluid Mech. 321, 395
[10] T. Krüger, F. Varnik, and D. Raabe, Comput. Math. Appl. 61, (1996).
3485 (2011). [26] S. K. Doddi and P. Bagchi, Phys. Rev. E 79, 046318 (2009).
[11] Y. Chen, RSC Adv. 4, 17908 (2014). [27] F. Guglietta, M. Behr, L. Biferale, G. Falcucci, and M.
[12] R. C. V. Coelho, D. P. F. Silva, A. M. R. Maschio, M. M. T. da Sbragaglia, Philos. Trans. R. Soc. A 379, 20200395 (2021).
Gama, and N. A. M. Araújo, Phys. Fluids 35, 013304 (2023). [28] F. Guglietta, M. Behr, L. Biferale, G. Falcucci, and M.
[13] T. Krüger, B. Kaoui, and J. Harting, J. Fluid Mech. 751, 725 Sbragaglia, Soft Matter 16, 6191 (2020).
(2014). [29] G. R. Lázaro, A. Hernández-Machado, and I. Pagonabarraga,
[14] G. Falcucci, G. Amati, P. Fanelli, V. K. Krastev, G. Polverino, Soft Matter 10, 7195 (2014).
M. Porfiri, and S. Succi, Nature (London) 595, 537 (2021). [30] G. R. Lázaro, A. Hernández-Machado, and I. Pagonabarraga,
[15] E. Brown and H. M. Jaeger, Phys. Rev. Lett. 103, 086001 Soft Matter 10, 7207 (2014).
(2009). [31] M. Foglino, A. N. Morozov, O. Henrich, and D. Marenduzzo,
[16] A. Fall, F. Bertrand, D. Hautemayou, C. Mezière, P. Phys. Rev. Lett. 119, 208002 (2017).
Moucheront, A. Lemaître, and G. Ovarlez, Phys. Rev. Lett. 114, [32] M. Foglino, A. N. Morozov, and D. Marenduzzo, Soft Matter
098301 (2015). 14, 9361 (2018).
[17] I. E. Zarraga, D. A. Hill, and D. T. Leighton, J. Rheol. 44, 185 [33] L. Fei, A. Scagliarini, K. H. Luo, and S. Succi, Soft Matter 16,
(2000). 651 (2020).

035106-8
EFFECT OF DROPLET DEFORMABILITY ON SHEAR … PHYSICAL REVIEW E 107, 035106 (2023)

[34] R. Benzi, S. Chibbaro, and S. Succi, Phys. Rev. Lett. 102, [47] R. C. V. Coelho, C. B. Moura, M. M. T. da Gama, and N. A. M.
026002 (2009). Araújo, Int. J. Numer. Methods Fluids 93, 2570 (2021).
[35] See Supplemental Material at http://link.aps.org/supplemental/ [48] L. Lanotte, J. Mauer, S. Mendez, D. A. Fedosov, J.-M.
10.1103/PhysRevE.107.035106 for further validation and the Fromental, V. Claveria, F. Nicoud, G. Gompper, and M.
effect of system size and packing fraction. Abkarian, Proc. Natl. Acad. Sci. USA 113, 13289 (2016).
[36] Z. Guo, C. Zheng, and B. Shi, Phys. Rev. E 65, 046308 (2002). [49] M. Kroupa, M. Soos, and J. Kosek, Phys. Chem. Chem. Phys.
[37] S. S. Chikatamarla and I. V. Karlin, Phys. Rev. E 79, 046701 19, 5979 (2017).
(2009). [50] C. D. Cwalina and N. J. Wagner, J. Rheol. 58, 949 (2014).
[38] T. Krüger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva, [51] A. Vananroye, P. Van Puyvelde, and P. Moldenaers, J. Rheol.
and E. M. Viggen, The Lattice Boltzmann Method (Springer, 51, 139 (2007).
Cham, 2017). [52] A. Vananroye, P. J. A. Janssen, P. D. Anderson, P. Van Puyvelde,
[39] R. Benzi, M. Sbragaglia, S. Succi, M. Bernaschi, and S. and P. Moldenaers, Phys. Fluids 20, 013101 (2008).
Chibbaro, J. Chem. Phys. 131, 104903 (2009). [53] H. Schlichting and K. Gersten, Boundary-Layer Theory, 9th ed.
[40] X. Shan, Phys. Rev. E 73, 047701 (2006). (Springer, Berlin, 2017), pp. 101–142.
[41] M. Sbragaglia, R. Benzi, L. Biferale, S. Succi, K. Sugiyama, [54] A. Koponen, M. Kataja, and J. Timonen, Phys. Rev. E 56, 3319
and F. Toschi, Phys. Rev. E 75, 026702 (2007). (1997).
[42] A. Montessori, G. Falcucci, M. L. Rocca, S. Ansumali, and S. [55] A. Guariguata, M. A. Pascall, M. W. Gilmer, A. K. Sum, E. D.
Succi, J. Stat. Phys. 161, 1404 (2015). Sloan, C. A. Koh, and D. T. Wu, Phys. Rev. E 86, 061311
[43] G. Falcucci, G. Bella, G. Shiatti, S. Chibbaro, M. Sbragaglia, (2012).
and S. Succi, Commun. Comput. Phys. 2, 1071 (2007). [56] P. Sajeesh, M. Doble, and A. K. Sen, Biomicrofluidics 8,
[44] B. Dollet, A. Scagliarini, and M. Sbragaglia, J. Fluid Mech. 766, 054112 (2014).
556 (2015). [57] M. Sbragaglia, R. Benzi, M. Bernaschi, and S. Succi, Soft
[45] G. Falcucci, S. Ubertini, C. Biscarini, S. D. Francesco, D. Matter 8, 10773 (2012).
Chiappini, S. Palpacelli, A. D. Maio, and S. Succi, Commun. [58] R. Benzi, M. Sbragaglia, A. Scagliarini, P. Perlekar, M.
Comput. Phys. 9, 269 (2011). Bernaschi, S. Succi, and F. Toschi, Soft Matter 11, 1271 (2015).
[46] C. Peng, L. F. Ayala, O. M. Ayala, and L.-P. Wang, Comput. [59] L. Fei, A. Scagliarini, A. Montessori, M. Lauricella, S. Succi,
Fluids 191, 104257 (2019). and K. H. Luo, Phys. Rev. Fluids 3, 104304 (2018).

035106-9