Numerical simulations of emulsions in shear flows

We present a modification of a recently developed volume of fluid method for multiphase problems, so that it can be used in conjunction with a fractional step-method and fast Poisson solver, and validate it with standard benchmark problems. We then consider emulsions of two-fluid systems and study their rheology in a plane Couette flow in the limit of vanishing inertia. We examine the dependency of the effective viscosity on the volume-fraction (from 10% to 30%) and the Capillary number (from 0.1 to 0.4) for the case of density and viscosity ratio 1. We show that the effective viscosity decreases with the deformation and the applied shear (shear-thinning) while exhibits a non-monotonic behavior with respect to the volume fraction. We report the appearance of a maximum in the effective viscosity curve and compare the results with those of suspensions of rigid and deformable particles and capsules. We show that the flow in the solvent is mostly a shear flow, while it is mostly rotational in the suspended phase; moreover this behavior tends to reverse as the volume fraction increases. Finally, we evaluate the contributions to the total shear stress of the viscous stresses in the two fluids and of the interfacial force between them.


Introduction
In the last decades, developments in colloidal science have proven to be crucial for fabrication of functional materials. Particles or droplets, with typical scales of micron, are manipulated to self-organize into controlled patterns with high precision, which then form basic building blocks for more complex structures (Sacanna and Pine, 2011). One important aspect to consider during the synthesis and assembly of innovative materials is the rheological behavior of the system (Mewis and Wagner, 2012). Indeed, the material properties will depend on the distribution of the dispersed phase and thus a more accurate control on the production process can help to generate materials of desired properties (Xia et al, 2000).
A great amount of work has been done in the past to study rigid and deformable particle suspensions (Freund, 2014;Takeishi et al, 2016;Alizad Banaei et al, 2017;Rosti et al, 2018;Rosti and Brandt, 2018). In his pioneering work, Einstein (1956) showed that, in the limit of vanishing inertia and for dilute rigid particle suspensions, the viscosity is a linear function of the particle volume fraction. Although Batchelor and Green (1972) and Batchelor (1977) added a second order correction, all existing analytical relations are not valid for moderately high concentrations and one needs to resort to empirical fits. One of the available empirical relations that provides a good description of the rheology at zero Reynolds number both in the high and low concentration limits is the Eilers fit (Ferrini et al, 1979;Zarraga et al, 2000;Singh and Nott, 2003;Kulkarni and Morris, 2008). Recently, inertia and deformation have been shown to introduce deviations from the behavior predicted by the different empirical fits, and these effect can be related to an increase and decrease of an effective volume fraction (Picano et al, 2013;Rosti et al, 2018).
Less attention has been given to emulsions, which are instead the object of the present work. Emulsions are biphasic liquid-liquid systems in which the phases are separated by a deformable interface subject to interfacial surface tension. These can be found in a variety of applications, ranging from advanced materials processing, waste treatment, enhanced oil recovery, food processing, and pharmaceutical manufacturing. Similarly to particle suspensions, it is often desirable to predict or manipulate the rheology and microstructure of emulsions, which in general exhibit highly varied rheological behaviors (Mason, 1999); however, there has been limited progress towards the creation of theoretical models that can reliably predict the rheology and microstructure of such flows (Loewenberg and Hinch, 1996;Loewenberg, 1998) and for years measurements of emulsion rheology were not quantitatively understood because of the complexity of the phenomenology involved and the difficulty to properly choose materials with controlled properties (Mason, 1999). Thus, numerical simulations can help to fill this gap, and indeed there has been considerable progress in the development of numerical simulations able to study such multiphase flows (Prosperetti and Tryggvason, 2009;Tryggvason et al, 2011).
Different techniques have been proposed to numerically tackle the problem at hand. The so-called front-tracking method is an Eulerian/Lagrangian method, used to simulate viscous, incompressible, immiscible two-fluid systems, first developed by Unverdi and Tryggvason (1992) and Tryggvason et al (2001). When dealing with moving and deformable boundaries, an alternative approach are the so-called front-capturing methods, which are fully Eulerian and handle topology changes automatically. A strong advantage of these methods is that they are easier to parallelize than their Lagrangian coun-terpart. Eulerian interface representations include essentially the volume of fluid (VOF) (Scardovelli and Zaleski, 1999) and level-set (LS) (Sussman et al, 1994;Sethian, 1999;Sethian and Smereka, 2003) methods. The volume of fluid method defines different fluids with a discontinuous color function, and its main advantage is the intrinsic mass conservation; however, it suffers from an inaccurate computation of the interface properties, such as normals and curvatures (Francois et al, 2006;Cummins et al, 2005). Differently from the volume of fluid, the level-set method prescribes the interface through a continuous function which usually takes the form of the signed distance to the interface. Thus, normals and curvatures can be readily and accurately computed, while mass loss or gain may occur. In this work we will employ the volume of fluid method.
In a conventional VOF method, the interface separating different fluids is piece-wisely reconstructed in each numerical cell by straight line segments, which are then used to calculate the numerical fluxes necessary to update the local volume of fluid function. This geometric reconstruction effectively eliminates the numerical diffusion that smears out the compactness of the transition layer of the interface. Different methodologies have been proposed to accurately recover the exact surface geometry from the discretized VOF function: the simple line interface calculation (SLIC) method (Noh and Woodward, 1976), the piecewise linear interface calculation (PLIC) (Youngs, 1982(Youngs, , 1984, the latter being further modified by several authors (Puckett et al, 1997;Rider and Kothe, 1998;Harvie and Fletcher, 2000;Aulisa et al, 2003;Pilliod Jr and Puckett, 2004). Another technique is the tangent of hyperbola for interface capturing (THINC) method (Xiao et al, 2005), which avoids the explicit geometric reconstruction by using a continuous sigmoid function rather than the Heaviside function, thus allowing a completely algebraic description of the interface; this enables the computation of the numerical flux partially analytically. An improvement was proposed by combining the original THINC method with the first-order upwind scheme in the so-called THINC/WLIC (THINC/weighted linear interface capturing) method (Yokoi, 2007). Recently, the method has been further developed in the multi-dimensional THINC (MTHINC) method where the fully multi-dimensional hyperbolic tangent function is used to reconstruct the interface (Ii et al, 2012). The numerical fluxes can be directly evaluated by integrating the hyperbolic tangent function and the normal vector, curvature and approximate delta function can be directly obtained from the derivatives of the function. Moreover, the scheme does not require the geometric reconstruction and a curved (quadratic) surface can be easily constructed as well.

Outline
In this work, we first present our numerical solver for multiphase incompressible flows and then employ it to study liquid-liquid systems (emulsions) in a plane Couette flow at low Reynolds number. The two fluids are Newtonian Fig. 1 Sketch of the channel geometry and coordinate system adopted in this study. and satisfy the full incompressible Navier-Stokes equations. We compare our results with those of suspensions of rigid and deformable particles and with capsules, consisting of a second fluid enclosed by a thin elastic membrane. In section 2, we first discuss the flow configuration and governing equations, and then present the numerical methodology used and its validation. The rheological study of the emulsions is presented in section 3, where we also discuss the role of the different non-dimensional parameters governing the flow. Finally, a summary of the main findings and some conclusions are drawn in section 4.

Formulation
We consider the flow of two incompressible viscous fluids, separated by an interface, in a channel with moving walls, i.e., in a plane Couette geometry. Figure 1 shows a sketch of the geometry and the Cartesian coordinate system, where x, y and z (x 1 , x 2 , and x 3 ) denote the streamwise, wall-normal and spanwise coordinates, while u, v and w (u 1 , u 2 , and u 3 ) denote the corresponding components of the velocity vector field. The lower and upper impermeable moving walls are located at y = −h and y = h, respectively, and move in opposite direction with constant streamwise velocity ±V w .
The two fluid motion is governed by the conservation of momentum and the incompressibility constraint, and the kinematic and dynamic interactions between the two fluid phases are determined by enforcing the continuity of the velocity and traction force at the interface between the two phases, i.e., where the suffixes f 1 and f 2 are used to indicate the two phases, σ ij denotes the Cauchy stress tensor, n i the normal vector at the interface, κ the interface curvature and σ the surface tension (assumed here to be constant).
To numerically solve the two-phase interaction problem at hand, we use the volume of fluid method following Ii et al (2012). We introduce an indicator (or color) function H to identify each fluid phase so that H = 1 in the region occupied by the fluid f 1 and H = 0 otherwise. Considering that the fluid is transported by the flow velocity, we update H in the Eulerian framework by Numerical simulations of emulsions in shear flows 5 the following advection equation written in divergence form: where u i is the local fluid velocity and φ the cell-averaged value of the indicator function.
Once φ is known, the two-fluid equations can be rewritten in the so called one-continuum formulation (Tryggvason et al, 2007), so that only one set of equations is solved over the whole domain. This is achieved by introducing a monolithic velocity vector field u i , defined everywhere and found by applying the volume averaging procedure (Takeuchi et al, 2010;Quintard and Whitaker, 1994). Thus, u i is governed by the following set of equations where ρ is the density, f i the surface tension force defined as f i = σκn i δ, being δ the delta function at the interface, and σ ij is the stress written in a mixture form, i.e., Note that, we have chosen φ to be the volume fraction of fluid 2, i.e., this is zero in the fluid 1, whereas φ = 1 in the fluid 2, and 0 ≤ φ ≤ 1 close to the interface. Both fluids are assumed to be Newtonian so that their stress tensors can be written as σ ij = −pδ ij + 2µD ij , where p is the pressure, δ ij the Kronecker delta, µ the dynamic viscosity and D ij the strain rate tensor (defined as D ij = (∂u i /∂x j + ∂u j /∂x i )/2). Finally, the mixture density ρ and dynamic viscosity µ are simply averaged in terms of the local φ: Note that, in order to solve equation (3) and equation (2), we need to determine the indicator function H, the normal vector n i and the curvature κ.

The MTHINC method
The indicator function H can be reconstructed in various ways; here, we use the multidimensional tangent of hyperbola for interface capturing (MTHINC) method, developed by Ii et al (2012), where a multidimensional hyperbolic tangent function is used as an approximated indicator function. In particular, the indicator function H is approximated as 6 Marco E. Rosti et al. where X, Y, Z ∈ [0, 1] is a centered local coordinate system defined in each cell, P is a three dimensional surface function, β a sharpness parameter and d a normalization parameter. The function P can be either a linear function (a plane) P (X, Y, Z) = a 100 X + a 010 Y + a 001 Z, or a quadratic function (a curved surface) The coefficients a l,m,n are determined algebraically by imposing the correct value of the three normal components n i and the six components of the Cartesian curvature tensor l ij = (∂n i /∂x j + ∂n j /∂x i ) /2 for the function P in each cell. Finally, the parameter d is found by enforcing the following constraint: The integration can be performed analytically in one direction, and numerically in the other two directions by the two-point Gaussian quadrature. The unit normal vector is defined as n i = m i /|∇φ|, being m i the gradient of the volume of fluid function, i.e., m i = ∂φ/∂x i . Here, we compute m i using the usual Youngs approach (Youngs, 1982(Youngs, , 1984, where first the values of the derivative at the cell corners are calculated, and then averaged to find the cellcenter value. Once the normal vector is known, the curvature κ can be easily found by taking the divergence of the normal vector, i.e., κ = −∂n i /∂x i , and the surface tension force f i can be computed by the continuum surface force (CSF) model (Brackbill et al, 1992), where the 1D approximate delta function δ is directly approximated by δ ≈ |∇φ|. Thus, we obtain

Numerical discretisation
The equation of motion are solved with an extensively validated in-house code (Picano et al, 2015;Rosti and Brandt, 2017;Rosti et al, 2018;Rosti and Brandt, 2018). The equations are solved on a staggered uniform grid with velocities located on the cell faces and all the other variables (pressure, stress and volume of fluid) at the cell centers. All the spatial derivatives are approximated with second-order centered finite differences, while the time integration is discussed hereafter. First, the volume of fluid function is updated in time from the time-step (n) to (n + 1) by solving equation (2), following the procedure proposed by Ii et al (2012). In particular, the time evolution of φ is calculated by evaluating the numerical fluxes sequentially in each direction, a robust and easy approach called directional splitting. Thus, equation (2) is discretised sequentially in the three Cartesian direction, i.e., where the subscript in parenthesis indicates the time iteration, with (n) and (n + 1) the old and new time steps, and ( * ), ( * * ) and ( * * * ) sub-iterations. Also, ∆t = t (n+1) − t (n) is the time-step, and f , g and h are the numerical fluxes defined later on. Note that, equation (11) is implicit in the function φ in each sub-step. Next, we solve an additional equation in order to ensure the divergence-free condition of the fully multi-dimensional operator (Puckett et al, 1997;Aulisa et al, 2003): Finally, we need to specify how the numerical fluxes are treated. These are defined as the space/time integration of the product of the velocity u i and the indicator function H, which is substituted by its approximate counterpart H, i.e., The temporal integration can be replaced by a spatial integration along the upwind path on the velocity field. For example, in the x-direction the upstream path is ∆x Similarly in the other two directions. Thus, the numerical fluxes can be computed as Similarly to equation (9), the numerical integration can be performed analytically in one direction, and numerically in the other two by the two-point Gaussian quadrature.
Once the volume of fluid function has been updated, i.e., φ n+1 is available, differently from what done by Ii et al (2012), the time integration of equation (3) is here performed with a fractional-step method (Kim and Moin, 1985) where the evolution equation is advanced in time with a second-order Adam-Bashforth scheme and a Fast Poisson Solver is used to enforce zero divergence of the velocity field. Due to the non-uniformity of the density, the Poisson equation used to enforce a divergence-free velocity field results in an equation with variable coefficients, i.e., where the pressure p and density ρ are evaluated at n + 1 and the velocity u i is the non-divergence free predicted velocity. In order to employ an efficient FFTbased pressure solver with constant coefficients (see also Ferrante, 2014, 2016), we use the following splitting of the pressure term (Dong and Shen, 2012): where ρ 0 is a constant density equal to the lowest density of the two phases, and p is an approximated pressure obtained by linear extrapolation, e.g., p = 2p n − p n−1 . With this splitting, the Poisson equation can be rewritten as Note that, also the correction step of the fractional-step method needs to be modified accordingly. Fig. 2 The 2D Zalesak's disk: simulation of a slotted disc undergoing solid body rotation (Zalesak, 1979). The black line denotes the exact initial solution, whereas blue and red the solutions obtained after one and five full revolutions. The two left figures are obtained with β = 1 and the two right ones with β = 2. In both sides, the results in the two panels are obtained with 33 and 66 grid points per diameter. Fig. 3 The 3D Zalesak's disk: the three panels represent the exact initial solution, and those after one full rotation with 33 and 66 grid points per diameter.

Code validation
The Zalesak's disk (Zalesak, 1979), i.e., a slotted disc undergoing solid body rotation, is a standard benchmark to validate numerical schemes for advection problems, since the initial shape should not deform under rigid body rotation. The set-up is the same described by Ii et al (2012) and the comparison of the initial shape (black line) and those after one (blue line) and five (red line) full rotations is shown in figure 2. We consider two grid resolutions here, 100 and 200 grid points per box size (being the disc of size 0.3), and two different values of the sharpness parameter β: 1 and 2. The deformed shape of the disk shows an overall good agreement with the initial one, with the comparison deteriorating when more rotations are performed. Better agreement is found on the finer grid and this is further slightly improved in the case with β = 2. As expected, the major differences are found on the sharp edges of the geometry, which are difficult to maintain undeformed. The test has been repeated in 3D and the results are shown in figure 3. Again, quite good agreement is found between the initial and final shapes, with the difference reducing with increasing resolution. Next, we study a deformed interface in a shear flow in order to evaluate the capability of the method to capture a heavily deformed and stretched interface (Rudman, 1997). We consider a circle of radius 0.2π, centered at [0.5π, −0.2(π + 1)] in a box of side 2π with the prescribed velocity field u 1 = sin(x) cos(y) and u 2 = − cos(x) sin(y) for t < 5π, and with the reversed velocity for t > 5π. First, the interface deforms due to the imposed velocity field reaching its maximum stretching at t = 5π; after the flow reversal, the interface deformation starts reducing and at time t = 10π the interface goes back to its initial position. Figure 4 shows the interface shape at time t = 0 (black line), at t = 5π (blue line) and at t = 10π (red line) for a mesh resolution of 100 (left panel), 150 (middle panel) and 200 (right panel) grid points in each direction. We can observe that, the tail of the interface is broken in smaller segments at its maximum stretching (t = 5π), and that the breakup reduces as the mesh is refined. In the reversed phase, all these interface segments merge and the shape is almost circular at the final instant (t = 10π), with some wiggles in the bottom part of the circle, yet reducing as the grid is refined.
We conclude this part on the numerical method performance by noting that the parasitic current often caused by numerical errors in the curvature estimation (Popinet and Zaleski, 1999;Torres and Brackbill, 2000) are effectively suppressed in the present methodology by the continuous VOF distribution with the curved surface reconstruction, as shown by Ii et al (2012).
Finally we present two different validations of our code where we compare our results with those reported in the literature: first, we consider a 3D deformable drop in a simple shear flow, and then 2D and 3D rising droplets in a buoyancy-driven flow.
Shear-driven flow -We consider a 3D unit spherical drop located at the center of a computational domain of size 8 × 4 × 8, with a resolution of 16 grid points for the drop diameter. The top and bottom boundaries move with opposite velocity ±U , giving a shear rateγ = U/2, while periodic boundary conditions are imposed in the streamwise x and spanwise z directions. The same density and viscosity are specified for the spherical drop and the surrounding fluid, the Reynolds number is fixed to 0.1 and three Capillary numbers are studied: 0.1, 0.2 to 0.3. Figure 5(left) shows the steady-state deformed shapes obtained by our numerical simulations (solid lines), which are in excellent agreement with the numerical results reported by Ii et al (2012), shown with points of the same color in the figure.
Buoyancy-driven flow -The method is further validated for a 2D droplet rising in a fluid as studied numerically by Hysing et al (2009). The domain is rectangular with size 1 × 2; no-slip boundary conditions are applied on the horizontal walls and free-slip boundary conditions on the sides. A round droplet with radius R is initially placed at the centerline of the channel at a distance of 0.5 from the bottom wall. The non-dimensional parameters governing the flow are the Reynolds number Re = ρ 1 U g 2R/µ 1 , the Eötvös number Eo = ρ 1 U 2 g 2R/σ, the Capillary number Ca = µ 1 U g /σ and the viscosity and density ratio k µ = µ 2 /µ 1 and k ρ = ρ 2 /ρ 1 , where we have chosen as reference length and velocity scales the drop diameter (2R) and the velocity U g = √ gL, being g the gravitational constant. The simulation is performed with the parameters corresponding to the benchmark test case 1 by Hysing et al (2009), i.e., Re = 35, Eo = 10, Ca = 0.2857, k µ = 10 and k ρ = 10. The comparison of the drop center of mass and rise velocity are reported in Figure 5(right); the nearly perfect agreement between our results and those taken from the literature thus verifies again our implementation. Finally, we extend the comparison to a 3D droplet rising in a fluid as studied experimentally by Legendre et al (2005), (Re = 259, Eo = 0.35, Ca = 0.0025, k µ = 1.56 and k ρ = 1.16). From our simulation we obtain a rising velocity of 8.5 cm/s, compared to 8.1 cm/s found by Legendre et al (2005), with an error of 5% well within the 10% uncertainty in the velocity measurement reported by those authors.

Results
We consider the plane Couette flow of two Newtonian fluids separated by an interface with interfacial tension σ. The density and viscosity of the two fluids are assumed to be equal (ρ 0 and µ 0 ), and the Reynolds number of the simulation is set equal to Re = ρ 0γ r 2 /µ 0 = 0.1, whereγ is the reference shear rate, so that we can consider inertial effects negligible. We denote the solvent as fluid 1 and the suspended phase as fluid 2. The total volume fraction Φ is defined as the volume average of the local volume fraction φ, i.e., Φ = φ . Hereafter, the double · indicates time and volume average while the single · average in time and in the homogeneous x and z directions. Four values of the total volume fraction Φ ≈ 0.0016, 0.1, 0.2, and 0.3 are considered, together with three values of the interfacial tension σ, resulting in the Capillary numbers Ca = µ 0γ /σ = 0.1, 0.2 and 0.4. Note that, the Capillary numbers considered in this study are below the critical value for droplet breakup under shear flow in similar conditions (Cristini et al, 2003;Caserta et al, 2008); indeed, Caserta et al (2008) and Caserta and Guido (2012) found that this system is mainly modified by deformation and coalescence, thus these effects will be investigated in our simulations. The numerical domain is a rectangular box of size 16r ×10r ×16r in the x, y, and z directions, discretised on a Cartesian uniform mesh with 16 grid points per radius r. No-slip boundary conditions are imposed on the solid walls, while periodic boundary conditions are enforced in the homogeneous x and z directions (see figure 1). All simulations start with a stationary flow and a random distribution of spherical droplets of fluid 2 with unit radius across the domain. Note that, the general set-up and most of the parameters used in this study are chosen as in Picano et al (2013) and Rosti and Brandt (2018) where suspensions of rigid and deformable spherical particles are examined. The independence of the results from the grid resolution and domain size were tested by performing two additional simulations for the case with high volume fraction and Capillary number, i.e., Φ = 0.3 and Ca = 0.4: i) a simulation with double grid points in all the directions; ii) a simulation with double domain size in the homogeneous directions x and z. The difference in the results was found to be less than 2%. We start the analysis of the two-fluid flow by showing the mean streamwise velocity profile of the solvent phase u f 1 in the left panel of figure 6. The blue, orange and red colors are used to distinguish the three different volume fractions Φ = 0.1, 0.2 and 0.3, respectively; the shaded area with the same color represents the spread of the data due to the different Capillary numbers studied in this work. In the Newtonian case, shown with a dashed black line, the velocity profile is linear, decreasing from V w at the wall, due to the no-slip condition, to zero at the center line, for symmetry. The wall-normal profile is not straight in the emulsion and the velocity decreases faster than in the Newtonian case close to wall, i.e., the wall-normal derivative of the velocity profile at the wall increases; also, the profile shows a local minimum around y ≈ 0.75h. These differences are enhanced for high values of volume fractions and Capillary numbers, and are strongly related to the distribution of the dispersed phase across the channel, reported in terms of average local volume fraction of the suspended fluid φ in the right panel of the same figure. Indeed, the suspended fluid has a non uniform distribution in the wall-normal direction, with a strong peak in the concentration around y ≈ 0.75h, (the position of the local minimum of velocity), which moves towards the wall for increasing volume fractions. From the figure we can also observe a significant variation with the Capillary number as the volume fraction increases; this result differs from what reported by Rosti and Brandt (2018) for deformable particles, where only minor variations where observed in the volume fraction distribution for different Capillary numbers.
The wall-normal derivative of the streamwise velocity at the wall can be used to estimate the effective viscosity of the non-Newtonian emulsion made of the two Newtonian fluids. We define the effective viscosity µ, normalized by the reference value µ 0 , as Figure 7(left) shows the effective viscosity µ as a function of the total volume fraction Φ for all the Capillary number Ca considered. We observe that the effective viscosity decreases with the Capillary number Ca and exhibits a nonmonotonic behavior with respect to the volume fraction: in particular, it first increases with the volume fraction up to Φ = 20% and then decreases. These data are obtained by time-averaging the results over a time of approximately 40γ −1 after a long transient of approximately 80γ −1 . Indeed, the system is reaching a quasi statistically steady state after this interval for all the cases we have investigated, as shown in figure 7(right) where we report the time history of the effective viscosity µ for three different volume fractions (Φ = 0.1, 0.2 and 0.3) and Capillary number Ca = 0.1. Note that, we have defined as quasi statistically steady state regime a condition where the global properties of the system (e.g., the effective viscosity) change by less than 3% over a time of 20γ −1 .
Next, we test whether an expression similar to that used for rigid and deformable objects can be used also for emulsions. To this end, we fit the effective viscosity of our simulations with the relation proposed by Batchelor and Green (1972), i.e., a second order extension of Einstein (1956) formula, where [µ] is the intrinsic viscosity (equal to 5/2 for rigid dilute particles) and B BG is a coefficient equal to 7.6 for non-Brownian spheres. In the two-fluids system studied here, we measure the intrinsic viscosity [µ] from the simulations with a single droplet in the computational domain (corresponding to the case Φ ≈ 0.0016), i.e., [µ] ≈ (µ − µ 0 ) / (µ 0 Φ). The coefficient B BG is kept as a fitting parameter, whose values are reported in table 1, together with the intrinsic viscosity extracted from our simulations. We see that the secondorder expression is applicable for values of the volume fraction up to Φ = 0.3 in the two-fluid system considered here, provided that both [µ] and B BG are modified to take into account the droplet deformation, i.e., they are assumed as functions of the Capillary number Ca. For rigid particle suspensions, this second-order relation is usually inaccurate for Φ 0.15, as the viscosity increases faster than a second order polynomial (Stickel and Powell, 2005); however, it has been shown recently to apply up to Φ ≈ 0.33 for suspensions of deformable particles . From the data in the table, we observe that in the case of emulsions both the intrinsic viscosity [µ] and the second-order coefficient B BG are decreasing with the Capillary number Ca. Also, B BG is negative for all the Capillary numbers Ca: the function µ = µ (Φ) has a maximum. This can be explained by considering the limiting behaviors for Φ = 0 and 1; indeed, the limit for Φ → 0 and Φ → 1 is µ = µ 0 , i.e., we recover the fluid viscosity, thus the function µ = µ (Φ) must show a maximum

Ca
[µ] B BG 0.1 1.7538 −2.2078 0.2 1.6615 −3.4116 0.4 1.4769 −3.6253. Table 1 The fitting parameter B BG in equation (19)  . The grey lines are the present results for a two-fluids system, the green and purple are results taken from the literature for deformable elastic particles  and capsules (Matsunaga et al, 2016), respectively, while the black dash-dotetd line represent the rigid particle limit (Picano et al, 2015).
for 0 < Φ < 1. These results are in good agreement with the experiments by Caserta et al (2006Caserta et al ( , 2008. Finally, we compare in figure 8 the effective viscosity µ of the present twofluids system with the results taken from the literature at similar volume fraction and Capillary number. In particular, the present results are shown in the figure with grey lines, while green lines represent the case of deformable hyper-elastic particles , the purple lines indicate the case of deformable capsules (Matsunaga et al, 2016) and the black line is the results for rigid particles (Picano et al, 2015). In the deformable cases reported (two-fluids system, capsules and elastic particles), results pertaining three different Capillary numbers are shown in the figure: Ca = 0.1 (dashed line), 0.2 (solid line) and 0.4 (dotted line). We clearly note that, all the deformable cases show smaller effective viscosity than in the presence of rigid particles, with the difference increasing with the volume fraction Φ and the Capillary number Ca. Also, the present results are those with lower effective viscosity, followed by capsules and finally by deformable particles whose results are the closest to the one for rigid particles. The reduction of the effective viscosity with respect to the case of rigid spheres indicates a reduced influence of the suspended phase on the flow of the carrier phase. This is due to the deformability which allows the intrusions to attain a shape which is able to reduce the obstruction to the flow and adapts to it.
We now study the behavior of the interface in the emulsions under investigations here. Figure 9 displays instantaneous configurations of the interface between the two fluids: the top row reports cases at Ca = 0.2 and increasing Φ, i.e., 0.1, 0.2 and 0.3, while the bottom row is for the dilute suspension Φ = 0.0016 and increasing Ca, i.e., 0.1, 0.2 and 0.4. We observe that, in the dilute case, the shape progressively changes from a sphere to an ellipsoid as the Capillary number Ca increases. However, in the case of an emulsion, the interface shape rapidly degenerates from that of an ellipsoid due to the merging process: indeed, we can observe the formation of elongated structures spanning large portions of the domain both in the streamwise and spanwise directions. To quantify the deformation, it is common in the literature to evaluate the so-called Taylor parameter where a and b are the semi-major and semi-minor axis of the inscribed ellipse passing through the center of the droplet in the x-y plane. However, due to the merging, only few interface shapes are truly ellipsoidal, and thus this measure is not appropriate for the present case. An alternative indicator of the deformation and merging of the droplets is reported in the left panel of figure 10, where the mean interface surface area S, normalized by its initial value S 0 , is depicted as a function of the total volume fraction Φ. We observe that S decreases (although not monotonically) with the volume fraction Φ and increases with the Capillary number Ca. In the dilute case (brown symbols) the surface is always greater than the initial spherical shape, as expected being the sphere the shape with lowest surface area to volume ratio. As Ca increases and the av- erage shape changes going from a sphere to an ellipsoid (see the bottom panels in figure 9) the surface area increases. Thus, at a fixed Capillary number and increasing volume fraction, the total surface area reduces due to the merging. We can therefore conclude the following: i) for the low Capillary number cases (Ca = 0.1, displayed with the symbol ×) the merging is the dominant effect because the droplets are only slightly deforming and most of the deviations from the initial configuration are due to the merging; ii) for the high Capillary number cases (Ca = 0.4, shown by ⊡) the deformation is dominant at low volume fractions when the surface area is increasing, but eventually, as Φ increases, the merging reduces the total surface area; iii) finally, the cases with intermediate Capillary number (Ca = 0.2, indicated with * ) present a mixed behavior, although at high Φ the merging always prevails, reducing the total area. In order to evaluate the different flow behaviors in the two fluids, we compute the so called flow topology parameter Q, successfully used by De Vita et al (2018) among others. The flow topology parameter is defined as where D 2 = D ij D ji and Ω 2 = Ω ij Ω ji , being Ω ij the rate of rotation tensor, Ω ij = (∂u i /∂x j − ∂u j /∂x i )/2. When Q = −1 the flow is purely rotational, regions with Q = 0 represent pure shear flow and those with Q = 1 elongational flow. The distribution of the flow topology parameter for the cases with Ca = 0.2 and different volume fractions Φ is reported in Figure 10(right). In particular, we show the probability density function (pdf) of Q in the two liquid phases separately. We observe that in the fluid 1 the flow is mostly a shear flow, as demonstrated by a sharp peak at Q = 0. On the other hand, the flow of fluid 2 shows a broad peak for Q < 0, meaning that the flow is more rotational. However, this peak displaces towards Q = 0 as the volume fraction is increased. This is caused by the increased merging at high Φ, which generates large structures (see figure 9) which spans larger and larger area of the domain.
To understand the rheological behavior of the emulsion, we consider the different contributions to the momentum transfer. The streamwise component of the mean momentum equations for a Couette flow with streamwise and spanwise homogeneous direction reads (Pope, 2001) where the ′ indicates the fluctuation with respect to the mean. Integrating the equation in the wall-normal direction y we obtain where we have defined the total shear stress σ 12 as the sum of the viscous and Reynolds stresses, with the addition of a contribution from the interfacial tension f i , i.e., Note that, although the integral of the interfacial tension over a full droplet is null, the wall-normal profile of its contribution is not, and it must thus be accounted for. Finally, the viscous and Reynolds stresses can be further decomposed into the contribution of the two fluids, see equation (4), giving the following expression: where the averages in the first two brackets are performed separately in the two phases. Each of these contribution have been averaged in time and in the homogeneous directions and displayed in figure 11(left) as function of the wallnormal distance y for the cases with Ca = 0.2 and volume fractions Φ = 0.1, 0.2 and 0.3 (top, middle and bottom panel in the left column). Note that, the stresses are normalized by the total wall value, thus they vary between 0 and 1 and that the Reynolds stress contributions are not shown being negligible in the inertialess limit. The viscous stress of fluid 1 is the only not null component at the wall, due to the fact that the fluid 2 has null concentration there (see figure 6(right)). At the lowest volume fraction shown here the fluid 1 viscous stress is the dominant contribution, being responsible for more than 70% of the total stress at  Fig. 11 Decomposition of the total shear stress σ 12 in its contributions (see equation (25)), normalized by the total wall value σ w 12 . The figures in the left column show the profile of the shear stress components as a function of the wall-normal distance for the case with Ca = 0.2 while those in the right column the percentage distribution of the shear stress contribution as a function of the Capillary number Ca. The top, middle and bottom row indicate the three volume fraction Φ = 0.1, 0.2 and 0.3, respectively. The shear stress contribution due to the interfacial tension, viscous stress of fluid 1 and viscous stress of fluid 2 are shown with the dash-dotted − · −, dashed − − and solid -lines in the right column and with the grey, blue and green bars in the histograms in the right column. each wall-normal location. It has a minimum value around y ≈ 0.75h corresponding to the location of maximum concentration of fluid 2 as previously discussed. As the volume fraction is increased, the fluid 1 viscous stress becomes smaller and smaller; this is compensated by an increase in the fluid 2 viscous stress contribution. Moreover, we observe that the interfacial tension term is almost uniform across the channel, and that its contribution increases with the volume fraction. In figure 11(right) we summarize the stress budgets analysis by showing the volume-averaged percentage contribution of all the non-zero components of the total shear stress, i.e., fluid 1 viscous stress (blue), fluid 2 viscous stress (green), and interfacial tension (grey). Each panel shows how the stress balance change with the Capillary number, and each panel corresponds to a different volume fraction (Φ = 0.1: top panel; 0.2: middle panel; 0.3: bottom panel). For every volume fraction, we observe that the percentage contribution of the viscous stress of the solvent and suspended fluids increases with the Capillary number; on the other hand, the contribution of the interfacial tension slightly decreases with the Capillary number (see figure 10(left)). These results are similar to the one reported for hyper-elastic particles by Rosti and Brandt (2018) who also found an increase of the particle contribution with the elasticity Capillary number (i.e., deformability).

Conclusion
We have implemented and validated a volume of fluid methodology based on the MTHINC method first proposed by Ii et al (2012) which can be used to study multiphase problems. A fully multidimensional hyperbolic tangent function is used to reconstruct the interface and the main advantages can be summarized as follows: i) the geometric reconstruction of the interface is not required; ii) both linear and quadratic surfaces can be easily constructed; iii) the continuous multidimensional hyperbolic tangent function allows the direct calculations of the numerical fluxes, derivatives and normal vectors; iv) the hyperbolic tangent function prevents the numerical diffusion that smears out the interface transition layer. Moreover, the implementation presented here allows the use of an efficient FFT-based fractional step method to solve the system of equations, and is based on the splitting of the pressure in a constant and a varying term, which makes the FFT solver applicable also when the density in the two phases differs.
We have studied the rheology of a system of two Newtonian fluids in a wall-bounded shear flow, i.e., plane Couette flow, at low Reynolds number such that inertial effects are negligible. The rheology of the emulsion is analyzed by discussing how the effective viscosity µ is affected by variations of the volume fraction Φ and Capillary number Ca. The effective viscosity µ is a non-linear function of both this parameters µ = µ (Φ, Ca) and the emulsion shows an effective viscosity lower than the one of suspensions of rigid particles, deformable elastic particles and capsules. Differently form these other cases, here the droplets can merge, especially at high volume fractions and high Capillary numbers. Also, the effective viscosity curve has a negative second derivative with respect to Φ, thus suggesting the presence of a maximum in the curve at higher volume fractions. The overall deformation of the emulsion has been quantified in terms of the interface surface area and the flow within the two phases studied by means of the flow topology parameter. In particular, we have shown that the flow of the suspending fluid is mainly a shear flow, while that of the dispersed fluid is more rotational.
Finally, we have analysed the contributions to the total shear stress of the two fluid phases and of the interfacial force and showed that the Reynolds stress contributions are negligible at this Reynolds number and that the viscous stress of the two fluids provide the dominant contribution. These have a non-uniform distribution across the channel, with the percentage contribution of both phases increasing with the volume fraction and Capillary number. An important contribution comes from the interfacial tension, whose effect increases with the volume fraction and decreases with the Capillary number.
The study presented in this work will be extended to consider different viscosity and density ratios between the two fluids, and to better understand the non-monotonic behavior of the effective viscosity and its effect on the rheological behavior of emulsions.