2.1.1. Governing Equations

The two-phase, incompressible fluid flow studied in this work is governed, at the "macro"-scale, by the Navier–Stokes equations, which in dimensionless form can be written in the space-time domain Ω × [0, *T*] as:

$$\begin{cases} \begin{aligned} & \frac{D\boldsymbol{\varrho} \cdot \boldsymbol{D} \mathbf{u}}{Dt} - \nabla \cdot \left[ \mu \left( \nabla \mathbf{u} + (\nabla \mathbf{u})^T \right) \right] + \nabla p = -\rho \mathbf{e}\_z \frac{R\boldsymbol{e}}{Fr^2} + \frac{c}{D\boldsymbol{e}} \nabla \cdot \mathbf{r}\_p + \frac{R\boldsymbol{e}}{\mathcal{W}\boldsymbol{e}} \kappa \delta \mathbf{r} (\boldsymbol{\varrho}) \mathbf{n}, \\ & \nabla \cdot \mathbf{u} = 0; \end{aligned} \\ & \mathbf{u} \left( \mathbf{x}, 0 \right) = \mathbf{u}\_0 \left( \mathbf{x} \right) \quad \forall \mathbf{x} \in \Omega, \\ & \mathbf{u} \left( \mathbf{x}, t \right) = \mathbf{0} \quad \text{on} \quad \delta D\_{\mathrm{no}\mathrm{-}\mathrm{lip}} \subset \delta \Omega, \forall t \in (0, T), \\ & \mathbf{u} \left( \mathbf{x}\_{\mathrm{imp}}, t \right) = (\mathbf{x}\_{\mathrm{ref}}, t) \quad \text{on} \quad \delta D\_{\mathrm{PBC}} \subset \delta \Omega, \forall t \in (0, T), \\ & \mathbf{u} \left( \mathbf{x}, t \right) \cdot \mathbf{n} = 0 \quad \text{and} \quad \mathbf{n} \cdot \boldsymbol{\tau}\_s \cdot \mathbf{t} = 0 \quad \text{on} \quad \delta D\_{\mathrm{free-slip}} = \delta \Omega \ ( \ \delta D\_{\mathrm{no-slip}} \cup \delta D\_{\mathrm{pbc}} ), \forall t \in (0, T). \end{aligned} \tag{1}$$

where *ρ* is the density, *μ* the viscosity, *D*()/*Dt* the total (convective) derivative operator, *p* the pressure, **v** the velocity, *<sup>τ</sup>p* the polymer ("extra-")stress tensor, Γ the interface of the drop, *κ* its curvature, *φ* the level set function (see Section 2.1.2), *Re* the Reynolds number, *Fr* the Froude number, *We* the Weber number, *c* the concentration parameter, *De* the Deborah number, *T* the final time, and *δ*Ω the boundary of the spatial domain Ω. The "no-slip", "free-slip", and "Periodic Boundary Conditions" (PBC), are implemented in Section 2.2.3 for the configurations explored in the Results and Discussion Section 3. Using a semi-Lagrangian approach, the total convective operator in Equation (1) is discretized along the "characteristic curves" of the flow without changing the Eulerian, underlying mesh. The resulting weak formulation can be written as a Stokes-like problem to be solved at each instant of time, with the following form:

$$\begin{cases} \frac{3Re}{2dt} \left( \boldsymbol{\rho}^\* \left( \boldsymbol{\phi}^n\_h \right) \mathbf{u}^n\_{h^\*} \boldsymbol{\psi}\_h \right) + \left( \boldsymbol{\mu}^\* \left( \boldsymbol{\phi}^n\_h \right) \nabla \mathbf{u}^n\_{h^\*} , \nabla \boldsymbol{\psi}\_h \right) - \left( p^n\_{h^\*} \nabla \cdot \boldsymbol{\psi}\_h \right) = \frac{2\mathcal{R}\varepsilon}{dt} \left( \boldsymbol{\rho}^\* \left( \boldsymbol{\phi}^n\_h \right) \mathbf{u}^{n-1}\_h , \boldsymbol{\psi}\_h \right) - \left( \boldsymbol{\rho}^\* \left( \boldsymbol{\rho}^\* \right) \mathbf{u}^n\_{h^\*} \right) \\ \quad - \frac{\mathcal{R}\varepsilon}{2dt} \left( \boldsymbol{\rho}^\* \left( \boldsymbol{\phi}^n\_h \right) \mathbf{u}^{n-2}\_h , \boldsymbol{\psi}\_h \right) - \frac{\mathcal{R}\varepsilon}{Fr^2} \left( \boldsymbol{\rho}^\* \left( \boldsymbol{\phi}^n\_h \right) \mathbf{e}\_z , \boldsymbol{\psi}\_h \right) + \\ \quad + \frac{c}{D\varepsilon} \left( \nabla \cdot \boldsymbol{\pi}^n\_{ph}, \boldsymbol{\psi}\_h \right) + \frac{\mathcal{R}\varepsilon}{\mathcal{W}\varepsilon} \left( \boldsymbol{\kappa}^n\_h \boldsymbol{\delta}\_{\Gamma\_{h/2}} \left( \boldsymbol{\phi}^n\_h \right) \mathbf{n}^n\_h, \boldsymbol{\psi}\_h \right), \forall \boldsymbol{\psi}\_h \in \mathbf{V}\_{h0}; \\ \quad \left( \nabla \cdot \boldsymbol{\operatorname{u}}^n\_h, q\_h \right) = 0, \; \forall q\_h \in Q\_h; \qquad \text{with } : (a, b) \equiv \int\_D \boldsymbol{ab} \, d\mathbf{x} \end{cases} \tag{2}$$

where the subscript *h* denotes the spatial discretization of the variables; the superscripts ∗ and *n* their non-dimensionality and temporal discretization at instant of time *n*, respectively; *dt* is the time step size used in the discretization of the full interval (0, *<sup>T</sup>*); *ψh* the basis functions of velocity vector space **V***h* and *qh* the basis functions of the pressure space *Qh*. The basis functions of the finite element spaces are chosen to be polynomials of degree {<sup>P</sup>2, P1, P2, <sup>P</sup>1} for velocity, pressure, the level set function, and the polymer stress tensor, respectively, thus satisfying the Ladyzhenskaya–Babuska–Brezzi condition. ˘

The spatial discretization of the domain Ω is accomplished using an unstructured mesh composed of simplices of average size *h*, and the surface tension effects are implemented making use of the Laplace–Beltrami operator, as detailed in [37].

In Equation (2), we collect the "macro-scale" physics involved in the FEM-based discretization of a multiphase flow of incompressible fluids subject to gravity, viscous, and surface tension forces, in which the interface is retrieved by a level set function *φ*, and the viscoelastic effects of the polymer solution are captured through the polymer stress tensor *<sup>τ</sup>p*. We now take a stochastic approach to the computation of *<sup>τ</sup>p* avoiding closed-form constitutive equations, performing instead Brownian dynamics simulations of polymer particles [7] modeled via the Hooke (equivalent to the Oldroyd-B constitutive equation) and FENE (Finitely Extensible Non-linear Elastic) kinetic models [51]. To compute *<sup>τ</sup>p*, we need to average over the internal configurations **Q** of each of the *Nd* dumbbells ("polymer particles"), which are scattered over the domain and follow the flow. Thus, we solve for each of the dumbbells the following stochastic (dimensionless) differential equation:

$$d\mathbf{Q} = \begin{cases} \left(\kappa \cdot \mathbf{Q} - \frac{1}{2D\varepsilon} \mathbf{Q}\right) dt + \frac{1}{\sqrt{D\varepsilon}} d\mathbf{W}, & \text{Hoke;}\\ \left(\kappa \cdot \mathbf{Q} - \frac{1}{2D\varepsilon} \frac{\mathbf{Q}}{1 - ||\mathbf{Q}||^{2}/b}\right) dt + \frac{1}{\sqrt{D\varepsilon}} d\mathbf{W}, & \text{FENE;} \end{cases} \tag{3}$$

a task that is accomplished by means of a (weak) second-order accurate predictor-corrector algorithm proposed by Öttinger [51]. In Equation (3), *κ* is the velocity gradient, **W** a stochastic Wiener process, the Deborah number *De* representing the ratio between the relaxation time of the polymer and a characteristic time of the flow, and *b* is the FENE extensibility parameter so that, as *b* → <sup>∞</sup>, Hooke →FENE. A variance-reduction technique can also be implemented [37,41], performing the average over the *Nd* dumbbells with uncorrelated noise that comprise each of the *Nens* ensembles, so that the *i*-th dumbbell of each ensemble is subject to the same (correlated) stochastic "kick", reducing the overall noise of the simulation. Finally, the polymer stress tensor is reconstructed and evaluated at the mesh nodes [52] using the ML-inspired techniques of Section 2.2.1.

### 2.1.2. Interface Capturing Technique

The interface Γ is captured as the zero isocontour of a level set function *φ*, which is transported by the characteristic curves of the flow through the total derivative operator:

$$\frac{D\phi}{Dt} = \frac{\partial\phi}{\partial t} + \mathbf{v} \cdot \nabla\phi = 0. \tag{4}$$

As we did with Equation (1), the convective terms in Equation (4) are dealt with using a semi-Lagrangian approach, as proposed first in [53]. The level set function *φ* is initialized as a signed-distance function to improve the regularity and behavior of the transported solution. However, as time goes on, discretization errors build up in the form of numerical artifacts that may eventually propagate into the interface and destroy it. To prevent this undesired effect, an Eikonal-based, redistancing procedure [54] is performed at each time step, so that we ensure that, within a band around the interface, the level set function preserves the signed-distance property. To improve mass conservation and shape definition, we add *Nmp* marker particles [18] that are passively advected by the flow, helping to correct the interface by defining local level set functions with a variable radius 0 ≤ *r*min ≤ *r<sup>n</sup> p* ≤ *r*max < *h*, accounting for lost resolution in the sub-grid scales (see [37] for details).
