Macro-Micro Coupled Simulations of Dilute Viscoelastic Fluids

1. 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. Common examples include liquid crystals, polymers, and colloids. A large number of research undertakings have been focused on the mathematical modeling of viscoelastic fluids, as a result computational efficiency is of particular importance for such modeling efforts to be useful in practice. However, modeling the behavior of viscoelastic materials requires significant computational power which clashes with the need for fast simulations. In recent years, GPU computing has emerged as a valuable tool to overcome this 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., [1] performed numerical simulations to find the equilibrium properties of polymer rings. Brito et al., [2] used GPUs to accelerate computations of smoothed particle hydrodynamics (SPH) for viscoelastic materials. Ingelsten et al., [3] extended their Lagrangian-Eulerian method for viscoelastic flow simulation to GPU computing. Yang et al., [4] simulated phase separation in polymer blends and polymer solutions using a GPU-based viscoelastic model. Bergamasco et al., [5] introduce a lattice-Boltzmann-based finite-volume formulation, where the macroscopic equations were solved by a finite-volume method and the microscopic equation by a lattice-Boltzmann method and implemented the solution to the microscopic equation in GPUs. Gutierrez-Castillo et al., [6] used GPU to accelerate their three-dimensional numerical simulations of viscoelastic fluids for both the Oldroyd-B and FENE-P models.

From a modeling perspective, it is well-known that the stress, $\mathit{\tau }$, of any Newtonian fluid is directly proportional to the imposed shear rate, $\mathit{\tau }=\eta \dot{\gamma }$, where the proportionally constant, $\eta$, is the viscosity of the material. Differences between materials will only affect the value of $\eta$, but not the functional form of the constitutive equation. Mathematical speaking this means that for Newtonian fluids there is a single, general constitutive equation capable of describing different types of Newtonian materials and 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 [7]. 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 [3,8,9,10,11,12].

Broadly speaking, the link between polymer fluids and external fields is relatively well-understood: 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, and 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 [10]. These approaches combine the solution of the conservation equations at the macro-scale with models derived from kinetic theory at the micro-scale. 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-grain representations of the microstructural dynamics. For more comprehensive reviews in these methods the reader is referred to [10,13,14,15,16].

Micro-macro methods are especially appealing when modifications to the microstructural dynamics need to be introduced without incurring in 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 the 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 [7]. Accordingly, these models had been extensively used to develop “an elementary but broad understanding of the relation between macro-molecular motions and rheological phenomena” [7]. 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 (Hooks law). These models, although physically unrealistic, has the advantage of being the only model in its class with an exact closure [7,17]. This closure leads to the well-known macroscopic constitutive equation called 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 example see Herrchen and Öttinger [18].

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 behavior and interactions of 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 since they are more computationally tractable than their micro alternatives. Parallel computing has surfaced as an alternative to close the gap between these two computational costs. In particular, GPU computers had 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 [19] 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. 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.

2. Governing Equations

2.1. Macroscopic Flow

At the macro-scale, 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,

$$\begin{array}{ccc}\hfill \tilde{\mathbf{\nabla }}·\tilde{\mathbf{u}}& =& 0,\hfill \\ \hfill \rho \left(\frac{\partial \tilde{\mathbf{u}}}{\partial \tilde{t}}+\tilde{\mathbf{u}}·\tilde{\mathbf{\nabla }}\tilde{\mathbf{u}}\right)& =& -\tilde{\mathbf{\nabla }}\tilde{p}+{\eta }_s{\tilde{\mathbf{\nabla }}}^2\tilde{\mathbf{u}}+\tilde{\mathbf{\nabla }}·\tilde{\mathit{\tau }},\hfill \end{array}$$

where $\tilde{\mathbf{u}}$ is the velocity vector, $\tilde{p}$ the pressure, $\tilde{\mathit{\tau }}$ the extra stress tensor, ${\eta }_s$ the solvent viscosity, and $\rho$ the density.

We scale the conservation equations as follows,

$$\begin{array}{c}\hfill t=\tilde{t}·\left(\frac{U}{L}\right),\mathbf{x}=\tilde{\mathbf{x}}·\left(\frac{1}{L}\right),\mathbf{u}=\tilde{\mathbf{u}}·\left(\frac{1}{U}\right),p=\tilde{p}·\left(\frac{1}{\rho U^2}\right),\mathit{\tau }=\tilde{\mathit{\tau }}·\left(\frac{\lambda }{{\eta }_p}\right).\end{array}$$

Here L is a characteristic macoscopic length; U is a characteristic velocity; ${\eta }_p$ is the zero shear rate polymeric viscosity; and $\lambda$ is the longest relaxation time. With this scaling the conservation equations become,

$$\begin{array}{ccc}\hfill \mathbf{\nabla }·\mathbf{u}& =& 0,\hfill \end{array}$$

(1a)

$$\begin{array}{ccc}\hfill \frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}·\mathbf{\nabla }\mathbf{u}& =& -\mathbf{\nabla }p+\frac{1}{Re}\left[\beta {\mathbf{\nabla }}^2\mathbf{u}+\frac{1-\beta }{De}\mathbf{\nabla }·\mathit{\tau }\right]\hfill \end{array}$$

(1b)

where

$$\begin{array}{c}\hfill Re=\frac{\rho UL}{\eta },\eta ={\eta }_s+{\eta }_p,\beta =\frac{{\eta }_s}{\eta },De=\lambda \frac{U}{L}.\end{array}$$

2.2. Constitutive Equations for Extra Stress

To solve the resulting system of conservation equations an expression for the extra stress tensor, $\mathit{\tau }$, is necessary. As discussed in the introduction, there are several ways in which the extra-stress can be modeled. Here we use micro-macro approach [13], where the coarse-grained approximation at the mesoscale are 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 to the two dumbbell models considered.

2.3. Dumbbell Models

The solvent is assumed to be an incompressible Newtonian fluid of viscosity ${\eta }_s$. The configuration of the dumbbell is described by the end-to-end connector vector $\mathbf{Q}={\mathbf{r}}_2-{\mathbf{r}}_1$ and the center-of-mass vector ${\mathbf{r}}_c=\frac{1}{2}\left({\mathbf{r}}_1+{\mathbf{r}}_2\right)$. Here we assume a homogeneous distribution, so that the only variable of interest is $\mathbf{Q}$.

The stochastic differential equation describing the evolution of the end-to-end vector is given by [7,20],

$$d{\tilde{\mathbf{Q}}}_t={\left(\tilde{\mathbf{\nabla }}\tilde{\mathbf{u}}\right)}^{\top }·{\tilde{\mathbf{Q}}}_t-\frac{2}{\zeta }F\left({\tilde{\mathbf{Q}}}_t\right)+\sqrt{\frac{4k_{B}T}{\zeta }}d{\mathbf{W}}_t.$$

(2)

Where $F(·)$ denotes the functional form of the spring force, ${\mathbf{W}}_t$ is a two-dimensional Wiener process, $\zeta$ 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 the Kramer’s relation [7,20],

$$\begin{array}{c}\hfill \tilde{\mathit{\tau }}=n〈F\left(\tilde{\mathbf{Q}}\right)\tilde{\mathbf{Q}}〉-n{k}_{B}T\mathbf{I},\end{array}$$

(3)

where n is the dumbbells number density, $n〈F\left(\mathbf{Q}\right)\mathbf{Q}〉$ is the contribution from the tension on the spring with spring law $F\left(\mathbf{Q}\right)$, and $n{k}_{B}T$ capture effects due to Brownian motion.

2.3.1. Hookean Dumbbells

For Hookean dumbbells a linear form is assumed for the spring force, $F\left(\mathbf{Q}\right)=H\mathbf{Q}$, where H is the spring constant. With this, Eqns. (2) and (3) become,

$$\begin{array}{ccc}\hfill d{\tilde{\mathbf{Q}}}_t& =& {\left(\tilde{\mathbf{\nabla }}\tilde{\mathbf{u}}\right)}^{\top }·{\tilde{\mathbf{Q}}}_t-\frac{2H}{\zeta }{\tilde{\mathbf{Q}}}_t+\sqrt{\frac{4k_{B}T}{\zeta }}d{\mathbf{W}}_t,\hfill \\ \hfill \tilde{\mathit{\tau }}& =& nH〈\tilde{\mathbf{Q}}\tilde{\mathbf{Q}}〉-n{k}_{B}T\mathbf{I}.\hfill \end{array}$$

These equations are scaled by the following characteristic microscopic length,

$$\begin{array}{ccc}\hfill x& =& \sqrt{\frac{k_{B}T}{H}},\hfill \end{array}$$

and all other scales are as before so that

$$\begin{array}{c}\hfill {\eta }_p=n{k}_{B}T\lambda ,\mathit{\tau }=\tilde{\mathit{\tau }}·\left(\frac{\lambda }{{\eta }_p}\right)=\tilde{\mathit{\tau }}·\left(\frac{1}{n{k}_{B}T}\right),\mathbf{Q}=\tilde{\mathbf{Q}}·\left(\sqrt{\frac{H}{k_{B}T}}\right).\end{array}$$

In addition, the longest relaxation time is defined as

$$\begin{array}{ccc}\hfill \lambda & =& \frac{\zeta }{4H}.\hfill \end{array}$$

The resulting non-dimensional equations are,

$$\begin{array}{ccc}\hfill d{\mathbf{Q}}_t& =& {\left(\mathbf{\nabla }\mathbf{u}\right)}^{\top }·{\mathbf{Q}}_t-\frac{1}{2De}{\mathbf{Q}}_t+\frac{1}{\sqrt{De}}d{\mathbf{W}}_t,\hfill \end{array}$$

(5a)

$$\begin{array}{ccc}\hfill \mathit{\tau }& =& 〈\mathbf{Q}\mathbf{Q}〉-\mathbf{I},\hfill \end{array}$$

(5b)

where the Deborah number is defined as before, $De=\lambda U/L$.

2.3.2. 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 non-linear connector force law [7],

$$\begin{array}{c}\hfill F\left(\tilde{\mathbf{Q}}\right)=\frac{H\tilde{\mathbf{Q}}}{1-{\left(\tilde{Q}/Q_{\max }\right)}^2},\end{array}$$

(6)

where $Q^2=\mathbf{Q}·\mathbf{Q}$ and $Q_{\max }$ is the maximum allowed extension of the dumbbell.

Introducing this spring law into Eqns. (2) and non-dimensionalizing as before gives

$$\begin{array}{ccc}\hfill d{\mathbf{Q}}_t& =& {\left(\mathbf{\nabla }\mathbf{u}\right)}^{\top }·{\mathbf{Q}}_t-\frac{1}{2De}\left(\frac{{\mathbf{Q}}_t}{1-Q^2/b}\right)+\frac{1}{\sqrt{De}}d{\mathbf{W}}_t,\hfill \end{array}$$

(7a)

$$\begin{array}{ccc}\hfill \mathit{\tau }& =& 〈\frac{\mathbf{Q}\mathbf{Q}}{1-Q^2/b}〉-\mathbf{I}.\hfill \end{array}$$

(7b)

where $b=H{Q}_{\max }^2/\left(k_{B}T\right)$.

2.4. Macro-constitutive equations

From the equation of motion of a dumbbell, Equation (2), the following evolution equation for the second moment of the end-to-end vector $\mathbf{Q}$ can be derived,

$$\begin{array}{c}\hfill {\left(〈\tilde{\mathbf{Q}}\tilde{\mathbf{Q}}〉\right)}_{\left(1\right)}+\frac{4}{\zeta }〈F\left(\tilde{\mathbf{Q}}\right)\tilde{\mathbf{Q}}〉=\frac{4k_{B}T}{\zeta }\mathbf{I},\end{array}$$

(8)

where ${(·)}_{\left(1\right)}$ is the upper convected derivative defined as

$$\begin{array}{c}\hfill {(·)}_{\left(1\right)}=\frac{\partial (·)}{\partial t}+\mathbf{u}·\mathbf{\nabla }(·)-(·)·\left(\mathbf{\nabla }\mathbf{u}\right)-{\left(\mathbf{\nabla }\mathbf{u}\right)}^{\top }·(·).\end{array}$$

2.4.1. Upper convected Maxwell Model

If the spring is linear (Hookean), $F\left(\mathbf{Q}\right)=H\mathbf{Q}$, the resulting constitute equation obtained from Equation (8) is known as the upper convected Maxwell (UCM) model,

$$\begin{array}{cc}& \tilde{\mathit{\tau }}=nH〈\tilde{\mathbf{Q}}\tilde{\mathbf{Q}}〉-n{k}_{B}T\mathbf{I},\\ & {\displaystyle {\left(〈\tilde{\mathbf{Q}}\tilde{\mathbf{Q}}〉\right)}_{\left(1\right)}+\frac{4H}{\zeta }〈\tilde{\mathbf{Q}}\tilde{\mathbf{Q}}〉=\frac{4k_{B}T}{\zeta }\mathbf{I}.}\end{array}$$

Using the same scaling as before leads to the following non-dimensional constitutive equations,

$$\begin{array}{cc}& \mathit{\tau }=〈\mathbf{Q}\mathbf{Q}〉-\mathbf{I},\end{array}$$

(9a)

$$\begin{array}{cc}& De{\left(〈\mathbf{Q}\mathbf{Q}〉\right)}_{\left(1\right)}+〈\mathbf{Q}\mathbf{Q}〉=\mathbf{I}.\end{array}$$

(9b)

2.4.2. FENE-P model

In a FENE model, with spring law given by Equation (6), the evolution equation, Equation (8), becomes

$$\begin{array}{c}\hfill {\left(〈\tilde{\mathbf{Q}}\tilde{\mathbf{Q}}〉\right)}_{\left(1\right)}+\frac{4H}{\zeta }〈\frac{H\tilde{\mathbf{Q}}\tilde{\mathbf{Q}}}{1-{(\tilde{Q}/Q_{\max })}^2}〉=\frac{4k_{B}T}{\zeta }\mathbf{I}.\end{array}$$

Different FENE models can be formulated based on the approximation used to calculate the second term in the equation above. In particular, the FENE-P model is obtained by the approximation,

$$\begin{array}{c}\hfill 〈\frac{H\tilde{\mathbf{Q}}\tilde{\mathbf{Q}}}{1-{(\tilde{Q}/Q_{\max })}^2}〉=\frac{〈\tilde{\mathbf{Q}}\tilde{\mathbf{Q}}〉}{1-〈{\tilde{Q}}^2〉/Q_{\max }^2},\end{array}$$

with this approximation the constitutive equations become,

$$\begin{array}{cc}& {\displaystyle \tilde{\mathit{\tau }}=\frac{nH}{1-〈{\tilde{Q}}^2〉/{\tilde{Q}}_{\max }^2}〈\tilde{\mathbf{Q}}\tilde{\mathbf{Q}}〉-n{k}_{B}T\mathbf{I},}\\ & & \\ & {\displaystyle {\left(〈\tilde{\mathbf{Q}}\tilde{\mathbf{Q}}〉\right)}_{\left(\tilde{1}\right)}+\frac{4H}{\zeta }\frac{〈\tilde{\mathbf{Q}}\tilde{\mathbf{Q}}〉}{1-〈{\tilde{Q}}^2〉/{\tilde{Q}}_{\max }^2}=\frac{4k_{B}T}{\zeta }\mathbf{I}.}\end{array}$$

Scaling as before we obtain,

$$\begin{array}{cc}& \mathit{\tau }=\frac{〈\mathbf{Q}\mathbf{Q}〉}{{\displaystyle 1-〈Q^2〉/b}}-\mathbf{I},\\ & {\displaystyle De{\left(〈\mathbf{Q}\mathbf{Q}〉\right)}_{\left(1\right)}+\frac{〈\mathbf{Q}\mathbf{Q}〉}{1-〈Q^2〉/b}=\mathbf{I}.}\end{array}$$

For convenience, one can define the non-dimensional configuration tensor, $\mathbf{A}$, as

$$\begin{array}{c}\hfill \mathbf{A}\equiv d\frac{〈\tilde{\mathbf{Q}}\tilde{\mathbf{Q}}〉}{{〈{\tilde{Q}}^2〉}_0},\end{array}$$

where ${〈{\tilde{Q}}^2〉}_0$ is the mean-square end-to-end spring length at equilibrium (absence of flow) and d is the dimensionality.

Note that if we scale the end-to-end vector as before, we have

$$\begin{array}{c}\hfill \mathbf{A}\equiv d\frac{\frac{k_{B}T}{H}〈\mathbf{Q}\mathbf{Q}〉}{{〈{\tilde{Q}}^2〉}_0}=d〈\mathbf{Q}\mathbf{Q}〉\left(\frac{k_{B}T}{H{〈{\tilde{Q}}^2〉}^2}\right)=d〈\mathbf{Q}\mathbf{Q}〉\left(\frac{1}{d}\right)=〈\mathbf{Q}\mathbf{Q}〉.\end{array}$$

In addition, trace$\left(\mathbf{A}\right)=d〈Q^2〉.$

The constitutive equation for $\mathbf{A}$ in two-dimensions is,

$$\begin{array}{cc}& {\displaystyle \mathit{\tau }=\frac{\mathbf{A}}{1-\mathrm{trace}\left(\mathbf{A}\right)/\left(2b\right)}-\mathbf{I}.}\end{array}$$

(10a)

$$\begin{array}{cc}& {\displaystyle De{\mathbf{A}}_{\left(1\right)}+\frac{\mathbf{A}}{1-\mathrm{trace}\left(\mathbf{A}\right)/\left(2b\right)}=\mathbf{I},}\end{array}$$

(10b)

2.5. Numerical Solutions

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 paralellizable, so that one can concurrently solve millions of equations. This approach eliminates the need for variance reduction approaches [21,22,23], since it is possible to resolve hundreds of thousands of realizations with very little computational time. The disadvantage of using GPUs is that the stochastic equations need to be fully decoupled; for instance, hydrodynamic interactions are harder to simulate using GPUs.

2.5.1. Numerical considerations for Hookean dumbbells

Equation (5a) will be solved for each dumbbell using a forward Euler scheme,

$$\begin{array}{c}\hfill \mathbf{Q}\left(t_{j+1}\right)=\mathbf{Q}\left(t_j\right)+\left[{\left(\mathbf{\nabla }\mathbf{u}\right)}^{\top }·\mathbf{Q}\left(t_j\right)-\frac{1}{2De}\mathbf{Q}\left(t_j\right)\right]\Delta t_j+\frac{1}{\sqrt{De}}\Delta {\mathbf{W}}_j.\end{array}$$

Boundary conditions for $\mathbf{Q}$ are imposed at the macro-scale through $\mathbf{u}$. Whenever necessary to calculate the stress gradient, ghost cells will be included in the macro-scale grid.

2.5.2. Numerical considerations for FENE dumbbells

Because of the non-linearity 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 [20],

$$\begin{array}{c}\hfill \overline{\mathbf{Q}}\left(t_{j+1}\right)=\mathbf{Q}\left(t_j\right)+\left[{\left(\mathbf{\nabla }\mathbf{u}\right)}^{\top }·\mathbf{Q}\left(t_j\right)-\frac{1}{2De}\left(\frac{\mathbf{Q}\left(t_j\right)}{1-Q{\left(t_j\right)}^2/b}\right)\right]\Delta t_j+\frac{1}{\sqrt{De}}\Delta {\mathbf{W}}_j,\end{array}$$

(11)

and,

$$\begin{array}{ccc}\hfill \left[1+\frac{1}{4De}\frac{\Delta t_j}{1-\frac{Q^2\left(t_{j+1}\right)}{b}}\right]\mathbf{Q}\left(t_{j+1}\right)& =& \mathbf{Q}\left(t_j\right)+\frac{1}{\sqrt{De}}\Delta {\mathbf{W}}_j+\hfill \\ \hfill & & \frac{1}{2}[{\left(\mathbf{\nabla }\mathbf{u}\right)}^{\top }\left(t_{j+1}\right)·\overline{\mathbf{Q}}\left(t_{j+1}\right)+{\left(\mathbf{\nabla }\mathbf{u}\right)}^{\top }\left(t_j\right)·\mathbf{Q}\left(t_j\right)\hfill \\ & & \left.-\frac{1}{2De}\left(\frac{\mathbf{Q}\left(t_j\right)}{1-\frac{Q^2\left(t_j\right)}{b}}\right)\right]\Delta t_j.\hfill \end{array}$$

(12)

Equation (12) can be written as

$$\begin{array}{c}\hfill Q_{j+1}^3-P{Q}_{j+1}^2-b\left(1+\frac{1}{4De}\Delta t\right)Q_{j+1}-bP=0,\end{array}$$

where P is the magnitude of the right-hand-side. This cubic equation has one real root in the interval $\left(0,\sqrt{b}\right)$[20,24]. The approach is to find the magnitude of $\mathbf{Q}\left(t_{j+1}\right)$ using the polynomial’s solution in this interval, then its direction is found by solving Equation (12) for each vector component.

Note that Equation (12) requires the velocity gradient tensor, ${\left(\mathbf{\nabla }\mathbf{u}\right)}^{\top }\left(t_{j+1}\right)$, at time $t+1$, which is unknown for non-viscometric flows. In our computations, $\mathbf{u}\left(t_{j+1}\right)$ is found by a second order extrapolation following [24],

$$\begin{array}{c}\hfill \mathbf{u}\left(t_{j+1}\right)=2\mathbf{u}\left(t_j\right)-\mathbf{u}\left(t_{j-1}\right).\end{array}$$

3. Solutions Under Viscometric Flows

3.1. Simple Shear

This flow field results from placing the fluid between two parallel plates, with the lower plate stationary and the upper plate moving a constant velocity $U_{\mathrm{top}}$. Under viscometric conditions the equations can be simplified by assuming a fully developed flow so that

$$\begin{array}{c}\hfill \tilde{\mathbf{u}}={\left[{\dot{\gamma }}_0\tilde{y},0\right]}^{\top },\frac{\partial \left(·\right)}{\partial x}=0,\end{array}$$

where ${\dot{\gamma }}_0=U_{\mathrm{top}}/L$, and L is the plates separation.

Using the same scaling as before simplifies the velocity as $\mathbf{u}={\left[y,0\right]}^{\top }$ and the velocity gradient tensor reduces to

$$\begin{array}{c}\hfill \mathbf{\nabla }\mathbf{u}=\left(\begin{array}{cc}0& 1\\ 0& 0\end{array}\right).\end{array}$$

3.2. Results

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 1 where results from both models in viscometric simple shear flow are compared.

Similarly, Figure 2 compares results from the FENE-P model with stochastic simulations of FENE dumbbells. The figure shows that the FENE-P approximation under predicts the nonlinear behavior of the shear stress, in agreement with previous studies [18].

Figure 3 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 [18].

The advantage of using dumbbell equations as opposed to closure approximations is that in addition to the average stress behavior, shown in Figure 1, Figure 2 and Figure 3, one can also capture the stretching and orientation of dumbbells under flow in a similar manner in which an experimentalist will analyze fluid behavior using scattering techniques [25]. Figure 4 shows the temporal evolution of the distribution of dumbbell lengths, $Q_x$ vs $Q_y$, as 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.

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 5. The figure shows that the nonlinear behavior differs from a 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,

$$\begin{array}{c}\hfill e=\sqrt{\frac{e_x^2+e_y^2}{2}}.\end{array}$$

(13)

where the axes of extension for the dumbbell `molecule’ are calculated as

$$\begin{array}{c}\hfill e_i=\sqrt{\frac{〈Q_i^2〉}{{〈Q_i^2〉}_0}},i=x,y\end{array}$$

here $<·{>}_0$ represents mean deformation at equilibrium.

Figure 6 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 axis of extension and analyze their behavior, as shown in Figure 7.

In Figure 7 we show the particulars of the two axis of extension in the nonlinear regime, $De=500$, together with their distributions at some selected times. In the scattering plots we have normalized the axis 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 \left(e_x\right)=10.54$, while at $t=1000:\max \left(e_x\right)=295.82$. This normalization was done to visualize the breaking of isotropy in the stretching as the system approaches the stress overshoot.

3.3. Elongational Flow

This flow field results from stretching the fluid in one direction at a rate ${\dot{\epsilon }}_0=d\epsilon /dt$, where $\epsilon$ is the strain. As before, we assume a fully developed flow so that the velocity field is given by,

$$\begin{array}{ccc}\hfill \tilde{\mathbf{u}}& =& {\left[{\displaystyle {\dot{\epsilon }}_0\tilde{x},-\frac{1}{2}{\dot{\epsilon }}_0\tilde{y},-\frac{1}{2}{\dot{\epsilon }}_0\tilde{z}}\right]}^{\top }\hfill \end{array}$$

Using the same scaling as before simplifies the velocity as $\mathbf{u}={\left[{\displaystyle x,-\frac{1}{2}y,-\frac{1}{2}z}\right]}^{\top }$ and the velocity gradient tensor is

$$\begin{array}{c}\hfill \mathbf{\nabla }\mathbf{u}=\left(\begin{array}{ccc}1& 0& 0\\ 0& {\displaystyle -\frac{1}{2}}& 0\\ 0& 0& {\displaystyle -\frac{1}{2}}\end{array}\right).\end{array}$$

Figure 8 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 [26].

4. Capillary Thinning

Transient elongational rheometry is used to obtain comparative information about elongational behavior of polymer solutions. 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 transmission of infectious diseases. Ardekani et al., [27] 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 optimal performance both under experimental and application settings. While most of the theoretical analysis is done using homogenized views of the microstructure, in the form of macroscopic constitutive equations, using methods that allow for a straightforward explicit integration of microstructural dynamics is appealing. Here our goal is to illustrate how coarse-grained mesoscopic simulations can be also 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 having high tensile strains. This instability, characteristic of fluids that have some elasticity, has become known as “beads-on-a-string” structure [28,29,30,31]. 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 [30].

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 [32]. 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 [33],

$$\begin{array}{ccc}\hfill & {\displaystyle \frac{\partial R^2}{\partial t}+\frac{\partial }{\partial z}\left(v{R}^2\right)=0}& \\ \hfill \end{array}$$

(14a)

$$\begin{array}{cc}\hfill & {\displaystyle \left(\frac{\partial v}{\partial t}+v\frac{\partial v}{\partial z}\right)=-\frac{\partial \kappa }{\partial z}+3Oh\frac{\beta }{R^2}\frac{\partial }{\partial z}\left(R^2\frac{\partial v}{\partial z}\right)+\frac{Oh}{De}\frac{\left(1-\beta \right)}{R^2}\frac{\partial }{\partial z}\left(R^2({\mathit{\tau }}_{zz}-{\mathit{\tau }}_{rr})\right),}\end{array}$$

(14b)

where

$$\begin{array}{c}\hfill {\displaystyle \kappa =\frac{1}{R{\left(1+R_z^2\right)}^{1/2}}-\frac{R_{zz}}{{\left(1+R_z^2\right)}^{3/2}}}\end{array}$$

is the curvature and $\mathit{\tau }$ 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, ${\displaystyle t_R=\sqrt{\rho R_0^3/\gamma },}$ as the macroscopic characteristic time scale; which gives the characteristic velocity as ${\displaystyle \sqrt{\gamma /\left(\rho R_0\right)}}$. Finally, the Ohnesorge number is defined as ${\displaystyle Oh={\eta }_0/\sqrt{\rho \gamma R_0}}$, 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 $\beta$, 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), Eqns. (14), is solved numerically via finite differencing [19]. 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,

$${\left(\frac{\partial f}{\partial z}\right)}_j=\frac{(1-g_j)(f_{j+1}-f_j)}{2\Delta z}+\frac{(1+g_j)(f_j-f_{j-1})}{2\Delta z},$$

where ${\displaystyle 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 [34]. Additionally, we have implemented variable step size, which greatly enhances the speed and efficiency of the solver. The numerical solver has been benchmarked against existing simulations of the Oldroyd-B model [27].

On 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. 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 [27] we set the initial shape of the jet as,

$$\begin{array}{c}\hfill R=R_0\left(1+0.01cos\left(kz\right)\right),\end{array}$$

(15)

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, $\mathbf{\nabla }\mathbf{u}$, at each grid point. Using these values, the dumbbell configurations are updated using Equation (7a) and the extra stress calculated with Equation (7b). The resulting values for the extra stress at each grid point are introduced in the macroscopic simulation through Equation (14b).

Below we show results 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 (13).

4.1. Deborah Number

Figure 9 shows the evolution of the jet morphology as a function of the Deborah number, $De$. As the elasticity of the fluid is increased −larger $De-$ the formation of 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 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 at the end leads to the formation of droplets.

4.2. Wavelength

The effect of the excitation wavenumber on the formation of droplets is well understood, see for example [27,35,36,37]. 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 10 the jet evolution for two different of 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.

4.3. Non-linear Parameter

Figure 11 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 similar trend 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 11, larger values of b suppress these secondary droplets. Furthermore, another main difference is observed at the pinch sites, with smaller b resulting in more pronounced pinching. This is a result of less elasticity, hence more viscous-like behavior.

5. Discussion

In this work we have considered multiscale simulations of polymeric liquids under homogeneous and non-homogeneous flow conditions. The Langevin equations describing the behavior of an elastic dumbbell are solved on GPUs, which allows for massively parallel simulations. Utilizing the dumbbell equations directly provides for more accurate information on the stretch and orientation of polymer chains than constitutive equations arising from solving for the second moment of the end-to-end vector. For example, the stretching and orientation of dumbbells are captured in a similar manner in which an experimentalist will analyze fluid behavior using scattering techniques.

Under simple shear and uniaxial extension, we showed that there is a significant difference between the linear and nonlinear flow regimes. In shear flow, 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 shows that `molecules’ are more oriented with respect to the flow in the high shear rate, 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 at the end 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 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 on 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 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 capabilities to other types of flow, a more careful coupling of the micro-macro approaches is required, with extra limitations arising from the use of GPUs. The results of this work show that coupling a GPU solver for the Langevin equations with an efficient macroscale solver for the flow field can lead to large scale multiscale simulations in complex flow fields.