Deformation of clean and surfactant-laden droplets in shear ﬂow

In this work we study the deformation of clean and surfactant-laden droplets in laminar shear-ﬂow. The simulations are based on Direct Numerical Simulation of the Navier–Stokes equations coupled with a Phase Field Method to describe interface topology and surfactant concentration. Simulations are performed considering both 2D (circular droplet) and 3D (spherical droplet) domains. First, we focus on clean droplets and we characterize the droplet shape and deformation. This enables us to deﬁne the range of parameters in which theoretical models well predict the results obtained from 2D and 3D simulations. Then, surfactant-laden droplets are considered; the main factors leading to larger droplet deformation are carefully described and quantiﬁed. Results obtained indicate that the average surface tension reduction and the accumulation of surfactant at the tips of the deformed droplet have a dominant role, while tangential stresses at the interface (Marangoni stresses) have a limited effect on the overall droplet deformation. Finally, the distribution of surfactant over the droplet surface is examined in relation to surface deformation and shear stress distribution.


Introduction
The deformation of a droplet in a simple shear flow is of fundamental relevance in a number of flowing system of industrial and biological interest. Possible applications include the formation and rheology of emulsions [7], emulsifying devices [18], polymer blending [9], oil recovery [19] and the study of red blood cells [36].
This problem was first tackled by Taylor [33,34], who developed an analytic formula able to predict the deformation of a droplet in shear flow. This formula, developed under the hypotheses of small deformations and negligible inertia effects, constitutes a simple tool for the calculation of droplet deformation. Within these hypotheses, the droplet steady-state deformation is a function of the capillary number, Ca (dimensionless number that expresses the ratio between viscous and surface tension forces) and of the viscosity ratio among the two phases. The capillary number highlights the two principal factors involved in droplet deformation: the shear rate, which tries to deform the droplet, and surface tension, which acts to restore the spherical shape of the droplet.
The simple configuration considered and the immediacy of the analytic formula made it a widely used tool for benchmark and validation of numerical methods and experimental facilities. Later, the accuracy of Taylor's formula was improved by Shapira and Haber [24], who introduced a correction for confinement effects. These sets of prediction tools are now a well assessed and easy to use benchmark for numerical methods and codes for multiphase flows [12,26,31,35,38,40].
Being a validation test, a fast and lightweight simulation is often preferable; for this reason twodimensional (2D) simulations are usually preferred to their three-dimensional (3D) counterparts. Clearly, 2D and 3D cases are substantially different: the formers are circular droplets (cylindrical when extended to 3D), while the latters are spherical droplets. Reducing a three-dimensional case to a two-dimensional one is indeed a drastic simplification, in which several effects are strongly affected (suppression of capillary instabilities, longer and stronger near-contact droplet interaction) but main effects are kept (development of high interfacial curvature regions, tip streaming) [41]. A first comparison between 2D and 3D simulations results was performed by Zhou and Pozrikidis [41] showing that 2D and 3D droplets exhibit a comparable behavior. Similarly, Tang et al. [32] arrived at analogous conclusions.
In this work we use Direct Numerical Simulation (DNS) of the Navier-Stokes equations coupled with a two-order-parameter phase field method (PFM) to describe the interface topology and the surfactant concentration [15,26,39]. We start by comparing the deformation of circular (2D) and spherical (3D) droplets and we extend previous works analyses [32,41] to investigate the limits of Taylor's formula. In particular, we trace back the similar deformation experienced by circular and spherical droplets at low Ca to the reduced shrinkage of the droplet in the direction normal to the velocity-velocity gradient plane. Then, we consider the effect of a soluble surfactant on the overall droplet deformation. The surfactant affects the droplet deformation introducing three additional effects: (1) surfactant reduces surface tension, thus leading to a lower average surface tension over the interface; (2) the external shear stresses accumulate surfactant at the droplet tips, further reducing the local surface tension (of the droplet tips); (3) the non-uniform surface tension over the interface generates stresses tangential to the interface (Marangoni). Simulations results allow us to quantify the relative contribution of each different effect in the overall droplet deformation. In particular, we found that the average surface tension reduction produced by the surfactant has a major contribution. As a consequence, rescaling the capillary number on the average surface tension (and thus considering the average surface tension decrease), droplet deformation can be well predicted by the Taylor analytic formula. Finally, we characterize the surfactant distribution over the interface and we compare qualitatively the different surface tension distribution (2D versus 3D) droplets and the generated Marangoni stresses.
The paper is organized as follows: first, the simulation framework employed is presented in Sect. 2, then the results obtained from clean and surfactantladen droplets, and the comparison with analytical predictions are presented and discussed in Sect. 3; conclusions are drawn in Sect. 4.

Methodology
The dynamics of a multiphase flow with a surfactant is modeled coupling pseudo-spectral-based solution of the Navier-Stokes (NS) equation with a Phase Field Method (PFM), used in a two-order-parameter formulation, to compute the interface dynamics and the surfactant concentration [15,26,39]. In the following, the governing equations of the two order parameters, phase field / and surfactant concentration w, will be derived and then coupled with continuity and Navier-Stokes (NS) equations to describe the hydrodynamics of the system.

Governing equations
We consider a ternary system composed of a soluble surfactant and two immiscible phases. In the framework of the phase field method, such a system can be described using two order parameters. A first-order parameter, the phase field /, is used to describe the interface. It is uniform in the bulk of the two phases (/ ¼ AE1) and it undergoes a smooth transition across the interface. An additional order parameter, w, is used to describe the surfactant concentration. This second order parameter is uniform in the bulk of the phases and reaches a maximum at the interface, where surfactant molecules preferentially accumulate. The phase field and the surfactant concentration are governed by two Cahn-Hilliard-like (CH) equations, which in dimensionless form are: In the above equations, u ¼ ðu; v; wÞ is the velocity vector, l / and l w are the two chemical potentials, M / and M w are the two mobilities (or Onsager coefficients) and Pe / and Pe w are the two Péclet numbers. The latter ones represent the ratio between convective and diffusive phenomena for the two order parameters. The expressions of the chemical potentials l / and l w are derived from a two-order-parameter Ginzburg-Landau free energy functional F ½/; r/; w. The functional is modeled as the sum of five different contributions: where X is the domain considered.
The first term, f 0 , is the ideal part of the free energy and describes the tendency of the system to separate into two pure fluids (/ ¼ AE1); this phobic behaviour is modelled with a double-well potential: The mixing energy term, f mix , accounts for the surface tension energy stored at the interface and is defined as: In the above expression the Cahn number, Ch, sets the thickness of the thin layer thickness. These two contributions are a function only of the phase field / and its gradient r/; their expressions match those adopted to describe a clean system (absence of surfactant) [21][22][23]. Surfactant is modeled with three additional contributions to the energy functional F ½/; r/; w. The first term is an entropy term, f w , and expresses the entropy decrease obtained when the surfactant is uniformly distributed in all the domain. Its expression is the following: This contribution bounds the value assumed by w to the range between w ¼ 0 (no surfactant) and w ¼ 1 (saturation of surfactant); the parameter Pi sets the surfactant diffusivity. The second term, f 1 , describes the accumulation of the surfactant at the interface; indeed surfactant molecules preferentially gather at the interface exposing their heads towards the water phase and their tails towards the other phase.
The last contribution, f E x , penalizes the presence of surfactant in the bulk of the two phases and it is defined as: The term f E x has a relevant contribution in the bulk of the two phases (/ ¼ AE1) and vanishes at the interface (/ ' 0). The parameter E x controls the surfactant solubility in the bulk of the two phases. The expressions of the chemical potentials are obtained by taking the variational derivative of the free energy functional with respect to / and w: In the phase field chemical potential, l / , the terms depending on the surfactant concentration have been removed. These terms induce an unphysical behavior of the interface [26,39] and, thus, to restore the correct interfacial behavior they have been neglected. Further details on this point can be found in Yun et al. [39] and Soligo et al. [26]. From the expressions of the chemical potentials, the equilibrium profiles of the two order parameters can be obtained. For the phase field /, the equilibrium profile is determined by the competition between f 0 and f mix . At the equilibrium, l / ¼ l eq / in the entire domain and thus: The phase field matches the values / ¼ AE1 in the bulk of the phases (x ! AE1) and undergoes a smooth transition following a hyperbolic tangent profile across the interfacial layer. Likewise, the surfactant equilibrium profile can be deduced from Eq. (10): at the equilibrium, the surfactant chemical potential is constant throughout the entire domain. The surfactant equilibrium profile results in: where the auxiliary variable w c is a function of the phase field solely: At the equilibrium, surfactant concentration matches the value w b in the bulk (/ ¼ AE1) and reaches its peak, w 0 ¼ wj /¼0 , at the interface. The peak value of w (at the equilibrium) depends on the surfactant bulk concentration, w b , and on the parameters E x and Pi. The mobilities are set to M / ¼ 1 [3] and to M w ðwÞ ¼ wð1 À wÞ respectively. Hence, the following expressions for the Cahn-Hilliard-like equations are obtained: Equations (14) and (15) describe the time evolution of phase field and surfactant concentration and are coupled with the Navier-Stokes (NS) equations to describe the hydrodynamics of the multiphase flow. In the most general case, this approach can handle nonmatched properties [8,21]. However, in this work, we focus on the effect of surfactant solely and we thus consider two phases with matched density (q ¼ q 1 ¼ q 2 ) and viscosity (g ¼ g 1 ¼ g 2 ). Therefore, continuity and Navier-Stokes equations can be written as follows: where u is the velocity field, p is pressure and the last term of the right-hand side is the interfacial term, which represents the surface tension forces. In the interfacial term, s c ¼ jr/j 2 I À r/ r/, is the Korteweg tensor [16] and f r ðwÞ is the equation of state for surface tension. The expression employed implicitly accounts for both normal (capillary) and tangential (Marangoni) components of the surface tension forces. In the Navier-Stokes equations, two dimensionless groups are present: the Reynolds number, Re, ratio between inertial and viscous forces and the Weber number, We, ratio between inertial and surface tension forces. In the definition of We, the surface tension of a clean interface (referred in the following as r 0 ) is used as reference. The action of the surfactant on surface tension is described using an Equation Of State (EOS) [4,6]; in this work we adopt a Langmuir EOS (Szyszkowski equation). In a dimensionless form this EOS reads: The elasticity number, b s , quantifies the strength of the surfactant action: for a fixed concentration, the higher is b s , the stronger is the surface tension reduction. The original Langmuir EOS is valid in the limit of moderate surfactant concentrations: experimental studies [6,30] showed that surface tension never reduces below roughly half of its clean value. This effect is not captured by the Langmuir EOS, which instead predicts an always decreasing surface tension for increasing surfactant concentrations. Therefore, this EOS well predicts the surface tension reduction up to a saturation concentration (w s ), which corresponds to the lowest value achievable by surface tension, f r ðwÞ ¼ 0:5.
To account for this feature, the Langmuir EOS has been modified to limit the surface tension decrease. The EOS adopted in this work is the following:

Numerical method
The governing Eqs. (14)(15)(16)(17) are solved in a closed channel geometry using a pseudo-spectral method [5,11,20]. In particular, the equations are discretized using Fourier series in the streamwise and spanwise directions (x and y) and Chebyshev polynomials along the wall-normal direction (z). The governing equations are advanced in time using an IMplicit-EXplicit (IMEX) scheme. The linear diffusive term of the equations is integrated using an implicit scheme, whereas the non-linear term is integrated using an explicit scheme. An Adams-Bashforth scheme is used for discretization of the non-linear terms of the Navier-Stokes equation, while a Crank-Nicolson scheme is used for the linear term. The two Cahn-Hilliard-like equations are time-discretized using again an Adams-Bashforth algorithm for the nonlinear terms, while the linear terms are discretized using an implicit Euler algorithm to improve the numerical stability [3,37]. The Navier-Stokes equation is rewritten and solved in the so-called velocityvorticity formulation. Instead of three 2nd order equations for each component of the velocity, a 4th order equation for the wall-normal component of the velocity and a 2nd order equation for the wall-normal component of the vorticity are obtained [13,27]. The phase field and the surfactant concentration transport equations, 4th and 2nd order respectively, are directly solved in the formulation presented above [25]. Specifically, the phase field is split into two 2nd order equations, while the surfactant transport equation can be directly solved (2nd order equation). Further details on the numerical method can be found in our recent work [26].

Boundary and initial conditions
A closed channel setup is employed to simulate the deformation of a single droplet in shear flow (see Fig. 1). Periodic boundary conditions are imposed on all variables in the x and y directions (streamwise and spanwise directions). For the flow field, no-slip boundary conditions are enforced at the two solid walls (z ¼ AEh). In particular, at the walls, the streamwise velocity is set equal to the top and bottom wall velocity, uðz ¼ AEhÞ ¼ AEu w ¼ AE1. For both phase field and surfactant concentration no-flux boundary conditions are enforced at the walls: The initial flow field is a linear profile along the z direction, uðzÞ ¼ z=h. The phase field is initialized so to obtain a spherical (circular for 2D simulations) droplet located at the channel center. The phase field has a uniform value in the droplet and carrier phase (/ ¼ AE1) and undergoes a smooth transition following its equilibrium profile across the interface, Eq. (11). The surfactant concentration is also initialized with its equilibrium profile, Eq. (12): in the bulk of the phases (/ ¼ AE1) it is equal to the surfactant bulk concentration, w b , while at the interface it reaches its maximum value.

Simulation setup
All simulations have been performed in a laminar shear flow configuration, Fig. 1. A single drop, with diameter d ¼ 0:8h, is initialized at the channel centre. The domain has dimensions L x Â L y Â L z ¼ 2ph Â ph Â 2h for the 3D simulations. When performing 2D simulations the domain is shrank in the y-direction and has dimensions L x Â L z ¼ 2ph Â 2h (red rectangle of Fig. 1). The computational domain is discretized with N x Â N y Â N z ¼ 512 Â 256 Â 513 grid points (3D simulations) and with N x Â N z ¼ 512 Â 513 (2D simulations). The grid spacing is uniform along x and y directions, while for the wall-normal direction we adopt Chebyshev collocation points. Once defined the grid, the Cahn number, parameter that controls the thickness of the interfacial layer can be set. In particular, the accurate description of the steep gradients at the interface requires a minimum of 5 grid points across the interface. To meet this requirement the Cahn number has been set to Ch = 0.025. The Péclet number for the phase field, Pe / , is determined based on the scaling Pe / ¼ 3=Ch ¼ 120. Concerning the surfactant parameters, the Péclet number has been set to Pe w ¼ 100, while the temperature dependent coefficient Pi and the solubility number E x has been set respectively to Pi ¼ 1:689 and E x ¼ 0:1. The elasticity number has been chosen within the range of a moderate strength surfactant, b s ¼ 0:5. Finally, in order to ensure creeping flow conditions and thus negligible inertial effects, the Reynolds number, Re ¼ qu w h=l, is set to 0.1. The simulations consider different capillary numbers, Ca ¼ lu w d=ðr 0 hÞ (ratio between viscous and surface tension) and surfactant bulk concentrations, w b (parameter controlling the total amount of surfactant present); please note that in defining the capillary number, the surface tension of the clean system, r 0 , is used as reference. We consider three different capillary numbers, starting from Ca ¼ 0:062 (higher surface tension) up to Ca ¼ 0:187 (lower surface tension), and three different surfactant loadings, starting from w b ¼ 0 (clean system) up to w b ¼ 0:02 (highest amount of surfactant in the system). A list of the simulation parameters is reported in Table 1; for each case, a 2D and a 3D simulation have been performed.
In addition, to investigate the effects introduced by the surfactant, two additional simulations (2D and 3D) at the highest capillary and at the intermediate surfactant loading (Ca ¼ 0:187 and w b ¼ 0:01, case L3) were performed neglecting the Marangoni stresses. These simulations are used to investigate their effect on the deformation of the droplet. To highlight the contribution of Marangoni stresses, the surface force term in the NS equations can be split in the normal (capillary) term and the tangential (Marangoni) term: 3 Results In the following, results obtained from the simulations will be presented and carefully discussed. First, we will focus on clean-droplets (C-series): the shape and the deformation of the droplet obtained from our 2D and 3D simulations will be compared against previous works [10,14,17,41] and with analytic predictions [24,34]. Then, we will consider surfactant-laden droplets (L-and H-series): the shape and the deformation of the surfactant-laden droplets will be compared against the clean cases and the analytic predictions [24,34]. The importance of the main factors which lead to a larger deformation of the droplet will be quantified and the surfactant distribution over the droplet interface will be characterized.
For each of these a clean and two surfactant-laden systems (w b ¼ 0:01 and w b ¼ 0:02) have been considered. Each combination of the parameters (capillary and surfactant bulk concentration) is simulated on a 2D and on a 3D domain

Clean droplets
We start our discussion considering the behavior of clean droplets under shear flow. This benchmark is often used as a validation tool for numerical codes. The droplet, initially spherical (circular in 2D simulations), is deformed by the imposed shear flow till a new steady-state shape is reached. The final shape is the result of the competition between the viscous forces, which try to deform the droplet, and the surface tension forces, which try to restore the spherical shape. The capillary number is the ratio between these two contributions, and, when a clean system is considered, is the main parameter that controls the final shape of the droplet. This final shape can be characterized by the length of the principal axes. Specifically, following the sketch reported in Fig. 2, we can identify the major axis of deformation, a, the minor axis of deformation, b, and the third axis, c, (only for the 3D simulations). Combining these lengths, the Taylor deformation parameter, D (ratio between the difference and the sum of major and minor axes) can be computed: In the limit of small deformations (and thus low capillary numbers) a steady-state deformation is always attained as the droplet never undergoes breakage. In these conditions, the value of the deformation parameter D can be also predicted via the analytic relation proposed by Taylor [34], which states that the droplet deformation parameter D is proportional to the capillary number (thus inversely proportional to surface tension): where the coefficient 35/32 is specific for the case of two phases with matched viscosity (as considered in this work). The relation can be modified to include also the confinement effects introduced by the two walls; following the work of Shapira and Haber [24], the final relation reads: with C SH being a numerical coefficient, C SH ¼ 5:6996 [24].

Droplet deformation
For each one of the clean cases considered (C1, C2 and C3), the time behavior of the deformation parameter D is computed and compared with the steady-state value predicted by the analytic relation of Taylor [24,34]. Figure 3a shows the time behavior of the deformation parameter D for the three capillary numbers analyzed.
Results from 2D simulations are reported with a dashed line while those from 3D simulations are reported with a solid line. The droplet is initially spherical (circular in the 2D simulations) and therefore D ¼ 0. Then, for all the cases considered here, the deformation parameter reaches a steady-state value after an initial transient. The time required to achieve this final configuration depends on the capillary number: for larger Ca (larger deformations), a longer time is required. From the results reported we can notice that for the range of capillary numbers here considered, there is an excellent agreement between the results obtained from 2D and 3D simulations. Even though similar deformations are achieved in the limit of low capillary numbers, 2D and 3D droplets deform in a different way: Figure 3b reports the length of the major and minor axes normalized by the initial droplet diameter, d 0 . As the capillary number is increased, the two axis start to differentiate between 2D and 3D cases: in particular, a clear difference can be Fig. 2 The final steady-state shape of the droplet can be characterized by the length of the three principal axes: the major axis of deformation, a, the minor axis of deformation, b, and the third axis, c. The latter one is aligned with the y direction and is thus not shown in this sketch. An additional reference frame (x 0 ; y 0 ; z 0 ) is defined, with the axes corresponding to the deformed droplet principal axes. The different regions of the droplet have been also highlighted for ease of reference: tips (green), belly (red) and sides (blue). (Color figure online) appreciated for the higher capillary number, Ca ¼ 0:187. It can be observed how 3D droplets elongate more than their 2D counterpart (higher a=d 0 ) but, at the same time, they undergo a lower compression along the minor axis (higher b=d 0 ). Thus, the longer a axis is balanced out by a longer b axis for 3D droplets, resulting in a similar drop deformation between 2D and 3D cases in the limit of low capillary numbers. As the capillary number is increased, e.g. Ca ¼ 0:187, it can be noticed how the major axis for 3D droplets elongates more than that of their 2D counterpart, while the minor axis is similar in both cases. This difference results in an increased deformation for 3D droplets with respect to 2D ones. This feature can also be appreciated by comparing the cross section of the deformed (steady-state deformation) 2D and 3D droplets, Figure 4. At the lowest capillary numbers (Ca ¼ 0:062 and Ca ¼ 0:125) the cross sections fall one on top of the other, while a clear difference can be appreciated in Fig. 4c, where the difference between the major axes is considerably larger than that between the minor axes. This observation agrees with the increased deformation obtained for 3D droplets.
To test the accuracy of the method, in Fig. 5, we compare the steady-state value of D obtained from our simulations with previous works. Our results are plotted as empty red squares (2D simulations) and as empty red diamonds (3D simulations). The analytic relation of Taylor [34], corrected by Shapira and Haber [24], is plotted with a black solid line. In addition, the results obtained by Zhou and Pozrikidis [41] (2D Boundary Integral Method), Guido and Villone [10] (experiments), Li et al. [17] (3D simulation) and Komrakova et al. [14] (3D simulation) are also reported. Comparing the different results, we can observe that an overall agreement among them is found. Specifically, our results are in very good agreement with those obtained from previous numerical simulations [14,17,41] and with the predictions of the analytic formula [24,34]. The results are also in fair agreement with previous experimental works [10]. However, the viscosity of the droplet used in the experiments is slightly larger and this leads to a smaller deformation and thus to slightly different results. Finally, we can observe that for the largest capillary number studied (Ca ¼ 0:187), the results obtained for 2D and 3D simulations start to deviate and the droplet deformation obtained from 3D simulation is slightly larger than that of 2D simulations. This behavior is in agreement with previous findings of Afkhami et al. [1], who performed an extensive comparison of 2D and 3D droplet deformation, showing that only at low capillary numbers a good agreement between 2D and 3D simulations is found.

Evolution of the droplet principal axes
To understand the origin of the agreement between 2D and 3D simulations at low capillary numbers and of the subsequent divergence at larger capillary numbers, we study the evolution of the shape of the droplet by examining the behavior of the three principal axes. We believe this is an important issue, as the deformation of a droplet in shear flow is a commonly used benchmark in the validation of numerical methods. Therefore, it is important to assess the validity of 2D simulation with respect to their 3D counterpart. In addition, this analysis will allow us to obtain further insights on the range of validity of Taylor analytical formula [34], which was obtained under the assumption that the shape of the deformed droplet is a prolate spheroid with major axis a and two minor equal axes b and c. Figure 6 shows the evolution over time of the axes length as calculated from our simulations. The major axis, a, is reported with solid lines, the minor axis, b, with dashed lines and the third axis, c, with dotted lines. The axes length is normalized by the initial droplet diameter, d 0 . The different colors refer to Ca ¼ 0:062 (black), Ca ¼ 0:125 (blue) and Ca ¼ 0:187 (red). Starting from the beginning of the simulation, t ¼ 0, the axis a elongates, the axis b shrinks, while the axis c does not change during this initial part of the simulation. When the axes a and b have almost reached their steady-state values, the axis c starts shrinking. Shrinkage/elongation magnitude increases with the capillary number, as the droplet becomes more deformable. With the assumption of a prolate spheroid, as in Taylor [34], axes b and c must be equal. This assumption leads to a larger integral of surface forces over the interface (on average the curvature is higher) and a stronger shear rate is needed to achieve the same deformation of an unconstrained droplet (axes a, b and c can vary independently). Hence, at larger capillary numbers, Taylor analytic formula underestimates the droplet deformation. This feature is indeed more pronounced at larger capillary numbers where results from experiments and simulations predict a higher droplet deformation. The observed shrinkage of the axis c plays a crucial role in the 2D simulations where it is constrained to be constant; being a two-dimensional case, no out-ofplane velocity can appear and thus no shrinkage is present. Due to this constraint, for a fixed Ca, the 2D circular droplet undergoes a lower deformation with respect to its 3D counterpart. Indeed, the shrinkage of the axis c during deformation enhances the droplet deformation. These findings are in excellent agreement with the behavior observed in the experiments by Guido and Villone [10].
The contribution of the third axis shrinkage to the overall deformation can be also graphically visualized. Figure 7 shows the spanwise velocity in a y 0 À z 0 plane (see Fig. 6 for further details on this reference frame) for the case C3 (Ca ¼ 0:187) and refers to time t ¼ 1:0. Observing the velocity map, the feeding from the sides of the droplet towards the core region is clear. These fluxes are a direct consequence of the shrinkage of the axis c and favor the droplet deformation.
From our simulations of clean droplets in simple shear flow, we confirmed that the third axis undergoes limited shrinkage at low capillary numbers, as found in experiments [10]. This limits the influence of this axis on the deformation of the droplet. Hence, a good agreement between simulations (2D and 3D), experiments and analytic predictions can be achieved. Increasing the capillary number (and thus the droplet deformation) analytic predictions and results from 2D simulations start to deviate from results obtained from 3D simulations and experiments.

Surfactant-laden droplets
We now move to the discussion of the results obtained from the simulations of the surfactant-laden droplets in shear flow. Compared with the clean droplet case, the final steady-state shape of the droplet will be influenced by three additional factors: (1) the surfactant decreases the average surface tension; (2) the surfactant accumulates on the tips of the deformed droplet producing non-uniform capillary forces, Fig. 8a; (3) a shear-induced inhomogeneous distribution of surfactant on the droplet interface gives rise to tangential stresses at the interface, Fig. 8b. Each of  Fig. 8 Panel a shows the surfactant distribution over the droplet interface. Due to the action of the external shear stress, surfactant accumulates at the droplet front and back, while the droplet belly is depleted of surfactant. Panel b depicts the dimensionless surface tension, rðwÞ=r 0 , over the droplet interface. The Marangoni stresses, which originate from the surface tension gradients, are reported using black unit length vectors. In both panels, results obtained from 2D cases are reported on the left while those obtained from 3D cases are reported on the right. The snapshots are taken at t ¼ 2:5 (steadystate deformation) and refers to the case H3 (w b ¼ 0:02, Ca ¼ 0:187) these effects contributes in a different way to the overall droplet deformation. In particular, the first (lower average surface tension) and second (more surfactant at the droplets tips) effects favor droplet deformation; conversely, the third effect (inhomogeneous distribution of surfactant and thus of surface tension) gives rise to tangential stresses that hinder the droplet deformation.
To investigate and quantify these effects, for each Ca considered before (cases C1, C2 and C3), two further cases with surfactant bulk concentrations w b ¼ 0:01 (cases L1, L2 and L3) and w b ¼ 0:02 (cases H1, H2 and H3) have been computed. In the following, we will first focus on the overall deformation; then the role played by the different effects will be analyzed and the surfactant distribution characterized.

Droplet deformation
To highlight the main differences between the behavior of clean and surfactant-laden droplets, we start considering the droplet deformation. Figure 9 shows the steady-state value of the deformation parameter D for the cases L1, L2 and L3, panel (a) and H1, H2 and H3, panel (b). With respect to the clean cases, whose results well match the analytic relation, we can immediately notice that the surfactant increases the deformation of the droplet. This effect becomes more pronounced increasing the surfactant bulk concentration and thus the total amount of surfactant present, cases H1, H2 and H3. The comparison between results obtained from 2D and 3D simulations exhibits a good agreement at low capillary numbers while, at larger Ca, results start to deviate. This difference becomes more pronounced increasing the surfactant bulk concentration (case H1, H2 and H3): the larger is the deformation experienced by the droplet, the larger is the difference between 2D and 3D simulations.
With the aim of quantifying the contribution to the droplet deformation produced by each of the three surfactant-induced effects, we can rescale the results by considering an effective capillary number, Ca e . This can be calculated using the actual value of the surface tension, hri, instead of the clean system surface tension, r 0 : The effective value of the surface tension is computed averaging the local value of the surface tension over the entire droplet interface. It is important to note that the non-uniform distribution of the surfactant over the interface produces in turn a non-uniform distribution of the surface tension.
Reporting the results using Ca e as reference parameter, the contribution of the average surface tension reduction on the overall deformation can be filtered out. In Fig. 10, the results obtained from the cases L1, L2, L3 (surfactant bulk concentration w b ¼ 0:01) and H1, H2, H3 (surfactant bulk concentration w b ¼ 0:02) are reported using the effective capillary number, Ca e . Interestingly, we can observe that employing Ca e as reference parameter, results are in good agreement with the analytic relation. A better agreement is obtained between simulations results and analytic predictions for the lowest surfactant bulk concentration (cases L1, L2 and L3), while for the highest surfactant bulk concentration (cases H1, H2 and H3) the agreement is slightly worse (especially considering the 3D simulations). This difference can be addressed to the different amount of surfactant available, which is lower for cases L1, L2 and L3. Therefore, the surfactant distribution (and surface tension) is more homogeneous and the droplet deformation is mainly determined by the average surface tension reduction. All these observations seem to confirm that for the range of parameters investigated here, the increased deformation experienced by the droplet is largely due to the average surface tension reduction.

Surfactant distribution over the droplet interface
To characterize the surfactant distribution over the droplet interface, we compute the joint Probability Density Function (PDF) of surfactant concentration at the interface and interface mean curvature, j. The mean curvature of the interface is the semi-sum of the two principal curvatures, j 1 and j 2 , and can be directly obtained from the phase field. In particular, it can be computed from the divergence of the normal to the interface, n, defined as [2,29]: and the curvature results in: In the above equations, the minus sign is needed to get the interface outward pointing normal (/ ¼ þ1 in the droplet and / ¼ À1 in the carrier fluid). For the 2D simulations (circular droplet), only one principal curvature is defined and thus the second principal curvature is j 2 ¼ 0. Figure 11 shows the joint PDF of the cases L1, L2 and L3, which refer to the surfactant bulk concentration w b ¼ 0:01. Results refer to 2D (left column) and 3D (right column) simulations, while capillary number increases from top (Ca ¼ 0:062) to bottom (Ca ¼ 0:187Þ. In the panels, a red vertical dashed line identifies the initial curvature, j 0 ¼ 2=d 0 (2D) and j 0 ¼ 4=d 0 (3D), while a red horizontal line identifies the initial surfactant concentration at the interface, w 0 . First, we can notice the effect of the capillary number on the droplet mean curvature. As the capillary number is increased, the droplet undergoes a stronger deformation, thus higher (at the tips) and lower (in the central region) values of the mean curvature can be found. Conversely, the surfactant range at the interface is not particularly affected by the capillary number. In addition, we can also observe the different range of curvature values sampled between 2D and 3D simulations: the second principal curvature for the circular droplet is always zero (cylindric surface) so the mean curvature for 2D cases is always about half that of their 3D counterpart. The analysis of the results highlights a clear common trend: higher surfactant concentrations are more likely to be found in high-curvature regions. However, as expected, the results obtained from 2D and 3D simulations exhibit a different behavior: for 2D simulations, a bimodal distribution is found while, for 3D simulations, a trimodal distribution is obtained.
The bimodal distribution of 2D simulations originates from the asymmetric distribution of surfactant with respect to the major axis, a: the imposed shear flow sweeps surfactant towards the back and the front of the droplet (see Fig. 2 for all the references to the droplet regions). This leads to regions with the same curvature that experience different surfactant concentrations. As the capillary number increases, panels (c) and (e), this asymmetric distribution becomes clearer and the two branches of the joint PDF part. This distribution can be better appreciated from the left part of Fig. 8a: surfactant accumulates at the droplet back (and front, not shown), while the belly is depleted of surfactant (dark red color). The left side of the joint PDF (low j, low w) corresponds to the droplet belly: low surfactant concentration in regions characterized by a lower curvature.
For the 3D simulations, a third branch appears. The third branch corresponds to the side area: at the sides, surfactant concentration is low while the interfacial curvature is relatively high, Fig. 8a. This additional branch becomes more evident at high capillary numbers, panels (d) and (f). Indeed, at low capillary numbers, the axes deformation is limited, thus a narrow range of curvature values is found and surfactant distribution is more homogeneous.
The results obtained from the cases H1, H2 and H3, surfactant bulk concentration w b ¼ 0:02, exhibit a similar behavior and thus they have not been reported. Overall, the findings confirm the tendency of the surfactant to accumulate in high curvature regions, which are also stagnation points [28]. However, the resulting surfactant distribution is not straightforward and it is affected by the flow condition found in the different regions of the droplet.

Effect of Marangoni stresses
In the previous section, we were able to show the preferential accumulation of the surfactant at the droplet tips (high curvature regions). This phenomenon, together with the average surface tension reduction, leads to an increase of the droplet deformation. The Marangoni stresses enter this picture with a negative contribution hindering the droplet deformation. Specifically, the Marangoni stresses follow the surface tension gradient: they are tangential to the interface and are directed from low surface tension regions (high surfactant concentration) towards high surface tension regions (low surfactant concentration), thus they are directed from the droplet tips towards the belly area, as can be appreciated from Fig. 8b, in which black arrows shows the direction of Marangoni stresses. To quantify their contribution to the overall droplet deformation, the case L3 has been recomputed neglecting the tangential stresses (see Eq. 22). The resulting values for the droplet deformation parameter D are reported in Table 2. The results indicate that the influence of the Marangoni stresses is almost negligible and the deviation is below 1% for both 2D and 3D simulations.

Conclusions
In this work, we studied the deformation of clean and surfactant-laden droplets in a laminar shear flow. The investigation is based on pseudo-spectral solution of the Navier-Stokes equations coupled with a Phase Field Method (PFM) to describe the interfacial phenomena (interfacial topology and surfactant concentration).
We first studied the behavior of clean droplets; the simulations outcomes enable us to characterize the shape and deformation of the droplet. This characterization allows us to define the range of parameters in which 2D and 3D simulations well match analytic predictions and experimental results. Specifically, at low capillary numbers (small deformations), results from 2D and 3D simulations are in excellent agreement and thus lightweight 2D simulations can be used as benchmark. By opposite, for larger capillary numbers (larger deformations), the shrinkage of the third axis of the droplet becomes considerable and results obtained from analytic predictions and 2D simulations deviate from those obtained from 3D simulations and experiments. Then, we focused on surfactant-laden droplets and we characterized the droplet deformation and the surfactant distribution over the droplet interface. The results indicate that the increased deformation of the surfactant-laden droplets is largely produced by the average surface tension reduction and by the accumulation of surfactant at the droplet tips. Conversely, the role played by the Marangoni stresses on the droplet deformation is negligible. This has a direct implication on the use of analytic predictions and interestingly, using an effective capillary that accounts for the average surface tension reduction, the results obtained well match the ones predicted by the analytic relation of Taylor [34]. In addition, a higher amount of surfactant (higher w b ) increases the strength of surfactant-induced effects, thus the range of agreement between 2D, 3D simulations and analytic predictions slightly reduces.