2.2.3. Boundary Conditions

The viscous and viscoelastic flows presented here for the analysis of drop deformation in polymer solutions make use of several types of boundary conditions: the "no-slip" boundary condition to ensure that the fluid adheres to the adjacent solid boundary; the "free-slip" boundary condition to neglect friction forces at the fluid-solid interface; and Periodic Boundary Conditions (PBCs) at the inlet and outlet of the domain to model sufficiently large domains.

In our FEM-based method, the "free-slip" condition can be applied naturally as a Neumann-type condition that is satisfied automatically by the weak formulation of the problem [7,37]. However, the implementation of PBCs in unstructured meshes for the configuration presented in Section 3.1 for drop deformation in shear flows is not straightforward. Following Pask et al. [71] and Sukumar and Pask [72], we carry out an efficient implementation of the PBCs in structured and unstructured finite element meshes using the idea of "reference" and "image" nodes at the inlet and outlet, respectively (see Figure 1). First, we obtain the connectivity lists for the assembled sub-matrices K\_00,K\_01,K\_10 of the block matrix K = [K00 K01; K10 K11]; then, we perform row operations (first) and column operations (later) to K00,K01,K10 and to the right-hand side F, as described in [72]; finally, we update the values of the (pressure,velocity) "image" nodes with the values obtained from the corresponding "reference" nodes. The addition of dumbbells and marker particles to deal with multiphase flows and non-Newtonian fluids adds complexity to the actual implementation, but not to the general concept.

As for the "no-slip" conditions, they translate into homogeneous or non-homogeneous Dirichlet-type conditions that can be applied efficiently using the symmetry-conserving procedure MatZeroRowsColumns to perform row and column operations on the PETSc-based assembled matrices K00,K01,K10 and right-hand side F.

**Figure 1.** Schematic of the computational domain for a drop deforming in shear flow.

### **3. Results and Discussion**

The main purpose of this paper is to explore the behavior of drops immersed in polymer solutions through a series of numerical experiments that try to highlight the effects of viscoelasticity in the shape of the drop, as well as in the emerging flow patterns. First, we focus on drop deformation under a steady, shear flow, offering results under two configurations: Newtonian drop in a Newtonian "matrix" (ambient fluid) and non-Newtonian drop in a Newtonian matrix. Then, we make use of our ML-enhanced, FEM-based method to investigate stronger viscoelasticity effects in flows where gravitational forces and surface tension effects can become relevant, showing the influence that the discretization employed and the numerical techniques utilized may have on the accurate reconstruction of the polymer stress tensor and on the quality of the discrete solution for the multiphase flow.

### *3.1. Drop Deformation in Steady, Shear Flow*

We start the exploration of viscoelastic effects on drop deformation by considering the problem of a drop immersed in another, ambient fluid and subject to a steady, shear flow. First, we show results for a Newtonian/Newtonian system; then, we introduce viscoelasticity in the system, comparing our results with others obtained using a constitutive (Oldroyd-B) model, equivalent to our stochastic, Hookean dumbbell approach. The configuration is represented in Figure 1: a drop of a viscous (real) fluid of radius *a* placed at the center of a domain [2*L* × <sup>2</sup>*<sup>H</sup>*], with *H* = 4*a* and *L* = 8*<sup>a</sup>*, experiences a steady, shear flow of rate *γ*˙ = *V*/*H* produced by the top and bottom lids moving at velocity *V* in opposite directions.

The inlet and outlet walls have Periodic Boundary Conditions (PBCs), the implementation of which is detailed in Section 2.2.3. At the top and bottom panels, the "no-slip" boundary condition is considered.

### 3.1.1. Newtonian Drop in a Newtonian Matrix

We place a Newtonian drop in a Newtonian ambient fluid as depicted in Figure 1 to study the shape of the interface for different values of the capillary number *Ca* = *We*/*Re* in terms of the "deformation parameter" *D*. This parameter is usually defined as *D* = (*L* − *B*) / (*L* + *<sup>B</sup>*), with *L* and *B* being the longest and shortest lengths from the center of the drop to the interface (corresponding also to the major and minor axes of the ellipse), respectively. We also choose a small value for the Reynolds number (*Re* = 0.1) so as to neglect strong inertial effects, comparing our results with those of Zhou and Pozrikidis [73], Yue et al. [20,22], and Afkhami et al. [74]. We employ an unstructured, uniform mesh of size *h* = 1/80, in a rectangular domain [0, 2] × [0, 1] with time step *dt* = 1/100 and *Nmp* = 2.5 × 10<sup>5</sup> marker particles. In Figure 2, we represent the steady value of the deformation parameter for increasing values of the capillarity number; in all cases, we obtain excellent agreemen<sup>t</sup> with the published results, especially with that by Afkhami et al. [74], deviating in the worst case scenario < 2.5% at *Ca* = 0.6 from the results offered by Zhou and Pozrikidis in [73]. The steady values of *D* show a linear behavior with *Ca* for small capillary numbers, then decreasing as *Ca* gets larger; as pointed out in [20], we expect the physical insights gained by this 2D exploration to be relevant to actual (three-dimensional) experiments.

**Figure 2.** Comparison of the steady-state deformation of a Newtonian drop in a Newtonian matrix between our results and those found in Zhou & Pozrikidis [73], Yue et al. [20,22], and Afkhami et al. [74].

We also collect in Figure 3 the evolution of the deformation parameter *D* as a function of the dimensionless time *t* ∗ = *tγ*˙ and compare it with the results of [74], for a wide range of capillary numbers from *Ca* = 0.1 to *Ca* = 0.6. Again, an excellent agreemen<sup>t</sup> is obtained throughout, with a maximum difference of less than 1.3% (for *Ca* = 0.6) between the two numerical methods at any given time.

**Figure 3.** Comparison of the evolution of deformation for a Newtonian drop in a Newtonian matrix of equal viscosity, for different values of capillarity *Ca*, between our results (solid line) and those found in Afkhami et al. [74] (dashed line).

Following Yue et al. [22], we now study the steady shape of the drop in polar coordinates, collecting the radius *r* for each orientation angle *θ* ∈ [0, *<sup>π</sup>*], with *θ* defined as the angle formed by the line connecting the center of gravity of the drop and a point in the free-surface, with the *x*-axis. Thus, we show in Figure 4 the shape in polar coordinates and the actual shape for a simulation with *Ca* = 0.1.

**Figure 4.** Comparison of the drop shape in (**a**) polar coordinates and (**b**) the actual shape, for *Ca* = 0.1, between our results (triangles) and those found in Yue et al. [22] (square markers).

We can observe in this Figure 4 very good agreemen<sup>t</sup> in the orientation angle of the droplet, which can be measured at ˆ *θ* ≈ 39.6◦, lower than 45◦ as numerical and experimental results confirm. However, there is a major discrepancy between the results of the two methods: the radii for each polar angle are smaller for the method of Yue et al. [22] when compared with ours (see left panel). To better understand what is happening here, we plot the actual shape of the drop according to both methods in the right panel of Figure 4. As can be observed, the reason for the discrepancy lies in the lack of mass conservation "mass-loss" incurred by the drop of [22], which loses around 5% of its initial mass at the end of the simulation, whereas we can ensure mass conservation up to 99.98% of the initial mass in this simulation.

### 3.1.2. Viscoelastic Drop in a Newtonian Matrix

After exploring the multiphase Newtonian/Newtonian system and comparing our results with those of the literature, we now proceed with a problem in which the viscoelastic effects are present. We focus on a viscoelastic drop modeled by the Hookean dumbbell model (equivalent to the Oldroyd-B constitutive equation), immersed in a Newtonian, ambient fluid, using the same configuration of the previous Newtonian/Newtonian analysis. The flow is impulsively started at *t* = 0 with shear rate *γ*˙ = 1, using a rather coarse mesh with 80 × 40 elements and *Np* = 10<sup>7</sup> dumbbells uniformly placed inside the droplet; the flow is continued until dimensionless time *t* ∗ = *tγ*˙ = 10, with a small time step to accurately solve the internal configurations of the dumbbells (*dt* = 1/200), taking *Nt* = 2000 time steps to finish each simulation; the number of marker particles to improve the definition of the interface was *Nmp* = 2.5 · 105. The ratio between viscous and inertial effects in the system, represented by the Reynolds number, are of utmost importance, in the sense that, if the method is not able to deal with extremely low *Re* ("creeping" or "Stokes" flows), the inertial effects become relevant when small time steps are used, thus affecting the history of the flow and, consequently, that of the polymer particles (dumbbells); it is for this reason that the value of *Re* = 0.1, deemed appropriate in the previous section for Newtonian flows, is now replaced by *Re* = 10−<sup>5</sup> to suppress any undesired effects that inertia may have on the computation of the steady-state values for the polymer solution. The remaining dimensionless parameters are chosen as those found in [21], with *Fr* → <sup>∞</sup>, *We* = <sup>10</sup>−6, so that the Capillary number *Ca* = 0.1; our concentration parameter *c*, according to the characteristic scales chosen in [21], corresponds to *c* = 1 − *β*, with *β* = 0.5 the retardation parameter of the Oldroyd-B fluid; the Deborah numbers studied are *De* = {0.25; 0.5; 1; <sup>2</sup>}; and the density and viscosity ratios between the drop and the outer ambient fluid (matrix) are taken as *ρ*2/*ρ*1 = *μ*2/*μ*1 = 1. We performed a set of simulations to obtain the evolution of the drop deformation and compared the results found in Figure 1 of [21].

Despite the rather coarse mesh used (equivalent to a uniform size *h* = 1/40), the not-so high number of dumbbells and the totally different approach taken by the two techniques compared—the diffuse-interface method combined with a phase-field approach ruled by Cahn–Hilliard dynamics to study the Newtonian/non-Newtonian problem in a unified way of [21] and our FEM-based, stochastic method—the results presented in Figure 5 are in remarkably good agreement, especially the steady-state values of the deformation parameter *D*. However, also the transient behavior shows a noteworthy resemblance, with the overshoot appearing for sufficiently high *De* values and the evolution of *De* = 2 being for *t* - 4 higher than those for *De* = 1. In any case, we notice the effect of the drop viscoelasticity as a means to reduce the deformation of the interface; plots of the actual shape of the interface (not included here for brevity's sake) show this same trend. Apart from the stochastic noise, the modification of the interface by the correction step of the marker particles and the mass conservation step add to the oscillatory behavior of *D* observed in Figure 5, which is explicitly computed from the discrete interface. Moreover, Figure 5 underscores the elastic effects of the drop, since an increase in the Deborah number *De* produces a larger ratio between the maximum and the steady-state value of the deformation parameter *D*, as a consequence of the longer time (slower response) that the polymer molecules take to recover from the applied shear strain.

We now compare the steady-state values of the polymer stress tensor obtained using the Oldroyd-B constitutive model of [22] and the Hookean dumbbell model. The results are presented in Figure 6. As we did for the Newtonian/Newtonian case, we represent the interface also in polar coordinates with *θ* ∈ [0, *<sup>π</sup>*]. The normal component of the polymer stress tensor *<sup>τ</sup>p*,*<sup>n</sup>* ≡ **n** · *<sup>τ</sup>p* · **n** (with **n** the outer normal at each point of the interface) is computed in Figure 6a.

**Figure 5.** Comparison of the evolution of deformation for a non-Newtonian (Oldroyd-B) drop in a Newtonian matrix of equal viscosity after the abrupt start of shear flow, for increasing values of the Deborah number *De*, between our results (solid line) and those found in Yue et al. [21] (dashed line).

**Figure 6.** Comparison of the (**a**) normal and (**b**) tangential stresses along the inner edge of the interface for a non-Newtonian (Oldroyd-B) drop in a Newtonian matrix of equal viscosity, for increasing values of the Deborah number, between our results (*De* = 0.5, green squares; *De* = 1, yellow triangles; *De* = 2, inverted red triangles) and those found in Yue et al. [22].

In spite of the largely dissimilar methods, we notice a very good agreemen<sup>t</sup> between them. We observe how the effects of the viscoelastic drop are not very strong on the normal component of the polymer stress tensor, though increasing the viscoelasticity by means of *De* reduces the maximum value, which is also moved towards smaller orientation angles, thus playing a part in diminishing drop deformation. Overall, the normal polymer stress component is weakened at the poles and at the equator of the drop as viscoelasticity increases, preventing drop deformation.

Figure 6b shows the tangential component of the polymer stress tensor *<sup>τ</sup>p*,*<sup>t</sup>* ≡ **t** · *<sup>τ</sup>p* · **t**, with **t** the unit vector tangent at all points to the interface and perpendicular to **n** (so that, in 2D, *tx* = *ny*, *ty* = −*nx*). Here, the polymer stresses are much larger than for the normal component, showing again an excellent agreemen<sup>t</sup> between our results and those in [22]. The larger *De* values intensify the tangential stresses at the equator, being responsible to a larger extent for the reduced deformation observed for the higher *De*. Despite small discrepancies, even with a coarse mesh and computationally inexpensive simulation, our results compare very well with those of the reference, capturing the minimum and maximum values of *<sup>τ</sup>p*,*t*, as well as the evolution of *<sup>τ</sup>p*,*<sup>t</sup>* along the interface.

### *3.2. Drop Deformation in Buoyancy-Driven Flow*

After studying the problem of drop deformation in unsteady, shear flow, in this section, we focus on following the behavior of a Newtonian drop immersed in an ambient, polymeric solution, rising due to buoyancy effects and the density ratio between the fluids [75]. By means of a series of numerical experiments, we gain insight into the deformation of the drop shape, characterized by the circularity ✁ *c*, the ratio between the perimeter of a circle whose area is equal to that of the drop and the perimeter of the drop. We also shed light on the response of the "extra-" stress tensor *<sup>τ</sup>p* to a varying degree of polymer concentration and relaxation time, when the drop undergoes large deformation and velocity gradients, highlighting the influence of different smoothness in the CSRBF. Finally, we observe the effects that viscoelasticity may have on the flow, with emphasis on the flow streamlines during the full, unsteady simulation.

Thus, we consider a dimensionless, 2D rectangular domain [0, 1] × [0, 2] in which a drop of radius *R* = 0.25 is placed at position (0.5, 0.5) inside an unstructured mesh of uniform size *h*. The outer, non-Newtonian polymer solution is modeled using the FENE and Hookean dumbbell models. The flow is further defined by the density and viscosity ratios *ρ*2/*ρ*1, *μ*2/*μ*1 and dimensionless numbers *Re*, *We*, *Fr*, *c*, *De*. The "no-slip" boundary condition is applied at the top and at the bottom, while we consider "free-slip" at the lateral boundaries. All simulations are carried until dimensionless time *t* = 3.
