Head-on collision of Newtonian drops in a viscoelastic medium

Depending on the application, one needs to either stabilize or destabilize the interfacial properties of an emulsion. An aspect of the dynamics that governs the stability of emulsions in general is the drainage time of the film that is formed when two drops collide. In this work, we study the effect that viscoelasticity of the matrix fluid has on this drainage time of two Newtonian drops that perform a head-on collision under an applied macroscopic extensional rate. For the modeling of the viscoelastic matrix material, the Giesekus model is chosen. A cylindrical coordinate system is applied with imposed axisymmetry and the resulting equations are solved using fully resolved numerical simulations employing a finite element discretization. Our results show that viscoelasticity reduces the drainage time, which is a combined effect of three different stages.


Introduction
The dynamics of drop interaction in flow has been extensively studied both experimentally (Orme 1997;Couder et al. 2005;Willis and Orme 2003) and theoretically (Eggers et al. 1999;Janssen and Anderson 2011) due to their vital importance in a diverse spectrum of multiphase systems. In some cases, for example in food and pharmaceutical industries, it is desired that the emulsion is stable for extended periods of time; i.e., the size distribution of drops should not change over time due to flow (Kogan and Garti 2006). However, in other applications, for instance in oil extraction, an unstable emulsion is preferred since it makes it easier to separate the drop phase from the suspending medium.
Active and passive manipulation of non-Brownian particles in fluids has been a topic of interest in microfluidics for several years now. The recent review by Lu et al. discusses particle manipulations in non-Newtonian microfluidics for various passive manipulations, including focusing, separation, washing and stretching of particles (Lu et al. 2017). Despite practical advantages of passive manipulation of particles (Kang et al. 2009), for several applications more control is required and an active manipulation is desired. Researchers have used different actuation mechanisms to displace and transport rigid and deformable particles like electrostatics (den Toonder et al. 2008;Khatavkar et al. 2007), magnetics (Zhang et al. 2019;Wang et al. 2013Wang et al. , 2015Wang et al. , 2016 and light (van Oosten et al. 2009).
Similar to industrial processes, also microfluidic multiphase systems have complex fluids as their components (Utracki 1989). The collision of drops is a delicate process which can be influenced by many different factors, such as the flow field type, the position and shape of the drops, the interfacial and viscoelastic properties of the phases. Despite all the previous work done, both experimentally and theoretically, the coalescence of drops in viscoelastic materials is still not very well understood.
Investigating the head-on collision of two drops can provide information, such as the drainage time, which is the basis of coalescence. Changes in the microstructure of blends, induced by coalescence and breakup, influence the size distribution of the drops and hence determine the final properties (Tucker and Moldenaers 2002). Information obtained from the evolution of the microstructure can This article is part of the topical collection "Particle motion in non-Newtonian microfluidics" guest edited by Xiangchun Xuan and Gaetano D'Avino.
1 3 87 Page 2 of 11 be used to derive models which characterize the droplet size and distribution macroscopically (Ferkl et al. 2016).
In this article, we investigate the effect of the viscoelasticity on the drainage time using numerical simulations in a regime relevant for microfluidics. A cylindrical coordinate system is used with imposed axisymmetry. For the modeling of the viscoelastic matrix material, the Giesekus model is used (Giesekus 1982). The Newtonian drops are suspended in the viscoelastic matrix and perform a headon collision under a macroscopically applied biaxial extensional flow. The problem is discretized using the finite element method.
The structure of the article is as follows. In Sect. 2, we define the problem and present the mathematical model that will be used for the balance and constitutive equations. In Sect. 3, the employed finite element method is presented. In Sect. 4, we show the validation of our method together with the obtained results. Finally, in Sect. 5, our main conclusions are summarized.

Problem definition
For the purpose of this work, two isolated drops having the same radius perform a head-on collision under biaxial extensional flow. The drops consist of a Newtonian fluid, whereas the matrix fluid is viscoelastic. To compute the viscoelastic stress, the Giesekus model is used. To reduce the computational cost, axisymmetry of the problem is assumed. The initial geometry of the colliding spherical drops is depicted in Fig. 1.

Governing equations
Since we are considering drops in a microfluidic environment, we assume that inertia can be neglected and that the fluid is incompressible. As a result, the momentum and the mass balance for both the drop and the matrix reduce to where u is the velocity vector and the stress tensor. The stress tensor is written as: where p is the pressure, I is the identity tensor and is the extra stress tensor. For the Newtonian drop, the extra stress tensor is given by with d the viscosity of the fluid in the drops and D = ( u + ( u) T )∕2 the symmetric part of the velocity gradient tensor. For the viscoelastic matrix fluid, there is a contribution from the Newtonian (solvent) part and a viscoelastic contribution p to the extra stress tensor: where s is the solvent viscosity. The viscoelastic extra stress p is governed by the Giesekus constitutive model, given by where c is the conformation tensor, is the relaxation time, is the mobility and ▽ () is the upper convected derivative: with D()∕Dt the material derivative. The polymer stress is given by herein the polymer modulus is G = p ∕ , where p is the polymeric viscosity. The zero-shear-rate viscosity of the matrix fluid is 0 = s + p .

Boundary and initial conditions
In this section, we define the boundary and initial conditions. It is assumed that the drop interface can be described by a sharp interface. This leads to a jump of the stress tensor at the boundary between the drops and the matrix, which can be expressed as: where the subscripts m and d indicate that the term is evaluated on the matrix and drop side, respectively, is the  On the outer boundaries Ω 1 , Ω 2 and Ω 3 (see Fig. 2), a macroscopic biaxial extensional velocity field is imposed, which in the axisymmetric case is given by where x is the position vector with (z, r) coordinates and ̇ is the extensional rate. Due to axisymmetry, symmetry conditions on the z-axis at r = 0 are needed. This yields the following equation: Boundary conditions at the inflow, i.e., boundaries Ω 1 and Ω 3 , for the conformation are needed. Assuming that the boundaries are sufficiently far away from the stagnation point and the drops, we compute the conformation tensor for a purely biaxial extensional flow without the presence of any drops and prescribe it as follows: where c e (ṫ) is the value of the conformation tensor for biaxial extensional flow without drops at time ṫ . Initial conditions for the conformation tensor are needed also, for which a stress-free state will be assumed: c(ṫ= 0) = I.

Interface tracking
Since we are dealing with sharp interfaces, it is needed to formulate appropriate equations that can track the motion of the interface. An interface can be described by a moving curvilinear coordinate system given by where represents the curvilinear coordinate and x is the mapping function that converts the curvilinear coordinates to spatial coordinates x . The velocity of the interface is given by and follows the material velocity u i.e., the interface moves in a Lagrangian way.

Numerical description
To solve the mathematical model described in Sect. 2, the finite element method is used. The equations are solved on moving boundary-fitted meshes, where the movement of the nodes is coupled with the flow problem.

Weak form
The weak form can be obtained by multiplying Eqs. (1) and (2) with test functions v, q : Find u, p such that in both the Newtonian drops and the matrix fluid. Using partial integration, Gauss' theorem, the interface condition Eq. (9) and the continuity of velocity Eq. (11), we obtain the following weak form: Find u and p such that using appropriate spaces for u , p, v and q, where S = S 1 ∪ S 2 . Now valid in the whole domain where in the Newtonian drops is given by Eq. (4) and in the matrix fluid by Eq. (5). During a single time step, the momentum and mass balance equations are solved first, where the polymer stress p is obtained by using a semi-implicit scheme as proposed in D'Avino and Hulsen (2010). With the updated velocity, the evolution equation for the conformation tensor is then solved. The mesh consists of both constant and moving boundaries. That means that the fluid is described on a mobile grid, which moves in a non-Lagrangian way. Hence, it is essential to use the arbitrary Lagrange-Euler (ALE) formulation, where the material derivative is rewritten as where the partial derivative with respect to time t is calculated at constant mesh coordinate x m and u m is the mesh velocity. For stabilization, we use the SUPG technique (Brooks and Hughes 1982), and together with the log-conformation approach (Hulsen et al. 2005;Fattal and Kupferman 2004), the equation for the evolution of the conformation tensor Eq. (6) in its weak form with test function d becomes: Find s such that where s = log c , is the SUPG parameter and u m is the mesh velocity. Further information about the log-conformation formulation and the expression for the function g can be found in Hulsen et al. (2005).

Discretization
The weak form is discretized using the finite element method employing a mesh of quadratic triangles. Quadratic interpolation (P 2 ) for the velocity and position and linear interpolation (P 1 ) for pressure and conformation are used. An additional Poisson equation is used for the movement of the mesh nodes which in its weak form using test function w reads: Find x such that where n is the gradient using the position at the beginning of each time step as a reference. The tracking of the interface can be done by solving two separate problems, usually by first predicting the position of the interface in the new (( n w) T , n x) = 0 for all w, time step and then solving the momentum and mass balance problem (Villone et al. 2014). In a decoupled approach, the time step would be limited by the mesh capillary time (Courant et al. 1967). That is why we choose to solve Eqs. (21), (22), (25), coupled with the Poisson problem which has the interface tracking Eq. (18) as a boundary condition. The implementation is similar as in Mitrias et al. (2017), but now the whole domain is perturbed. In this way, our scheme is not limited by the characteristic mesh size time limit.
Regarding the time discretization, a semi-implicit stress formulation is implemented. That is not necessary for the case that will be studied in this work since a nonzero solvent viscosity is present, but the model has been derived such that s = 0 could be chosen. Thus, Eq. (6) using an explicit Euler scheme similar to D'Avino and Hulsen (2010) a second-order time integration (BDF2) is used (Ascher and Petzold 1998). Furthermore, the mesh velocity u m is computed by a BDF2 formula and updated during iterations. Although a first-order scheme is used for the viscoelastic stress tensor in the momentum balance, the actual error is still O( t 2 ) . After u n+1 , p n+1 , x n+1 have been obtained, s n+1 is computed from the discretized Eq. (24) using BDF2 and thus we can calculate n+1 p = G(exp s n+1 − I). When the moving mesh becomes distorted, remeshing is performed. All the variables that need time integration, i.e., the position of the interfaces and conformation field at previous time steps, are projected onto the new mesh (Hu et al. 2001). To be able to compute the mesh velocity, the mesh coordinates at previous time steps are projected as well (Jaensson et al. 2015). The mesh is generated using Gmsh (Geuzaine and Remacle 2009).
The components of the conformation tensor form multiple systems of equations with the same matrix but a varying right-hand side. We choose to solve these systems using a direct solver provided by the HSL library (HSL 2013), which allows for the reuse of the LU-decomposition (Lay and Stade 2005). The system of equations formed by the momentum balance, mass balance and Poisson equations can grow fast due to the refinement of the approaching drop interfaces. This system is solved using a GMRES iterative solver from the Sparskit library (Saad 2001), with a modified preconditioner as explained in Jaensson et al. (2018). While the mesh is moving, the connectivity of the nodes does not change until remeshing is needed and the structure of the system matrix remains the same. Thus, it is possible to reuse the same preconditioner for multiple time steps, which significantly reduces the computational time. To reduce the size of the preconditioner, the graph that represents the system matrix is renumbered using MeTiS (Karypis and Kumar 1998).

Results
In the following sections, all the results will be presented in dimensionless form. There are several dimensionless groups that characterize different aspects of our problem, which are where Wi is the Weissenberg number, Ca is the capillary number, D is the normalized initial drop interface distance d 0 with the drop diameter 2R, is the ratio of the drop viscosity and the zero-shear viscosity of the viscoelastic fluid, is the ratio of the solvent viscosity and zero-shear viscosity of the viscoelastic fluid and is the mobility parameter of the Giesekus model. For all the simulations that consist of a viscoelastic matrix fluid, the viscosity ratio = 0.5 , the mobility parameter = 0.2 and the viscosity ratio = 1 will be fixed. The initial distance of the drop interfaces D = 3 is chosen such that there is enough time to build up viscoelastic stresses in the layer of fluid in between the drops. Note that these stresses require a scaled time of at least ṫ= Wi to develop. If the drops are initially placed close to each other, the role of viscoelastic effects would be insignificant. The size of the domain is 20 times the radius R both in axial and radial directions, so that the imposed boundary conditions for the conformation tensor do not affect the dynamics close to the drops. The schematic in Fig. 3 shows the interface distance at the center h cent and the minimum interface distance h min of the drops. When studying the film drainage, these are the important parameters that govern the drainage dynamics.

Convergence test
In order to verify our numerical model, spatial and time convergence is investigated. For the constant parameters, we choose Ca  Table 1. The relative error is defined as where h ref cent (ṫ) is the reference interface distance at time ṫ and h cent (ṫ) the interface distance for the corresponding n e at time ṫ . In Fig. 4, the error e r at ṫ= 0.01 of the interface distance at the center of the drops for the cases of n e = 40, 45, 60 and 75 which is relative to the one with n e = 500 is plotted. It can be seen that the error converges giving a maximum relative error of 10 −7 . We now chose the number of elements on the drop equator to be n e = 60 , and we perform simulations with a varying time step. In Fig. 5, the relative error at ṫ= 0.1 of time steps ṫ= 10 −3 , 5 × 10 −4 , 2.5 × 10 −4 and 10 −4 relative to the case of a time step ṫ= 10 −6 is plotted. The numerical method shows convergence, and in the remainder of this article, we are going to use n e = 60 and ṫ= 10 −3 , which is a compromise between accuracy and numerical efficiency.
To be able to correctly capture the flow in the thin layer that forms between the drop interfaces, the drop interfaces are refined so that there are enough elements in the fluid layer. In Figs. 6, 7, 8, the dimensionless first-normal stress difference of the polymer stress, N 1 ∕̃ , where N 1 = p,11 − p,22 , is shown using 4, 8 and 16 elements between the drops. Here, ̃= (̇e ff )̇e ff is the total shear stress, where (̇e ff ) is the steady-state viscosity function of the viscoelastic model. Furthermore, the effective strain rate may be expressed as the magnitude of the rate of deformation tensor ̇e ff = √ 2D ∶ D . Eight elements are shown to be enough to accurately describe the polymer stress in the thin fluid layer.

Validation
We commence by validating our code by comparing to results obtained using the boundary integral method (BIM) of Janssen  et al. (2006,2008,2011), where both the matrix and the drop consist of a Newtonian fluid. The parameters in this case are Ca = 0.05 , = 1 and D = 0.5. In Fig. 9, it can be seen that there is a good agreement between our model and that of Janssen et al. (2006). The dimple shape as depicted in Fig. 3 is observed only at the later stages of the process. From this comparison, we can also conclude that our computational domain is sufficiently large, because of the good agreement with the boundary integral method which assumes a velocity field imposed at infinity.

Effect of viscoelasticity on the film drainage
Although there is limited literature regarding the film drainage during head-on collision of drops suspended in a viscoelastic matrix, it has been reported that the viscoelasticity promotes coalescence; i.e., it speeds up the film drainage (Yue et al. 2005). Yue et al. used a two-dimensional implementation of the diffuse-interface formulation to study the effect of viscoelasticity on the film drainage. As they mentioned, the observed speedup is a result of two stages with different effects. Initially, the viscoelastic stresses are negligible and thus the drops are allowed to move faster than the Newtonian case. In the second stage, the viscoelastic stresses are fully developed and the movement of the drops slows down. The first stage dominates the dynamics, and it is the one that is responsible for the speedup.
In our simulations, apart from the first stage that contributes to the speedup, we also observed that the complex nature of the flow between the two drops can lead to faster drainage depending on the Ca number. Noticeably, a third stage that contributes to the speedup can be identified, when the drop interfaces are close to each other. At that moment, the local strain rates are minor and thus the viscoelastic stresses of the fluid inside the film layer are small. To better understand the dynamics during this process, in Fig. 10 we plot the viscoelastic stress magnitude √ p ∶ p for the case of Ca = 0.05 and Wi = 0.2 at different times. Initially, in Fig. 10a, it can be seen that at time ṫ= 0.3 the stresses are low, since we start from a stress-free condition. When the drops approach each other, the stresses start to build up as it can be seen in Fig. 10b at time ṫ= 1.5 . Finally, in the last stage where the drops start to slow down, the viscoelastic stresses in the thin layer between the drop interfaces relax as shown in Fig. 10c.
Thus, the dynamics between the drops are governed not only by the global Weissenberg number which is defined by the macroscopic imposed strain rates, but also by the local strain rates that the moving interfaces enforce. We therefore define a local Weissenberg number as: herein ̇c ent is the local strain rate imposed by the drop interfaces given by where u cent is the magnitude of the velocity of the approaching interfaces at the centerline. In Fig. 11, the increase in the Wi cent is plotted over time, for the case of Wi = 0.4 and a range of Ca . It can be seen that the local Weissenberg number Wi cent is increasing, as the drops are approaching each other. A maximum value is found for each Ca which increases with decreasing Ca . This is expected since stiffer drops, i.e., small Ca , deform less and thus impose higher strains on the matrix fluid. For the case of the smallest Ca = 0.005 that we studied, the maximum value of Wi cent is more than ten times larger than the macroscopically applied strain rate. Finally, at the last stage where the drops slow down, the value of Wi cent significantly decreases.
The speedup of the film drainage can be also identified by looking at the local Wi cent that the drop interfaces impose for different Wi numbers and a constant Ca . As shown in Fig. 12, for Ca = 0.1 the scaled Wi cent with the global Wi increases with increasing Wi , meaning that the drop interfaces move faster toward each other and thus enforcing higher local strains. The small bump that is observed after ṫ= 3.75 is related to the drop interfaces stop moving toward each other locally and start to divert, forming the dimple shape as shown in Fig. 3. Looking at the velocity profiles in the thin film between the drops, it can be seen that higher Wi leads to larger velocities as expected (Fig. 13).
The effect that viscoelasticity has on the drainage time is summarized in Figs. 14,15,16,17,18. For drops with a large Ca , this effect is significant, as it can be seen in Fig. 14    not observe any considerable differences as shown in Figs. 17 and 18. It is worth noting that for the cases where Ca < 0.05 , due to the increasing computational costs, we could not reach the point where the dimple shape of the interface starts to form and the interface distances might start to differ. One of the reasons why the viscoelastic effect gets smaller with decreasing Ca is the drastic increase in terms of Wi cent as shown in Fig. 11. This increase would also lead to increased viscoelastic stresses and thus contribute to the aforementioned slowdown of the film drainage for higher Wi . Similar effects have been reported in the literature by Dreher et al. (1999) where they presented experimental results for different sizes of droplets showing that larger drops reduce drainage time, whereas smaller drops had a slower film drainage than the Newtonian case.

Conclusions
We have presented direct numerical simulation of Newtonian drops in viscoelastic matrices under biaxial extensional flow using the Giesekus model and imposed axisymmetry. The model was validated for the case of a Newtonian matrix with the implementation of Janssen et al. (2006). The influence of viscoelasticity on the drainage time was investigated for a range of Wi numbers to study the effect of viscoelasticity. It was shown that viscoelasticity reduces the drainage time, which is a result of three different stages. In the first stage, the viscoelastic stresses are negligible and thus the drops are moving faster than the Newtonian case. During the second stage, the stresses start to build up, but for Ca numbers above 0.01 they were not sufficient in order to slow down the drops. On the contrary, the drops move faster with increasing Wi number. For small Ca, there was no significant effect of viscoelasticity up to the interface distances that we were able to study. Finally, in the third stage, when the drop interfaces are close to each other, the local strain rates are small and thus the viscoelastic stresses of the fluid inside the film layer are negligible. It is worth noting that for smaller Ca numbers and interface distances than the ones presented in this work, different effects might be present.
Funding The research leading to these results has received funding from the European Commission under the grant agreement number 604271 (Project acronym: MoDeNa; call identifier: FP7-NMP-2013-SMALL-7).
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creat iveco mmons .org/licen ses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.