Multiscale Modeling of Viscoelastic Fluids

1. Introduction

Although you may be unfamiliar with the term viscoelasticity, viscoelastic fluids (VEFs) are ubiquitous in our daily lives. The shampoo you used this morning, the salad dressing you ate yesterday, and about every biological fluid in your body; all of them are VEFs. As their name implies, viscoelastic materials exhibit both viscous and elastic behaviors. Viscous behavior is related to a fluid’s resistance to flow. The higher the viscosity, the more the fluid resists motion. Honey, for example, has a high viscosity, while water has a low viscosity. Conversely, elasticity pertains to reversible deformations, such as the snapback of a rubber band.

Under applied deformations, VEFs display an instantaneous pure elastic response, followed by a time-dependent mechanical response and energy dissipation, characteristic of viscous liquids. The differences from one VEF to another comes from the relative timescales of these elastic and viscous responses. This “duality” in the behavior of VEFs plays a critical role in their applications. For example, paints can be thin enough to be applied with a brush, yet thick enough to stay on the wall. And although mayonnaise appears semisolid in a jar, it can be easily spread on bread.

The study of VEFs falls within the field of rheology. Rheology investigates how materials deform or flow during and after a load is applied. Measuring rheological properties is pertinent to all materials, from liquids such as water, polymers, and protein solutions to semisolids such as gels and creams and to solid polymers such as resins. Within rheology, at the most basic level, fluids can be divided into Newtonian and non-Newtonian according to their response to flow. From a modeling point of view, all Newtonian fluids are described by the well-known Navier–Stokes equations Bat99. This set of equations works well on systems in which the flow does not alter the dynamics of individual constituents. In contrast, for non-Newtonian fluids, applied fields can alter the local microstructure. Hence, there is no single set of equations that can comprehensively describe all non-Newtonian materials. What these fluids have in common is that their properties emerge from the collective behavior of many microstructural components.

The field of VEFs offers many opportunities for mathematical exploration, particularly in the development of new constitutive models and numerical techniques. The absence of a unified equation for VEFs and their wide range of applications in industrial and biological processes has resulted in extensive research activity in this field. However, there are still many issues that need to be addressed. On the numerical side, key challenges involve the loss of accuracy and convergence of numerical methods due to nonlinearities in the constitutive equations. Another challenge is the change of type of the partial differential equations (PDEs), which sometimes leads to a loss of well-posedness. Additionally, in certain cases, fluid flows result in rapid changes of the solutions in specific regions, making it necessary to implement adaptive mesh techniques. For reviews on this area, the reader is referred to OP02Keu04AOP21.

Fractional calculus is another rapidly growing area in the field of VEFs. It provides a more detailed understanding of the memory effect through the use of fractional derivatives Mai22. The behavior of VEFs falls on a spectrum between that of fully elastic solids and fully viscous fluids. Fractional models provide a unified framework to understand this entire range by varying the order of the fractional derivative. This allows fractional models to achieve comparable accuracy to classical models, but with fewer parameters, making data fitting less complex.

In addition, the field of VEFs offers numerous opportunities for the mathematical analysis of existing constitutive equations. To understand the dynamics of VEFs, it is crucial to evaluate the stability of these systems and understand how the interplay between viscous and elastic properties influences their behavior. However, there is still much to be explored regarding how different flow conditions affect the solutions of these models. Only a limited subset of these equations has been thoroughly studied in terms of their existence and uniqueness RT21. Investigating the connections between VEFs constitutive models and broader concepts in dynamical systems could uncover previously unexplored complexities.

Perturbation analysis is another valuable field commonly used to gain insight and address unresolved issues pertaining to VEFs. Singular perturbation methods are especially well-suited for examining VEFs because they can effectively manage the wide range of temporal and spatial scales inherent in the dynamics of these materials. By analyzing extreme cases of high and low elasticity, we can gain valuable insights into the underlying dynamics, circumventing the need to solve complex systems of equations.

This review aims to introduce the reader to basic concepts related to VEFs, focusing mainly on constitutive modeling. We note that each class of VEFs has unique properties, posing different mathematical challenges. Providing a complete overview of every class of VEFs and their mathematical representations is beyond the scope of this introductory review. Instead, we discuss some commonalities and focus on a specific class of VEFs. When suitable, the reader will be referred to more in-depth reviews and textbooks.

2. Multiscale Modeling of VEFs

From a mathematical perspective, modeling challenges of VEFs arise from the need to describe the dynamics resulting from the complex interactions among microscopic constituents and how such interactions dictate material properties and functions at the macroscale (Fig. 1). In short, there is a need for continuous communication across multiple scales of time and space. Ideally, a constitutive equation for VEFs should provide sufficient insight into the microscopic changes that lead to a given macroscopic response.

Figure 1.

VEFs exist at different scales because their macroscopic responses depend on the dynamics of their microstructural components.

This coupling of micro and macro scales is particularly crucial for some materials. For example, one might like to study biochemical changes in the mucus network and their impact on the mucociliary clearance process, or perhaps investigate how changes in the viscoelastic properties of the cytoplasm affect cellular function. In these situations, modeling platforms that consider dynamics at the microscale are better suited.

To model the microstructure explicitly, we consider a VEF to be composed of two major components: a solvent and some dynamic network at the microscale. We consider the system to be a continuum. This means it can be described by classical mechanics and its state under a deformation is determined by the fundamental hydrodynamic fields of density, momenta, and energy, Bat99BCAH87,

$$\begin{gather*} \text{Conservation of mass:} \qquad \frac{D\rho }{Dt} + \rho \boldsymbol{\nabla } \cdot \mathbf{u} = 0,\\[3.0pt] \text{Conservation of momentum:} \qquad \frac{D \left( \rho \,\mathbf{u}\right)}{Dt} = \boldsymbol{\nabla } \cdot \boldsymbol{\sigma } + \rho \mathbf{f}_{\mathrm{b}},\\[3.0pt] \text{Conservation of energy:} \qquad \rho \frac{De}{Dt} = \boldsymbol{\sigma }:\boldsymbol{\nabla }\mathbf{u} + \boldsymbol{\nabla } \cdot \left( \kappa \boldsymbol{\nabla } T\right). \end{gather*}$$

Here $\rho$ represents the material density, $\mathbf{u}$ the velocity field, $\boldsymbol{\sigma }$ the Cauchy stress tensor, $\mathbf{f}_{\mathrm{b}}$ body forces, $e$ the internal energy per unit mass, $\kappa$ the thermal conductivity and $T$ the temperature. And, the material derivative is defined as,

$$\begin{equation*} \frac{D\left( \cdot \right)}{Dt} = \frac{\partial \left( \cdot \right)}{\partial t} + \mathbf{u} \cdot \boldsymbol{\nabla } \left( \cdot \right). \end{equation*}$$

If the flow is incompressible ($\rho =$ constant), isothermal ($T=$ constant), and in the absence of body forces ($\mathbf{f}_{\mathrm{b}} = 0$), the conservation equations become,

$$\begin{align} \boldsymbol{\nabla } \cdot \mathbf{u} &= 0,\tag{1a}\\ \rho \frac{D\mathbf{u}}{Dt} &=\boldsymbol{\nabla } \cdot \boldsymbol{\sigma }. \tag{1b} \end{align}$$

The Cauchy stress tensor can be decomposed into isotropic and extra stress components,

$$\begin{equation*} \boldsymbol{\sigma } = -p\,\boldsymbol{\delta } + \boldsymbol{\tau }, \end{equation*}$$

where $\boldsymbol{\delta }$ is the identity tensor. At equilibrium, the isotropic component is the thermodynamic pressure, $p$, while the extra stress tensor, $\boldsymbol{\tau }$, vanishes Gra18. Under these considerations, the resulting conservation equations are,

$$\begin{align} \boldsymbol{\nabla } \cdot \mathbf{u} &= 0,{\tag{2a}}\\ \rho \frac{D\mathbf{u}}{Dt} &=-\boldsymbol{\nabla } p + \boldsymbol{\nabla } \cdot \boldsymbol{\tau }. {\tag{2b}} \end{align}$$

For a viscous or Newtonian fluid, the stress is directly proportional to the strain rate, $\dot{\boldsymbol{\gamma }} = \boldsymbol{\nabla } \mathbf{u} + \left( \boldsymbol{\nabla } \mathbf{u}\right)^{\intercal }$, so that $\boldsymbol{\tau } = \eta \, \dot{\boldsymbol{\gamma }}$. In this case $\boldsymbol{\tau }$ is known as a viscous stress and $\eta$ is the fluid’s viscosity. We note that in some references, the strain rate tensor is defined as $\mathbf{D} =\frac{1}{2}\dot{\boldsymbol{\gamma }}$, and $\boldsymbol{\tau } = 2\,\eta \,\mathbf{D}$. Here we follow the notation proposed in BCAH87 and use $\dot{\boldsymbol{\gamma }}$. The resulting conservation of momentum equation is known as the incompressible Navier–Stokes equation,

$$\begin{equation} \rho \frac{D\mathbf{u}}{Dt} = -\boldsymbol{\nabla } p + \eta \boldsymbol{\Delta } \mathbf{u}. {\tag{3}} \end{equation}$$

Inspection of 3 shows that differences between materials will only affect the value of $\eta$, but not the functional form of the governing equations. Then, by varying the viscosity constant, the same numerical algorithms can describe different materials, as long as they behave as Newtonian fluids.

Coming back to 2a–2b, in the case of VEFs there is a need to introduce an extra term in stress, $\boldsymbol{\tau }_p$, which accounts for the contributions from the microstructure, $\boldsymbol{\tau } = \eta \,\dot{\boldsymbol{\gamma }} + \boldsymbol{\tau }_p$. In our discussion of 3 we noted that for Newtonian fluids there is a single, general constitutive equation capable of describing different Newtonian materials, so that 3 suffices to describe a large class of fluids. In contrast, there is not a single all-encompassing constitutive equation for viscoelastic materials. The details of $\boldsymbol{\tau }_p$, particularly the constitutive equation that defines it, will differ depending on the material or the specific flow conditions. Thus, the formulation of constitutive equations describing different viscoelastic materials is a prolific area of research BCAH87Lar99.

Broadly speaking, we understand relatively well the connection between VEFs and external flow fields: the dynamics of the underlying microstructure directly control the behavior of a VEF under these deformation fields. Therefore, when dealing with constitutive equations for VEFs, challenges arise at three main levels:

•

At the level of their derivation, since different materials require different mathematical descriptions of their microstructure;

•

At the level of their numerical simulation, since the resulting equations are of different class and each has its unique numerical challenges;

•

At the level of their mathematical treatment, since researchers have only established the existence and uniqueness of the solution of VEF constitutive equations in only a few cases (See for example Ren85LM00RT21).

The mathematical representation of the microstructure, the length scale at which this representation will be rendered, and the numerical methods used to solve the resulting constitutive equations will all depend on the particulars of the fluid and the chosen mathematical description. Accordingly, many studies have tackled one aspect or another and even combinations of them. However, a full description is out of the scope of this review. Moving forward, we focus our attention on one type of models which originates from kinetic theory and a specific category of VEFs: polymeric fluids.

2.1. Coarse-grained representation of the microstructure

To better understand how we can develop mathematical models of VEFs based on representations of the microstructure, here we will focus on polymeric fluids. In these fluids, polymer chains compose the microstructure. The underlying dynamics driving the material’s response to deformation result from both individual chain configurations and interchain dynamics. Among others, these include coiling and uncoiling processes, hindered motion due to physical entanglements, hydrodynamic effects caused by the presence of other molecules, and in some cases, physical cross-linking between the polymer chains.

Within the context of kinetic models and polymeric fluids, various approaches have developed to describe the coupling between macroscopic responses and microstructure dynamics. One family of models is 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. These models are based on molecular physics by considering the interaction between individual polymer chains and the surrounding fluid. The extra stress arising from the polymer molecules, i.e., the microstructure, depends strongly on their spatial configuration, the most important features being their orientation and their extension.

The simplest of these models considers only two beads and it is known as the elastic dumbbell model. The configuration of each dumbbell is fully specified by its stretching and orientation. Although it is widely recognized that a dumbbell is too simple to be able to describe any complicated dynamics in polymeric molecules, it is also well known that stretching and orientation alone suffice to give a qualitative description of steady-state rheological properties and flows with slow characteristic timescales BCAH87. Accordingly, these models had been extensively used to develop “an elementary but broad understanding of the relation between macro-molecular motions and rheological phenomena” BCAH87.

2.2. Dumbbell models

To model a given microstructure this class of models uses a coarse-grained approximation at the mesoscale consisting of noninteracting elastic dumbbells. The VEF system will be described by the dynamics of these dumbbells in a solvent. 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)$.

To keep it short, this review only discusses homogeneous flows. This means that the velocity gradient is assumed to be the same everywhere. We refer the reader to BAB91 for a discussion on how to introduce spatial variations into the dumbbell equations.

To capture the dynamics of each dumbbell, we start from a balance of forces at the inertia-less limit. That is we use Newton’s second law, $m \vec{\mathbf{a}} = \sum \mathbf{F}$, but assume the mass is negligible BCAH87Ött96,

$$\begin{align} 0 &= \underbrace{- \zeta \left[ d \mathbf{r}_{i}(t) + \left( \boldsymbol{\nabla } \mathbf{u}\right)^{\intercal } \cdot \mathbf{r}_{i}(t) \right]}_{\text{Drag force}} - \underbrace{\mathbf{F}_i(t)}_{\text{Spring force}} \\[3.0pt] &\quad {} + \underbrace{\sqrt {4 k_BT\, \zeta }\, d\mathbf{W}_{i}(t)}_{\text{Thermal noise}} \qquad i =1,2. {\tag{4}} \end{align}$$

Here $\zeta$ is the drag coefficient, $\mathbf{u}$ is the fluid velocity, $\mathbf{F}(t)$ denotes a functional form of the spring force, and $k_B T$ is the thermal energy. $\mathbf{W}(t)$ denotes a Wiener process where each component of the vector $\mathbf{W}(t)$ is a random number drawn from a normal distribution with zero mean and variance equal to 1.

Equation 4 shows that changes in the configuration of the dumbbell, namely its orientation and extension, are the result of three competing forces. The drag force, imposed by the solvent molecules onto the beads, has the tendency of aligning the dumbbells with the macroscopic flow. The thermal or Brownian force tends to randomize their configuration. And, the spring force tends to bring both beads together, which counteracts the stretching effects of the drag and thermal forces.

Rearranging 4, gives the following two stochastic differential equations (SDEs) describing the evolution of each bead’s position,

$$\begin{align} d\, \mathbf{r}_{1}(t) &= \left( \boldsymbol{\nabla } \mathbf{u}\right)^{\intercal } \cdot \mathbf{r}_{1}(t) - \frac{\mathbf{F}_{1}(t)}{\zeta } + \sqrt {\frac{4k_BT}{\zeta }} \,d\mathbf{W}_{1}(t), {\tag{5a}}\\ d\, \mathbf{r}_{2}(t) &= \left( \boldsymbol{\nabla } \mathbf{u}\right)^{\intercal } \cdot \mathbf{r}_{2}(t) - \frac{\mathbf{F}_{2}(t)}{\zeta } + \sqrt {\frac{4k_BT}{\zeta }} \,d\mathbf{W}_{2}(t). {\tag{5b}} \end{align}$$

Since the spatial location of the dumbbells is not relevant under the homogeneous flow assumption, the only variable of interest is the end-to-end vector, $\mathbf{Q}$. By subtracting 5a from 5b, we obtain the SDE describing the evolution of $\mathbf{Q}$,

$$\begin{equation} d {\mathbf{Q}}_t = \left( {\boldsymbol{\nabla }} {\mathbf{u}} \right)^{\intercal } \cdot {\mathbf{Q}}_t - \frac{2}{\zeta }\mathbf{F}\left( {\mathbf{Q}}_t\right) + \sqrt {\frac{4k_BT}{\zeta }}\, d\mathbf{W}_t. {\tag{6}} \end{equation}$$

The only remaining task is to establish the form of the spring law, denoted as $\mathbf{F}(\mathbf{Q})$, where $\mathbf{F}_1 = -\mathbf{F}_2 = \mathbf{F}$.

2.2.1. Hookean dumbbells

For Hookean dumbbells $\mathbf{F}(\mathbf{Q}) = H\mathbf{Q}$, where $H$ is the spring constant. With this, 6 becomes,

$$\begin{equation} d {\mathbf{Q}}_t = \left( {\boldsymbol{\nabla }} {\mathbf{u}} \right)^{\intercal } \cdot {\mathbf{Q}}_t - \frac{2H}{\zeta } {\mathbf{Q}}_t+ \sqrt {\frac{4k_BT}{\zeta }}\, d\mathbf{W}_t. {\tag{7}} \end{equation}$$

For convenience, we will make these equations non-dimensional. To couple these equations with the conservation equations, we will scale the time using the macroscopic timescale, $t^*$. In addition, we introduce the following characteristic microscopic time and length scales, respectively,

$$\begin{equation*} \lambda = \frac{\zeta }{4H}, \qquad L_m = \sqrt {\frac{k_B\,T}{H}}. \end{equation*}$$

The nondimensional variables are then given by,

$$\begin{equation*} \tilde{\mathbf{Q}} = {\mathbf{Q}} \cdot \left( \sqrt {\frac{H}{k_BT}}\right), \quad \tilde{t} = \frac{t}{t^*}. \end{equation*}$$

Scaling 7 and dropping the tildes gives,

$$\begin{equation} d\mathbf{Q}_t = \left( \boldsymbol{\nabla } \mathbf{u} \right)^{\intercal } \cdot \mathbf{Q}_t - \frac{1}{2\,De}\mathbf{Q}_t+ \sqrt {\frac{1}{De}}\,d\mathbf{W}_t, {\tag{8}} \end{equation}$$

where the nondimensional group $De = \lambda /t^*$ is the so-called Deborah number. Since it is the ratio of micro-to-macro characteristic timescales, this dimensionless group compares how long it takes for a material to adapt to deformations relative to the process’s characteristic timescale.

Note that, we could have chosen a different macroscopic timescale, namely, $L/U$, where $L$ and $U$ are, respectively, characteristic macroscopic length and velocity. Here, the resulting nondimensional group, $Wi = \lambda \,U/L$, is called the Weissenberg number, and it represents the ratio of elastic to viscous forces. For many applications $De =Wi$ and it is very common to confuse these two nondimensional groups. To better understand the difference between $De$ and $Wi$, we recommend reading Poo12.

2.2.2. FENE dumbbells

The linear spring law used in the Hookean dumbbell model is unphysical, since it allows the end-to-end vector, $\mathbf{Q}$, to stretch without limit. One modification of this law uses finitely extensible nonlinear elastic (FENE) springs laws. FENE-type models are derived by introducing Warner’s force law BCAH87,

$$\begin{equation} F \left( {\mathbf{Q}}\right) = \frac{H {\mathbf{Q}}}{1- \left( {Q}/Q_{\mathrm{max}}\right)^2}, {\tag{9}} \end{equation}$$

where $Q^2 = \left|\mathbf{Q}\right|^2 = \mathbf{Q} \cdot \mathbf{Q}$ and $Q_{\mathrm{max}}$ is the maximum allowed extension of the dumbbell.

Introducing this spring law into 6 and nondimensionalizing as before gives

$$\begin{equation} d\mathbf{Q}_t = \left( \boldsymbol{\nabla } \mathbf{u} \right)^{\intercal } \cdot \mathbf{Q}_t - \frac{1}{2De}\left( \frac{\mathbf{Q}_t}{1- Q^2/b}\right) + \frac{1}{\sqrt {De}}d\mathbf{W}_t, {\tag{10}} \end{equation}$$

where $b = H\,Q_{\mathrm{max}}^2/(k_BT)$.

Finally, we note that, in this dumbbell formulation, we only considered the most basic dynamics. Still, dumbbell models can integrate other physical assumptions. For instance, certain models account for a nonisotropic drag BD83 while others include breaking and reforming dynamics VMC07.

Before we move forward, we note that dumbbell models are not the only, nor the most prominent, class of models used to describe polymeric systems. We have chosen this particular class of models for specific reasons. Among them, the Hookean dumbbell model, although physically unrealistic, stands out as the only model in its class with an exact closure BCAH87BAB91. This allows for a more coherent explanation of the transition through the different length scales that are discussed in this review. We acknowledge that dumbbell models have their limitations, especially when dealing with complex flows. These drawbacks primarily stem from their oversimplification of several key aspects in polymer physics. For instance, they do not account for entanglements, excluded volume effects, nor internal chain dynamics. For a deeper understanding of more physically relevant models, we highly recommend the book by Doi and Edwards DEE88. For readers interested in the analytical and numerical solutions to SDEs-based models, the book by Öttinger is another excellent starting point Ött96.

2.3. Fokker–Planck representation

In the previous section, we discussed how SDEs can describe the evolution of the end-to-end vector $\mathbf{Q}_t$. This implies that $\mathbf{Q}_t$ is a stochastic or random process. To fully capture the system’s dynamics, it is necessary to solve thousands of SDEs. Now, if instead of “following” individual realizations of this process, i.e., solving 8 or 10, we decide to “follow” ensembles of realizations, we would need a different set of equations. This new set of equations should describe the same system, but instead of using stochastic variables, it uses deterministic variables that fluctuate because of stochasticity. To accomplish this we use the probability density function (PDF), $\psi \left( \mathbf{Q}, t\right)$, which describes the probability of finding dumbbells with configurations in the interval $\left( \mathbf{Q}, \mathbf{Q} + d\mathbf{Q}\right)$ at time $t$. Risken’s book Ris84 excellently explains the differences between these two representations, which we summarize in Fig. 2.

Figure 2.

Levels of description of a system using Langevin and Fokker–Planck Equations. Figure adapted from Ris84.

Let $\Psi \left( \mathbf{r},\mathbf{Q},t\right)$ represent the configuration number density function, so that the number density of dumbbells with end-to-end vector $\mathbf{Q}$ and center of mass at position $\mathbf{r}$, at time $t$ is given by,

$$\begin{equation*} n\left( \mathbf{r},t\right) = \int \Psi \, d\mathbf{Q}. \end{equation*}$$

In the homogeneous case the spatial dependence can be neglected, so that $\Psi = n\, \psi \left(\mathbf{Q},t\right)$.

When dealing with PDFs, one can do an expansion with similar flavor as the Taylor expansion taught in calculus. This expansion is known as the Kramers–Moyal expansion Ris84. If the expansion is truncated after the second term, the resulting equation is called a Fokker–Planck equation, also known as the forward Kolmogorov equation. A brief summary of how Brownian dynamics can be described by Langevin and their corresponding Fokker–Planck equations is given in MDV20, but for a more comprehensive treatment, see Ris84.

The general form of a Fokker–Planck equation on the variable $\mathbf{Q}$ is

$$\begin{align} \frac{\partial \psi }{\partial t} = -\frac{\partial }{\partial \mathbf{Q}} \left[ A\left( \mathbf{Q},t\right)\, \psi \left(\mathbf{Q},t\right)\right] + \frac{1}{2} \frac{\partial ^2}{\partial \mathbf{Q}^2} \left[ \mathbf{B}\left( \mathbf{Q},t\right) \, \mathbf{B}^{\intercal }\left( \mathbf{Q},t\right)\, \psi \left( \mathbf{Q},t\right)\right], \tag{11} \end{align}$$

which corresponds to the Langevin equation,

$$\begin{equation} d\mathbf{Q}_t = A\left( \mathbf{Q}_t,t\right) + \mathbf{B}\left( \mathbf{Q}_t,t\right)\,d\mathbf{W}_t. {\tag{12}} \end{equation}$$

Thus, we can use 12 together with 8 or 10 to obtain the Fokker–Planck equations corresponding to the Hookean and FENE models.

Hookean dumbbells

$$\begin{align} \frac{\partial \psi }{\partial t} = -\frac{\partial }{\partial \mathbf{Q}} \left[ \left( \left( \boldsymbol{\nabla } \mathbf{u} \right)^{\intercal } \cdot \mathbf{Q} - \frac{1}{2De}\,\mathbf{Q}\right)\, \psi \right] + \frac{1}{De} \frac{\partial ^2\psi }{\partial \mathbf{Q}^2} {\tag{13}} \end{align}$$

FENE dumbbells

$$\begin{align} \frac{\partial \psi }{\partial t} = -\frac{\partial }{\partial \mathbf{Q}} \left[ \left( \left( \boldsymbol{\nabla } \mathbf{u} \right)^{\intercal } \cdot \mathbf{Q}\vphantom {\frac{1}{2De}} \right.\right. - \left. \left. \frac{1}{2De}\left( \frac{\mathbf{Q}}{1- Q^2/b}\right)\right)\, \psi \right] +\frac{1}{De} \frac{\partial ^2\psi }{\partial \mathbf{Q}^2}. {\tag{14}} \end{align}$$

Just like it was the case for the Langevin equations, we can relate each term in the Fokker–Planck equations with their physical counterpart. Since the derivatives are on $\mathbf{Q}$, this means these terms depend on the dumbbell’s configuration. In 13–14 the first term on the right-hand side describes transport by the macroscopic flow. The second term relates to the spring force, whose effect is to concentrate the distribution function about $\mathbf{Q} = 0$. The last term represents diffusion in configuration space; its effect is to “spread out” the distribution function.

Just as with the previous section, we will avoid going into details about solution strategies for these Fokker–Planck equations in order to maintain brevity and accessibility in this review. However, the book by Risken Ris84 is a good starting point for those interested in this subject.

2.4. Macroscopic representation

So far, we have described two ways in which we can capture the dynamics of a microstructure comprising chains represented by elastic dumbbells. The first representation considers the equation of motion of individual dumbbells given by SDEs. The second representation considers ensembles of dumbbells and gives the evolution equation for their PDFs as partial differential equations (PDEs). The next logical step is to discuss how to couple these dynamics with the macroscopic flow field described by the conservation of mass and momentum equations.

As discussed in previous sections, the equations for the conservation of mass and momentum of an isothermal, incompressible, viscoelastic fluid, in the absence of external forces, are given by

$$\begin{align} \boldsymbol{\nabla } \cdot \mathbf{u} &= 0,\tag{15a}\\ \rho \frac{D\mathbf{u}}{Dt} &=-\boldsymbol{\nabla } p + \eta _s \boldsymbol{\Delta } \mathbf{u} + \boldsymbol{\nabla } \cdot \boldsymbol{\tau }_p. \tag{15b} \end{align}$$

Here we have used the notation $\eta _s$ to denote the viscosity of the Newtonian component of the VEF, i.e., the solvent. We will scale this equations as,

$$\begin{equation*} \tilde{t} = \frac{t}{t^*},\quad \tilde{\mathbf{x}} = \frac{\mathbf{x}}{U t^*}, \quad \tilde{p} = \frac{p}{\rho U^2}, \quad \tilde{\mathbf{u}} = \frac{\mathbf{u}}{U}, \quad \tilde{\boldsymbol{\tau }}_p = \frac{\boldsymbol{\tau }_p}{nk_BT}, \end{equation*}$$

where $t^*$ and $U$ are characteristic time and velocity, $n$ is the dumbbells’ number density, and $k_BT$ is the thermal energy. Dropping the tildes, gives

$$\begin{align} \boldsymbol{\nabla } \cdot \mathbf{u} &= 0,{\tag{16a}}\\ \frac{D\mathbf{u}}{Dt} &= -\boldsymbol{\nabla } p + \frac{\beta }{Re}\, \boldsymbol{\Delta } \mathbf{u} + \frac{1-\beta }{Re\,De}\,\boldsymbol{\nabla } \cdot \boldsymbol{\tau }_p. {\tag{16b}} \end{align}$$

Here, $\beta = \eta _s/\eta _0$ is the ratio of the solvent to the total viscosity, where $\eta _0 = \eta _s + \eta _p$ is the zero-shear rate viscosity of the fluid, with $\eta _p = nk_BT\,\lambda$ being the polymer contribution to the viscosity. And, $Re = \rho U^2t^*/\eta _0$ is the Reynolds number, which is the ratio of inertial to viscous forces.

Together with a constitutive equation for $\boldsymbol{\tau }_p$, 16 give the mathematical description, in time and space, of the resulting flow field of a VEF. However, these equations are applicable only when considering the fluid as a continuous medium. This means that we need to find an expression for $\boldsymbol{\tau }_p$ that can “translate” the dynamics at the microscopic and mesoscopic levels, described in previous sections, to the continuum level.

As dumbbells move about the solvent fluid, there is a drag force imposed on the beads by the fluid’s velocity, $\mathbf{u}$. This drag on the beads causes an extra stress on the solvent, which in turns changes its velocity. This exchange between the dumbbells and the fluid depends on the momentum transferred between the beads in each dumbbell, which is modulated by the connector vector $\mathbf{Q}$. Because of this dependence on $\mathbf{Q}$, momentum transfer perpendicular to the flow exerts additional viscous forces, i.e., resistance to flow. While momentum transfer parallel to the flow gives elastic properties to the fluid.

In order to understand how single dumbbells contribute to the stress, let’s start by examining the significance of each entry in the stress tensor. In three-dimensional space, the stress tensor is a $3 \times 3$ matrix. If we consider a piece of fluid as a cube, then the $i,j$ entry corresponds to the stress resulting from a force in the $\hat{i}$ direction imposed in the face with a normal vector in the $\hat{j}$ direction; see Fig. 3.

Figure 3.

Components of the stress tensor, where $\tau _{ij}$ are stresses resulting from forces in the $i$-direction, acting in the face of the volume with normal vector in the $j$-direction. In general, $\boldsymbol{\tau }$ is a symmetric tensor; see proof for symmetry in Appendix A1 of OP02.

The contribution to the stress from a single dumbbell is then given by Kramer’s relation BCAH87Ött96,

$$\begin{equation*} \boldsymbol{\mathcal{T}} = H f(\mathbf{Q})\otimes \mathbf{Q} - k_BT\boldsymbol{\delta }, \end{equation*}$$

where $F(\mathbf{Q}) = Hf(\mathbf{Q})$ is the spring force, $\boldsymbol{\delta }$ the identity matrix, and $\otimes$ indicates the tensor product of vectors. And, we can connect $\mathbf{Q}$, at the microscopic or mesoscopic scales, to the macroscopic stress using ensemble averages. In the configuration space represented by $\mathbf{Q}$, the ensemble average of any function $g(\mathbf{Q})$ is given by

$$\begin{equation} \left< g(\mathbf{Q})\right> = \int g(\mathbf{Q})\,\Psi (\mathbf{r} , \mathbf{Q}, t)\, d\mathbf{Q}. {\tag{17}} \end{equation}$$

We use 17 to find the stress resulting from the collective dynamics all of dumbbells as

$$\begin{align*} \boldsymbol{\tau }_p = \left<\boldsymbol{\mathcal{T}} \right> &= n\int H f(\mathbf{Q})\,\mathbf{Q}\,\psi \, d\mathbf{Q} - nk_BT \boldsymbol{\delta },\\ &= nH\left< f( {\mathbf{Q}})\, {\mathbf{Q}}\right> - nk_BT \boldsymbol{\delta }, \end{align*}$$

where $n\left< F(\mathbf{Q})\, \mathbf{Q}\right>$ is the contribution from the tension on the spring with spring law $F(\mathbf{Q})$, and $nk_BT$ capture effects due to Brownian motion.

•

For Hookean springs, $f(\mathbf{Q}) = \mathbf{Q}$ and we obtain,$$\begin{equation} \boldsymbol{\tau }_p = \left< {\mathbf{Q}} {\mathbf{Q}}\right> - \boldsymbol{\delta }. {\tag{18}} \end{equation}$$

•

For FENE springs $f(\mathbf{Q}) = \mathbf{Q}/(1-\left( Q/Q_{\mathrm{max}}\right)^2)$, so that,$$\begin{equation} \boldsymbol{\tau }_p = \left< \frac{{\mathbf{Q}} {\mathbf{Q}}}{1 - {Q}^2/b}\right> - \boldsymbol{\delta }. {\tag{19}} \end{equation}$$

Here, we used the same nondimensionalization as before.

To find a constitutive equation for $\boldsymbol{\tau }_p$, we start with the general form of the Fokker–Planck equation for elastic dumbbells,

$$\begin{equation*} \frac{\partial \psi }{\partial t} = -\frac{\partial }{\partial \mathbf{Q}} \left[ \left( \left( \boldsymbol{\nabla } \mathbf{u} \right)^{\intercal } \cdot \mathbf{Q} - \frac{f(\mathbf{Q})}{De}\right)\, \psi \right] + \frac{1}{De} \frac{\partial ^2\psi }{\partial \mathbf{Q}^2}. \end{equation*}$$

Recall that in this representation we assume the spring force is of the form $\mathbf{F}(\mathbf{Q}) = Hf(\mathbf{Q})$. To find an expression for $\left<\mathbf{Q}\mathbf{Q}\right>$, we multiply 20 throughout by $\mathbf{Q}\mathbf{Q}$, and integrate over the configuration space. We use the divergence theorem and the fact that $\psi \to \infty$ as $\mathbf{Q} \to \infty$ to obtain BCAH87

$$\begin{equation*} \frac{\partial \left<\mathbf{Q}\mathbf{Q}\right>}{\partial t}= \left( \boldsymbol{\nabla } \mathbf{u} \right)^{\intercal } \cdot \left<\mathbf{Q}\mathbf{Q}\right> + \left<\mathbf{Q}\mathbf{Q}\right> \cdot \left( \boldsymbol{\nabla } \mathbf{u}\right) - \frac{1}{De} \left<\mathbf{Q} f(\mathbf{Q})\right> + \frac{1}{{De}} \boldsymbol{\delta }. \end{equation*}$$

If we define the upper convected derivative as,

$$\begin{equation*} \left( \cdot \right)_{(1)} = \frac{\partial \left( \cdot \right)}{\partial t} - \left( \boldsymbol{\nabla } \mathbf{u} \right)^{\intercal } \cdot \left( \cdot \right) -\left( \cdot \right) \cdot \left( \boldsymbol{\nabla } \mathbf{u} \right), \end{equation*}$$

we arrive at,

$$\begin{equation} \left<\mathbf{Q}\mathbf{Q}\right>_{(1)} = - \frac{1}{De} \left<\mathbf{Q} f(\mathbf{Q})\right> + \frac{1}{{De}} \boldsymbol{\delta }. {\tag{20}} \end{equation}$$

In the next sections we show how we can use 20 to formulate constitutive equations for $\boldsymbol{\tau }_p$ corresponding the Hookean and FENE dumbbells.

2.4.1. Hookean dumbbells

In this case $f(\mathbf{Q}) = \mathbf{Q}$, so that

$$\begin{equation*} \left<\mathbf{Q}\mathbf{Q}\right>_{(1)} = - \frac{1}{De} \left<\mathbf{Q} \mathbf{Q}\right> + \frac{1}{{De}} \boldsymbol{\delta }. \end{equation*}$$

Using 18 we obtain the constitutive equation for the extra stress tensor of Hookean dumbbells,

$$\begin{equation} De\,\boldsymbol{\tau }_{p,(1)} + \boldsymbol{\tau }_p = De\, \dot{\boldsymbol{\gamma }}. {\tag{21}} \end{equation}$$

This is the well-known Upper Convected Maxwell (UCM) model BCAH87. We should emphasize that the UCM model provides an exact closure to the Hookean dumbbell model because 21 can be directly obtained from 8. As we will see below, this is not the case for the FENE dumbbells.

Finally, if instead of considering a constitutive equation for only $\boldsymbol{\tau }_p$, we consider the total extra stress, $\boldsymbol{\tau } = \eta _s \boldsymbol{\Delta }\mathbf{u} + \boldsymbol{\tau }_p$, the constitutive equation for $\boldsymbol{\tau }$ is known as the Oldroyd-B Model. For a discussion of the mathematical considerations and challenges arising from the description of VEF using this model, see RT21.

2.4.2. FENE dumbbells

For FENE dumbbells, 20 gives,

$$\begin{equation*} \left<\mathbf{Q}\mathbf{Q}\right>_{(1)} =- \frac{1}{De} \left<\frac{\mathbf{Q} \mathbf{Q}}{1-Q^2/b}\right> + \frac{1}{{De}} \boldsymbol{\delta }. \end{equation*}$$

Because of the nonlinear term, it is not possible to obtain a close-form constitutive equation of $\boldsymbol{\tau }_p$ for FENE dumbbells. Instead, several closures have been suggested to allow the formulation of macroscopic constitutive equations. Here we will discuss the so-called Peterlin approximation, which results in the FENE-P model BCAH87. Other FENE closures are discussed in DLY05.

Peterlin proposed a separate average of the numerator and denominator of the spring law BCAH87,

$$\begin{align} \left<\mathbf{Q}\mathbf{Q}\right>_{(1)} &= - \frac{1}{De} \frac{\left<\mathbf{Q} \mathbf{Q}\right>}{1-\left<Q^2\right>/b} + \frac{1}{De}\boldsymbol{\delta },\tag{22}\\ \boldsymbol{\tau }_p & = \frac{\left<{\mathbf{Q}} {\mathbf{Q}}\right>}{1-\left<Q^2\right>/b} - \boldsymbol{\delta }. {\tag{23}} \end{align}$$

In this way, instead of restricting the length of individual dumbbells to be less than $Q_{\mathrm{max}}$, the Peterlin’s approximation relaxes the restriction where only the average dumbbell length needs to be less than the prescribed maximum extension. Individual dumbbell lengths can thus exceed $Q_{\mathrm{max}}$ as long as the average stays within bounds.

For convenience, we define a nondimensional configuration tensor, $\mathbf{A}$, as

$$\begin{equation} \mathbf{A} \equiv d\,\frac{\left< {\mathbf{Q}} {\mathbf{Q}}\right>}{\left< {Q}^2\right>_0}, \tag{24} \end{equation}$$

where $\left< {Q}^2\right>_0$ is the mean-square end-to-end spring length at equilibrium (absence of flow) and $d=3$ is the dimensionality BCAH87. Note that if we scale the end-to-end vector as before, we have a nondimensional conformation tensor $\mathbf{A} =3 \left< \mathbf{Q} \mathbf{Q} \right>$. The constitutive equation for $\mathbf{A}$ in three-dimensions is then found as,

$$\begin{align} & {\boldsymbol{\tau }}_p = \frac{ \mathbf{A}}{1 - \operatorname {trace}\left(\mathbf{A}\right)/(3b)} - \boldsymbol{\delta }.\tag{25a}\\[3.0pt] & De \mathbf{A}_{({1})} + \frac{\mathbf{A}}{1 - \operatorname {trace}\left(\mathbf{A}\right)/(3b)} = \boldsymbol{\delta }. \tag{25b} \end{align}$$

To elucidate the trade-offs involved in the Peterlin approximation, in Fig. 4 we show solutions of 8, 10, 21, and 25 under simple shear flow. We compare the macroscopic closures given by the UCM and FENE-P models with their corresponding stochastic counterparts: the Hookean and FENE models. Since the UCM model is an exact closure of the Hookean dumbbell model, solutions to 8 and 21 will always agree with each other, as shown in Fig. 4(A). On the other hand, the approximation used in the FENE-P formulation leads to deviation between solutions of 10 and 25. Fig. 4(B) shows that these differences are more noticeable at higher deformation rates, when dumbbells are close to their maximum extension.

Figure 4.

Solutions in simple shear flow for macroscopic models (solid lines) compared to stochastic simulations (dashed lines). For simple shear flow, the velocity is prescribed as $\mathbf{u} = \left[ \dot{\gamma _0}\,y, 0,0\right]^{\intercal }$. Defining $De = \lambda \,\dot{\gamma }_0$ gives the nondimensional velocity as $\mathbf{u} = \left[ y, 0,0\right]^{\intercal }$. (A) Hookean dumbbell vs. UCM model. (B) FENE dumbbell vs. FENE-P model.

As mentioned above, FENE-P is not the only closure proposed for the FENE model. The mathematical analysis of various closures for the FENE model remains an active area of research, involving the exploration of different approximations and/or higher order moments, e.g., DLY05.

3. Conclusions

This review aims to introduce the reader to the fundamental aspects of mathematical modeling of viscoelastic fluids (VEFs). The most important factor being that the underlying microstructure of VEFs is what determines their properties at the macroscale. This microstructure can comprise polymer molecules, colloidal particles, emulsion drops, etc. The common characteristic is that these structures are larger than the solvent molecules and, as a result, bring about additional stresses to the system. This results in two components of the extra stress, one from the Newtonian component (solvent) and the other from the microstructure. We discussed three levels of description used in modeling this microstructure. The first description focuses on the evolution of individual dumbbells using Langevin-type SDEs. The second description is concerned with the evolution of the PDF of the dumbbell’s configuration through Fokker–Planck equations. The third description provides information at the macroscopic level using PDEs.

Among the different representations, macroscopic constitutive models offer a higher level of computational feasibility, enabling us to find solutions in complex flows or geometries. However, by using these models, one must make a compromise regarding how accurately we can describe the underlying molecular physics. Similarly, although computationally expensive, Langevin or Fokker–Planck descriptions are more amenable to incorporating additional degrees of freedom. This allows us to develop models that better capture the intricacies of physical processes. Hence, in finding the appropriate level of description from a mathematical standpoint, it is important to strike a balance between the complexity of molecular information and the computational costs involved.