1 Introduction
Emulsions are a biphasic liquid–liquid system in which the two fluids are partially or totally immiscible. They are present in many biological and industrial applications such as waste treatment, oil recovery and pharmaceutical manufacturing. Emulsions are relevant also in the field of colloidal science, where the accuracy and the control of the production process of functional materials rely on a knowledge of the complex microstructure of the suspension (Xia et al. Reference Xia, Gates, Yin and Lu2000).
The study of the rheology of suspensions can be traced back to the pioneering work of Einstein (Reference Einstein1906, Reference Einstein1911) who found that the effective viscosity
$\unicode[STIX]{x1D707}_{e}$
of a dilute suspension of rigid spheres increases linearly with the volume fraction
$\unicode[STIX]{x1D719}$
as
$\unicode[STIX]{x1D707}_{e}=\unicode[STIX]{x1D707}_{m}(1+5\unicode[STIX]{x1D719}/2)$
, where
$\unicode[STIX]{x1D707}_{m}$
is the viscosity of the matrix phase. Batchelor & Green (Reference Batchelor and Green1972) extended this relation including a term
$O(\unicode[STIX]{x1D719}^{2})$
, obtaining
$\unicode[STIX]{x1D707}_{e}=\unicode[STIX]{x1D707}_{m}(1+5\unicode[STIX]{x1D719}/2+5\unicode[STIX]{x1D719}^{2})$
. These two analytical expression can well predict the suspension effective viscosity only in the dilute regime, at higher volume fraction no closed theory exist and we still rely on empirical relations. Eilers (Reference Eilers1941) proposed the following expression for emulsions
$\unicode[STIX]{x1D707}_{e}=\unicode[STIX]{x1D707}_{m}[1+(5\unicode[STIX]{x1D719}/4)/(1-\unicode[STIX]{x1D719}/\unicode[STIX]{x1D719}_{c})]$
, which gives good approximation of the effective viscosity for low and high volume fractions. Here
$\unicode[STIX]{x1D719}_{c}$
represent the maximum packing volume fraction, which depends on the properties of the suspended phase e.g. shape and deformability. Taylor (Reference Taylor1932) was the first to account for the deformation of the particles by introducing the viscosity ratio
$\unicode[STIX]{x1D706}$
, defined as the disperse phase over the matrix phase viscosity, into the Einstein formula
$\unicode[STIX]{x1D707}_{e}=\unicode[STIX]{x1D707}_{m}[1+2.5\unicode[STIX]{x1D719}(\unicode[STIX]{x1D706}+0.4)/(\unicode[STIX]{x1D706}+1)]$
. Using perturbation analysis, Schowalter, Chaffey & Brenner (Reference Schowalter, Chaffey and Brenner1968) were the first to find that for a suspension of deformable particles the first normal stress difference
$N_{1}$
is positive and the second normal stress difference
$N_{2}$
is negative, and later Choi & Schowalter (Reference Choi and Schowalter1975) introduced a correction with
$O(\unicode[STIX]{x1D719}^{2})$
. Pal (Reference Pal2002, Reference Pal2003) derived expressions for the effective viscosity of infinitely dilute and concentrated emulsions using the effective medium approach. In general, all these previous relations are limited to the dilute regime for inertialess flows, are based on pair interaction between particles and typically assume limiting cases to model surface tensions effects (either going to zero or infinite). To overcome these limitations numerical simulations, providing access to all details of velocity and stresses in the system, play therefore an important role.
As concerns suspensions of rigid particles, many numerical investigations have been performed in the past years. For these systems the only relevant parameters are the solid volume fraction and the particle Reynolds number. Kulkarni & Morris (Reference Kulkarni and Morris2008) studied non-Brownian suspensions at finite particle Reynolds number, up to 5, using a lattice Boltzmann method. These authors found that inertia increases the particle contribution to the effective viscosity and breaks the fore–aft symmetry of the pair distribution function, see also Picano et al. (Reference Picano, Breugem, Mitra and Brandt2013) for an analysis of the relative particle motion in the presence of inertia. Mari et al. (Reference Mari, Seto, Morris and Denn2014) reported shear thickening behaviour when a suspension of particles changes from the contactless to contact dominated regime and relate this effect to jamming. For a detailed review on suspensions of rigid particles, the reader is refereed to the recent perspective by Guazzelli & Pouliquen (Reference Guazzelli and Pouliquen2018).
When considering deformable particles, emulsions being one relevant example, the scenario is further complicated by the occurrence of deformation, coalescence and breakup. Zhou & Pozrikidis (Reference Zhou and Pozrikidis1993) simulated numerically two-dimensional emulsions employing a boundary element method (BEM), which later has been extended to three-dimensional flows by Loewenberg & Hinch (Reference Loewenberg and Hinch1996). Srivastava, Malipeddi & Sarkar (Reference Srivastava, Malipeddi and Sarkar2016) investigated the inertial effects on emulsions using a front tracking method. These authors reported an inversion of the normal stress difference for increasing Reynolds number. Matsunaga et al. (Reference Matsunaga, Imai, Yamaguchi and Ishikawa2015), also using a BEM method, investigated the rheology of a suspension of capsules up to 40 % of volume fraction and found that corrections of order
$O(\unicode[STIX]{x1D719}^{3})$
to the effective viscosity are negligible with respect to the rigid spheres. Rosti & Brandt (Reference Rosti and Brandt2018) and Rosti, Brandt & Mitra (Reference Rosti, Brandt and Mitra2018) simulated suspensions of deformable particles described as a viscous hyperelastic solid material and showed that the effective viscosity of the suspension can be estimated by the Eilers formula, valid for rigid particles, when computing an effective volume fraction based on the mean deformation of the particles. As mentioned before, all these studies do not account for coalescence or breakup, which however play key roles in the rheological behaviour of emulsions. The aim of this work is therefore to show how coalescence affects the rheology of emulsions.
The macroscopic properties of emulsions in shear flow strongly depend on their microstructure, mainly droplets size and distribution. The dynamics of the disperse phase is dictated by the interface deformation and the collision rate. In a dilute emulsion the interface deforms only under the action of the viscous stress exerted by the external flow. The relative importance between this and the interfacial stress, which tends to keep the drop spherical, is known as the capillary number (
$Ca$
). When
$Ca$
exceeds a critical value, droplets do not reach a steady state shape but break into two or more fragments in order to recover the balance between viscous and interfacial stresses. When the volume fraction increases, i.e. in the non-dilute regime, the flow can induce collisions between two or more drops. The outcome of a collision depends on the interaction force between the droplets. Chesters (Reference Chesters1991) proposed a conceptual framework to describe the complex dynamic of two colliding drops. The interaction can be thought of as the combined action of an external flow, responsible of the frequency, force and duration of the collision, and an internal flow which accounts for the deformation of the interface and the drainage of the thin liquid film between the two particles. When the collision duration is larger than the drainage time, droplets coalesce whereas in the opposite case they repel. In the first case the emulsion is attractive, in the latter case it is repulsive. Many experimental studies have been conducted to describe size evolution and deformation of droplets in microconfined shear flow (Guido & Simeone Reference Guido and Simeone1998; Sibillo et al.
Reference Sibillo, Pasquariello, Simeone, Cristini and Guido2006; Vananroye, Van Puyvelde & Moldenaers Reference Vananroye, Van Puyvelde and Moldenaers2006). Interestingly, Caserta & Guido (Reference Caserta and Guido2012) have shown that the rheological curve, effective viscosity versus volume fraction, exhibits negative curvature for emulsions, contrary to the behaviour of particles (indeed the term
$O(\unicode[STIX]{x1D719}^{2})$
in the Eilers formula is positive). In this work the authors reported also a phase separation for viscosity ratios smaller than one in large domains: they observed droplet rich and droplet poor regions oriented in the flow direction and alternated in the vorticity direction (vorticity banding). These results suggest that coalescence plays a key role in the rheological behaviour of emulsions and can affect their dynamics under shear flow.
Performing numerical simulations of emulsions in shear flow is a challenging problem due to the large scale separation between the external flow and the internal flow. The gap between the plates can be order of
$\unicode[STIX]{x03BC}\text{m}$
while the thin liquid film where the van der Waals force acts is order of
$nm$
(Chesters Reference Chesters1991). For this reason most of the numerical works on emulsions in the literature neglect coalescence. Only a few works that include coalescence are available: Shardt, Derksen & Mitra (Reference Shardt, Derksen and Mitra2013) investigated collision between two droplets in dilute emulsions using a lattice Boltzmann method; Rosti, De Vita & Brandt (Reference Rosti, De Vita and Brandt2019) present a numerical model for simulations of emulsions at moderate concentrations using a volume of fluid (VoF) approach and present results for different capillary numbers. Additionally, the presence of surfactants (surface active agents) can significantly affect the interface dynamics in terms of droplet deformability and distribution evolution. The modelling of this mechanisms in a numerical simulation is not trivial (Khatri & Tornberg Reference Khatri and Tornberg2011; Soligo, Roccon & Soldati Reference Soligo, Roccon and Soldati2019). An alternative to the previously mentioned methods is to employ mesoscale methods, which consider also thermal effects, as for example in Sega et al. (Reference Sega, Sbragaglia, Kantorovich and Ivanov2013). In this study, we complete the study of the rheology of emulsions by including the effect of the viscosity ratio and, more importantly, the role of coalescence/breakup on the macroscopic behaviour. The use of a VoF technique allows us to take into account the coalescence of droplets, while we model repulsive short-range interactions.
The paper is structured as follow: in § 2 we describe the governing equations, the numerical method and the set-up we adopt; in § 3 we present a fully Eulerian framework to introduce a collision force to the solver which is able to delay or prevent the coalescence between droplets; in § 4 we discuss the results in term of the effective viscosity, stress budget and normal stresses for emulsions with different coalescence rates. Finally we summarise the main findings and draw some conclusions in § 5.
2 Numerical method and set-up
In this study we investigate the rheology of emulsions in shear flow simulating a liquid–liquid biphasic system in a simple Couette flow. The flow dynamics is governed by the incompressible Navier–Stokes equations which in the one-fluid formulation for multiphase flows read
To track in time the position of the interface we employ a VoF technique based on the multi-dimensional tangent of hyperbola interface capturing (MTHINC) method (Ii et al.
Reference Ii, Sugiyama, Takeuchi, Takagi and Matsumoto2012). To identify the two fluids we define a colour function
${\mathcal{H}}(\boldsymbol{x},t)$
so that
${\mathcal{H}}=1$
in one fluid and
${\mathcal{H}}=0$
in the other. The VoF function
${\mathcal{T}}(\boldsymbol{x},t)$
is defined as the cell-average value of the colour function
and represents the volume fraction in each cell of the domain. The VoF function is advected by the flow field as
The key point of the MTHINC method is to approximate the colour function with a tangent of hyperbola
where
$(X,Y,Z)\in [0,1]$
is a local coordinate system,
$\unicode[STIX]{x1D6FD}$
a sharpness parameter,
$d$
a normalisation factor and
$P$
a three-dimensional surface function which can be either linear (plane) or quadratic (curved surface) with no additional cost. This discretisation allows us to solve the fluxes of (2.3) by integration of the approximated colour function in each computational cell. The material properties of the two fluids are linked to the VoF function
${\mathcal{T}}$
as follows:
where the subscript 1 stands for the disperse phase, the subscript 0 for the carrier fluid and
${\mathcal{T}}$
is equal to 1 in the disperse phase and 0 in the carrier fluid. Finally the surface tension force is approximated using the continuum surface force (CSF) approach (Brackbill, Kothe & Zemach Reference Brackbill, Kothe and Zemach1992)
See Rosti et al. (Reference Rosti, De Vita and Brandt2019) for a detailed description and validation of the code employed in this work.
Sketch of the computational domain and coordinate system used.

2.1 Flow configuration
A random monodisperse distribution of droplets, with radius
$r$
, is initialised in a bi-periodic system (
$x$
being the streamwise direction and
$y$
the spanwise) delimited by two walls (
$z$
direction) moving in opposite directions with velocity
$U$
at a constant distance
$L_{z}$
. Choosing as reference length the initial radius of the droplets,
$r$
, the size of the domain is
$16r\times 16r\times 10r$
(
$x$
,
$y$
and
$z$
, respectively), as sketched in figure 1. This configuration has been widely adopted in the literature for the study of the rheology of suspensions of rigid and deformable particles and emulsions (see Picano et al.
Reference Picano, Breugem, Mitra and Brandt2013; Rosti & Brandt Reference Rosti and Brandt2018; Rosti et al.
Reference Rosti, De Vita and Brandt2019). The flow is governed by four non-dimensional parameters, namely the Reynolds number
$Re$
, the capillary number
$Ca$
, the viscosity ratio
$\unicode[STIX]{x1D706}$
and the volume fraction
$\unicode[STIX]{x1D719}$
, with
The applied shear rate is equal to
$\dot{\unicode[STIX]{x1D6FE}}=2U/L_{z}$
and, in order to keep the inertial effect negligible, is chosen to give a Reynolds number equal to 0.1. The interface tension coefficient
$\unicode[STIX]{x1D70E}$
is chosen so that the capillary number based on the initial radius varies between 0.05 and 0.2. The simulations are performed with fixed Reynolds number, for different viscosity ratios,
$\unicode[STIX]{x1D706}=[0.01,0.1,1]$
and for four volume fractions,
$\unicode[STIX]{x1D719}=[0.00164,0.1,0.2,0.3]$
, the first one corresponding to one single droplet initially located at the centre of the domain. In all the cases the density of the two fluids is the same, thus there are no buoyancy effects. All the simulations have been performed with a mesh of
$256\times 256\times 160$
cells, which corresponds to a grid size of
$\unicode[STIX]{x1D6E5}=r/16$
. The independence of the results from the grid size, measured by the suspension effective viscosity and surface area evolution, was verified by performing simulations with half and double resolution.
3 Collision force model
The outcome of a collision between two droplets is function of hydrodynamic forces, capillary forces and geometric parameters, such as radius and relative distance. Droplets coalesce when the capillary number is smaller than a critical value
$Ca_{c}$
whereas they slide away when the capillary number is larger than
$Ca_{c}$
(Shardt et al.
Reference Shardt, Derksen and Mitra2013). To investigate the effect of coalescence on the rheology of emulsions we introduce a repulsive force with the aim of completely inhibiting the merging of droplets. From a physical point of view, this force can be seen as a model for the presence of surfactants in the emulsion, which alter the value of the interfacial tension coefficient and change the capillary number. The action of surfactants occurs on the nanoscale, which cannot be resolved on our mesh and therefore needs to be modelled. To prevent coalescence one should impose a very strong interfacial tension coefficient which would make the simulations expensive by requiring a very small time step. Because we are only interested in inhibiting coalescence, we introduce an Eulerian repulsive force of the same form as a lubrication force (Bolotnov et al.
Reference Bolotnov, Jansen, Drew, Oberai, Lahey and Podowski2011)
where
$a$
and
$b$
are two coefficients and
$\unicode[STIX]{x1D713}$
is the signed distance from the interface. Operationally, we reconstruct every time step the distance function from the interface (defined by the VoF function
${\mathcal{T}}$
) as follows: first we initialise the distance function
$\unicode[STIX]{x1D713}_{0}$
following Albadawi et al. (Reference Albadawi, Donoghue, Robinson, Murray and Delauré2013)
with
$\unicode[STIX]{x1D6E5}$
the grid spacing; then the distance function is propagated in a region of thickness
$3\unicode[STIX]{x1D6E5}$
from the interface solving the redistancing equation
where
$S(\unicode[STIX]{x1D713}_{0})$
is the sign function and
$\unicode[STIX]{x1D70F}$
is an artificial time; the equation is marched in time with step
$\unicode[STIX]{x0394}\unicode[STIX]{x1D70F}=0.1\unicode[STIX]{x1D6E5}$
. We solve the previous equation using the second-order redistancing algorithm of Russo & Smereka (Reference Russo and Smereka2000). At the beginning of the simulation, every droplet is tagged with an integer number
${\mathcal{I}}$
, going from 1 to the total number of droplets. This index is used to determine where to apply the collision force, as sketched in figure 2. The grey area in the figure represents the region, of thickness
$3\unicode[STIX]{x1D6E5}$
, where the distance function is propagated and the red region is the area where the force is applied. This region is found by checking, inside the grey area, if in a stencil
$7\unicode[STIX]{x1D6E5}\times 7\unicode[STIX]{x1D6E5}\times 7\unicode[STIX]{x1D6E5}$
there are two different values of
${\mathcal{I}}$
. If this condition is verified, the interaction force is computed following (3.1). Two additional indexes are used to identify the walls and model collisions with them, also to avoid mass losses at the boundary. Every time step, the droplet index
${\mathcal{I}}$
is updated based on the new position of
${\mathcal{T}}$
: if
${\mathcal{T}}<0.5$
we set
${\mathcal{I}}=0$
otherwise we search in a stencil
$3\unicode[STIX]{x1D6E5}\times 3\unicode[STIX]{x1D6E5}\times 3\unicode[STIX]{x1D6E5}$
the non-zero value of
${\mathcal{I}}$
. To avoid instabilities resulting from the application of the force directly on the surface, the force is applied only if
$\unicode[STIX]{x1D713}$
is larger than
$\sqrt{2}\unicode[STIX]{x1D6E5}/2$
in two dimensions and
$\sqrt{3}\unicode[STIX]{x1D6E5}/2$
in three dimensions. The algorithm currently implemented to advance in time the index function does not handle topological changes, i.e. breakup or coalescence, hence we need to ensure that none of the droplets coalesces or breaks during the simulations. By trial and error, we found that the minimum values of the coefficients in (3.1) that always prevent coalescence, for the given non-dimensional parameters of the flow, are
$a=55$
and
$b=3.5$
. These values have been used to obtain the results presented in the next section.
Sketch of the collision model: the grey area represents the region of propagation of the distance function
$\unicode[STIX]{x1D713}$
, the external contour corresponding to
$\unicode[STIX]{x1D713}=3\unicode[STIX]{x1D6E5}$
. The red region represents the area where the force is applied.

Time evolution of the interface at three different time instants; from left to right:
$t=10$
,
$t=15$
and
$t=20$
. Panels (a–c) zero force; panels (d–f)
$a=13.75$
and
$b=0.875$
; panels (g–i)
$a=27.5$
and
$b=1.75$
; panels (j–l)
$a=55$
and
$b=3.5$
.

Before moving to the next section, we show the effect of the repulsive force (3.1) on the interaction between two droplets in a simple shear flow and the sensitivity to the coefficients
$a$
and
$b$
. We consider two drops of equal initial radius
$r=1$
, the length of the domain is
$12r$
and the height is
$6r$
. The two droplets are placed in the centre of the domain with a horizontal shift
$\unicode[STIX]{x0394}x=3$
(half of the channel height) and with a vertical shift
$\unicode[STIX]{x0394}y=1.5$
. The particle Reynolds number
$Re=0.1$
, the capillary number
$Ca=0.1$
and the viscosity and density ratios equal to unity. Note that this set-up is similar to that adopted in other numerical and experimental studies (Guido & Simeone Reference Guido and Simeone1998; Shardt et al.
Reference Shardt, Derksen and Mitra2013). Figure 3 shows the interface location at three different time instants,
$t=10$
,
$t=15$
and
$t=20$
, from the left to the right. Each row of the figure, from top to bottom, corresponds to a different case: first row (panels a–c) corresponds to the case of zero force, when droplets are free to coalesce (this is the case with maximum drainage velocity); second row (panels d–f) corresponds to the case with coefficients
$a=13.75$
and
$b=0.875$
; third row (panels g–i) corresponds to the case with coefficients
$a=27.5$
and
$b=1.75$
; last row (panels j–l) corresponds to the case with coefficients
$a=55$
and
$b=3.5$
, the same as in the results section. As soon as the droplets come closer to each other, they start to flatten and to deform. As a consequence of the applied shear and of the interaction of the particles with the local flow field they start to rotate. For the adopted capillary number and initial displacement and in the absence of any collision force, the two drops coalesce at a certain instant,
$t\approx 16$
. If the applied force is small (second row) the time needed to actually merge is larger and the coalescence is delayed. This is an interesting result because by tuning the value of the force it is possible to change the time scale of the drainage of the film between the droplets. Finally, panels (j–l) show the case of a fully repulsive force, completely preventing coalescence.
Figure 4 shows the time history (from the left to the right in a curvilinear system) of the vertical displacement
$\unicode[STIX]{x0394}y$
over the horizontal displacement
$\unicode[STIX]{x0394}x$
for the four cases of figure 3. When the droplets coalesce both the horizontal and vertical displacements decrease quickly, as illustrated by the red and green curves in figure 4. In the opposite case, after the collision, droplets continue to move away in the horizontal direction with a steady state value of the vertical displacement. Although the force in the fourth case is twice as large as that in the third case, the final vertical displacement differs only by approximately 3 %, showing that the force does not affect significantly the dynamics of the droplets after the collision. To check the effect of the grid size we also performed one simulation without the collision force and with double resolution, reported with a black line in the same figure. The final value of the vertical displacement for the cases with the force is approximately 6 % different from the case with double resolution.
Plot of relative distance between droplets
$\unicode[STIX]{x0394}y$
versus
$\unicode[STIX]{x0394}x$
: case with no force (○); case with
$a=13.75$
and
$b=0.875$
(♢); case with
$a=27.5$
and
$b=1.75$
(▫); case with
$a=55$
and
$b=3.5$
(▵). Black dots (●) correspond to the time instant in figure 3. The black line corresponds to the simulation without a collision force and with double resolution.

4 Results
We start our analysis by first considering the cases without the repulsive force given by (3.1). Although the two fluids considered are both Newtonian, emulsions can exhibit non-Newtonian behaviour such as shear thinning, normal stress differences or viscoelasticity (Foudazi et al.
Reference Foudazi, Qavi, Masalova and Malkin2015). A first parameter used to characterise the rheological behaviour of a suspension is the effective shear viscosity of the system,
$\unicode[STIX]{x1D707}_{e}$
. The wall-normal gradient of the streamwise velocity gradient at the walls can be used to evaluate the effective viscosity of an emulsion as
$$\begin{eqnarray}\unicode[STIX]{x1D707}_{e}=\frac{\displaystyle \left\langle \unicode[STIX]{x1D707}_{0}\frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}z}\biggm\vert_{w}\right\rangle }{\dot{\unicode[STIX]{x1D6FE}}},\end{eqnarray}$$
where the symbol
$\langle \rangle$
indicates that the quantity has been averaged in the two homogeneous directions and in time. Following Yang, Krishnan & Shaqfeh (Reference Yang, Krishnan and Shaqfeh2016), it can be shown that the wall-normal velocity gradient at the wall is equivalent to the bulk shear stress. Here, we show how to include the interfacial tension contribution inside the shear stress balance. At steady state the streamwise momentum equation, after averaging in the homogeneous direction, reduces to
where
$f_{1}$
represents the streamwise component of the interfacial tension force. This equation can be rewritten as
where we define
${\mathcal{G}}(z)=\int _{z}\langle f_{1}\rangle \,\text{d}z$
the integral of
$\langle f_{1}\rangle$
in the wall-normal direction. To compute this integral we approximate the relation
$\unicode[STIX]{x2202}{\mathcal{G}}/\unicode[STIX]{x2202}z=\langle f_{1}\rangle$
with a finite difference scheme and invert it. To determine the constant of integration we impose
${\mathcal{G}}=0$
at the walls. We can now integrate equation (4.3) from the wall,
$z=0$
, to a generic section
$z=\unicode[STIX]{x1D701}$
The first term on the left-hand side is zero because the velocity is zero at the wall and the last term on the left-hand side is zero due to the imposed boundary condition, therefore
Because of the arbitrary choice of the section
$\unicode[STIX]{x1D701}$
, we can average the right-hand side in the wall-normal direction
or similarly
where the symbol
$\langle \langle \rangle \rangle$
indicates that the quantity has been averaged in the whole domain. Providing that the first term on the right-hand side is negligible, we have proved that the wall-normal shear stress is a measure of the bulk shear stress in the domain, given by the sum of the viscous part and the interfacial tension contribution. Before proceeding to quantify the shear stress components we need to underline two aspects: (i) the interfacial term
${\mathcal{G}}$
at the wall is zero only if there are no droplets reaching the wall and this is always verified in our simulations, as shown in the following section; (ii) in order to consider our solution at steady state we average in time over several shear units (typically between
$20$
and 40) and check if the mean value of the shear rate varies less than a few per cent, therefore all the quantities have to be considered also as averaged in time.
(a) Stress distribution along the wall-normal direction
$z$
(i.e. the direction of the velocity gradient) for the case with
$Ca=0.1$
,
$\unicode[STIX]{x1D719}=0.1$
and
$\unicode[STIX]{x1D706}=1.0$
: disperse phase viscous stress (▫); outer fluid viscous stress (○);
${\mathcal{G}}$
(▵); advection (♢). (b) Histogram of the stress components for the case with
$Ca=0.1$
,
$\unicode[STIX]{x1D719}=0.1$
and three viscosity ratios: interfacial force (solid blue); outer fluid viscous stress (green dense net); disperse phase viscous stress (red sparse net); the stresses are normalised with the wall shear stress.

To illustrate the outcome and indications of such an analysis, we display in figure 5(a) the profile of the stress components over the wall-normal direction
$z$
, that is the direction of the velocity gradient, obtained from (4.7). As already anticipated, the advection term is small (order
$10^{-6}$
) and for this reason will be neglected in the discussion hereinafter. The main contribution to the total shear stress is given by the viscous stress of the outer fluid, which is also the only non-zero term at the walls, and the sum of all the components is constant along the vertical direction. This confirms that the shear stress at the wall is equal to the bulk shear stress and that the effective viscosity can be evaluated using the wall-normal gradient of the streamwise velocity at the wall, see (4.1). In (b) of the same figure, we report the histogram of the shear stress components for three different values of the viscosity ratio. Here, we note that the contribution due to the viscous stress of the dispersed phase (red sparse net) is almost negligible for
$\unicode[STIX]{x1D706}$
lower than 1, while the interfacial tension (in blue) remains constant, as the volume fraction is constant.
Effective viscosity
$\unicode[STIX]{x1D707}_{e}$
as function of the volume fraction:
$\unicode[STIX]{x1D706}=1.0$
and different
$Ca$
(a);
$Ca=0.1$
and different
$\unicode[STIX]{x1D706}$
(b). The effective viscosity is normalised with the outer fluid viscosity.

We next examine in detail the global suspension viscosity for the emulsions with coalescence. The effective viscosity for different volume fractions, capillary numbers and viscosity ratios is shown in figure 6, together with the experiments in Caserta & Guido (Reference Caserta and Guido2012). The rheological curve exhibits a negative curvature, unlike rigid and deformable particles (Guazzelli & Pouliquen Reference Guazzelli and Pouliquen2018; Rosti & Brandt Reference Rosti and Brandt2018). In panel (a) of the figure, we display the results for emulsions with a viscosity ratio equal to one: the effective viscosity is always greater than the outer fluid viscosity and exhibits a maximum at approximately
$\unicode[STIX]{x1D719}=0.2$
. As
$Ca$
decreases, droplets become less deformable, leading to an increase of the effective viscosity, whereas, for the higher capillary case, the decrease of effective viscosity at high volume fraction is more prominent. When the disperse phase viscosity decreases (panel b of the same figure 6) the peak moves to lower volume fractions, approximately
$\unicode[STIX]{x1D719}=0.1$
; in addition, the effective viscosity of the emulsion becomes smaller than 1 at larger volume fraction, when increasing the volume of the low viscosity fluid. To check the effect of the droplet to system size ratio, we have also performed simulations with droplets of half the initial radius, shown in black in figure 6(b): also in this case the curvature of the effective viscosity is negative. The differences in the values of the effective viscosity are a consequence of the different capillary number in the case of smaller droplets with the same interfacial tension.
Caserta & Guido (Reference Caserta and Guido2012) performed experiments with emulsions in shear flow and reported a phase separation associated with the negative curvature of the curve of the effective viscosity versus the volume fraction. This separation results in bands that are aligned in the direction of the flow and alternated in the vorticity direction; the process of banding has been observed only for viscosity ratios smaller than 1. We can observe that in the case of
$\unicode[STIX]{x1D706}=1$
the concavity of the viscosity curve is only marginally negative, while for
$\unicode[STIX]{x1D706}<1$
the sign of the curvature is clearly negative both experimentally and numerically. For this reason the banding phenomenon was visible only in emulsions with
$\unicode[STIX]{x1D706}<1$
. Because the characteristic width of the bands is of the order of the gap between the walls we cannot reproduce with the present simulations this separation of phases. For this reason we compare results only for the case with
$\unicode[STIX]{x1D706}=1$
, when the distribution of the dispersed phase in homogeneous in space. The numerical results in figure 6 are qualitatively in good agreement with the experimental results, also reported in the same figure. It is difficult to perform an exact comparison between the numerical simulations and the experimental results due to the uncertainty on the initial size of the droplets. The initial size distribution in the experiment is a consequence of the application of a strong pre-shear, higher than the one corresponding to the critical capillary number for breakup, to fragment the droplets and remove any effect of the initial configuration. By knowing the value of the pre-shear
$(40~\text{s}^{-1})$
and assuming that the droplets have a monodisperse radius distribution corresponding to the critical capillary number of 0.5 (for
$\unicode[STIX]{x1D706}=1$
), we can estimate the initial value of the radius. With this estimated radius the
$Ca$
for the data in figure 6 corresponding to the experiments is 0.12. Therefore, the small mismatch between our initial
$Ca$
and the estimate of the experiment can be most likely explained by the coalescence efficiency, which is close to unity in our simulations.
It is worth mentioning that shear banding is a phase separation process observed only for viscosity ratios smaller than unity and in large domains, in the vorticity direction, and over a very long time, i.e. more than 1000 shear units. In this case the dispersed phase distribution is not homogenous and there is a reduction in the effective viscosity of the system (Caserta & Guido Reference Caserta and Guido2012). It is common in the literature to distinguish between the constitutive curve, defined as the effective viscosity as a function of the volume fraction for a homogeneous microstructure and shear rate, and the flow curve effective viscosity versus volume fraction measured in the rheometer which can include vorticity banding. In this work we force the dispersed phase to be approximatively homogenous due to lateral confinement, hence the curve reported here should be considered as the constitutive curve of the emulsion. This should be compared to experimental measures run at short time (see Caserta & Guido Reference Caserta and Guido2012) when the bands have not yet formed. The effect of the vorticity banding on the effective viscosity will be an object of future works.
4.1 Rheology of repulsive emulsions
In this section we study the emulsion behaviour when suppressing the coalescence. This is equivalent to assuming that the characteristic drainage time tends to infinity; hence, the collision dynamics is faster than the coalescence time scale and the droplets never merge. We have therefore performed simulations for the same parameters as in the results discussed in the previous section, now with the collision force given by (3.1). The effective shear viscosity obtained when varying the viscosity ratio and the volume fraction
$\unicode[STIX]{x1D719}$
is reported in figure 7.
When the coalescence is prohibited, the curve of viscosity versus concentration exhibits a positive curvature, as for the case of rigid and deformable particles. In the same figure, we display also the Eilers fit, valid for solid particles, and the fit from Rosti & Brandt (Reference Rosti and Brandt2018) for deformable particles. In particular, Rosti et al. (Reference Rosti, Brandt and Mitra2018) have shown that it is possible to estimate the effective viscosity of a suspension of deformable particles with the Eilers formula by computing an effective volume fraction based on the mean deformation of the particles (see appendix A for more details). The results in figure 7 demonstrate that, in the absence of coalescence, emulsions behave as suspensions of deformable particles.
The change of sign in the curvature can be explained by examining the stress components, as shown previously for droplets coalescing. In this case an additional force, the collision force, needs to be included in the stress budget, which is treated in the same way as the interfacial tension force,
with
${\mathcal{C}}=\int _{z}\langle F_{c}\rangle \,\text{d}z$
the wall-normal integral of the streamwise component of the collision force.
Effective viscosity
$\unicode[STIX]{x1D707}_{e}$
for the cases with collision force and
$Ca=0.1$
. The line denoted as ‘deformable particles’ corresponds to the fit proposed in Rosti & Brandt (Reference Rosti and Brandt2018) for initially spherical viscoelastic particles with the same capillary number based on the material shear elastic module.

(a) Wall-normal distribution (
$z$
) of the stresses for the case with
$Ca=0.1$
,
$\unicode[STIX]{x1D719}=0.3$
,
$\unicode[STIX]{x1D706}=1.0$
and a collision force: disperse phase viscous stress (▫); carrier fluid viscous stress (○); interfacial force (▵); collision force (♢). The stresses are normalised with the wall shear stress. (b) Histogram of the stress components for the case with
$Ca=0.1$
,
$\unicode[STIX]{x1D719}=0.3$
with (left column) and without (right column) the collision force: interfacial force (solid blue); carrier fluid viscous stress (green dense net); disperse phase viscous stress (red sparse net); collision force (purple oblique bars). The stresses are normalised with the single-phase wall shear stress.

(a) Time history of the total surface for different volume fractions and
$\unicode[STIX]{x1D706}=0.01$
. Cases without collision force:
$\unicode[STIX]{x1D719}=0.1$
(▫),
$\unicode[STIX]{x1D719}=0.2$
(○),
$\unicode[STIX]{x1D719}=0.3$
(▵); cases with collision force: same colour and solid symbols. (b) Comparison of the wall-normal distribution of the average volume fraction
$\langle \unicode[STIX]{x1D719}\rangle$
for the cases with
$\unicode[STIX]{x1D719}=0.3$
and
$Ca=0.1$
: without collision force (▫); with collision force (○).

The stress budget for the case with
$Ca=0.1$
,
$\unicode[STIX]{x1D719}=0.3$
and
$\unicode[STIX]{x1D706}=1$
is shown in figure 8(a). The interfacial tension contribution is now of the same order of magnitude as the outer fluid viscous stress, whereas in the absence of a collision force it is approximately two times smaller, as shown in (b). By comparing the stress budget for the cases with the collision force and without the collision force (figure 8
b), we notice that the increase in effective viscosity is mostly due to the interfacial tension term, with a small contribution due to the collision force, which is approximately 10 % of the total stress. Without the collision force, coalescence decreases the total interface area, which leads to a reduction of the energy associated with the surface tension. This is clearly demonstrated by computing the time history of the total surface. In figure 9(a) the total surface, normalised with the initial value, is displayed for different volume fractions at the viscosity ratio
$\unicode[STIX]{x1D706}=0.01$
for the simulations with and without the collision force. The latter cases exhibit a strong decrease of the total surface before reaching a regime configuration, associated with a reduction of the interfacial area of up to 80 % of the initial value for the highest volume fraction considered. This statistically steady state is reached faster for the cases with higher volume fractions. On the contrary, when coalescence is prohibited the total surface area slightly increases, owing to the deformation of the individual droplets. This explains the increased surface tension contribution to the stress budget observed in the absence of coalescence.
Figure 9(b) depicts the wall-normal distribution of the average (in the homogeneous
$x$
and
$y$
directions) local volume fraction
$\langle \unicode[STIX]{x1D719}\rangle$
for the cases with
$\unicode[STIX]{x1D719}=0.3$
and
$Ca=0.1$
. Here, we note a significant increase of the volume fraction close to the walls for the cases with the collision force, which can be associated with the increased effective viscosity. On the other hand, in the cases with coalescence, the larger droplets are localised more towards the centre of the domain. This result is in agreement with previous experimental observations (Hudson Reference Hudson2003; Caserta et al.
Reference Caserta, Sabetta, Simeone and Guido2005) where the migration was attributed to the combined action of wall migration and shear-induced diffusion of drops. This is also confirmed in figure 10 where we report a snapshot of the droplet distribution for the cases in figure 9(b).
Instantaneous droplet distribution for emulsions with volume fraction
$\unicode[STIX]{x1D719}=0.3$
and capillary number
$Ca=0.1$
extracted from the simulation without a collision force (a) and with a collision force (b).

A measure of the coalescence is the droplet size distribution in the system. To study this, we evaluate the volume distribution in the whole domain at steady state and draw the cumulative distribution in figure 11 for the cases without the collision force. On the horizontal axis we report the value of the equivalent diameter computed with the volume of each dispersed droplet, normalised with the initial value, and on the vertical axis the fraction of the total volume occupied by droplets of size smaller and equal to the corresponding abscissa. At volume fraction
$\unicode[STIX]{x1D719}=0.1$
the maximum value of the diameter is approximately 3.5 the initial one and, by reducing the viscosity ratio, the distribution shows a reduction of droplets with small diameter due to the increase of coalescence. By increasing the volume fraction droplets merge into larger structures which have a maximum equivalent diameter of approximately five times the initial one. In the simulations with the collision force, the mean spherical equivalent diameter is the same for all the droplets and equal to the initial value, depicted by the black line in figure 11.
Cumulative volume distribution as a function of the equivalent diameter, normalised with the initial value, for the case with
$\unicode[STIX]{x1D719}=0.1$
and different
$\unicode[STIX]{x1D706}$
(a) and for the cases with different volume fractions and
$\unicode[STIX]{x1D706}=1$
. The black line represents the initial distribution, which is a Heaviside on the initial value of the diameter.

PDF of the flow topology parameter for simulations without a repulsive force among the droplets (a) and with (b): solid lines represent the outer fluid domain and dashed lines the dispersed phase. Colours refer to different volume fractions:
$\unicode[STIX]{x1D719}=0.1$
red,
$\unicode[STIX]{x1D719}=0.2$
blue,
$\unicode[STIX]{x1D719}=0.3$
black.

4.2 Flow topology and normal stress
Normal stresses can arise in emulsions when sheared and provide evidence for the non-Newtonian behaviour of the system (Loewenberg & Hinch Reference Loewenberg and Hinch1996; Pal Reference Pal2011; Srivastava et al. Reference Srivastava, Malipeddi and Sarkar2016). To evaluate the effect of the coalescence and of the volume fraction on the extensibility of the flow we compute the flow topology parameter (De Vita et al. Reference De Vita, Rosti, Izbassarov, Du, Tammisola, Hormozi and Brandt2018; Rosti et al. Reference Rosti, De Vita and Brandt2019) defined as
where
$D^{2}=\unicode[STIX]{x1D60B}_{ij}\unicode[STIX]{x1D60B}_{ij}$
and
$\unicode[STIX]{x1D6FA}^{2}=\unicode[STIX]{x1D6FA}_{ij}\unicode[STIX]{x1D6FA}_{ij}$
with
$\unicode[STIX]{x1D734}=(\unicode[STIX]{x1D735}\boldsymbol{u}^{T}-\unicode[STIX]{x1D735}\boldsymbol{u})/2$
the rate of rotation tensor. When
$Q=-1$
the flow is purely rotational, whereas regions with
$Q=0$
represent pure shear flow and those with
$Q=1$
elongational flow. In figure 12 we display the probability distribution function (PDF) of
$Q$
for different volume fractions and for simulations with and without the collision force. We observe that in the presence of coalescence the outer fluid has a peak of the PDF in correspondence with
$Q=0$
, which implies that the flow is almost pure shear flow as in the case of single-phase Couette flow. The dispersed phase, instead, shows a peak for negative values of
$Q$
(rotational flow), which moves towards zero on increasing the volume fraction. When the collision force is applied, we note a strong reduction of the peaks in the outer fluid distribution, compensated by an increase of the underlying area of the curve for positive values of
$Q$
. This means that, when the coalescence is prohibited, we find in the external flow regions not only shear flow but also extensional flow, which is related to the velocities in the gaps between the droplets. For these cases we also observe that the dynamics of the dispersed phase flow is substantially unaffected by the volume fraction and has always components of rotational flow (inside the deforming droplets).
We conclude our analysis by computing the normal stress difference following Batchelor’s formulation (Batchelor Reference Batchelor1970; Srivastava et al.
Reference Srivastava, Malipeddi and Sarkar2016). Here, however, we cannot follow the same methodology used for the shear stress, because it is not possible to obtain the integration constant needed to determine the function
${\mathcal{G}}$
when integrating in the homogeneous directions. Following Batchelor’s formulation, the bulk stress
$\langle \langle \unicode[STIX]{x1D749}\rangle \rangle$
can be decomposed into the sum of the stress due to the outer fluid, in the case of zero volume fraction,
$\langle \langle \unicode[STIX]{x1D749}^{0}\rangle \rangle$
and that arising from the presence of the disperse phase
$\langle \langle \unicode[STIX]{x1D749}^{1}\rangle \rangle$
The last term on the right-hand side can be further decomposed as the sum of three terms: the first due to the viscosity difference
$\langle \langle \unicode[STIX]{x1D749}^{\unicode[STIX]{x1D707}}\rangle \rangle$
, the second to the interfacial tension
$\langle \langle \unicode[STIX]{x1D749}^{\unicode[STIX]{x1D70E}}\rangle \rangle$
and the third to the perturbation in the velocity field induced by the presence of the droplets
$\langle \langle \unicode[STIX]{x1D749}^{ptb}\rangle \rangle$
$$\begin{eqnarray}\displaystyle \langle \langle \unicode[STIX]{x1D749}^{1}\rangle \rangle & = & \displaystyle \langle \langle \unicode[STIX]{x1D749}^{\unicode[STIX]{x1D707}}\rangle \rangle +\langle \langle \unicode[STIX]{x1D749}^{\unicode[STIX]{x1D70E}}\rangle \rangle +\langle \langle \unicode[STIX]{x1D749}^{ptb}\rangle \rangle \nonumber\\ \displaystyle & = & \displaystyle \frac{\unicode[STIX]{x1D707}_{1}-\unicode[STIX]{x1D707}_{0}}{{\mathcal{V}}}\int _{{\mathcal{S}}}(\boldsymbol{u}\boldsymbol{n}+\boldsymbol{n}\boldsymbol{u})\,\text{d}{\mathcal{S}}-\frac{\unicode[STIX]{x1D70E}}{{\mathcal{V}}}\int _{{\mathcal{S}}}\left(\boldsymbol{n}\boldsymbol{n}-\frac{\boldsymbol{I}}{3}\right)\,\text{d}{\mathcal{S}}-\frac{1}{{\mathcal{V}}}\int _{{\mathcal{V}}}(\unicode[STIX]{x1D70C}\boldsymbol{u}^{\prime }\boldsymbol{u}^{\prime })\,\text{d}{\mathcal{V}},\qquad\end{eqnarray}$$
where
${\mathcal{S}}$
is the total surface of the droplets and
${\mathcal{V}}$
is the averaging volume, i.e. the computational domain. The normalised normal stress differences are then defined as
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle N_{1}=\frac{\langle \langle \unicode[STIX]{x1D70F}_{xx}^{1}\rangle \rangle -\langle \langle \unicode[STIX]{x1D70F}_{zz}^{1}\rangle \rangle }{\unicode[STIX]{x1D707}^{0}\dot{\unicode[STIX]{x1D6FE}}}\\ \displaystyle N_{2}=\frac{\langle \langle \unicode[STIX]{x1D70F}_{zz}^{1}\rangle \rangle -\langle \langle \unicode[STIX]{x1D70F}_{yy}^{1}\rangle \rangle }{\unicode[STIX]{x1D707}^{0}\dot{\unicode[STIX]{x1D6FE}}}.\end{array}\right\}\end{eqnarray}$$
When
$N_{1}$
is greater then zero and
$\unicode[STIX]{x1D70F}_{xx}^{1}>\unicode[STIX]{x1D70F}_{zz}^{1}$
droplets are elongated in the direction of the flow and compressed in the wall-normal direction. The opposite configuration corresponds to negative values of
$N_{1}$
.
In figure 13, we report the values of
$N_{1}$
and
$N_{2}$
computed from our simulations. First, we observe that the magnitude of the normal stress difference is larger for the cases without the collision force. This can be explained by noting that, when the collision force is turned on, the droplets exhibit smaller deformation due to the higher packing inside the domain (see figure 10). The magnitude of the normal stress differences decreases with the viscosity ratio and increases with the volume fraction. The data in the figure also reveal that the second normal difference is always negative whereas the first normal difference is positive for
$\unicode[STIX]{x1D706}=1$
and
$\unicode[STIX]{x1D706}=0.1$
and becomes negative for smaller viscosity ratios. Additionally, the ratio of
$N_{1}$
and
$N_{2}$
is almost constant with the volume fraction, showing that they have the same dependency on
$\unicode[STIX]{x1D719}$
.
Normal stress differences for the cases with a repulsive force among the droplets (open symbols and dashed lines) and without (solid symbol and solid lines). Colours represent the different viscosity ratios under investigation:
$\unicode[STIX]{x1D706}=1$
red,
$\unicode[STIX]{x1D706}=0.1$
blue,
$\unicode[STIX]{x1D706}=0.01$
black.

5 Conclusion
We have studied the rheology of emulsions under shear flow in the dilute and moderate concentration regimes by means of numerical simulations. The multiphase flow dynamics is governed by the Navier–Stokes equations while the interface is tracked in time employing a VoF technique. This approach naturally allows us to take into account the coalescence of droplets, which significantly affects the rheological behaviour of emulsions. To single out the effect of coalescence, we have therefore developed an Eulerian collision model which allows us to delay or fully prohibit the merging of droplets. We have reported here simulations at different capillary numbers, viscosity ratios and volume fractions with a constant Reynolds number, small enough to assume inertial effects to be negligible. In this work we have focused on the characterisation of the constitutive curves of emulsions by simulating approximatively homogeneous distributions of the dispersed phase and neglecting the vorticity banding observed in experiments at viscosity ratios smaller than unity.
We show that the curvature of the rheological curve (effective viscosity versus volume fraction) is negative when coalescence is allowed whereas it changes sign when we introduce the collision force, which prevents merging. The decrease of the effective viscosity in the former case is a consequence of the reduction of the total surface of the system, which in turn reduces the contribution to the stress tensor due to the interface tension stress. When the coalescence is prohibited, this term is responsible for approximately half of the total effective viscosity and emulsions behave similarly to suspensions of deformable particles.
In the case of coalescence, we observe the formation of relatively large droplets, migrating towards the channel centre. Moreover, by examining the probability distribution function of the flow topology parameter, we show that the external flow is mostly shear while the dispersed phase exhibits some rotational flow regions, which become smaller when increasing the volume fraction (on average large droplets). On the other hand, when we introduce the collision force, the dominant presence of shear flow in the outer fluid decreases and regions of extensional flows emerge. The dispersed phase, in this case characterised by smaller deformed droplets in the flow, is mostly unaffected by the volume fraction and shares equally distributed regions of shear and rotational flow. Furthermore, analysis of the droplet size distribution at steady state reveals that in the presence of coalescence the mean equivalent diameter of the droplets in the emulsions increases by up to three times the initial value and is a function of both the volume fraction and viscosity ratio. For simulations with the collision force, instead, every droplet has an equivalent diameter equal to the initial one.
To characterise the viscoelastic system behaviour, we show that the first normal stress difference is positive and the second negative, as in suspensions of capsules and deformable particles. The magnitude of the normal stress difference always increases with the volume fraction; also noteworthy is that we find an inversion of the first normal difference for small values of the viscosity ratio.
To conclude, it is worth noting once more that we have investigated two limiting cases, one with coalescence efficiency tending to unity and one with efficiency approaching zero. In a real scenario, the coalescence efficiency is likely to have intermediate values and the behaviour of emulsions can therefore differ from the limiting cases considered here. Future works may therefore be devoted to handling both coalescence and collisions in order to simulate different types of emulsions and to study the formation of the vorticity banding.
Acknowledgements
The work is supported by the Microflusa project, the European Research Council grant no. ERC-2013-CoG-616186, TRITOS. The Microflusa project receives funding from the European Union Horizon 2020 research and innovation program under grant agreement no. 664823. The authors acknowledge computer time provided by SNIC (Swedish National Infrastructure for Computing).
Appendix A
In this appendix, we compare our results with expressions available in the literature for suspensions of particles and emulsions. Pal (Reference Pal2003) proposed the following equation for the viscosity of a concentrated emulsion
with
$\unicode[STIX]{x1D707}_{r}=\unicode[STIX]{x1D707}_{e}/\unicode[STIX]{x1D707}_{0}$
,
$\unicode[STIX]{x1D719}_{m}$
the maximum packing volume fraction and
Effective viscosity versus volume fraction. (a) Simulations with
$\unicode[STIX]{x1D706}=1$
and different capillary numbers: solid lines represent equation (A 1), dashed lines with symbols represent our numerical simulations. Colours represent different capillary numbers:
$0.05$
red,
$0.1$
green,
$0.2$
blue. Black diamonds represent experiment by Caserta & Guido (Reference Caserta and Guido2012). (b) Cases with
$Ca=0.1$
and different
$\unicode[STIX]{x1D706}$
: solid coloured lines represent equation (A 1), dashed lines with symbols represent our numerical simulations, black solid line represents data from Rosti et al. (Reference Rosti, Brandt and Mitra2018) and black dashed line represents equation (A 3). Colours represent different
$\unicode[STIX]{x1D706}$
:
$1$
red,
$0.1$
green,
$0.01$
blue.

Next, see figure 14(b), we compare (A 1) with results for the cases with the collision force, Eilers formula
$$\begin{eqnarray}\unicode[STIX]{x1D707}_{r}=\left[1+\frac{1.25\unicode[STIX]{x1D719}}{\displaystyle 1-\frac{\unicode[STIX]{x1D719}}{\unicode[STIX]{x1D719}_{m}}}\right]^{2}\end{eqnarray}$$
and the data from Rosti & Brandt (Reference Rosti and Brandt2018) for deformable particles of a viscous hyperelastic material. Despite the differences in the specific values, the trend is similar for all curves, with (A 1) underpredicting the effective viscosity. It is worth noticing that, if the parameter
$\unicode[STIX]{x1D719}_{m}$
is used as a fitting parameter, it is possible to get a better agreement between our results and equation (A 1), in particular using
$\unicode[STIX]{x1D719}_{m}\approx 0.4$
.
Appendix B
Rosti & Brandt (Reference Rosti and Brandt2018) showed that the effective viscosity of a suspension of deformable particles can be predicted using the Eilers formula (A 3) by modifying the volume fraction to take into account the particle deformation. This can be measured by the Taylor parameter
where
$a$
and
$b$
are the semi-major and semi-minor axes of the inscribed ellipse. The authors showed that this fit provides good predictions for different capillary numbers and viscosity ratios, both for fluid-filled capsules and red blood cells. Here, we wish to show that the same scaling can be applied also for the emulsions considered in this study when coalescence is inhibited. We therefore evaluate the effective viscosity
$\unicode[STIX]{x1D719}_{e}$
based on spheres of radius equal to the semi-minor axis
$b$
of the inscribed ellipsoid
where
$N$
is the number of droplets in the computational box of volume
${\mathcal{V}}$
. We then depict in figure 15 the effective viscosity of the emulsions as a function of the average volume fraction
$\unicode[STIX]{x1D719}$
and of the effective volume fraction
$\unicode[STIX]{x1D719}_{e}$
. The figure shows that, indeed, the effective viscosity of the emulsion successfully collapses onto the Eilers formula, with an error of approximately 6 % only for the highest volume fraction. It is worth noticing that this scaling works only for the cases with the collision force; indeed, in the case of coalescence, the shape of the interface is not well approximated by an ellipsoid, hence the Taylor parameter is not a good measure of the droplet deformation. Additionally, the curvature of the curve representing the effective viscosity versus the volume fraction is negative and thus it is not possible to collapse the results on the Eilers formula, which is instead monotonic. These results further prove that, in the absence of coalesce, emulsions behave similarly to suspensions of deformable particles and capsules and that coalescence should be considered as properly describing the rheological response of emulsions.
Effective viscosity as a function of the volume fraction
$\unicode[STIX]{x1D719}$
(▫) and of the effective volume fraction
$\unicode[STIX]{x1D719}_{e}$
(○). The black solid line represents the Eilers formula (A 3).















