3.2.1. Convergence Results

We turn our attention to the evolution of the drop shape during the full simulation, as the mesh size is decreased from *h* = 1/40 down to *h* = 1/320, and the number of ensembles scattered over the domain changes, using two different kinetic models, Hooke, and FENE, and different degrees of viscoelasticity, (*c* = 1, *De* = 1) and (*c* = 5, *De* = 3). The other relevant parameters in this case are: *Re* = 35, *We* = 10, *Fr* = 1, *ρ*2/*ρ*1 = 10−<sup>1</sup> = *μ*2/*μ*1.

In Figure 7, we present the results for convergence under mesh refinement, for the FENE and Hooke models: for the case with weaker viscoelastic effects (*c* = 1, *De* = 1), the number of ensembles is *Nens* = 5000, while the number of dumbbells inside each ensemble is *Nd* = 150, 000 for both stochastic models. When the polymer concentration is increased (*c* = 5, *De* = 3), we scatter more ensembles in the ambient (polymeric) fluid *Nens* = 5000 to capture the more dramatic effects in shape deformation, reducing the number of dumbbells within each ensemble to *Nd* = 15, 000 to keep the computational demands at a similar level. The results obtained here compare very well with those provided in [7] for purely uncorrelated dumbbells, showing excellent convergence with mesh refinement for both kinetic models and degrees of polymer concentration. Only in the coarser mesh, *h* = 1/40, do we observe a slightly different behavior to that produced by the other unstructured meshes. This pattern is more evident when using the Hooke instead of the FENE dumbbell model, due to the larger stresses produced by the former model during drop deformation: the minimum value of the circularity decreases for both *c* = 1 and *c* = 5, and we observe a lower local maximum at *t* - 2.2 for *c* = 5, *De* = 3, with the drop showing a lack of deformation at the latter stages of the simulation (*t* ≥ 2.75).

**Figure 7.** Evolution of the circularity of a drop rising due to buoyancy effects in a polymer solution, using the FENE and Hookean kinetic models, variable mesh size *h* = {1/40, 1/80, 1/160, 1/320}, and different degrees of viscoelasticity: *c* = 1, *De* = 1; and *c* = 5, *De* = 3.

Figure 8 shows the evolution of the circularity when the number of ensembles *Nens* and dumbbells per ensemble *Nd* changes, keeping the product *Nens* × *Nd* = 750 × 10<sup>6</sup> constant so as to maintain the computational cost; the mesh size is *h* = 1/320. For the case with lower polymer concentration and viscoelastic effects, the convergence results are remarkable for both models. When the concentration parameter and Deborah number are increased to *c* = 5, *De* = 3, we observe a small discrepancy in the circularity values, proving that a certain number of ensembles are required to reconstruct, in an accurate way, the stronger polymer stresses. However, there seems to be an optimal value of *Nens* and *Nd* for a given computational cost, with this value of *Nens* being larger for Hooke than for the FENE model. In any case, such a modest diverging behavior is only observed at or after the minimum circularity (*c* = 1) or local maximum circularity (*c* = 5) are attained.

**Figure 8.** Evolution of the circularity of a drop rising due to buoyancy effects in a polymer solution, using the FENE and Hookean kinetic models, variable number of ensembles (*Nens* × *Nd* = 750 × 106), and different degrees of viscoelasticity: *c* = 1, *De* = 1; and *c* = 5, *De* = 3.

### 3.2.2. Impact of CSRBF smoothness on the polymer stress tensor

Next, we investigate the influence that the degree of CSRBF smoothness has on the flow. In particular, we are interested in observing qualitative and quantitative deviations in the shear component and normal stress difference of the polymer stress tensor *<sup>τ</sup>p*, in a case with strong viscoelastic effects (*c* = 5, *De* = 3), when Wendland's CSRBFs, *ϕ*3,*k*, ∀*j* = 0, .., 3, are used. The reason behind our pursuing this investigation is twofold: since the "extra-"stress is responsible for the viscoelastic effects, it is mandatory to represent it as accurately as possible in simulations that aim to highlight the non-Newtonian behavior of the polymer solution; at the same time, this should be

accomplished as efficiently as possible, opting for the CSRBFs that provide results nearly as accurate as the more computationally expensive alternatives.

In the following simulations, we take an unstructured mesh of size *h* = 1/320, keeping the same values of the dimensionless parameters found in the previous section. The number of ensembles for the FENE model is *Nens* = 30, 000, and the number of dumbbells within each ensemble *Nd* = 25, 000; for the numerical experiments using the Hookean dumbbell model, the parameters are *Nens* = 75, 000 and *Nd* = 10, 000. Thus, Figsures 9–12 show the shear component *<sup>τ</sup>p*,<sup>12</sup> and the normal stress difference *<sup>τ</sup>p*,<sup>11</sup> − *<sup>τ</sup>p*,22, at (dimensionless) instants of time *t* = 2 and *t* = 3; for the ease of representation, the Hookean stresses are scaled down a factor of 1.81 from the results obtained with the FENE model.

**Figure 9.** Representation of the shear component of the polymer stress tensor, *<sup>τ</sup>p*,12, for the FENE model (*b* = 50) with *c* = 5, *De* = 3, at dimensionless instants of time *t* = 2 and *t* = 3, for variable types of CSRBF. From left to right: *ϕ*3,0, *ϕ*3,1, *ϕ*3,2, *ϕ*3,3. Top row: *t* = 2; bottom row: *t* = 3.

For the shear component *<sup>τ</sup>p*,<sup>12</sup> of the polymer stress tensor using the FENE model (Figure 9), all CSRBFs provide very smooth solutions at *t* = 2. At the final time of the simulation *t* = 3, when larger stresses are produced, the Wendland *ϕ*3,0 CSRBF shows small oscillations and a certain roughness that may translate into an insufficiently resolved velocity field and an overall loss of symmetry in the drop shape. Wendland *ϕ*3,1 offers much better accuracy, as well as an improved spatial resolution of the stresses, to such an extent that increasing the smoothness of the CSRBF (*ϕ*3,0, *ϕ*3,3) offers almost no distinct enhancement.

**Figure 10.** Representation of the shear component of the polymer stress tensor, *<sup>τ</sup>p*,12, for the Hookean dumbbell model (Oldroyd-B) with *c* = 5, *De* = 3, at dimensionless instants of time *t* = 2 and *t* = 3, for variable types of CSRBF. From left to right: *ϕ*3,0, *ϕ*3,1, *ϕ*3,2, *ϕ*3,3. Top row: *t* = 2; bottom row: *t* = 3.

When the Hookean dumbbell model is used (Figure 10), *<sup>τ</sup>p*,<sup>12</sup> shows much larger values throughout the whole simulation. The influence of the type of CSRBF becomes evident at the later stages, when extreme and very localized stresses are observed. The less smooth *ϕ*3,0 is unable to provide a symmetric solution, while the other Wendland functions offer much better results in terms of symmetry.

**Figure 11.** Representation of the normal stress difference of the polymer stress tensor, *<sup>τ</sup>p*,<sup>11</sup> − *<sup>τ</sup>p*,22, for the FENE model (*b* = 50) with *c* = 5, *De* = 3, at dimensionless instants of time *t* = 2 and *t* = 3, for variable types of CSRBF. From left to right: *ϕ*3,0, *ϕ*3,1, *ϕ*3,2, *ϕ*3,3. Top row: *t* = 2; bottom row: *t* = 3.

If we now inspect the results for the normal stress difference, *<sup>τ</sup>p*,<sup>11</sup> − *<sup>τ</sup>p*,22, we observe that when the FENE model is used (Figure 11), remarkable smoothness is retrieved by all the CSRBFs explored, but for *ϕ*3,0, which somewhat underperforms at the last instants of the simulation, when the viscoelastic effects have shown up and the drop attains maximum deformation.

**Figure 12.** Representation of the normal stress difference of the polymer stress tensor, *<sup>τ</sup>p*,<sup>11</sup> − *<sup>τ</sup>p*,22, for the Hookean dumbbell model (Oldroyd-B) with *c* = 5, *De* = 3, at dimensionless instants of time *t* = 2 and *t* = 3, for variable types of CSRBF. From left to right: *ϕ*3,0, *ϕ*3,1, *ϕ*3,2, *ϕ*3,3. Top row: *t* = 2; bottom row: *t* = 3.

Figure 12 represents the normal stress difference, *<sup>τ</sup>p*,<sup>11</sup> − *<sup>τ</sup>p*,22, for the Hooke model. In this case, a higher degree of smoothness has a positive effect, improving the symmetry of the solution, to a certain extent. However, unless the spatial resolution of the velocity field and the velocity gradients is improved, the number of ensembles distributed over the domain is large enough, and the dumbbells per ensemble ensure a sufficiently accurate computation of the "extra-"stress at the ensembles, there is no point in further increasing the quality of the CSRBF for polymer stress reconstruction, since the extremely large, very localized stresses cannot be accurately computed unless the previous conditions are met. On the other hand, all previous conditions can be satisfied so long as the necessary computational resources are available. As a result of all the previous investigations, we find Wendland's *ϕ*3,1 CSRBF as a good choice in terms of accuracy, smoothness, and computational demands, used in this work unless otherwise noted.

### 3.2.3. Flow Pattern under Increasing Viscoelastic Effects

Finally, we explore the flow pattern, streamlines, and isocontours of the polymer stress tensor, in a series of simulations under increasing polymer concentration and relaxation times, using the FENE dumbbell model with extensibility parameter *b* = 50; see Figures 13–16. For all the subsequent numerical experiments, a uniform, unstructured mesh of size *h* = 1/320 is used, with the number of dumbbells per ensemble *Nd* and the number of ensembles *Nens* being collected in Table 1. Notice that, for stronger viscoelastic effects, a larger number of ensembles are scattered throughout the domain, reducing the computationally available number of dumbbells per ensemble, which in turn translates into a larger stochastic noise. Two different density and viscosity ratios are explored to underscore the effects caused by the polymer solution on the shape deformation and flow pattern of the drop (*ρ*2/*ρ*1 = 10−<sup>1</sup> = *μ*2/*μ*1, *Re* = 35, *We* = 10, *Fr* = 1) and of the bubble (*ρ*2/*ρ*1 = <sup>10</sup>−3, *μ*2/*μ*1 = <sup>10</sup>−2, *Re* = 35, *We* = 125, *Fr* = 1).

**Table 1.** Each entry in the table shows *Nd* / *Nens* for increasing viscoelastic effects and different density ratios. *Nd* is the number of dumbbells per ensemble and *Nens* the number of ensembles scattered in the domain.


In Figure 13, we show the final shape and filled isocontours of the normal stress difference of the polymer stress tensor for increasing concentration of the polymer solution and longer relaxation times for the drop (density ratio of 10), while Figure 14 presents the behavior of the shear component under the same circumstances, for the bubble immersed in the polymeric media (density ratio 1000). For the drop and normal stress difference, the results present an excellent degree of symmetry and smoothness, despite the underlying stochastic procedure, even for the strongest viscoelastic effects at (*c* = 5, *De* = <sup>3</sup>), which produce extreme values at the top and bottom of the rising drop.

The bubble depicted in Figure 14 with the shear stress isocontours shows how the harder computational demands of this simulation translate into a modest reduction of the overall symmetry along with more dramatic changes of the bubble shape. Nevertheless, the smoothness of the discrete solution obtained for the polymer stress tensor, even when the stronger polymer concentration is used, is such that it is possible to retrieve the maximum values at the wake of the bubble, with a high degree of accuracy and symmetry. Figures 13 and 14 also evince elastic effects, with heightened drop deformation caused by larger values of the Deborah number (i.e., relaxation time of the polymer) and polymer concentration, while viscous and surface tension forces are kept at the same level.

**Figure 13.** Drop shape and isocontours of normal stress difference of the polymer *<sup>τ</sup>p*,<sup>11</sup> − *<sup>τ</sup>p*,<sup>22</sup> for the FENE model (*b* = 50) with *ρ*2/*ρ*1 = 10−<sup>1</sup> = *μ*2/*μ*1, at *t* = 3. From left to right, *c* = 1, *De* = 1; *c* = 3, *De* = 1; *c* = 5, *De* = 3; and *c* = 9, *De* = 5.

Next, we focus on the streamlines of the flow under increasing viscoelastic effects. In Figure 15, we represent the flow pattern for the drop (moderate density and viscosity ratios), with Figure 15a presenting a snapshot at time *t* = 2, while Figure 15b shows the pattern at the final instant of time, *t* = 3. For these high values of polymer concentration and relaxation times, at *t* = 3, we observe the downwards velocities that characterize the "negative wake" effect [23,76].

For large density and viscosity ratios, Figure 16 depicts the behavior of the deforming bubble at *t* = 3. At the higher levels of polymer concentration and relaxation times, the area in which downwards velocities are observed, comprising the "negative wake", is much larger than in the case of the deforming drop.

**Figure 14.** Drop shape and shear component of the polymer stress tensor *<sup>τ</sup>p*,<sup>12</sup> for the FENE model (*b* = 50) with *ρ*2/*ρ*1 = <sup>10</sup>−3, *μ*2/*μ*1 = <sup>10</sup>−2, at *t* = 3. From left to right, *c* = 1, *De* = 1; *c* = 3, *De* = 1; *c* = 5, *De* = 3; and *c* = 9, *De* = 5.

**Figure 15.** Drop shape and streamlines for the FENE model (*b* = 50), with *ρ*2/*ρ*1 = 10−<sup>1</sup> = *μ*2/*μ*1 at different instants of time: (**a**) *t* = 2 and (**b**) *t* = 3. For each panel, the left figure (yellow line) shows the streamlines for *c* = 5, *De* = 3; the right figure (red line), those for *c* = 9, *De* = 5.

**Figure 16.** Drop shape and streamlines for the FENE model (*b* = 50) with *ρ*2/*ρ*1 = <sup>10</sup>−3, *μ*2/*μ*1 = <sup>10</sup>−2, at *t* = 3. From left to right, *c* = 1, *De* = 1 (blue line); *c* = 3, *De* = 1 (green line); *c* = 5, *De* = 3 (yellow line); and *c* = 9, *De* = 5 (red line).
