**1. Introduction**

Bubble and drop dynamics in non-Newtonian fluids are a topic of undeniable interest within the community [1], owing largely to the number of real-world situations that may benefit from a comprehensive knowledge of the underlying physics: from drop formation mechanisms [2], to biomedical equipment involving droplet manipulation [3] or engineering devices in which breakup plays a central role [4,5], from droplet impact on liquid surfaces [6] to the study of drop dynamics within polymer gels and solutions [7,8]; a deeper understanding of this type of multiphase flows not only would improve existing manufacturing processes, but also encourage the development of new applications and spur breakthroughs in scientific research.

To study the multiphase flow of polymeric liquids, one should choose an appropriate discretization method capable of providing an accurate description of the interface. Mesh-based methods [9] offer such representation of the interface, either in an explicit (interface-tracking) or in an implicit (interface-capturing) form. Among the former, front-tracking and Arbitrary Lagrangian-Eulerian (ALE) schemes [10,11] display excellent performance in terms of mass conservation and shape preservation; however, remapping techniques are often found necessary [12] under extreme deformation of the moving interface. In contrast, interface-capturing methods follow a Eulerian approach, with the 'Volume-Of-Fluid' (VOF) method positioning itself as one of the most popular techniques within this context [13,14], showing good conservation properties, but requiring additional tools for handling geometrical quantities derived from the interface [15]. As an alternative, Level Set (LS) methods [16,17] capture the interface as the zero isocontour of a certain scalar function, which is advected by the flow, with noticeable mass loss and shape degradation if excessive diffusion is introduced during the advection stage. Despite these shortcomings, the LS method is widely used for interface problems undergoing dramatic deformation and topological changes and can be readily enhanced via "hybrid" schemes such as the Particle Level Set (PLS) method [18] that are able to ease, to a large extent, many of its drawbacks. Whatever the multiphase technique, the correct implementation of appropriate Boundary Conditions (BC) for the configuration of interest is a topic that requires careful consideration of the chosen discretization method [19], as it may have an enormous influence on the actual shape of the interface.

A common approach to handling the polymer interaction with the flow at a macroscopic scale is to represent the polymer contribution to the stress tensor by means of a closed-form, "constitutive" equation; thus, e.g., the work of Yue et al. [20–22] on drop deformation and complex two-phase flow using a diffuse-interface method and constitutive modeling; in [23], Pillapakkam employed an LS method to study rising bubbles in viscoelastic media, while Foteinopoulou and Laso [24] used a Phan–Thien-–Tanner model together with an elliptic mesh-deformation algorithm to investigate bubble oscillation; Castillo et al. proposed an LS method with a pressure-enriched FE space to study the two-fluid flow problem along with a Giesekus model for the polymeric liquid [25]; Fraggedakis et al. [26] characterized the critical volume of a bubble rising in a viscoelastic fluid using an FEM-based method and the exponential Phan–Thien and Tanner model; using a coupled LS-VOF ("VOSET") method, Wang et al. [27] studied drag reduction in cavity flow; Xie et al. [28] focused on droplet oscillation under a Maxwell model using lattice Boltzmann techniques. In contrast to constitutive modeling, the "micro-macro" approach [29] tackles the polymer-flow interaction using stochastic and Brownian Dynamics (BD) simulations [30–33] to retrieve the polymer stress tensor from the internal configurations of the polymer particles advected by the flow. Taking the CONNFFESSITapproach of Laso and Öttinger [34], Cormenzana and co-workers [35] and later Grande et al. [36] successfully handled free surface flows of polymer solutions, while Prieto [7,37] conducted multiphase simulations in viscoelastic fluids using a variance-reduced, stochastic implementation of a "micro-macro" method [38]. Further study of multiphase non-Newtonian flows was carried out by Bajaj et al. [39] and Xu et al. [40] using the Brownian Configuration Field (BCF) method of Hulsen et al. [41].

At present, there is a growing interest in Machine Learning (ML) [42,43] in the context of polymer simulation: Doblies et al. [44] employed ML as a means of predicting mechanical properties, while Jackson and collaborators [45] focused on the optoelectronic properties of conjugated polymers. It was also used by Kopal and co-workers in the prediction of the viscoelastic behavior of elastomer systems [46] using Radial Basis Functions (RBFs) or as a data-driven classification method to determine polymer/solvent compatibility by Chandrasekaran et al. [47]. Fluid dynamics of multiphase flow also benefits from this recent focus on ML: Ma and collaborators used statistical learning in [48] for bubbly systems; Ladický et al. accelerated multiphase simulations in GPUs using a data-driven approach with regression forests [49]; and Gibou and co-workers employed deep learning techniques with sharp interface methods in [50]. Hence, it seems unquestionable that ML and data-driven simulations are already having a tremendous impact within the scientific community, and the future possibilities seem nearly endless.

The main purpose of this paper is to highlight the effects caused by the non-Newtonian behavior of the polymer solution in a multiphase flow system, performing a series of numerical simulations by means of a computationally efficient, accurate, and robust numerical method based on a finite element discretization of the governing equations that uses ML-inspired techniques for the reconstruction of the polymer stress tensor, a feature that proves to be of the greatest importance for the accurate characterization of the viscoelastic flow. Thus, after presenting this Introduction, we describe the governing equations along with the aspects of the computational implementation in the Materials and Methods Section 2. Then, we move to Section 3, where we present results for drop deformation under steady, shear flows and in situations where gravitational effects drive the dynamics of the flow. Finally, Section 4 offers some conclusions and future lines of work.

### **2. Materials and Methods**

In this section, we describe the main ideas, mathematical background, and computational implementation of the FEM-based, ML-enhanced method used to perform the series of numerical experiments carried out in Section 3 to gain insight into the impact that viscoelasticity may have on drop deformation. We start with the finite element discretization, then moving on to the computational implementation.

### *2.1. Finite Element Discretization*

The problem of drop deformation in polymer solutions can be tackled using different methods. In this work, we use an FEM-based discretization of the "macro-scale" equations governing the fluid flow, employing a stochastic modeling of the stresses arising from the polymer fluid to account for the "micro"-scale. Finally, we use a Particle Level Set (PLS) method to capture the interface of the deforming drop.
