Macro–Micro-Coupled Simulations of Dilute Viscoelastic Fluids

: Modeling the ﬂow of polymer solutions requires knowledge at various length and time scales. The macroscopic behavior is described by the overall velocity, pressure, and stress. The polymeric contribution to the stress requires knowledge of the evolution of polymer chains. In this work, we use a microscopic model, the ﬁnitely extensible nonlinear elastic (FENE) model, to capture the polymer’s behavior. The beneﬁt of using microscopic models is that they remain faithful to the polymer dynamics without information loss via averaging. Their downside is the computational cost incurred in solving the thousands to millions of differential equations describing the microstructure. Here, we describe a multiscale ﬂow solver that utilizes GPUs for massively parallel, efﬁcient simulations. We compare and contrast the microscopic model with its macroscopic counterpart under various ﬂow conditions. In particular, signiﬁcant differences are observed under nonlinear ﬂow conditions, where the polymers become highly stretched and oriented.


Introduction
Viscoelastic fluids are materials that exhibit complex behaviors due to their unique composition and structure.Uses of these materials can be found in many industrial, biomedical, and basic science applications [1][2][3][4][5][6].Common examples include liquid crystals, polymers, and colloids.In practice, modeling the behavior of viscoelastic materials requires significant computational power, which clashes with the need for fast simulations.An increasing number of research studies focus on the mathematical modeling of viscoelastic fluids at different scales.To allow such modeling efforts to be useful in practice, there is a critical need for more efficient computational methods.In recent years, GPU computing has emerged as a valuable tool to overcome these challenges.Here, we describe numerical algorithms aimed at using GPU computing in viscoelastic flow simulations.
GPU computing refers to the use of graphics processing units (GPUs) to perform numerical computations.GPU computing has become an important tool for the study of complex fluids due to its ability to handle large amounts of data and perform computations in parallel.Because GPUs can handle massive computational demands, they can provide results in a fraction of the time of traditional CPU-based methods.As computing power continues to increase, GPU-based simulations can become increasingly important for understanding the behavior of complex materials and complex flows.GPUs have been previously used in the modeling of polymeric systems.Reith et al. [7] performed numerical simulations to find the equilibrium properties of polymer rings.Brito et al. [8] used GPUs to accelerate computations of smoothed particle hydrodynamics (SPH) for viscoelastic materials.Ingelsten et al. [9] extended their Lagrangian-Eulerian method for viscoelastic flow simulations to GPU computing.Yang et al. [10] simulated phase separation in polymer blends and polymer solutions using a GPU-based viscoelastic model.Bergamasco et al. [11] introduced a lattice-Boltzmann-based finite-volume formulation, where the macroscopic equations were solved by a finite-volume method and the microscopic equation was solved by a lattice-Boltzmann method.These authors implemented the solution to the microscopic equation in GPUs.Gutierrez-Castillo et al. [12] used GPUs to accelerate their three-dimensional numerical simulations of viscoelastic fluids for both the Oldroyd-B and FENE-P models.Schram and Barkema [13] simulated ring polymer melts via GPU computing using an elastic lattice polymer model.Behbahani et al. [14] investigated high-molecular-weight entangled polymer melts using HOOMD-blue, a GPU accelerated simulation software.In this paper, we use GPU computing to couple dynamics at the macroscale with those at the mesoscale, described by Langevin-type equations.
From a modeling perspective, it is well-known that the stress, τ, of any Newtonian fluid is directly proportional to the imposed shear rate, τ = η γ, where the proportionality constant η is the viscosity of the material.Differences between materials will only affect the value of η but not the functional form of the constitutive equation.Mathematically speaking, this means that for Newtonian fluids, there is a single, general constitutive equation capable of describing different types of Newtonian materials.In principle, this means that the same numerical algorithms can be used for any material by only changing the value of the viscosity constant.In contrast, there is not a single all-encompassing constitutive equation for viscoelastic materials.This is because, in this case, the evolution of the stress as a function of the applied strain can vary from material to material or even for the same material under different flow conditions.Accordingly, the formulation of constitutive equations describing different types of viscoelastic materials is a prolific area of research [15].It is not surprising, then, that because of the lack of a single constitutive equation and due to their ubiquity in many industrial and biological processes, numerical methods for solving the flow of these viscoelastic materials under flows fields of different levels of complexity have been a very active research field.Some extensive reviews in this subject can be found in [9,[16][17][18][19][20].
Broadly speaking, the link between polymer fluids and external fields is relatively wellunderstood: how a polymeric fluid behaves under a given field is directly governed by the dynamics of their corresponding underlying microstructure.For polymeric fluids, this microstructure comprises polymer chains.The underlying dynamics include coiling and uncoiling processes, hindered motion due to physical entanglements or hydrodynamic effects caused by the presence of other polymer molecules in their vicinity, and in some cases, physical cross-linking between chains.In this paper, we focus on a subclass of the available computational techniques: the so-called micro-macro approaches [18].These approaches combine the solution of the conservation equations at the macroscale with models derived from kinetic theory at the microscale.The advantage of using kinetic models is that instead of a homogenized view of the rheology of the material through a macroscopic constitutive equation, the viscoelastic contribution to the stress is obtained using coarse-grained representations of the microstructural dynamics.For more comprehensive reviews on these methods, the reader is referred to [18,[21][22][23][24][25].
Micro-macro methods are especially appealing when modifications to the microstructural dynamics need to be introduced without incurring large numerical changes or extreme computational costs.This is because the evolution of the microscopic description is somehow decoupled from the macroscopic flow solvers.This has the potential of making the model highly versatile and applicable to a range of deformation conditions.In addition, micro-macro methods facilitate the investigation of the effects of various parameters to gain insights into the underlying physics of a given phenomenon and lead to an optimization of a range of industrial processes.On the other hand, micro-macro methods are more computationally demanding than their corresponding macroscopic counterparts.In practice, one must choose between models that are computationally tractable but perhaps do not offer a direct link between the flow-induced development of the micro-structure and the applied flow field or be limited to simple geometries and viscometric flows when using macro-micro methods.
Within the context of microscopic models, different models have been developed to describe the coupling between macroscopic response and microstructure dynamics in polymeric fluids.One family of models are the so-called bead-spring models.These models use molecular coarse-graining to describe the behavior of polymer chains represented as beads connected by massless springs.The simplest of these models considers only two beads, and it is known as the elastic dumbbell model.Although it is widely recognized that a single dumbbell is too simple to be able to describe any complicated dynamics in polymeric systems, it is also known that stretching and orientation alone suffice to give a qualitative description of steady-state rheological properties and flows with slow characteristic time scales [15].Accordingly, these models have been extensively used to develop "an elementary but broad understanding of the relation between macro-molecular motions and rheological phenomena" [15].One of the most common dumbbell models is the Hookean dumbbell model, where the spring connecting the two beads is assumed to follow a linear elastic law ( Hooke's law).This model, although physically unrealistic, has the advantage of being the only model in its class with an exact closure [15,26]; this closure leads to the well-known macroscopic constitutive equation called the Oldroyd-B model.If a different spring law is assumed, the transition between micro and macro models is not directly possible.For instance, the family of models derived from using a finite extensible (FENE) spring force can have many types of closure depending on which approximation is made in the closure step.Many studies have investigated the effects of the different FENE closures; for an example, see Herrchen and Öttinger [27].
In short, the use of dumbbell models allows for the simulation of complex fluids at a molecular level, which, although somewhat unrealistic, can still provide valuable insight into the microstructural responses in these materials.However, simulating complex fluids using dumbbell models requires large amounts of computational power, which has hindered the ability to conduct large-scale simulations until recently.That is why, although some information is lost by averaging the effect of the dumbbells in a single differential equation, macroscopic models are still more widely used in complex flow.The main reason for this is that PDE models are more computationally tractable than their SDE counterparts.Parallel computing has surfaced as an alternative to close the gap between these two computational costs.In particular, GPU computing emerged in past decades as a cheap option for performing parallel computing.The only drawback when dealing with GPUs is that memory transfer between the CPU and the GPU is very costly, which is why these devices perform better on systems that are trivially parallelizable.
Here, we demonstrate the use of GPU computing using non-interacting dumbbell models, which are by design trivially parallelizable.First, we explore two viscometric flows: simple shear and uniaxial extensional flow.Next, we show the coupling between an existing macro solver [28] and the GPU dumbbell models to simulate capillary thinning.Capillary thinning is a phenomenon where a liquid bridge formed between two solid surfaces undergoes thinning due to surface tension forces until it breaks.This is a common phenomenon in a variety of industrial processes, such as coating, printing, and microfabrication [1,3,5].The use of simulations in understanding capillary thinning can provide insight into the underlying physics and aid in the design of improved processes.An advantage of using dumbbell models is that they enable the investigation of the effects of various microstructural parameters on capillary thinning.

Governing Equations 2.1. Macroscopic Flow Equations
At the macroscale, considering an unsteady, incompressible viscoelastic flow and in the absence of external body forces, the flow is governed by the conservation of mass and the conservation of momentum equations where ũ is the velocity vector, p is the pressure, τ is the extra stress tensor, η s is the solvent viscosity, and ρ is the density.We scale the conservation equations as follows: Here, L is a characteristic macoscopic length, U is a characteristic velocity, η p is the zero shear rate polymeric viscosity, and λ is the longest relaxation time.With this scaling, the conservation equations become Several nondimensional groups arise from this scaling: • The Reynolds number, Re = ρUL η , which is the ratio of inertial to viscous forces.

•
The ratio of the solvent to total zero-shear viscosity: , which is the ratio of the fluid relaxation time to a characteristic time for the macroscopic flow.Note that in the case of viscometric flow, in which U/L = γ0 , this is also referred to as the Weissenberg number.

Constitutive Equations for Extra Stress
To solve the resulting system of conservation equations, an expression for the extra stress tensor, τ, is necessary.As discussed in the introduction, there are several ways in which the extra stress can be modeled.Here we use a micro-macro approach [21], where the coarse-grained approximation at the mesoscale utilizes non-interacting elastic dumbbells.Below, we introduce the evolution equations of the dumbbell models used in this study.For completeness, we also introduce the macroscopic closure equations corresponding to each of the two dumbbell models considered.

Microscopic Description-Dumbbell Models
The solvent is assumed to be an incompressible Newtonian fluid of viscosity η s .The configuration of the dumbbell, shown in Figure 1, is described by the end-to-end connector vector Q = r 2 − r 1 and the center-of-mass vector r c = 1 2 (r 1 + r 2 ).Here, we assume a homogeneous distribution, so that the only variable of =interest is Q.The stochastic differential equation describing the evolution of the end-to-end vector is given by [15,29], where F(•) denotes the functional form of the spring force, W t is a two-dimensional Wiener process, ζ is the drag coefficient acting on each bead, and k B T is the thermal energy.
The viscoelastic contribution to the stress is given by Kramer's relation [15,29], where n is the dumbbells' number density, n F(Q) Q is the contribution of the tension on the spring with the spring law F(Q), and nk B T captures effects due to Brownian motion.

Hookean dumbbells:
For Hookean dumbbells, a linear form is assumed for the spring force, F(Q) = HQ, where H is the spring constant.With this, Equations ( 5) and ( 6) become These equations are scaled using Equation (3) and the characteristic microscopic length, With these, In addition, the longest relaxation time is defined as The resulting non-dimensional equations are then where the Deborah number is defined as before, De = λU/L.

FENE Dumbbells
If instead of a linear spring connector a finitely extensible spring is considered, we obtain FENE-type models.These models are derived by introducing Warner's finitely extensible nonlinear connector force law [15], where Introducing this spring law into Equation ( 5) and non-dimensionalizing as before gives where the nonlinear FENE parameters are b = H Q 2 max /(k B T).

Macro-Constitutive Equations-Closure Models
In this section, we introduce macroscopic closures for the two dumbbell models discussed in Section 2.2.1.We introduce these closures here and use them in Sections 3.1 and 3.2 to compare their results with those of their corresponding microscopic dumbbell model.The main objective is to show how much information is lost at the microscale using closure approximations.These shortcomings of the macroscopic closures highlight the importance of micro-macro models to properly capture the underlying dynamics that govern the flow in complex fluids.
Recall that the extra stress is found as a function of the second moment of the end-toend vector, Q.From the equation of motion of a dumbbell, Equation ( 5), the following evolution equation for the second moment of the end-to-end vector Q can be derived: where (•) (1) is the upper convected derivative, defined as •

Upper convected Maxwell model:
If the spring is linear (Hookean), where F(Q) = HQ, the resulting constitute equation obtained from Equation ( 14) is known as the upper convected Maxwell (UCM) model Using the same scaling as before leads to the following non-dimensional constitutive equations: •

FENE-P model
In a FENE model, with the spring law given by Equation ( 12), the evolution equation, Equation ( 14), becomes Different FENE models can be formulated based on the approximation used to calculate the second term in Equation (18).Here, we focus on the FENE-P model, which is obtained by the approximation With this approximation, the constitutive equations become Scaling as before, we obtain For convenience, one can define the non-dimensional configuration tensor, A, as where Q2 0 is the mean-square end-to-end spring length at equilibrium (absence of flow) and d is the dimensionality [15].Note that if we scale the end-to-end vector as before, we have In addition, trace The constitutive equation for A in two dimensions is

Numerical Considerations
In this work, we use graphics processing units (GPU) to perform calculations at the mesoscale.This approach takes advantage of the fact that the stochastic differential equations for non-interacting dumbbells are trivially parallelizable, so that one can concurrently solve millions of equations.This approach eliminates the need for variance reduction approaches [30][31][32] since it is possible to resolve hundreds of thousands of realizations with very little computational time.The main computational scheme is summarized in Figure A1, particular considerations for the Hookean and FENE dumbbells are described below.
One thing to keep in mind when using GPUs is that this approach works best when the dynamics of the individual dumbbells are trivially parallelizable.That is why, for instance, modeling hydrodynamics comes with a much higher computational costs using GPUs.In this study, we limit ourselves to equations of motion that are trivially parallelizable.

•
Parameter values: Unless otherwise noted, the parameter values used in the simulations are given in Table 1.Equation (11a) will be solved for each dumbbell using a forward Euler scheme, Boundary conditions for Q are imposed at the macroscale through u.Whenever necessary to calculate the stress gradient, ghost cells will be included in the macroscale grid.

FENE Dumbbells
Because of the nonlinearity introduced in the evolution equation, a semi-implicit, first-order algorithm is used instead of the forward Euler scheme.Briefly, a two-step operator splitting procedure is implemented, where the two steps are given by Equations ( 26) and ( 27) [29] as follows: and Equation ( 27) can be written as where P is the magnitude of the right-hand side of Equation ( 27).This cubic equation has one real root in the interval 0, √ b [29,35].The approach is to find the magnitude of Q t j+1 using the polynomial's solution in this interval; then, its direction is found by solving Equation ( 27) for each vector component.Note that Equation ( 27) requires the velocity gradient tensor, (∇u) t j+1 , at time t + 1, which is unknown for non-viscometric flows.In our computations, u t j+1 is found by a second-order extrapolation following [35],

Simple Shear
This flow field results from placing the fluid between two parallel plates, with the lower plate stationary and the upper plate moving at a constant velocity U top .Under viscometric conditions, the equations can be simplified by assuming a fully developed flow so that where γ0 = U top /L and L is the plates' separation.
Using the same scaling as before simplifies the velocity to u = [y, 0] , and the velocity gradient tensor reduces to All initial conditions correspond to equilibrium conditions.For dumbbell models, this means the end-to-end vectors are drawn from a normal distribution.For the closure models, this means the extra stress is zero.
The UCM model is an exact closure of the Hookean dumbbell equation, and as expected results from these two models will always agree, this is confirmed in Figure 2, where results from both models in viscometric simple shear flow are compared.11)) and UCM model (Equation ( 17)).The first normal stress difference is calculated as N 1 = τ yy − τ xx .All variables are non-dimensionalized using the scaling given in Equation ( 3) and parameter values given in Table 1.
Similarly, Figure 3 compares results from the FENE-P model with stochastic simulations of FENE dumbbells.The figure shows that the FENE-P approximation underpredicts the nonlinear behavior of the shear stress, in agreement with previous studies [27].24)).All variables are non-dimensionalized using the scaling given in Equation ( 3) and parameter values given in Table 1.
Figure 4 shows that as the nonlinearity, represented by the parameter b, is increased, the discrepancies between the micro and macro representations become more pronounced.For a more in-depth comparison between the FENE model and different FENE closures, the reader is referred to [27].
The advantage of using dumbbell equations, as opposed to closure approximations, is that in addition to the average stress behavior, shown in Figures 2-4, one can also capture the stretching and orientation of dumbbells under flow.In this way, the model predictions can be compared to data obtained from scattering techniques [36].Figure 5 shows the temporal evolution of the distribution of dumbbell lengths, Q x vs. Q y , as a function of shear flow both in the linear regime, De = 10, and the nonlinear regime, De = 500.In the linear regime, dumbbells are stretched and oriented in the direction of the flow, but the distribution of lengths remains Gaussian.In the nonlinear regime, Gaussianity is lost around the same time the overshoot in the shear stress is observed.3) and parameter values given in Table 1.
We can also examine the probability distribution function of the dumbbell orientation as a function of the Deborah number; these results are shown in Figure 6.The figure shows that the nonlinear behavior differs from linear behavior in two main ways.First, the 'molecules' are more oriented with respect to the flow, as proven by the fact that the mean orientation angle at steady state for De = 10 is not quite zero, while for De = 500, the mean angle is zero.Second, narrower distributions for De = 500 indicate that more 'molecules' are aligned with the flow at steady state compared with molecules in the linear regime.A final checkpoint for our simulations is the expansion ratio, defined as where the axes of extension for the dumbbell 'molecule' are calculated as where < • > 0 represents the mean deformation at equilibrium.Figure 7 shows the expansion ratio for different values of De.Compared to the information gathered from the shear stress curves, the expansion ratio does not offer any new insight into the flow dynamics.However, our simulations at the mesoscale allow us to take a closer look at the individual axes of extension and analyze their behavior, as shown in Figure 8.In Figure 8, we show the particulars of the two axes of extension in the nonlinear regime, De = 500, together with their distributions at some selected times.In the scattering plots, we have normalized the axes with respect to their maximum values.However, we should note that the maximum values vary widely for different times.For instance, at t = 1 : max(e x ) = 10.54, while at t = 1000 : max(e x ) = 295.82.This normalization was carried out to visualize the breaking of isotropy in the stretching as the system approaches the stress overshoot.In the nonlinear regime, t > 100, the molecules stretch more in the direction of the flow field, i.e., the x-direction.As the stretching approaches the maximum dumbbell length (b = 300), the stretching is 'accumulated' at the end of the spectrum.

Elongational Flow
This flow field results from stretching the fluid in one direction at a rate ε0 = dε/dt, where ε is the strain.As before, we assume a fully developed flow so that the velocity field is given by ũ Using the same scaling as before simplifies the velocity as u = x, − 1 2 y, − 1 2 z and the velocity gradient tensor as As before, we assume that initially, the system is at equilibrium, as discussed in Section 3.1.
Figure 9 shows the comparison of the macro and stochastic simulations.As expected, the FENE-P closure deviates from the FENE dumbbell model, and this deviation increases in the nonlinear regime.For a complete discussion of FENE models under rapid elongational flow, the reader is referred to [37].24)) with stochastic simulations (Equation ( 13)) under viscometric elongational flow.

Capillary Thinning
Transient elongational rheometry is used to obtain comparative information about the elongational behavior of polymer solutions.The investigation of the elongational behavior of complex fluids is relevant to applications defined by spraying and jetting, such as microfluidic devices, printing, fiber spinning, and understanding how droplet emissions from respiratory events affect the transmission of infectious diseases [1][2][3][4][5][6].Ardekani et al. [34] showed that three conditions should be satisfied for an elongational experiment to be successful: an elastocapillary balance must be established, the range of diameters of the liquid jets should be experimentally accessible, and the formation of secondary droplets along the thread must be suppressed.Clearly, understanding the underlying microstructural dynamics that govern one or more of these conditions is important to evaluate the optimal performance both under experimental and application settings.While most of the theoretical analysis is carried out using homogenized views of the microstructure in the form of macroscopic constitutive equations, using methods that allow for the straightforward explicit integration of microstructural dynamics is appealing.Here, our goal is to illustrate how coarse-grained mesoscopic simulations can also be used to gain insight into the development of droplets in a viscoelastic fluid jet.We note that the objective of this publication is to demonstrate the feasibility of the methods, and we take advantage of symmetry and one-dimensional flow fields to simplify many of the computations.More in-depth studies will be left for future investigations.
In capillary thinning experiments, a piece of fluid is stretched between two separating surfaces, with surface tension being the main driver of the fluid motion.When a Newtonian fluid is stretched, its radius decreases due to surface tension until it breaks into many droplets.For viscoelastic materials, the underlying microstructure can give rise to different configurations.In polymeric materials, these dynamics involve nonlinear effects derived from the finite extensibility of polymer chains within the fluid.As the fluid jet starts to form droplets, the polymer inside the drops remains relaxed, while polymer chains in the filaments between drops become strongly oriented and stretched, essentially leading to the formation of bridges with high tensile strains.This instability, characteristic of fluids that have some elasticity, has become known as a "beads-on-a-string" structure [38][39][40][41].Moreover, under certain flow regimes, when there is a balance of capillary, viscous, and inertial effects, additional structures can appear.For a more in-depth review of the formation of bead-on-a-string structures, see [40].
To solve this flow, a thin film formulation is used where only the leading-order approximation in an expansion on the jet radius, R(z, t), is considered [42].This essentially assumes that the flow inside the fluid jet is one-dimensional so that the only relevant variables are the radius and the axial velocity, v(z, t).In this case, the non-dimensional conservation of mass and of momentum equations are simplified as [43], where 3/2 (37) is the curvature and τ is the polymer contribution to the stress.The equations have been non-dimensionalized by taking the initial radius of the jet, R 0 , as the characteristic macroscopic length scale and the Rayleigh time, t R = ρR 3 0 /γ, as the macroscopic characteristic time scale; this gives the characteristic velocity as γ/(ρR 0 ).Finally, the Ohnesorge number is defined as Oh = η 0 / ργR 0 , and all other dimensionless groups are defined as before.
With this scaling, the dynamics of the axisymmetric slender jet are governed by three non-dimensional groups: the Ohnesorge number, Oh, which is the ratio of surface tension to inertial forces; the parameter β, which is the ratio between the solvent and fluid viscosities; and the Deborah number, De, which measures the ratio of the material's characteristic time to the time scale of the flow.An additional parameter captures the nonlinear viscoelastic behavior of the polymer chains, which, in the case of FENE models, is the finite extensibility parameter, b, defined above.
The nonlinear system of partial differential equations (PDEs), Equation (36), is solved numerically via finite differencing [28].All spatial terms, except convection, are discretized using second-order central differencing.For the convective terms, a velocity-weighted upwind discretization scheme is used.In particular, the derivative is interpolated between pure central differencing and pure upwind.The result is that in regions of slow flow, a central differencing is used; however, as the velocity increases, the weighting moves closer to a pure upwind differencing.For the convective terms where the derivative depends on the velocity, the first derivative is calculated as follows: where g j = v j /v max .We use the implicit Euler method for time discretization combined with time filters, which increase the order of accuracy from first to second order [44].Additionally, we have implemented variable step sizes, which greatly enhances the speed and efficiency of the solver.The numerical solver has been benchmarked against existing simulations of the Oldroyd-B model [34].
In the spatial domain, periodic boundary conditions are used for the velocity, the radius, and the polymer configurations.The initial condition considers a homogeneous base state corresponding to a uniform cylindrical shape and molecules at equilibrium conditions.To this base state a small amplitude perturbation is added.We note that the eigenvalues of the linearized system resulting from expanding the governing equations to first order in the amplitude of these perturbations determines how fast the perturbations grow towards a necked state or whether they decay to leave a uniform filament.Following [34], we set the initial shape of the jet as where k is the wavenumber.At each time step, the system of PDEs is solved to find the velocity field.The field is then used to find the velocity gradient, ∇u, at each grid point.Using these values, the dumbbell configurations are updated using Equation (13a) and the extra stress calculated with Equation (13b).The resulting values for the extra stress at each grid point are introduced in the macroscopic simulation through Equation (36b).
Below, we show the results obtained when the Deborah number, the wavelength values, and the maximum extensibility of the FENE dumbbells are varied.As a metric for microstructural organization, we use the expansion ratio defined in Equation (32).

Deborah Number
Figure 10 shows the evolution of the jet morphology as a function of the Deborah number, De.As the elasticity of the fluid is increased-a larger De-the formation of droplets is inhibited.This inhibition is related to the higher tensile stress in the filament prior to the formation of the droplets, represented by a larger expansion ratio in the figures.This, in turn, hinders the reorganization of the strands that would result in a lower expansion ratio in the middle of the jet, which is what ultimately leads to the formation of droplets.

Wavelength
The effect of the excitation wavenumber on the formation of droplets is well understood; see, for example, [34,[45][46][47].Briefly, there are three regimes: (i) for very small wavelengths, the jets are characterized by traveling capillary waves; (ii) in the middle regime, the jet develops single and multiple drops; and (iii) at large wavelength numbers, no formation of droplets is observed.Our simulations show the same behavior.In Figure 11 the jet evolutions for two different wavenumbers are shown.The expansion ratio shows the same trends discussed above, where the main stretching is observed at the bridges.In addition, this figure shows that smaller droplets result in a higher level of stretching in other regions of the filament.

Nonlinear FENE Parameter
Figure 12 shows the evolution of droplets for two values of the nonlinear parameter b.Recall that b controls the maximum stretch allowed in the dumbbells.It is this stretch that gives polymer solutions their elasticity.In this manner, filament formation has a similar trend when varying b and De, i.e., as maximum extensibility is increased, satellite bead formation is suppressed.Under conditions where satellite beads do form, smaller, secondary droplets between the satellite and main beads may appear.As shown in Figure 12, larger values of b suppress these secondary droplets.Furthermore, another main difference is observed at the pinch sites, with smaller values of b resulting in more pronounced pinching.This is a result of less elasticity, hence more viscous-like behavior.

Discussion
In this work, we have considered multiscale simulations of dilute polymeric solutions under homogeneous and non-homogeneous flow conditions.At the microscale, Langevin equations describing the behavior of elastic dumbbells are solved using GPUs, which allows for massively parallel simulations.Considering dumbbell equations provides more accurate information on the stretch and orientation of polymer chains, as opposed to macroscopic constitutive equations.In this way, the stretching and orientation of dumbbells are captured in a similar manner in which an experimentalist will analyze fluid behavior using scattering techniques.
Under viscometric conditions of simple shear and uniaxial extension, we showed that there is a significant difference between predictions from macroscopic and mesoscopic models.These differences are exacerbated under nonlinear flow conditions.In the linear regime, dumbbells are stretched and oriented in the direction of the flow, but the distribution of lengths remains Gaussian.In the nonlinear regime, Gaussianity is lost around the same time the overshoot in the shear stress is observed.Furthermore, the PDF of dumbbell orientation shows that 'molecules' are more oriented with respect to the flow in the nonlinear regime compared with molecules in the linear regime.
Under capillary thinning, we investigated the effects of the Deborah number, wavenumber, and maximum extensibility on the formation of satellite drops.As the elasticity of the fluid (De) is increased, the formation of satellite droplets is inhibited.This inhibition is related to a higher tensile stress in the filament prior to the formation of the droplets, represented by a larger expansion ratio.This, in turn, hinders the reorganization of the strands that would result in a lower expansion ratio in the middle of the jet, which is what ultimately leads to the formation of droplets.The finite extensibility plays a similar role.By increasing the extensibility, the polymer becomes more elastic, suppressing droplet formation.For less elastic materials, significant pinching can be observed where the bead and string meet, as one might expect from a viscous filament.The wavenumber of perturbation, although not a material feature, plays a significant role in filament evolution.The expansion ratio shows the same trends as elastic effects, where the main stretching is observed at the bridges.In addition, smaller droplets result in a higher level of stretching in other regions of the filament.
The objective of this work is to further examine and reveal the power and utility of microscopic descriptions of macromolecules.One of the primary concerns in using micro-macro multiscale simulations is computational efficiency.In homogeneous flow fields, configurations are easy to follow around due to the absence of spatial resolution.We extended the simulation capabilities to an inhomogeneous flow field, namely onedimensional capillary thinning.The symmetry and thin film assumptions resulted in the problem being reduced to one spatial dimension.This allowed the decoupling of the dumbbell stretching in two dimensions, which simplifies computational time and effort.The goal now is to consider more general flow fields.In order to extend the solver's capabilities to other types of flow, a more careful coupling of the micro and macro

Figure 1 .
Figure 1.Dumbbell configuration given by the coordinates of each dumbbell, r 1 and r 2 , or its center of mass, r c and end-to-end vector, Q.

Figure 2 .
Figure 2. Extra stress results in viscometric simple shear flow for Hookean dumbbells (Equation (11)) and UCM model (Equation (17)).The first normal stress difference is calculated as N 1 = τ yy − τ xx .All variables are non-dimensionalized using the scaling given in Equation (3) and parameter values given in Table1.

Figure 3 .
Figure 3. Extra stress in viscometric simple shear flow for FENE models compared to the FENE-P closure (Equation (24)).All variables are non-dimensionalized using the scaling given in Equation (3) and parameter values given in Table1.

Figure 4 .
Figure 4. Stochastic simulations for the FENE model with different values of the nonlinear parameter b.All variables are non-dimensionalized using the scaling given in Equation (3) and parameter values given in Table1.

Figure 5 .
Figure 5. Scatter plots of Q x and Q y in simple shear flow for two values of the Deborah number De. Note: the scale of the scatter plots is not the same across the different plots.

Figure 6 .
Figure 6.Probability distribution function (PDF) of dumbbells' orientation as a function of time and Deborah number (De).

Figure 7 .
Figure 7. Expansion ratio as a function of time and Deborah number.The ratio is calculated using Equation (32).

Figure 8 .
Figure 8. Axes of extension as a function of time and De = 500.Distribution plots at the bottom have been normalized by their corresponding maximum values.Between 0 and 100, the system is in the linear regime, and the stretching in the xand y-directions are proportional to each other.In the nonlinear regime, t > 100, the molecules stretch more in the direction of the flow field, i.e., the x-direction.As the stretching approaches the maximum dumbbell length (b = 300), the stretching is 'accumulated' at the end of the spectrum.

Figure 10 .
Figure 10.Filament stretching evolution coupled with microstructure evolution.In this model, the microstructure is represented by non-interacting dumbbells, and their conformation is given by the expansion ratio in Equation (32).Other parameters are β = 0.27, Oh = 0.04, k = 0.8, and b = 10 5 .