**4. Conclusions**

This paper presents a series of numerical experiments in an effort to gain insight into the impact that viscoelasticity may have on drop deformation, under a number of multiphase flow configurations. The numerical method employed is based on an FEM-based spatial discretization, using a semi-Lagrangian approach to deal with the convective terms, a kinetic modeling of the polymer contribution to the stress tensor, and ML-inspired techniques for building, over the whole domain, each of the components of the extra-stress tensor, which effectively decouples the "microscopic" and "macroscopic" scales. Hence, the process of mesh refinement is no longer hampered by the low number of polymer particles at a certain computational cell, allowing us to use very refined meshes and achieve excellent smoothness and accuracy in kinetic-based, complex flow simulations.

The results on drop deformation under shear flow agree extremely well with those in the literature produced by an equivalent constitutive (Oldroyd-B) model and indicate that drop viscoelasticity prevents deformation to some degree, underscoring the efficiency and viability of the stochastic approach in multiphase flows. For more demanding situations in buoyancy-driven flow, we obtain remarkable mesh (macro) and ensemble (micro) convergence for both the FENE and Hooke models; as for the degree of CSRBFs smoothness required to accurately build the "extra-"stress tensor, the numerical study points to a minimum degree of *C*<sup>2</sup> smoothness guaranteed by *ϕ*3,1, with higher smoothness providing small benefits compared to the additional computational cost. Finally, we observe dramatic changes of flow pattern, including viscoelastic ("negative-wake") effects and extremely large values of shear and normal stress differences close to the drop interface, as the polymer concentration and relaxation times increase.

In future investigations of multiphase flow of polymeric liquids, we will try to address improved resolutions and problem configurations taking advantage of this ML-enhanced method to combine isotropic and anisotropic mesh adaptation [77] with a kinetic modeling approach.

**Funding:** This research was funded by "Ministerio de Ciencia, Innovación y Universidades" under Grant PGC2018-097565-B-I00.

**Conflicts of Interest:** The author declares no conflict of interest.

**Abbreviations** The following abbreviations are used in this manuscript:

